常微分方程數值計算方法

03-06

Runge-Kutta-Fehlbrg Method

該方法的出現是為了解決 $Runge-Kutta-Method$ 的缺陷。因為每次計算都要比較 ${h}/{2}和{h}$ 兩次步長產生的結果差。如果結果不好，就要縮小步長，這又要重新開始兩次計算。所以 $Runge-Kutta-Fehlbrg-Method$ 出現了，又名RKF45方法。

方法思路：

1. 每一次計算，都使用連個不同的計算方法求解，得出兩個值

2. 如果兩次求解結果接近，則接受求解結果。
3. 如果兩次求解結果的誤差超過了指定的容許值，則減小步長。如果誤差超過了有效位數，則增加步長。

方法步驟：

1. 每一步要使用下面六個值： $left.egin{aligned} k1 &= hf(t_{k}, y_{k})\ k2 &= hf(t_{k}+dfrac{1}{4}h, y_{k}+dfrac{1}{4}k_{1})\ k3 &= hf(t_{k}+dfrac{3}{8}h, y_{k}+dfrac{3}{32}k_{1}+dfrac{9}{32}k_{2})\ k4 &= hf(t_{k}+dfrac{12}{13}h, y_{k}+dfrac{1932}{2197}k_{1}-dfrac{7200}{2197}k_{2}+dfrac{7296}{2197}k_{3})\ k5 &= hf(t_{k}+h, y_{k}+dfrac{439}{216}k_{1}-8k_{2}+dfrac{3680}{513}k_{3}-dfrac{845}{4104}k_{4})\ k6 &= hf(t_{k}+dfrac{1}{2}h, y_{k}-dfrac{8}{27}k_{1}+2k_{2}-dfrac{3544}{2565}k_{3}+dfrac{1859}{4104}k_{5}-dfrac{11}{40}k_{4})\ end{aligned} ight}$
2. 用4階Runge-Kutta Method求出一個近似解： $y_{k+1}=y_{k}+dfrac{25}{216}k_{1}+dfrac{1408}{2565}k_{3}+dfrac{2197}{4101}k_{4}-dfrac{1}{5}k_{5}$

3. 用5階Runge-Kutta Method求出一個更好的近似解：
$z_{k+1}=y_{k}+dfrac{16}{135}k_{1}+dfrac{6656}{12825}k_{3}+dfrac{28561}{56430}k_{4}-dfrac{9}{50}k_{5}+dfrac{2}{55}k_{6}$
4. 最佳步長 $sh$ 可以通過當前步長值乘以標量 $s$ 得到，標量 $s$ 為: $=Big(dfrac{ m tol h}{2|z_{k+1}-y_{k+1}}Big)^{1/4}approx 0.84Big(dfrac{ m tol h}{2|z_{k+1}-y_{k+1}}Big)^{1/4}$

計算程序(待修改):

1. Runge-Kutta四階經典演算法。

def f(y): return 1 + y ** 2def rk4(a, b, y_init,M): h = (b - a)/M y = y_init for i in range(M): k1 = f(y) k2 = f(y + h*k1/2) k3 = f(y + h*k2/2) k4 = f(y + h*k3) y = y + h*(k1 + 2*k2 + 2*k3 + k4)/6 print("iteration time {1:2d}: y is {0:.6f}".format(y, i+1))c = rk4(0, 1.4, 0, 14)列印結果:iteration time 1: y is 0.100335iteration time 2: y is 0.202710iteration time 3: y is 0.309336iteration time 4: y is 0.422793iteration time 5: y is 0.546302iteration time 6: y is 0.684137iteration time 7: y is 0.842289iteration time 8: y is 1.029639iteration time 9: y is 1.260159iteration time 10: y is 1.557406iteration time 11: y is 1.964747iteration time 12: y is 2.572072iteration time 13: y is 3.601563iteration time 14: y is 5.791975

計算結果與參考文獻 $^{[1]}$ P347表9.11一樣。改程序實際上應該調整的，沒有自變數 $m{t}$ ，可以設想自變數 $m{t} = 1$ 來考慮。

2. Runge Kutta Fehlbrg方法

def fsolve(t, y): return 1 + y ** 2def rkfMethod(a, b, y_init,M): h = (b - a)/M y = y_init t = a #-------k coefficient---------- a1 = 1; b1 = 1 a2 = 1/4; b2 = 1/4 a3 = 3/8; b3 = 3/32; c3 = 9/32 a4 = 12/13; b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197 a5 = 1; b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104 a6 = 1/2; b6 = -8/27; c6 = 2 ; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40 #--------y coefficient--------- r1 = 16/135; r3 = 6656/12825; r4 = 28561/56430; r5 = -9/50; r6 = 2/55 for i in range(M): k1 = fsolve(t, y) k2 = fsolve(t + a2*h, y + b2*k1*h) k3 = fsolve(t + a3*h, y + b3*k1*h + c3*k2*h) k4 = fsolve(t + a4*h, y + b4*k1*h + c4*k2*h + d4*k3*h) k5 = fsolve(t + a5*h, y + b5*k1*h + c5*k2*h + d5*k3*h + e5*k4*h) k6 = fsolve(t + a6*h, y + b6*k1*h + c6*k2*h + d6*k3*h + e6*k4*h + f6*k5*h) #------------------------------ y = y + h*(r1*k1 + r3*k3 + r4*k4 + r5*k5 + r6*k6) print("iteration time {1:2d}: y is {0:.6f}".format(y, i+1)) return y列印結果：rkfMethod(0, 1.4, 0, 14)iteration time 1: y is 0.100335iteration time 2: y is 0.202710iteration time 3: y is 0.309336iteration time 4: y is 0.422793iteration time 5: y is 0.546303iteration time 6: y is 0.684137iteration time 7: y is 0.842288iteration time 8: y is 1.029639iteration time 9: y is 1.260159iteration time 10: y is 1.557409iteration time 11: y is 1.964762iteration time 12: y is 2.572157iteration time 13: y is 3.602127iteration time 14: y is 5.798128

參考文獻：

John.H.Mathews. 數值方法(Matlab)第四版.北京：機械工業出版社. 2017.02
Ramin S. Esfandiari. Numerical Methods for Engineers and Scientists Using MATLAB. CRC Press. 2017