數值常微分方程(IVP)

Runge-Kutta Method

Runge-Kutta method 經常被簡稱為RK2 與RK4 分別代表兩種不同收斂速度的數值積分公式。RK2 的公式如下:
\(k_1 = hf(t_n,y_n)\)

\(k_2 = hf(t_n+h/2,y_n+\frac{k_1}{2})\)

\(y_{n+1} = y_n + k_2 +O(h^3) \)

RK2 的2 代表這個方法是二階(truncation error 是 \(O(h^3)\))的誤差,證明的方法如下:

首先,我們從最原始的題目 \(y'(t) = f(t,y)\) 開始,得到

\(y^{(2)}(t) = f_t(t,y) + f_y(t,y)y'(t)\)

將原始題目代入後得到

\(y^{(2)}(t) = f_t(t,y) + f_y(t,y)f(t,y)\)

有了這個式子和原方程式,\(y(t+h)\) 對 \(y(t)\) 的泰勒展開式可以寫成

(RK式一):\(y(t+h) = y(t) + hf(t,y) + \frac{1}{2}h^2 f_t(t,y) + \frac{1}{2}h^2 f_y(t,y) f(t,y) + O(h^3)\)

我們對應 RK2 的形式把 \(y(t+h)\) 的疊代式寫成

(RK式二):\(y(t+h) = y(t) + A hf_0 + B hf_1\)

的形式,其中 \(A\), \(B\) 是未知的常數,\(f_0 = f(t,y)\),\(f_1 = f(t+P h, y+Q hf_0)\),其中 \(P\), \(Q\) 是另外兩個要決定的常數。

我們把 \(f_1\) 對 \(f(t,y)\) 的泰勒展式寫開,得到

\(f_1 = f(t,y) + P h f_t(t,y) + Q h f_y(t,y)f(t,y) + O(h^2)\)

我們把上式以及 \(f_0\) 代入(RK式二)中,經過整理得到

(RK式三):\(y(t+h) = y(t) + (A+B) h f(t,y) + B P h^2f_t(t,y) + B Q h^2f_y(t,y)f(t,y) + O(h^3) \)

比較(RK式三)與(RK式一)我們得到

- \(A+B=1\)
- \(BP = \frac{1}{2}\)
- \(BQ = \frac{1}{2}\)


四個未知數,三個方程式,我們可以有很多種 \((A,B,P,Q)\) 的組合滿足這三個方程式。為了讓 RK2 的公式變得簡潔,我們選擇讓 \(A = 0\),因此得到 \(B = 1\),\(P = \frac{1}{2}\) 以及 \(Q = \frac{1}{2}\),這就是 RK2 的由來。

更高階的 Runge-Kutta Method 是 RK4,意思是 truncation error 是\(O(h^5)\)。他的公式如下:

\(k_1 = hf(t_n,y_n)\)

\(k_2 = hf(t_n+h/2,y_n+K_1/2)\)

\(k_3 = hf(t_n+h/2,y_n+k_2/2)\)

\(k_4 = hf(t_n+h,y_n+k_3)\)

\(y_{n+1} = y_n + \frac{k_1 + 2k_2 + 2k_3 + k_4}{6} + O(h^5)\)

要證明這個方法的 truncation error 是 \(O(h^5)\),其證明方式與 RK2 類似,我們令

\(y_{n+1} = y_n + w_1k_1 + w_2k_2 + w_3k_3 + w_4k_4\)

並且
\(k_1 = hf(t_n,y_n)\)

\(k_2 = hf(t_n+a_1 h,y_n+b_1k_1)\)

\(k_3 = hf(t_n+a_2 h,y_n+b_2 k_1 + b_3 k_2)\)

\(k_4 = hf(t_n+a_3 h,y_n + b_4 k_1 +b_5 k_2+ b_6 k_3)\)

同樣地,我們把每一項泰勒展開式都寫到 \(h^5\) 項,整理後就可以得知為何 RK4 會有 \(\frac{1}{6}(1,2,2,1)\) 以及其他的對應係數了。 RK2 與 RK4 的Python 程式如下:

def ode_RK2(f,t0,tf,y0=0,n=100):
    """
    First order ODE (y' = f(t,y)) Solver using RK2 method
    
    t0: initial value of independent variable
    tf: final value of independent variable
    y0: initial value of dependent variable
    n : number of steps
    f : function of f(t,y)

    The return is list of y(t)
    """

    t = np.linspace(t0,tf,n)
    y = list([y0])
    for i in range(n-1):
        h = t[i+1]-t[i]
        k1 = h*f(t[i],y[-1])
        k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
        y.append(y[-1]+k2)

    y = np.array(y)
    
    return y

def ode_RK4(f,t0,tf,y0=0,n=100):
    """
    First order ODE (y' = f(t,y)) Solver using RK4 method
    
    t0: initial value of independent variable
    tf: final value of independent variable
    y0: initial value of dependent variable
    n : number of steps
    f : function of f(t,y)

    The return is list of y(t)
    """

    t = np.linspace(t0,tf,n)
    y = list([y0])
    for i in range(n-1):
        h = t[i+1]-t[i]
        k1 = h*f(t[i],y[-1])
        k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
        k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0)
        k4 = h*f(t[i]+h,y[-1]+k3)
        y.append(y[-1]+(k1+2*k2+2*k3+k4)/6)

    y = np.array(y)
    
    return y

 

學了上述的數值ODE方法,對於IVP 我們可以比較一下這些方法的差異。我們給一個知道真解的例子,來比較每一個數值方法會造成多大的誤差。

>>> t0,tf,y0,N = 0,10,0,50
>>> def dfun(t,y):
.... return -y+1

>>> def tfun(t):
.... return 1-np.exp(-t)

>>> t = np.linspace(t0,tf,N)
>>> yt = tfun(t)
>>> y_euler = ode_Euler(lambda t,y:dfun(t,y),t0,tf,y0,N)
>>> y_heun = ode_Heun(lambda t,y:dfun(t,y),t0,tf,y0,N)
>>> y_RK2 =ode_RK2(lambda t,y:dfun(t,y),t0,tf,y0,N)
>>> y_RK4 = ode_RK4(lambda t,y:dfun(t,y),t0,tf,y0,N)
>>> from pylab import *
>>> plot(t,yt,'r')
[<matplotlib.lines.Line2D object at 0x1c48c10>]
>>> plot(t,y_euler,'b')
[<matplotlib.lines.Line2D object at 0x5c74cf0>]
>>> plot(t,y_heun,'g')
[<matplotlib.lines.Line2D object at 0x5c74c70>]
>>> plot(t,y_RK2,'r-.')
[<matplotlib.lines.Line2D object at 0x5c74c10>]
>>> plot(t,y_RK4,'k')
>>> fig1 = figure(1)
>>> fig1.show()
>>> fig1.savefig('ode_compare.png')
>>> plot(t,y_heun-yt,'r')
[<matplotlib.lines.Line2D object at 0x7f963d0>]
>>> plot(t,y_RK2-yt,'b')
[<matplotlib.lines.Line2D object at 0x7f96590>]
>>> plot(t,y_RK4-yt,'k')
>>> fig2 = figure(2)
>>> fig2.show()
>>> fig2.savefig('ode_err_compare.png')

程式中的 dfun() 是定義 \(y'(t) = f(t,y)\) 的 \(f(t,y)\),tfun() 是真解,其餘代號都顯而易見。我們將這些數值解和真解畫在同一張圖上,我們可以看見只有Euler Method 特別例外,可見 Euler Method 只能當作入門的教材,在真實計算上並不適用。其餘的數值解幾乎和真解貼平無法分辨。我們惕除 Euler 的解,單畫誤差,我們看見同樣是二階的 Huen's Method 和 RK2 貼齊了,另外的 RK4 因為太小所以看不見。這裡要說得是,RK4 是被廣泛使用的,在不太要求精準度的情況下,RK4 幾乎是多數的數學軟體預設的解 IVP 方法。