エクセルで描く

　Lorenz方程式を描く 　　

               dx/dt = -px+py
               dy/dt = -xz+rx-y
               dz/dt = xy-bz　　　　

　(例）p=10　r=30　b=2 (初期値：x=1　y=2　z=1)の場合


Sub Lorez()
Dim k0(2), k1(2), k2(2), k3(2)

dt = 0.01
tmax = 50

x = 1
y = 2
z = 1
p = 10
r = 30
b = 2

For t = 0 To tmax Step dt
    k0(0) = dt * f1(t, x, y, z, p)
    k0(1) = dt * f2(t, x, y, z, r)
    k0(2) = dt * f3(t, x, y, z, b)

    k1(0) = dt * f1(t + dt / 2, x + k0(0) / 2, y + k0(1) / 2, z + k0(2) / 2, p)
    k1(1) = dt * f2(t + dt / 2, x + k0(0) / 2, y + k0(1) / 2, z + k0(2) / 2, r)
    k1(2) = dt * f3(t + dt / 2, x + k0(0) / 2, y + k0(1) / 2, z + k0(2) / 2, b)

    k2(0) = dt * f1(t + dt / 2, x + k1(0) / 2, y + k1(1) / 2, z + k1(2) / 2, p)
    k2(1) = dt * f2(t + dt / 2, x + k1(0) / 2, y + k1(1) / 2, z + k1(2) / 2, r)
    k2(2) = dt * f3(t + dt / 2, x + k1(0) / 2, y + k1(1) / 2, z + k1(2) / 2, b)

    k3(0) = dt * f1(t + dt, x + k2(0), y + k2(1), z + k2(2), p)
    k3(1) = dt * f2(t + dt, x + k2(0), y + k2(1), z + k2(2), r)
    k3(2) = dt * f3(t + dt, x + k2(0), y + k2(1), z + k2(2), b)

    x = x + (k0(0) + 2 * k1(0) + 2 * k2(0) + k3(0)) / 6
    y = y + (k0(1) + 2 * k1(1) + 2 * k2(1) + k3(1)) / 6
    z = z + (k0(2) + 2 * k1(2) + 2 * k2(2) + k3(2)) / 6

    Worksheets(1).Cells(6 + t / dt, 3) = x
    Worksheets(1).Cells(6 + t / dt, 4) = y
    Worksheets(1).Cells(6 + t / dt, 5) = z
Next t

End Sub

Function f1(t, x, y, z, p)
f1 = -p * x + p * y
End Function

Function f2(t, x, y, z, r)
f2 = -x * z + r * x - y
End Function

Function f3(t, x, y, z, b)
f3 = x * y - b * z
End Function

３Ｄ描画する(GraphR ver2.32にて)

Math TOP