數值偏微分方程

Parabolic PDE

接下來我們用一維的熱傳導方程來介紹Parabolic PDE,我們令 \(u(x,t)\) 是熱在一維的 \(x\) 空間上隨時間 \(t\) 的分佈,方程如下:

\(A \frac{\partial^{2}u(x,t)}{\partial x^{2}} = \frac{\partial u(x,t)}{\partial t}\),

其中 \(0\leq x\leq x_f\), \(0\leq t \leq T\)

並且有下列邊界條件:

\(u(0,t) = b_0(t)\), \(u(x_f,t) = b_f(t)\),以及 \(u(x,0) =i_0(x)\).


The Explicit Forward Euler Method

我們一樣使用finite difference method 把 \([0,x_f]\) 區間切割成 \(M\) 等分,把 \([0,T]\) 時間區間切割成 \(N\) 等分。因此,我們有 \(\Delta x = \frac{x_f}{M}\) 以及 \(\Delta t = \frac{T}{N}\)。並且我們用central difference approximation 去取代 \(u(x,t)\) 在空間上的二次偏微分,用forward difference approximation 取代對時間的偏微分。我們得到下列差分方程式:

\(A\frac{u_{i+1}^k-2u_i^{k}+u_{i-1}^{k}}{\Delta x^{2}} = \frac{u_{i}^{k+1}-u_{i}^{k}}{\Delta t}\)

這裡的 \(u_{i}^{k}\) 代表 \(u(x_{i},t_{k})\) 因此 \(u_{i}^{k+1}\) 的時間比其他項更新,我們可以將上式寫成下列的疊代式:

(PB 式一) \(u_{i}^{k+1} = r(u_{i+1}^{k}+u_{i-1}^{k})+(1-2r)u_{i}^{k}\) ,其中 \(r = A\frac{\Delta t}{\Delta x^{2}}\) ,\(i = 1,2,..., M-1\)

和先前的Laplace's equation 不一樣,因為不知道T 時間的狀態,我們隨著時間的增加,若先前的計算有誤差,不知道這誤差會怎麼被累積下去。這和解常微分方程的Euler 方法一樣,我們擔心數值解會有不穩定的情形發生。若我們希望這個疊代式穩定,(PB式一)會對應一個線性系統 \(Au^{k+1} = u^{k}\),理所當然地矩陣 \(A\) 的特徵值要小於1。我們除了把對應的矩陣寫開並且計算其特徵值之外,還可以用下列方法決定適合的 \(\Delta t\) 和 \(\Delta x\)。

我們把 \(u_{i}^{k}\) 用傅立葉基底表示,假設

\(u_{i}^{k} = \lambda^{k}e^{j\frac{i\pi}{P}}\) 其中 \(j\) 是 \(\sqrt{-1}\) 的意思,\(P\) 是任意非零的整數
代入(PB式一),我們得到

\(\lambda = 1-2r(1-cos(\frac{\pi}{P}))\)

當我們希望 \(u_{i}^{k}\) 不會隨著 \(k\) 增加而爆掉時,我們希望 \(|\lambda|\leq 1\),因此我們有

\(r=A\frac{\Delta t}{\Delta x^{2}} \leq \frac{1}{2}\)

用Explicit Forward Euler Method 解 heat equation 的Python 程式如下,若要解一般的Parabolic PDE 對應的函數與公式讀者要自行調整。

