単振り子とルンゲ・タック法

※VBAカテにするかどうか迷ったのですが、微分方程式の数値解の話ですのでこちらで質問させてください。

Excel VBA で単振り子の厳密解をルンゲ・タック法で解いてみたのですが、初期条件とした与えた角度（単振動で言えば振幅）が時間が経過すると共に大きくなっていきます。どこに問題があるのでしょうか？
　ここの画像は小さいので（いい加減改善してほしい）
http://imepic.jp/20220714/862800
でご確認願います。

Sub SiPendulum()
　　Worksheets("単振り子").Activate

　　Dim K1 As Double, K2 As Double, K3 As Double, K4 As Double
　　Dim L1 As Double, L2 As Double, L3 As Double, L4 As Double
　　Dim t As Double, x As Double, v As Double
　　Dim h As Double, Dx As Double, Dv As Double
　　Dim i As Long

　　
　　Cells(5, 5).Value = "t"
　　Cells(5, 6).Value = "θ"
　　Cells(5, 7).Value = "ω"
　 ' d^2θ/dt^2 = -(g/L)sin(θ) ------>　dω/dt = -(g/L)sin(θ).
　 '　　　　　　　　　　　　　　　　　　 dθ/dt = ω
　 ' 変数名を変える θ→x　　ω→v
　
　 ' d^2x/dt^2 = -(g/L)sin(x) ------>　dv/dt = -(g/L)sin(x).
　 '　　　　　　　　　　　　　　　　　 dx/dt = v

　　h = 0.05
　　t = 0
　　x = Cells(5, 11).Value '初期位相
　　v = 0
　　
　　Cells(6, 5).Value = t
　　Cells(6, 6).Value = x
　　Cells(6, 7).Value = v

　　For i = 1 To 81
　　　　K1 = h * Px_t(t, x, v) 'dx/dt
　　　　L1 = h * Pv_t(t, x, v) 'dv/dt
　　　　K2 = h * Px_t(t + h / 2, x + K1 / 2, v + L1 / 2)
　　　　L2 = h * Pv_t(t + h / 2, x + L1 / 2, v + L1 / 2)
　　　　K3 = h * Px_t(t + h / 2, x + K2 / 2, v + L1 / 2)
　　　　L3 = h * Pv_t(t + h / 2, x + L2 / 2, v + L1 / 2)
　　　　K4 = h * Px_t(t + h, x + K3, v + L3)
　　　　L4 = h * Pv_t(t + h, x + L3, v + L3)
　　　　Dx = (K1 + 2 * K2 + 2 * K3 + K4) / 6
　　　　Dv = (L1 + 2 * L2 + 2 * L3 + L4) / 6
　　　　t = t + h
　　　　x = x + Dx
　　　　v = v + Dv
　　　　Cells(6 + i, 5).Value = t
　　　　Cells(6 + i, 6).Value = x
　　　　Cells(6 + i, 7).Value = v
　　Next i
End Sub

'dx/dt = v
Function Px_t(t, x, v)
　　Px_t = v
End Function

'dv/dt = -(g/L)sin(x)
Function Pv_t(t, x, v)
　　Worksheets("単振り子").Activate
　　Dim GL As Double
　　GL = Cells(6, 4).Value
　　Pv_t = -GL * Sin(x)
End Function

No.1 ベストアンサー 回答者： Tacosan 回答日時： 2022/07/15 00:50

古典的な 4段 4次のルンゲ・クッタ法だと思うけど, L2, K3, L3, L4 の計算は間違ってない? 例えば, L2 の式

L2 = h * Pv_t(t + h / 2, x + L1 / 2, v + L1 / 2)
は x+L1/2 でいいの?

あと t の関数 x(t) について dx/dt が陽に現れない微分方程式
d^2x/dt^2 = f(t, x)
だと同じ 4段 4次でももうちょっと精度を上げることができたはず.

この回答へのお礼

> L2 = h * Pv_t(t + h / 2, x + L1 / 2, v + L1 / 2)
> は x+L1/2 でいいの?
　あちゃ～、お恥ずかしい（笑）。

> あと t の関数 x(t) について dx/dt が陽に現れない微分方程式
> d^2x/dt^2 = f(t, x)
> だと同じ 4段 4次でももうちょっと精度を上げることができたはず.
　知りませんでした。感謝いたします。

お礼日時：2022/07/15 02:49