エクセルを用いたルンゲクッタフェールベルグ法(Runge-Kutta-Fehlberg method：RKF45)

エクセルを用いたルンゲクッタフェールベルグ法(Runge-Kutta-Fehlberg method：RKF45)

Excelを用いた科学技術計算が第2版になりました　30年10月！

amazon：Excelで操る! ここまでできる科学技術計算第2版

楽天： Excelで操る! ここまでできる科学技術計算 [ 神足　史人 ]

　　　　　　　　　図１　ルンゲクッタフェールベルグ法で解いた例題のグラフ

　常微分方程式のポピュラーな解法は４段４次のルンゲクッタ法である。ここでは、より高性能である６段５次のルンゲクッタフェールベルグを説明する。Runge-Kutta-Fehlberg法は誤差をチェックし、時間刻み幅を変更しながら計算するため、４次のルンゲクッタ法と比較して非常に精度が良い方法である。

ルンゲクッタフェ－ルベルグ法の説明

　次の１階1元微分方程式で説明する。

　dy/dt = F( t , y)

　Fは時間tと変数yの関数である。

初期値　t=t0　の時 y=y0 とすると、Δt後の6個の変数値は下式となる。

　　k1 = Δt・F( t0, y0)

　　k2 = Δt・F( t0+Δt/4, y0+k1/4)

　　k3 = Δt・F( t0+3Δt/8, y0+3k1/32+9k2/32)

　　k4 = Δt・F( t0+12Δt/13, y0+1932k1/2197-7200k2/2197+7296k3/2197)

　　k5 = Δt・F( t0+Δt, y0+439k1/216-8k2+3680k3/513-845k4/4104)

　　k6 = Δt・F( t0+Δt/2, y0-8k1/27+2k2-3544k3/2565+1859k4/4104-11k5/40)

Δt後の変数値y1は4次の近似とし、次式となる。

　　y1 = y0 + 25k1/216 + 1408k3/2565 + 2197k4/4104 - k5/5

また、5次の近似は次式となる。

　　y1' = y0 + 16k1/135 + 6656k3/12825 + 28561k4/56430 - 9k5/50 + 2k6/55

５次の近似式と４次の近似式の差を誤差(R)とする。

　　R = |y1' - y1| = |k1/360 - 128k3/4275 - 2197k4/75240 + k5/50 + 2k6/55|

　誤差の判定は次式でｓをもとめ、s>=１の場合、近似値が許容誤差内と判定し次の時間ステップに進む。そうでない場合、Δt=sΔtに変更し計算をやりなおす。

　　s = [ εΔt/(2R)]^1/4

　　　　　　ε：許容誤差=0.00001

　例題のVBAプログラムでは、s>0.1の制限を設定している。