アップデート 2003.9.6

単振り子のマクロ

---

Sub pendulum()
' 二階微分方程式 Runge-Kutta Method
' 単振り子
' x'=y
' y'=-omega*omega*sin(x)
[B10].Select 'データを書き出す最初のセルの名前
steps = Range("steps") 'stepsと名前の付いたセルからデータを取り込む
t_end = Range("t_end") 'stendと名前の付いたセルからデータを取り込む
x = Range("x0") 'x0と名前の付いたセルからデータを取り込む
y = Range("y0") 'y0と名前の付いたセルからデータを取り込む
omega = Range("omega") 'omegaと名前の付いたセルからデータを取り込む
w2 = omega ^ 2 ' 定数値の計算
h = t_end / steps '刻み幅の計算
t = 0 '時刻の初期値
For i = 1 To steps '
ActiveCell.Offset(i - 1, 0).Value = t '時刻tを書き出す
ActiveCell.Offset(i - 1, 1).Value = x '変位xを書き出す
ActiveCell.Offset(i - 1, 2).Value = y '速度yを書き出す
Call rk2(h, t, x, y, w2) 'ルンゲ・クッタ法のサブルーチンを呼び出す
Next
End Sub
Sub rk2(h, t, x, y, w2) 'ルンゲ・クッタ法のサブルーチン
k1 = h * f(t, x, y) '関数fの値の計算
l1 = h * g(t, x, y, w2) '関数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, w2)
x2 = x + k2 / 2
y2 = y + l2 / 2
k3 = h * f(t1, x2, y2)
l3 = h * g(t1, x2, y2, w2)
t = t + h
x3 = x + k3
y3 = y + l3
k4 = h * f(t, x3, y3)
l4 = h * g(t, x3, y3, w2)
x = x + (k1 + 2 * (k2 + k3) + k4) / 6
y = y + (l1 + 2 * (l2 + l3) + l4) / 6
End Sub
Function f(t, x, y) '関数f
f = y
End Function
Function g(t, x, y, w2) '関数g
g = -w2 * Sin(x)
End Function