def pde_heat_exp(a,xf,T,Fx,M,N):
    """
    solve a u_xx = u_t for 0<=x<=xf, 0<=t<=T by explicit method
    u(x,0) = it0(x)
    u(0,t) = bx0(t)
    u(xf,t) = bxf(t)
    M : the number of subintervals along x axis
    N : the number of subintervals along t axis
    Fx : the list of function, Fx = [it0,bx0,bxf]
    
    For example:
    Fx = [lambda x:np.sin(np.pi*x), lambda t: 0.0*t, lambda t: 0.0*t]
    u, x, t = pde_heat_exp(1.0,1.0,0.1,Fx,20,100)
    Notice that r = a*dt/dx^2 must < 1/2 for stablility
    """

    dx, dt = float(xf)/M, float(T)/N
    x = np.linspace(0,xf,M+1)
    t = np.linspace(0,T,N+1)
    
    if len(Fx)<>3:
        print "Fx = [it0,bx0,bxf]"
        return
    it0, bx0, bxf = Fx[0], Fx[1], Fx[2]

    u = np.zeros([N+1,M+1])
    X = np.zeros([N+1,M+1]) #the mesh grids of x
    Y = np.zeros([N+1,M+1]) #the mesh grids of t

    for i in range(N+1):
        for j in range(M+1):
            X[i,j] = x[j]
            Y[i,j] = t[i]
            
    #set the boundary condition
    u[0,:] = it0(x)
    u[:,0] = bx0(t)
    u[:,-1] = bxf(t)

    r = a*dt/(dx**2)
    r1 = 1 - 2*r
    
    for k in range(N):
        for i in range(1,M):
            u[k+1,i] = r*(u[k,i+1]+u[k,i-1])+r1*u[k,i]
    
    return u, X, Y

 

因為是3D 的圖形,所以output 的 \(x\) 與 \(t\) 我們調整成和u 一樣大小,方便繪圖。要顯示我們的計算結果,Python 的3D 繪圖方法如下:

 

>>> Fx = [lambda x:np.sin(np.pi*x), lambda t: 0.0*t, lambda t: 0.0*t]
>>> u, x, t = pde_heat_exp(1.0,1.0,0.1,Fx,20,100)
>>> import pylab as p
>>> import mpl_toolkits.mplot3d.axes3d as p3
>>> fig = p.figure()
>>> ax = p3.Axes3D(fig)
>>> ax.plot_wireframe(u,x,t)
<mpl_toolkits.mplot3d.art3d.Line3DCollection object at 0x1d500d0>
>>> p.show()

 

The Implicit Backward Euler Method

若我們把差分方程式的右端改成backward 差分,我們得到

\(A\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{\Delta x^{2}} = \frac{u_{i}^{k}-u_{i}^{k-1}}{\Delta t}\)

上式整理過後成為

(PB 式二) \(-r(u_{i-1}^{k}+(1+2r)u_{i}^{k}-ru_{i+1}^{k}) = u_{i}^{k-1}\) ,其中 \(r = A\frac{\Delta t}{\Delta x^{2}}\) ,\(i = 1,2,..., M-1\)

這時我們無法對每一個 \(i\) 位置的 \(u\) 進行單純地迭代得到答案,而是解一個三對角的矩陣。若我們解的是Dirichlet 形式的邊界條件問題,我們有 \(u_{0}^{k}\) 和 \(u_{M}^{k}\)。若我們解的是Neumann 形式的邊界條件問題,我們使用 \(\frac{u_{1}^{k}-u_{-1}^{k}}{2\Delta x} = b'_{0}(k)\),這樣 \(u_{-1}^{k}\) 這一項才能被 \(b'_{0}(k)\) 表示出來。雖然解矩陣是比較麻煩的方式,但若我們套用先前的傅立葉分析方法將 \(u_{i}^{k} = \lambda^{k}e^{j\frac{i\pi}{P}}\)帶入(PB 式二),經過整理算式後我們發現\(|\lambda|\) 自動小於 \(1\),這表示不用特別調整 \(\Delta x\) 和 \(\Delta t\) 的大小,我們用Implicit Backward Euler Method 可以得到比較穩定的解。

下列是Python 的heat equation 用 implicit backward Euler Method 的程式:

