數值常微分方程(BVP)

 

Finite Difference Method


Finite Difference Method 的概念是把討論的區間 \([t_0,t_f]\) 切割成 \(N\) 個均勻片段,然後把 \(y_i = y(x_i)\) 當作未知變數,利用 \(\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}\) 代替 \(y^{(2)}(t_i)\) 以及 \(\frac{y_{i+1}-y_{i-1}}{2h}\) 代替 \(y'(t_i)\) ,帶入要解的 BVP 方程式中,會成為一個\(y_i\) 的線性方程組。例如 BVP 方程為

\(y^{(2)}(t) + a_1(t)y'(t) + a_0(t)y(t) = u(t)\), with \(y(t_0) = y_0\), \(y(t_f) = y_f\)

則對應的方程組為

\((2-ha_{1,i})y_{i-1} + (-4+2h^{2}a_{0,i})y_i + (2+ha_{1,i})y_{i+1} = 2h^2u_i\)

其中 \(a_{1,i} = a_1(t_i)\), \(a_{0,i} =a_0(t_i)\), \(u_i = u(t_i)\),\(i = 1,\dots,N-1\)。這個系統剛好是解一個三對角的矩陣。對於 \(i=1\) 及\(i=N-1\) ,因為兩端點 \(y_0\) 與 \(y_f\) 的值是給定的,所以會把已知項搬到右邊。這裡我們又看到一個解三對角矩陣的例子,並且當 \(N\) 很大時,對應的矩陣也跟著變大,這樣的矩陣都很容易有一個特性,就是隨著矩陣的變大,condition number 也越便越糟。因此,如何精準地解一個三對角矩陣,自然就變成數值分析很重要的一類課題。

Finite difference method 的Python 程式如下:

def ode_FD2(F,t0,tf,y0,yf,n=100):
    """
    The solver of Second order Boundary Value Problem
        y'' + a_1(t)y'+a_0(t)y = u(t)  by finite difference method

    F   : the list of function [a_1(t), a_0(t), u(t)]
    t0  : initial value of independent variable
    tf  : final value of independent variable
    y0  : initial value of dependent variable
    yf  : final value of dependent variable
    n   : number of time steps

    The return is a list of y(t)
    """
    import numpy.linalg as lg
    h = float(tf-t0)/n
    t = np.linspace(t0,tf,n+1)
    h2 = 2*h**2
    a1, a0, u = F[0], F[1], F[2]

    A = np.zeros([n-1,n-1])
    c1 = a1(t[1:-1])
    c2 = a0(t[1:-1])
    b = u(t[1:-1])

    ha = h*c1[0]
    #the first equation
    A[0,0] = -4+h2*c2[0] 
    A[0,1] = 2+ha
    b[0] = b[0] + (ha-2)*y0

    #the middle equations
    for m in range(1,n-2):
        ha = h*c1[m]
        A[m,m-1] = 2-ha
        A[m,m] = -4+h2*c2[m]
        A[m,m+1] = 2+ha

    #the last equation
    ha = c1[-1]
    A[-1,-2] = 2-ha
    A[-1,-1] = -4+h2*c2[-1]
    b[-1] = b[-1]-(ha+2)*yf

    #solve the tridigonal linear system
    x = lg.solve(A,b)
    y = [y0]
    y.extend(list(x))
    y.append(yf)

    return np.array(y)

 

 

使用這個函式要很小心,特別注意我們定義 \(F\) 的方法是一個函數所構成的串列,並且我們解的二階常微分方程為

\(y^{(2)} + a_1(t)y'+a_0(t)y = u(t)\)

這種形式。假設我們要解的題目為

\(y^{(2)}(t) + \frac{2}{t}y'(t) - \frac{2}{t^2}y(t) = 0\) with \(y(1) = 5\) , \(y(2) = 3\)

那麼使用 ode_FD2() 的計算方式如下:

 

>>> def a1(t):
return 2.0/t

>>> def a0(t):
return 2.0/(t**2)

>>> def u(t):
return 0.0*t

>>> F = [a1,a0,u]
>>> y = ode_FD2(F,1,2,5,3,100)

 

這裡要特別注意,因為要設定 \(u(t) = 0\) ,若 \(u(t)\) 中的 return 後面直接接 \(0\),則程式無法把向量形的 \(t\)  對應到向量形的 output,一定要讓 \(t\) 出現在 return 之後,output 的形式才會自動變成向量型。