數值常微分方程(IVP)

Predictor-Corrector Method

Adams-Bashforth-Moulton Method

我們知道 RK4 是一個廣泛被使用的方法,這個方法的好處是 truncation error 的階數較高,我們預期他會產生的誤差比較小,但縱使 RK4 是一個可能不錯的方法,但我們也擔心,誤差會隨著計算次數的增加,一直把誤差堆疊到後面的 \(y_n\) 項。我們也知道,一路走來我們的分析工具都是使用泰勒展開式,如果我們有

\(\{ (t_{k-3},f_{k-3}), (t_{k-2},f_{k-2}), (t_{k-1},f_{k-1}), (t_{k},f_{k}) \}\)

我們也可以利用這些值來製作出一個三次的多項式來逼近函數 \(f(t)\),用這個多項式\(l_3(t)\) 代替未來要計算的 \([t_k,t_{k+1}]\) 的 \(f(t)\),我們有 \(y'(t) = f(t,y) \simeq l_3(t)\),得到

\(p_{k+1} = y_k + \int_{tk}^{t_{k+1}} l_3(t) dt = y_k + \frac{h}{24}(-9f_{k-3}+37f_{k-2}-59f_{k-1}+55f_{k})\)

當我們有了新的 \(y_{k+1}\) 值,我們重新利用

\(\{ (t_{k-2},f_{k-2}), (t_{k-1},f_{k-1}), (t_{k},f_{k}), , (t_{k+1},f(t_{k+1},p_{k+1})) \}\)

這四個點重算一次逼近函數 \(f(t)\) 的多項式 \(l'_3(t)\),得到 Corrector 來取代 \(y_{k+1}\)

\(y_{k+1} = y_k + \int_{tk}^{t_{k+1}} l_3(t) dt = y_k + \frac{h}{24}(f_{k-2}-5f_{k-1}+19f_{k}+9f_{k+1})\)

要分析這 Predictor 和 Corrector 逼近的階數,我們一樣用老招泰勒展開式

\(y_{k+1} = y_k + hf_k + \frac{h^2}{2}\left( (\frac{-1}{3}f_{k-3}+\frac{3}{2}f_{k-2}-3f_{k-1}+\frac{11}{6}f_k)/h +\frac{h^3}{4}f_{k}^{(4)}+...\right)\\

+\frac{h^3}{3!} \left((-f_{k-3}+4f_{k-2}-5f_{k-1}+2f_{k})/h^2 +\frac{11}{12}h^2 f_{k}^{(4)} + ...\right)\\

+\frac{h^4}{4!} \left((-f_{k-3}+3f_{k-2}-3f_{k-1}+f_{k})/h^3 +\frac{3}{2}hf_{k}^{(4)}+...\right) + \frac{h^5}{120} f_k^{(4)}+...\\

= y_k + \frac{h}{24}( -9f_{k-3} + 37f_{k-2} -59 f_{k-1} +55 f_{k})) + \frac{251}{720} h^5 f_{k}^{(4)} + ...\\

\simeq p_{k+1} + \frac{251}{720}h^{5} f_{k}^{(4)}(\xi)\)

(ABM 式一)

上式的第一行 \(h^2\) 項 裡的 \(f_{k-3}\) 到 \(f_k\) 的對應係數 \(\{-1/3, 3/2, -3, 11/6\}\) 是讓 \((-f_{k-3}+4f_{k-2}-5f_{k-1}+2f_{k})/h^2\) 成為 \(f'_k\) 的逼近,並且truncation error 是 \(h^3\) 。同理, \((-f_{k-3}+4f_{k-2}-5f_{k-1}+2f_{k})/h^2\) 是 \(f_k^{(2)}\) 的逼近,truncation error 是 \(h^2\); \((-f_{k-3}+3f_{k-2}-3f_{k-1}+f_{k})/h^3\) 是 \(f_k^{(3)}\) 的逼近,truncation error 是 \(h^1\)。如此,整理後的係數,就剛好是Predictor 的係數,這說明 \(p_{k+1}\) 的global 誤差是四階的。類似的證法我們得到

\(y_{k+1} = y_k + hf_{k+1} - \frac{h^2}{2} \left( (\frac{-1}{3}f_{k-2}+\frac{3}{2}f_{k-1}-3f_{k}+\frac{11}{6}f_{k+1})/h +\frac{h^3}{4}f_{k+1}^{(4)}+...\right)\\

+\frac{h^3}{3!} \left( (-f_{k-2}+4f_{k-1}-5f_{k}+2f_{k+1})/h^2 +\frac{11}{12}h^2 f_{k+1}^{(4)} + ...\right)\\

+\frac{h^4}{4!} \left( (-f_{k-2}+3f_{k-1}-3f_{k}+f_{k+1})/h^{3} +\frac{3}{2}hf_{k+1}^{(4)}+...) + \frac{h^5}{120} f_{k+1}^{(4)}+...\right)\\

= y_{k} + \frac{h}{24}( f_{k-2} - 5f_{k-1} + 19 f_{k} + 9 f_{k+1}) - \frac{19}{720} h^{5} f_{k+1}^{(4)} + ...\\

\simeq c_{k+1} + \frac{19}{720}h^{5} f_{k}^{(4)}(\xi)\)

(ABM式二)

這裡顯示 Corrector 也和 Predictor 一樣有四階的 global 誤差。所以,ABM Method 應該和 RK4 有接近的數值結果。讀者可能會有這樣的疑問,既然都是四階的誤差,那麼我們何必大費周章地介紹ABM 方法?回顧一下在數值積分與數值微分時,truncation error 的階數與我們真正在計算時的誤差控制是如何的關聯?在數值積分與數值微分計算時,雖然我們不知道truncation error 項的帶頭係數,但是我們可以利用不同的step size 來估計這一個帶頭係數,進一步我們可以知道數值解與真解的誤差。然而在解IVP 類型的ODE 時,如果我們把step size 做改變,例如變為先前的一半,然後呢?我們有什麼?第一次以 \(h\) 計算的 \(y_1\) 要和第二次以 \(h/2\) 計算的\(y_2\) 做比較,而這兩個對應的泰勒展開式會多一個未知的第二次以\(h/2\) 計算\(y_1\) 中的truncation error 項的係數。所以,我們還是無法估計到底在 \(y(t+h)\) 這個位置的 truncation error 帶頭係數是多少?但是,在ABM 方法中,我們卻可以在沒有真解的情況下估計這個誤差項。

帶有修正程序的ABM Method

從(ABM 式一)與(ABM 式二)中我們得到

\(E_{p,k+1} = y_{k+1} - p_{k+1} \simeq \frac{251}{720}h^5f_k^{(4)} \simeq \frac{251}{720}Kh^5\\

E_{c,k+1} = y_{k+1} - c_{k+1} \simeq \frac{-19}{720}h^5f_k^{(4)} \simeq \frac{-19}{720}Kh^5\)

這裡我們欠缺的是 \(K\) 到底是多少,這樣我們才知道我們計算的 \(p_{k+1}\) 或是\(c_{k+1}\) 到底有多靠近真解。在沒有 \(y_{k+1}\) 的情況下,我們有

\(E_{p,k+1}-E_{c,k+1} = c_{k+1}-p_{k+1} \simeq \frac{270}{720} Kh^5 = \frac{270}{251}E_{p,k+1} = \frac{-270}{19} E_{c,k+1}\)

這裡有一件巧妙的事,就是 \(c_{k+1}\) 與 \(p_{k+1}\) 是我們知道的,\(h\) 也是有的,所以通過上述的式子我們可以估計 \(K\) 的值,然後就知道對應的 \(E_{p,k+1}\) 與 \(E_{c,k+1}\) 是多少?

再來,既然我們知道真解 \(y_{k+1}\) 與 Predictor \(p_{k+1}\) 及 Corrector \(c_{k+1}\) 的誤差可能是多少?自然地,我們把這可能的誤差再修正一次,所得到的解,就應該有很大的機會比 RK4 表現優良。因此,修正後的 ABM 方法為

- Predictor: \(p_{k+1} = y_k + \frac{h}{24}( -9f_{k-3} + 37f_{k-2} -59 f_{k-1} +55 f_{k})\)
- Modifier: \(m_{k+1} = p_{k+1} + \frac{251}{270}(c_k-p_k)\)
- Corrector: \(c_{k+1} = y_k + \frac{h}{24}( f_{k-2} - 5f_{k-1} + 19 f_k + 9 f(t_{k+1},m_{k+1}))\)
- Modification: \(y_{k+1} = c_{k+1}-\frac{19}{270}(c_{k+1}-p_{k+1})\)


帶有修正的 ABM 方法的Python 程式如下 :

def ode_ABM(f,t0,tf,y0=0,n=100):
    """
    First order ODE (y' = f(t,y)) Solver using ABM method with modification
    
    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)
    tfrk4 = t[3]
    y = ode_RK4(f,t0,tfrk4,y0,4)
    y = list(y)
    p = y[-1]
    c = y[-1]
    for i in range(3,n-1):
        h = t[i+1]-t[i]
        p1 = y[-1] + (h/24)*(-9*f(t[i-3],y[-4])+37*f(t[i-2],y[-3])-59*f(t[i-1],y[-2])+55*f(t[i],y[-1]))
        m = p1 + (251/270.0)*(c-p)
        c1 = y[-1] + (h/24)*(f(t[i-2],y[-3])-5*f(t[i-1],y[-2])+19*f(t[i],y[-1])+9*f(t[i+1],m)) 
        y.append(c1-(19/270.0)*(c1-p1))
        c, p = c1, p1
        
    y = np.array(y)

 

仿造先前的範例,我們比較RK4 和ABM 的誤差,程式如下:

>>> 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_RK4 = ode_RK4(lambda t,y:dfun(t,y),t0,tf,y0,N)
>>> y_ABM = ode_ABM(lambda t,y:dfun(t,y),t0,tf,y0,N)
>>> from pylab import *
>>> fig1 = figure(1)
>>> plot(t,y_RK4-yt,'r')
[<matplotlib.lines.Line2D object at 0x1c48c50>]
>>> plot(t,y_ABM-yt,'b')
[<matplotlib.lines.Line2D object at 0x1c48c10>]
>>> fig1.savefig('RK4_ABM_error.png')

我們從圖形的結果發現,ABM 真的又比RK4 優秀一些。