數值常微分方程(BVP)

angry bird 

 

Shooting Method 的精神就和玩憤怒鳥遊戲一樣,一開始彈弓的位置是 \(y_0\),使用者可以調整發射的角度 \(y'_0\),然後發射(當我們有 \(y_0\) 和 \(y'_0\) 時,我們就可以用常用的 IVP 的解法,求 \(y_i\))。如果沒有打到 \(y_f\),就調整 \(y'_0\) 直到打中為止。說到這裡,很自然地隱含著這類 BVP 問題至少是二次的常微分方程。若是一次的常微分方程,\(y'_0\) 是被起始值決定,根本沒有空間讓我們調整。也因此這類問題的 IVP 解法也是高階的常微分方程解法。這是常微分方程課程裡的常識,利用變數變換,把高階的常微分方程變成一階向量形的常微分方程。例如

\(y^{(2)}(t) = f(t,y(t),y'(t))\)

我們令 \(y_1(t) = y(t)\),\(y_2(t) = y'(t)\),則

\(y'_1(t) = y_2(t)\\

y'_2(t) = f(t,y_1(t),y_2(t))\)

我們用向量 \(Y(t) = (y_1(t), y_2(t))^T\) 把上述方程組,看成 \(Y'(t) = f(t,Y(t))\)。直覺上當我們的 \(y(t)\) 變成向量型式的 \(Y(t)\), \(f(t,y)\) 也變成向量型式時,對應的 Python code 需要做一些調整。但是非常剛好的是,我們先前使用的 Python code 可以完全不用更動,唯獨要注意的是 \(f\) 的 output 和\(y_0\) 要換成 numpy.ndarray 型別。假設我們要計算的微分方程式為

\(y^{(2)}(t) = 2y^2(t) + 4ty(t)y'(t)\)

我們令 \(y_1 = y\), \(y_2 = y'\) 得到

\(y'_1(t) = y_2(t)\\

y'_2(t) = (2y_1(t)+4ty_2(t))y_1(t)\)

在 Python 裡我們定義 \(F(t,y_1,y_2)\) 如下:

>>> def F(t,Y):
.... return np.array([Y[1],Y[0]*(2*Y[0]+4*t*Y[1])])

意思是Y[0] 是 \(y(t)\),Y[1] 是 \(y'(t)\),輸出時用 numpy 的函式 array() 更換型別。這樣就可以正常套用原本的 ode_RK4()。而 BVP 的 Shooting method 對應的 Python 程式如下:

def ode_shooting(F,t0,tf,y0,yf,n=100,tol=10**(-8),kmax=10):
    """
    The solver of Boundary Value Problem (y'' = f(t,y,y')) by shooting method

    F   : the vector type function
    t0  : initial value of independent variable
    tf  : final value of independent variable
    y0  : initial value of dependent variable
    yf  : final value of dependent variable
    tol : the tolerance 
    n   : number of time steps
    kmax: the maximal number of shootings  

    The return is a list of y(t)
    """
    yp0 = float(yf-y0)/(tf-t0) #The initial guess of y'(0)
    Y0 = np.array([y0,yp0])
    Y = ode_RK4(F,t0,tf,Y0,n)
    err = Y[-1][0]-yf #The 1st mismatching
    Y0[1] = yp0-0.1*np.sign(err)
    
    for k in range(kmax-1):
        Y = ode_RK4(F,t0,tf,Y0,n)
        err1 = Y[-1][0]-yf
        ddy = Y0[1]-yp0 #difference between successive derivatives
        if abs(err1)<tol or abs(ddy)<tol:
            y = list()
            for x in Y:
                y.append(x[0])
            return np.array(y), Y0[1], err

        deddy = (err1-err)/ddy #The gradient of mismatching error
        yp0 = Y0[1]
        err = err1
        Y0[1] = Y0[1] -err/deddy #modify by secant method
    y = list()
    for x in Y:
        y.append(x[0])
    return np.array(y), Y0[1], err

 

從這個程式看到,要執行向量形式的 RK4,我們輸入的\(Y_0\) 也是用 numpy 的 array() 函式更改其型別,這樣就可以做向量的加法,避免程式的錯誤。程式中預設的 tolerance 是 \(10^{-8}\),並且我們用 Secant method 來進行對 \(y'(t_0)\) 的修正。

對於上列問題,假設起始值是 \(y(t_0)=\frac{1}{4}\),邊界值 \(y(t_f)=\frac{1}{3}\),時間區間為 \([0,1]\),計算方法如下:

>>> def F(t,Y):
return np.array([Y[1],Y[0]*(2*Y[0]+4*t*Y[1])])

>>> y,yp,err = ode_shooting(F,0,1,1/4.0,1/3.0,100)
>>> len(y)
100
>>> yp
9.9651168851077772e-11
>>> err
-1.076104994002236e-08