アップデート 2009.7.6

ローレンツの方程式のマクロ

Sub Lorenz() ' 三元一階微分方程式
'Runge-Kutta Method' Lorenz方程式
'x'=a(y-x)
'y'=-bx-y-xz
'z'=-xy-cz
[B10].Select 'データを書き出す最初のセルの場所
steps = Range("steps") 'stepsと名前の付いたセルからデータを取り込む,
Iteration = Range("iteration") 'iterationと名前の付いたセルからデータを取り込む
a = Range("ca") 'caと名前の付いたセルからデータを取り込む
b = Range("cb") 'cbと名前の付いたセルからデータを取り込む
c = Range("cc") 'ccと名前の付いたセルからデータを取り込む
x = Range("x0") 'x0と名前の付いたセルからデータを取り込む
y = Range("y0") 'y0と名前の付いたセルからデータを取り込む
z = Range("z0") 'z0と名前の付いたセルからデータを取り込む
stwps = Range("steps") 'stepsと名前の付いたセルからデータを取り込む
t = 0 '時刻の初期値
For i = 1 To Iteration '繰り返し
Call rk2(steps, t, x, y, z, a, b, c) 'ルンゲ･クッタ法のサブルーチンを呼び出す
ActiveCell.Offset(i - 1, 3).Value = t '時刻tを書き出す
ActiveCell.Offset(i - 1, 4).Value = x '変位xを書き出す
ActiveCell.Offset(i - 1, 5).Value = y '変位yを書き出す
ActiveCell.Offset(i - 1, 6).Value = z '変位zを書き出す
Next i
End Sub

Sub rk2(h, t, x, y, z, a, b, c) 'ルンゲ･クッタ法のサブルーチン
k1 = h * f1(t, x, y, z, a, b, c) '関数fの値の計算
l1 = h * f2(t, x, y, z, a, b, c) '関数gの値の計算
m1 = h * f3(t, x, y, z, a, b, c) '関数hの値の計算
t1 = t + h / 2
x1 = x + k1 / 2
y1 = y + l1 / 2
z1 = z + m1 / 2
k2 = h * f1(t1, x1, y1, z1, a, b, c)
l2 = h * f2(t1, x1, y1, z1, a, b, c)
m2 = h * f3(t1, x1, y1, z1, a, b, c)
x2 = x + k2 / 2
y2 = y + l2 / 2
z2 = z + m2 / 2
k3 = h * f1(t1, x2, y2, z2, a, b, c)
l3 = h * f2(t1, x2, y2, z2, a, b, c)
m3 = h * f3(t1, x2, y2, z2, a, b, c)
t = t + h
x3 = x + k3
y3 = y + l3
z3 = z + m3
k4 = h * f1(t, x3, y3, z3, a, b, c)
l4 = h * f2(t, x3, y3, z3, a, b, c)
m4 = h * f3(t, x3, y3, z3, a, b, c)
x = x + (k1 + 2 * (k2 + k3) + k4) / 6
y = y + (l1 + 2 * (l2 + l3) + l4) / 6
z = z + (m1 + 2 * (m2 + m3) + m4) / 6
End Sub

Function f1(t, x, y, z, a, b, c) '関数x'=a(y-z)
f1 = a * (y - x)
End Function

Function f2(t, x, y, z, a, b, c) '関数=bx-y-xz
f2 = b * x - y - x * z
End Function

Function f3(t, x, y, z, a, b, c) '関数z'=xy-cz
f3 = x * y - c * z
End Function