ローレンツアトラクタをルンゲクッタ法とオイラー法でそれぞれ実装してみました.

ルンゲクッタ法でローレンツアトラクタ

In [57]:
# 関数をそれぞれ次にように定義する.
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
In [62]:
import numpy as np

n = 10000

#刻み幅
h = 0.1

# 時刻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
In [63]:
# ルンゲクッタ法
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,s3,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,s3,s3)
        m3 = h * f3(s1,s2,s3)
        s1 = y1[i] + k3/2
        s2 = y2[i] + l3/2
        s3 = y3[i] + m3/2

        k4 = h * f1(s1,s2)
        l4 = h * f2(s1,s3,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
In [64]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)

ax.plot(y1,y2,y3)
Out[64]:
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x1158b1128>]

オイラー法でローレンツアトラクタ

In [68]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)

h = 0.01
n =10000

y1 = np.zeros(n)
y2 = np.zeros(n)
y3 = np.zeros(n)

#初期値
y1[0] = 3
y2[0] = 10
y3[0] = 45

for i in range(n):
    if i < n-1:
        y1[i+1] = y1[i] + h*f1(y1[i],y2[i])
        y2[i+1] = y2[i] + h*f2(y1[i],y2[i],y3[i])
        y3[i+1] = y3[i] + h*f3(y1[i],y2[i],y3[i])
ax.plot(y1,y2,y3)
plt.show()

ローレンツアトラクタをルンゲクッタ法とオイラー法でそれぞれ実装してみました.」への2件のコメント

  1. kmさんありがとうございます.確かにそうですね.大変失礼いたしました.

    新しいコードは次のようになります.また刻み幅を変えると滑らかに描画できました,

    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)

  2. ルンゲクッタ法の部分で
    > l3 = h * f2(s1,s3,s3)
    がs3ダブっているのと
    > s1 = y1[i] + k3/2
    > s2 = y2[i] + l3/2
    > s3 = y3[i] + m3/2
    が1/2で割ってしまっているところが間違っていますね。

コメントを残す

メールアドレスが公開されることはありません。

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください