數值常微分方程(IVP)

歐拉法(Euler Method)

Euler Method 不是一個好方法,它的精度也差,但幾乎提到數值微分方程,都一定會介紹到這個方法。原因是它很簡單,並且從介紹 Euler Method 中可以學到基本的數值微分方程的精神。舉例來說,若我們有一個微分方程式為

\(y'(t)+a y(t) = r\), \(y(0) = y_0\)

這個問題所對應的真解是

\(y(t) = (y_0-\frac{r}{a})e^{-at}+\frac{r}{a}\)

這類形式的微分方程式,可以用Laplace 轉換來求得真解。當我們知道真解時,便能方便地與數值解作比較。我們如果把 \(y'(t)\) 這一項改用 \(\frac{y(t+h)-y(t)}{h}\) 來代替,我們得到

\(y(t+h) = (1-ah)y(t) + hr\)

用一個差分式來代替微分值,沒有意外會有誤差。當然,原來的微分方程式變成新的迭代式子所代出來的 \(y(t+h)\) 值當然也會有誤差。這類的誤差不是由電腦的浮點數儲存限制所造成的,而是由算法本身所造成的。我們從上面的迭代式得知,當我們知道 \(y(0)\) 的值,我們可以推得 \(y(h)\) 的值,當我們知道 \(y(h)\) 的值,我們又可以推得\(y(2h)\) 的值,依此類推,我們可以得到每一個 \(y(nh)\) 的值。我們也可以預見,誤差也會由第一次的迭代,累積到第二次,第三次,到第 \(n\) 次的迭代。

到底誤差是怎麼被累積的呢?我們用泰勒展示,把 \(y(t+h)\) 寫成對 \(y(t)\) 的展開,就可以很明顯的看出來。其實迭代式帶著一個 \(O(h^2)\) 的誤差。這裡的\(O(h^2)\) 詳細寫下來是 \(\frac{f^{(2)}(\xi)}{2!}\) ,其中 \(\xi\) 是介於 \((t,t+h)\) 的一個常數。因此,若 \(t\) 附近的 \(f^{(2)}\) 很靠近 \(0\) ,那們誤差被放大的情況會好一點,反之,誤差會很嚴重。我們只是從 \(y_0\) 往前推進一步,這裡我們把 truncation error 寫入迭代式中,我們看到走這一步會產生 \(O(h^2)\) 的誤差,我們稱 Euler method 是一階的誤差。

 

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

    y = np.array(y)
    
    return y

 

使用這個Euler ODE Solver,我們要先定義一個函數 \(f(t,y)\) 作為 \(y'(t)\),然後帶入此函式。用法如下:

>>> from numpy import *
>>> def f(t,y):
.... return 4-3*y
....
>>>
>>> y = ode_Euler(lambda x,y: f(x,y),0,2,1,100)
>>> from pylab import *
>>> t = linspace(0,2,100)
>>> plot(t,y)
>>> fig1 = figure(1)
>>> fig1.show()
>>> fig1.savefig('ode.png')

要將一個函數引入另外一個函數,我們使用lambda :這個指令,讓Python 程式可以分辨哪些變數是屬於外部函數。要看看到底計算出來的東西長甚麼樣子,我們使用pylab 內的套件繪圖。t = linspace(0,2,100) 會建立一個建立一個numpy.narray,起始值是 0,結束於 2,長度剛好是 100 的陣列。fig1 = figure(1) 會指定 fig1 為圖形的型別,並且指明是第一張圖。因為我們知道這個範例的真解為何,所以我們可以將真解與數值解同時畫在一張圖上。範例如下:

>>> yt = (1.0-4.0/3)*exp(-3*t)+4.0/3
>>> plot(t,y,'r',t,yt,'b')
[<matplotlib.lines.Line2D object at 0xa6a9b0c>, <matplotlib.lines.Line2D object
at 0xa6b67cc>]
>>> fig2 = figure(2)
>>> plot(t,y,'r',t,yt,'b')
[<matplotlib.lines.Line2D object at 0xa6d658c>, <matplotlib.lines.Line2D object
at 0xa8962cc>]
>>> fig2.show()
>>> fig2.savefig('ode2.png')

其中 plot() 中的參數使用和 octave 類似,'r' 代表紅色,'b' 代表藍色。