アップデート 2001.3.8

ダフィングの方程式のマクロ

---

Sub duffing()
' 二階微分方程式 Runge-Kutta Method
' Duffingの方程式
' x'=y
' y'=-a*y-x^3-b*cos(t)
[B10].Select 'データを書き出す最初のセルの場所
Steps = Range("steps") 'stepsと名前の付いたセルからデータを取り込む
Period = Range("period") 'periodと名前の付いたセルからデータを取り込む
Iteration = Range("iteration") 'iterationと名前の付いたセルからデータを取り込む
x = Range("x0") 'x0と名前の付いたセルからデータを取り込む
y = Range("y0") 'y0と名前の付いたセルからデータを取り込む
a = Range("a") 'aと名前の付いたセルからデータを取り込む
b = Range("b") 'bと名前の付いたセルからデータを取り込む
h = Period / Steps '刻み幅の計算
For i = 1 To Iteration
t = 0 '時刻の初期値
For j = 1 To Steps '
Call rk2(h, t, x, y, a, b) 'ルンゲ･クッタ法のサブルーチンを呼び出す
Next j
ActiveCell.Offset(i - 1, 4).Value = x '変位xを書き出す
ActiveCell.Offset(i - 1, 5).Value = y '速度yを書き出す
Next i
　
kaisu =2
t = 0 '時刻の初期値
For j = 1 To steps * kaisu '
ActiveCell.Offset(j - 1, 0).Value = t '時刻tを書き出す
ActiveCell.Offset(j - 1, 1).Value = x '変位xを書き出す
ActiveCell.Offset(j - 1, 2).Value = y '速度yを書き出す
Call rk2(h, t, x, y, a, b) 'ルンゲ･クッタ法のサブルーチンを呼び出す
Next j
End Sub
Sub rk2(h, t, x, y, a, b) 'ルンゲ･クッタ法のサブルーチン
k1 = h * f(t, x, y) '関数fの値の計算
l1 = h * g(t, x, y, a, b) '関数gの値の計算
t1 = t + h / 2
x1 = x + k1 / 2
y1 = y + l1 / 2
k2 = h * f(t1, x1, y1)
l2 = h * g(t1, x1, y1, a, b)
x2 = x + k2 / 2
y2 = y + l2 / 2
k3 = h * f(t1, x2, y2)
l3 = h * g(t1, x2, y2, a, b)
t = t + h
x3 = x + k3
y3 = y + l3
k4 = h * f(t, x3, y3)
l4 = h * g(t, x3, y3, a, b)
x = x + (k1 + 2 * (k2 + k3) + k4) / 6
y = y + (l1 + 2 * (l2 + l3) + l4) / 6
End Sub
Function f(t, x, y) '関数x'=f
f = y
End Function
Function g(t, x, y, a, b) '関数y'=g
g = -a * y - x * x * x + b * Cos(t)
End Function

---