def pde_heat_imp(a,xf,T,Fx,M,N):
    """
    solve a u_xx = u_t for 0<=x<=xf, 0<=t<=T by explicit method
    u(x,0) = it0(x)
    u(0,t) = bx0(t)
    u(xf,t) = bxf(t)
    M : the number of subintervals along x axis
    N : the number of subintervals along t axis
    Fx : the list of function, Fx = [it0,bx0,bxf]
    
    For example:
    Fx = [lambda x:np.sin(np.pi*x), lambda t: 0.0*t, lambda t: 0.0*t]
    u, x, t = pde_heat_imp(1.0,1.0,0.1,Fx,20,100)
    """

    dx, dt = float(xf)/M, float(T)/N
    x = np.linspace(0,xf,M+1)
    t = np.linspace(0,T,N+1)
    
    if len(Fx)<>3:
        print "Fx = [it0,bx0,bxf]"
        return
    it0, bx0, bxf = Fx[0], Fx[1], Fx[2]

    u = np.zeros([N+1,M+1])
    X = np.zeros([N+1,M+1]) #the mesh grids of x
    Y = np.zeros([N+1,M+1]) #the mesh grids of t

    for i in range(N+1):
        for j in range(M+1):
            X[i,j] = x[j]
            Y[i,j] = t[i]
            
    #set the boundary condition
    u[0,:] = it0(x)
    u[:,0] = bx0(t)
    u[:,-1] = bxf(t)

    r = a*dt/(dx**2)
    r2 = 1 + 2*r

    #set the linear system
    A = np.zeros([M-1,M-1])
    for i in range(M-1):
        A[i,i] = r2
        if i>1:
            A[i-1,i] = -r
            A[i,i-1] = -r
    #compute each time step        
    for k in range(N):
        b = u[k,1:-1]
        b[0] += r*u[k,0]
        b[-1] += r*u[k,-1]
        u[k+1,1:-1] = lg.solve(A,b)
    
    return u, X, Y

 

這個程式的使用方法與ode_heat_exp() 類似,讀者可以試著調整M, N 的值,比較兩者的差別。

 

The Crank-Nicholson Method

回顧一下implicit backward Euler method 的遞迴公式

\(A\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{\Delta x^{2}} = \frac{u_{i}^{k}-u_{i}^{k-1}}{\Delta t}\)

