アップデート 2001.2.14

太陽と地球の運動のマクロ

---

Sub earch()
' 二階微分方程式 Runge-Kutta Method
' 太陽を回る地球の運動
' r"-a^2/r^3=-b/r^2
' θ'=a/r^2
[B13].Select 'データを書き出す最初のセルの場所
n = Range("n") 'nと名前の付いたセルから繰り返し回数データを取り込む
h = Range("h") 'hと名前の付いたセルから刻み幅データを取り込む
a = Range("a") 'aと名前の付いたセルから係数データを取り込む
b = Range("b") 'bと名前の付いたセルから係数データを取り込む
r1 = Range("initr1") 'initr1と名前の付いたセルからrの初期値データを取り込む
r2 = Range("initr2") 'initr2と名前の付いたセルからr'の初期値データを取り込む
r3 = Range("inits0") 'inits0と名前の付いたセルからθの初期値データを取り込む
t = 0 '時刻の初期値
For i = 1 To n + 1
ActiveCell.Offset(i, 0).Value = t '時刻tを書き出す
ActiveCell.Offset(i, 1).Value = r1 '距離の変位rを書き出す
ActiveCell.Offset(i, 2).Value = r2 '距離の速度r'を書き出す
ActiveCell.Offset(i, 3).Value = r3 '角θを書き出す
ActiveCell.Offset(i, 4).Value = r1 * Cos(r3) 'x座標を書き出す
ActiveCell.Offset(i, 5).Value = r1 * Sin(r3) 'y座標を書き出す
Call rk2(h, t, r1, r2, r3, a, b) 'ルンゲ・クッタ法のサブルーチンを呼び出す
Next
End Sub
Sub rk2(h, t, r1, r2, r3, a, b) 'ルンゲ・クッタ法のサブルーチン
k1 = h * f1(h, t, r1, r2, r3, a, b)
l1 = h * f2(h, t, r1, r2, r3, a, b)
m1 = h * f3(h, t, r1, r2, r3, a, b)
t1 = t + h / 2
s1 = r1 + k1 / 2
s2 = r2 + l1 / 2
s3 = r3 + m1 / 2
k2 = h * f1(h, t1, s1, s2, s3, a, b)
l2 = h * f2(h, t1, s1, s2, s3, a, b)
m2 = h * f3(h, t1, s1, s2, s3, a, b)
s1 = r1 + k2 / 2
s2 = r2 + l2 / 2
s3 = r3 + m2 / 2
k3 = h * f1(h, t1, s1, s2, s3, a, b)
l3 = h * f2(h, t1, s1, s2, s3, a, b)
m3 = h * f3(h, t1, s1, s2, s3, a, b)
t = t + h
s1 = r1 + k3
s2 = r2 + l3
s3 = r3 + m3
k4 = h * f1(h, t1, s1, s2, s3, a, b)
l4 = h * f2(h, t1, s1, s2, s3, a, b)
m4 = h * f3(h, t1, s1, s2, s3, a, b)
r1 = r1 + (k1 + 2 * (k2 + k3) + k4) / 6
r2 = r2 + (l1 + 2 * (l2 + l3) + l4) / 6
r3 = r3 + (m1 + 2 * (m2 + m3) + m4) / 6
End Sub
Function f1(h, t1, s1, s2, s3, a, b)
f1 = s2
End Function
Function f2(h, t1, s1, s2, s3, a, b)
f2 = a ^ 2 / s1 ^ 3 - b / s1 ^ 2
End Function
Function f3(h, t1, s1, s2, s3, a, b)
f3 = a / s1 ^ 2
End Function