數值常微分方程(IVP)

 

Heun's Method

Heun's Method 又稱為 Trapezoidal Method,作法是將 \(y'(t) = f(t,y)\) 兩邊積分,變成

\(y(t)|_{t_k}^{t_{k+1}} = y(t_{k+1})-y(t_{k}) = \int_{t_k}^{t_{k+1}} f(t,y) dt\)

因此,我們有下列的遞迴式

\(y(t_{k+1}) =y(t_k) + \int_{t_k}^{t_{k+1}} f(t,y) dt\) , 其中 \(y(t_0) = y_0\)

若我們用常數函數 \(f(t_k,y(t_k))\) 來逼近這個小區間的函數值,那麼這個方法就退化成 Euler's method。若我們用 trapezoidal rule 來計算這個小區間的積分值,則遞迴式變成

\(y(t_{k+1}) =y(t_k) + \frac{h}{2}(f(t_k,y_k)+f(t_{k+1},y_{k+1}))\)
然而上式的右邊有一個 \(y_{k+1}\) 還是未知,因此我們用 \(y_k + h f(t_k,y_k)\) 來取代 \(y_{k+1}\),因此新的遞迴式成為

\(y(t_{k+1}) =y(t_k) + \frac{h}{2}(f(t_k,y_k)+f(t_{k+1},y_{k}+hf(t_k,y_k)))\)

同樣地,我們若把 truncation error 寫入迭代式,我們知道 Heun's method 的誤差是 \(O(h^3)\),所以我們稱 Heun's Method 的 Global 誤差是二階。Python 的 ODE Heun's Method 程式如下:

def ode_Heun(f,t0,tf,y0=0,n=100):
    """
    First order ODE (y' = f(t,y)) Solver using Heun's 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]
        fk = f(t[i],y[-1])
        fk1 = f(t[i+1],y[-1]+h*fk)
        y.append(y[-1]+h*(fk+fk1)/2.0)

    y = np.array(y)
    
    return y