請注意一下,左端點的差分公式是對 \(u(x_{i},t_{k})\) 做展開,而右端的差分公式是對 \(u(x_{i},(t_{k}+t_{k-1})\) 這一點展開,這裡存在時間上的不一致,就會增加更多的truncation error。Crank-Nicholson Method 就是為了改進時間上不一致而有的方法,他的差分方程式如下:

\(\frac{A}{2}(\frac{u_{i+1}^{k+1}-2u_{i}^{k+1}+u_{i-1}^{k+1}}{\Delta x^{2}}+\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{\Delta x^{2}}) = \frac{u_{i}^{k+1}-u_{i}^{k}}{\Delta t}\)

經過這樣的調整,等號兩邊的位置與時間都保持一致,簡化後的差分式如下:

\(-ru_{i+1}^{k+1}+2(1+r)u_{i}^{k+1}-ru_{i-1}^{k+1} = ru_{i+1}^{k}+2(1-r)u_{i}^{k}+ru_{i-1}^{k}\) 其中 \(r = A\frac{\Delta t}{\Delta x^2}\)

仿照同樣的分析穩定性的方法,我們得到對應的\(|\lambda|\) 也是自動小於等於 \(1\) 。說明這個方法也是一個穩定的解法,並且我們預期它的誤差會比implicit backward Euler method 來的好。對應的Crank-Nicholson 程式如下:

def pde_heat_CN(a,xf,T,Fx,M,N):
    """
    solve a u_xx = u_t for 0<=x<=xf, 0<=t<=T by explicit method
    u(x,0) = it0(x)
    u(0,t) = bx0(t)
    u(xf,t) = bxf(t)
    M : the number of subintervals along x axis
    N : the number of subintervals along t axis
    Fx : the list of function, Fx = [it0,bx0,bxf]
    
    For example:
    Fx = [lambda x:np.sin(np.pi*x), lambda t: 0.0*t, lambda t: 0.0*t]
    u, x, t = pde_heat_CN(1.0,1.0,0.1,Fx,20,100)
    """

    dx, dt = float(xf)/M, float(T)/N
    x = np.linspace(0,xf,M+1)
    t = np.linspace(0,T,N+1)
    
    if len(Fx)<>3:
        print "Fx = [it0,bx0,bxf]"
        return
    it0, bx0, bxf = Fx[0], Fx[1], Fx[2]

    u = np.zeros([N+1,M+1])
    X = np.zeros([N+1,M+1]) #the mesh grids of x
    Y = np.zeros([N+1,M+1]) #the mesh grids of t

    for i in range(N+1):
        for j in range(M+1):
            X[i,j] = x[j]
            Y[i,j] = t[i]
            
    #set the boundary condition
    u[0,:] = it0(x)
    u[:,0] = bx0(t)
    u[:,-1] = bxf(t)

    r = a*dt/(dx**2)
    r1 = 2*(1+r)
    r2 = 2*(1-r)

    #set the linear system
    A = np.zeros([M-1,M-1])
    for i in range(M-1):
        A[i,i] = r1
        if i>1:
            A[i-1,i] = -r
            A[i,i-1] = -r
    #compute each time step        
    for k in range(N):
        b = r*(u[k,0:M-1]+u[k,2:M+1])+r2*u[k,1:M]
        b[0] += r*(u[k,0]+u[k+1,0])
        b[-1] += r*(u[k,-1]+u[k+1,-1])
        u[k+1,1:-1] = lg.solve(A,b)
    
    return u, X, Y

 

我們用下列heat equation 來驗證上述三種數值分法的優劣。考慮parabolic PDE

\(\frac{\partial^{2} u(x,t)}{\partial x^{2}} = \frac{\partial u(x,t)}{\partial t}\), for \(0\leq x \leq 1\), \(0 \leq t \leq 0.1\)

且邊界條件為

\(u(x,0) = sin(\pi x)\), \(u(0,t) = 0\), \(u(1,t) = 1\)

這個方程的真解為

\(u(x,t) = sin(\pi x)e^{-\pi ^{2}t}\)

我們用下列程式來觀察這三種方法的比較。

from glophynum import *
import numpy.linalg as lg

a, xf, T = 1.0, 1.0, 0.1
M, N = 25, 100
Fx = [lambda x:np.sin(np.pi*x), lambda t: 0.0*t, lambda t: 0.0*t]
u_exp, x, t = pde_heat_exp(a,xf,T,Fx,M,N)
u_imp, x, t = pde_heat_imp(a,xf,T,Fx,M,N)
u_cn, x, t = pde_heat_CN(a,xf,T,Fx,M,N)

#evalue the true solution
def u(x,t):
    return np.sin(np.pi*x)*np.exp(-t*np.pi**2)

m,n = x.shape
u_true = np.zeros([m,n])
for i in range(m):
    for j in range(n):
        u_true[i][j] = u(x[i][j],t[i][j])

err_exp = lg.norm(u_true-u_exp)
err_imp = lg.norm(u_true-u_imp)
err_cn = lg.norm(u_true-u_cn)

print 'Error of exp method is: ', err_exp
print 'Error of imp method is: ', err_imp
print 'Error of CN method is: ', err_cn

import pylab as p
import mpl_toolkits.mplot3d.axes3d as p3
fig1 = p.figure(1)
fig2 = p.figure(2)
fig3 = p.figure(3)
ax1 = p3.Axes3D(fig1)
ax2 = p3.Axes3D(fig2)
ax3 = p3.Axes3D(fig3)
ax1.plot_wireframe(u_exp,x,t)
ax2.plot_wireframe(u_imp,x,t)
ax3.plot_wireframe(u_cn,x,t)
p.show()

這個例子,explicit method 是爆掉的。我們也看到,再同樣的條件下,Crank-Nicholson 方法是最好的。