import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
%matplotlib inline
# 関数をそれぞれ次にように定義する.
def f1(y1, y2):
return -10 * y1+ 10 * y2
def f2(y1, y2, y3):
return -y1*y3+28*y1-y2
def f3(y1, y2, y3):
return y1*y2- (8/3)*y3
n = 10000
#刻み幅
h = 0.005
# 時刻t
t = np.zeros(n)
t[1] = 50
y1 = np.zeros(n)
y2 = np.zeros(n)
y3 = np.zeros(n)
#初期値
y1[1] = 10
y2[1] = 29
y3[1] = 4
# ルンゲクッタ法
for i in range(n) :
i= i + 1
if i < n-1 :
k1 = h * f1(y1[i],y2[i])
l1 = h * f2(y1[i],y2[i],y3[i])
m1 = h * f3(y1[i],y2[i],y3[i])
s1 = y1[i] + k1/2
s2 = y2[i] + l1/2
s3 = y3[i] + m1/2
k2 = h * f1(s1,s2)
l2 = h * f2(s1,s2,s3)
m2 = h * f3(s1,s2,s3)
s1 = y1[i] + k2/2
s2 = y2[i] + l2/2
s3 = y3[i] + m2/2
k3 = h * f1(s1,s2)
l3 = h * f2(s1,s2,s3)
m3 = h * f3(s1,s2,s3)
s1 = y1[i] + k3
s2 = y2[i] + l3
s3 = y3[i] + m3
k4 = h * f1(s1,s2)
l4 = h * f2(s1,s2,s3)
m4 = h * f3(s1,s2,s3)
y1[i +1] = y1[i] + (k1 + 2*k2 + 2*k3 + k4)/6
y2[i +1] = y2[i] + (l1 + 2*l2 + 2*l3 +l4) / 6
y3[i +1] = y3[i] + (m1 + 2*m2 + 2*m3 + m4) /6
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(y1,y2,y3)
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x113bced30>]
using Plots
# 関数をそれぞれ次にように定義する.
function f1(y1, y2)
return -10 * y1+ 10 * y2
end
function f2(y1, y2, y3)
return -y1*y3+28*y1-y2
end
function f3(y1, y2, y3)
return y1*y2- (8/3)*y3
end
n = 10000
#刻み幅
h = 0.005
y1 = Float64[]
y2 = Float64[]
y3 = Float64[]
#初期値
push!(y1, 3)
push!(y2, 10)
push!(y3, 45)
for i in 1:n
if i < n
push!(y1, y1[i] + h*f1(y1[i],y2[i]))
push!(y2, y2[i] + h*f2(y1[i],y2[i],y3[i]))
push!(y3, y3[i] + h*f3(y1[i],y2[i],y3[i]))
end
end
plot(y1,y2,y3)
using Plots
# 関数をそれぞれ次にように定義する.
function f1(y1, y2)
return -10 * y1+ 10 * y2
end
function f2(y1, y2, y3)
return -y1*y3+28*y1-y2
end
function f3(y1, y2, y3)
return y1*y2- (8/3)*y3
end
n = 10000
#刻み幅
h = 0.005
y1 = Float64[]
y2 = Float64[]
y3 = Float64[]
#初期値
push!(y1, 3)
push!(y2, 10)
push!(y3, 45)
for i in 1:n
if i < n
k1 = h*f1(y1[i],y2[i])
l1 = h*f2(y1[i],y2[i],y3[i])
m1 = h*f3(y1[i],y2[i],y3[i])
k2 = h*f1(y1[i]+h,y2[i]+k1)
l2 = h*f2(y1[i]+h,y2[i]+k1,y3[i])
m2 = h*f3(y1[i]+h,y2[i]+k1,y3[i])
push!(y1, y1[i] + 0.5*(k1+k2))
push!(y2, y2[i] + 0.5*(l1+l2))
push!(y3, y3[i] + 0.5*(m1+m2))
end
end
plot(y1,y2,y3)