數值偏微分方程

Hyperbolic PDE


著名的Hyperbolic PDE 是波動方程(Wave equation),對於位置 \(x\) 與時間 \(t\) 的振幅函數 \(u(x,t)\) 滿足下列形式

\(A\frac{\partial^{2} u(x,t)}{\partial x^{2}} = \frac{\partial^{2} u(x,t)}{\partial t^{2}}\) for \(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)\) 以及 \(\frac{\partial u}{\partial t}|_{t=0} (x,0) = i'_{0}(x)\)

 

The Explicit Central Difference Method

 

因為等號兩端都是二次偏導數,理所當然地我們使用Central Difference 來代替二次偏導數,得到的差分方程如下:

\(A(\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{\Delta x^{2}}) = \frac{u_{i}^{k+1}-2u_{i}^{k}+u_{i}^{k-1}}{\Delta t^{2}}\) 其中 \(\Delta x = \frac{x_{f}}{M}\), \(\Delta t = \frac{T}{N}\)

經過整理後得到

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

由於我們沒有 \(u_{i}^{-1}\) 的值,所以我們無法開始進行疊代的工作。我們利用 \(\frac{u_{i}^{1}-u_{i}^{-1}}{2\Delta t}\) 來逼近 \(i'_{0}(x_{i})\)。對於 \(k=1\) 的情況,經過整理後我們有

\(u_{i}^{1} = \frac{1}{2}r(u_{i+1}^{0}+u_{i-1}^{0}) + (1-r)u_{i}^{0} + i'_{0}(x_{i})\Delta t\)

一旦有了 \(k=1\) 以及 \(k=0\) 的值,我們就可以開始進行疊代。和先前分析穩定性的方式一樣,我們計算出 \(r\) 的穩定範圍為小於等於 \(1\)。對應的Python 程式如下:

def pde_wave_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_t(x,0) = i1t0(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,i1t0,bx0,bxf]
    
    For example:
    Fx = [lambda x:x*(1-x), lambda x: 0.0*x, lambda t: 0.0*t, lambda t: 0.0*t]
    u, x, t = pde_wave_exp(1.0,1.0,2,Fx,20,50)
    Notice that r = a*dt^2/dx^2 must <=1  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)<>4:
        print "Fx = [it0,i1t0,bx0,bxf]"
        return
    it0, i1t0, bx0, bxf = Fx[0], Fx[1], Fx[2],Fx[3]

    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*dt/(dx**2)
    r1 = r/2
    r2 = 2*(1-r)

    #set the u(x,t_{1})
    u[1,1:-1] = r1*u[0,0:-2]+(1-r)*u[0,1:-1]+r1*u[0,2:]+dt*i1t0(x[1:-1])
    
    for k in range(2,N+1):
        u[k,1:-1] = r*u[k-1,0:-2]+r2*u[k-1,1:-1]+r*u[k-1,2:]-u[k-2,1:-1]
        
    return u, X, Y

 

我們執行下列程式用動態繪圖觀察我們計算的結果。

>>> from glophynum import *
>>> Fx = [lambda x:x*(1-x), lambda x: 0.0*x, lambda t: 0.0*t, lambda t: 0.0*t]
>>> u,x,t = pde_wave_exp(1.0,1.0,2,Fx,20,50)
>>> import pylab as p
>>> import time
>>> u.shape
(51, 21)
>>> p.ion()
>>> p.hold(False)
>>> for i in range(51):
... ll = p.plot(x[i,:],u[i,:])
... ax = p.axis([0,1,-0.3,0.3])
... p.draw()
... time.sleep(1)
...

因為我們要畫動態的圖,所以我們在p.plot() 指令前給一個變數,讓Python 不要在shell 裡面顯示一些我們不太關注的顯示項目。axis() 函式是強迫繪圖的範圍,這樣才不會因為Python 自動調整軸的級距,誤導我們對圖形的認識。其他指令,請參考Glophy 首頁上的 Python 入門。

Two-Dimensional Hyperbolic PDE

這一節我們考慮二維的波動方程u(x,y,t) :

\(A(\frac{\partial^{2} u(x,y,t)}{\partial x^{2}}+\frac{\partial^{2} u(x,y,t)}{\partial y^{2}}) = \frac{\partial^{2} u(x,y,t)}{\partial t^{2}}\) for \(0 \leq x \leq x_{f}\), \(0 \leq y \leq y_{f}\), \(0 \leq t \leq T\)

其邊界條件為
\( u(0,y,t) = b_{x_0}(y,t)\), \(u(x_{f},y,t) = b_{x_f}(y,t)\), \(u(x,0,t) = b_{y_0}(x,t)\) 及 \(u(x,y_{f},t) = b_{y_f}(x,t)\)

並且給予起始條件

\(u(x,y,0) = i_{0}(x,y)\) 以及 \(\frac{\partial u}{\partial t}|_{t=0}(x,y,0) = i'_{0}(x,y)\)

和一維的波動方程一樣,我們用central difference 方法來替代偏導數值,得到

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

其中 \(\Delta x = \frac{x_{f}}{M_{x}}\), \(\Delta y = \frac{y_{f}}{M_{y}}\), \(\Delta t = \frac{T}{N}\)

經過整理後得到

\(u_{i,j}^{k+1} = r_{x}(u_{i,j+1}^{k}+u_{i,j-1}^{k}) +2(1-r_{x}-r_{y})u_{i,j}^{k} + r_{y}(u_{i+1,j}^{k}+u_{i-1,j}^{k})-u_{i,j}^{k-1}\)

其中 \(r_{x} = A\frac{\Delta t^{2}}{\Delta x^{2}}\), \(r_{y} = A\frac{\Delta t^{2}}{\Delta y^{2}}\)

和一維的起始情況一樣,因為我們沒有時間 \(k=-1\) 的值,所以我們用central difference 配合 \(k=0\) 時的偏導數,把 \(u_{i,j}^{-1}\) 用偏導數值替換。經過整理得到

\(u_{i,j}^{1} = \frac{1}{2}\{r_{x}(u_{i,j+1}^{0}+u_{i,j-1}^{0}) + r_{y}(u_{i+1,j}^{0}+u_{i-1,j}^{0})\}+2(1-r_{x}-r_{y})u_{i,j}^{0} + i'_{0}(x_{j},y_{i})\Delta t\)

這個差分式的穩定條件為

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

這個章節特別提到二維的情形只是讓讀者有一個認識,一維時我們在定義域上做切割,若不是explicit 的形式,對應的矩陣會隨著切割定義域的細緻度而增加,可想而知,若二維的情形我們也用解矩陣的方式處理,那麼矩陣會有多大?再來,當我們處理真實世界的3D 問題時,我們需要多少的計算資源與技巧才能顯現3D 的即時結果。

2D 的波動方程Python 程式如下:

def ode_wave2(a,D,T,Fx,Mx,My,N):
    """
    solve a(u_xx + u_yy) = u_tt for D[0]<=x<=D[1], D[2]<=y<=D[3],
    0<=t<=T with
    Initial condition: u(x,y,0) = Fx[0](x,y),
                       u_t(x,y,0) = Fx[1](x,y)
    Boundary condition: u(0,y,t) = Fx[2](y,t)
                        u(x_f,y,t) = Fx[3](y,t)
                        u(x,0,t) = Fx[4](x,t)
                        u(x,y_f,t) = Fx[5](x,t)
    Mx: the number of subintervals along x axis
    My: the number of subintervals along y axis
    N : the number of subintervals along t axis

    For example:
    D = [0,2,0,2]
    Fx = [lambda x,y:0.1*np.sin(np.pi*x)*np.sin(np.pi*y/2),
          lambda x,y:0.0*(x+y),
          lambda y,t:0.0*(y+t),
          lambda y,t:0.0*(y+t),
          lambda x,t:0.0*(x+t),
          lambda x,t:0.0*(x+t)]
          
    u,X,Y,t = ode_wave2(0.25,D,2,Fx,40,40,40)
    """
    dx = float(D[1]-D[0])/Mx
    dy = float(D[3]-D[2])/My
    x = np.linspace(D[0],D[1],Mx+1)
    y = np.linspace(D[2],D[3],My+1)
    dt = float(T)/N
    t = np.linspace(0,T,N+1)

    it0,i1t0,u0y,ufy,ux0,uxf = Fx[0],Fx[1],Fx[2],Fx[3],Fx[4],Fx[5]

    
    #Initialization
    u = np.zeros([N+1,My+1,Mx+1]) #u(x_j,y_i,t_k) = u[k,i,j], k = 0,1,2,...,N
    ut = np.zeros([My+1,Mx+1])
    for j in range(1,Mx):
        for i in range(1,My):
            u[0,i,j] = it0(x[j],y[i])
            ut[i,j] = i1t0(x[j],y[i])
    adt2 = a*dt*dt
    rx, ry = adt2/(dx*dx), adt2/(dy*dy)
    rxy1 = 1-rx-ry
    rxy2 = rxy1*2

    for k in range(1,N+1):
        #set boundary condition
        u[k,0,:] = u0y(y,t[k])
        u[k,-1,:] = ufy(y,t[k])
        u[k,:,0] = ux0(x,t[k])
        u[k,:,-1] = uxf(x,t[k])

        if k==1:
            for i in range(1,My):
                for j in range(1,Mx):
                    u[k,i,j] = 0.5*(rx*(u[k-1,i,j-1]+u[k-1,i,j+1])+ry*(u[k-1,i-1,j]+u[k-1,i+1,j]))+rxy1*u[k-1,i,j]+ut[i,j]*dt
        else:
            for i in range(1,My):
                for j in range(1,Mx):
                    u[k,i,j] = rx*(u[k-1,i,j+1]+u[k-1,i,j-1])+rxy2*u[k-1,i,j]+ry*(u[k-1,i+1,j]+u[k-1,i-1,j])-u[k-2,i,j]
     
    m,n = My+1, Mx+1
    X = np.zeros([m,n])
    Y = np.zeros([m,n])
    for i in range(m):
        for j in range(n):
            X[i,j] = x[j]
            Y[i,j] = y[i]

    return u,X,Y,t

 

若我們解下列波動方程

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

給予起始條件

\(u(x,y,0) = 0.1*sin(\pi x) sin(0.5 \pi y)\),

\(\frac{\partial u(x,y,t)}{\partial t} = 0\), for \(t = 0\)

其餘邊界條件都是 \(0\) 函數。我們用下列Python 程式計算。

>>> D = [0,2,0,2]
>>> Fx = [lambda x,y:0.1*np.sin(np.pi*x)*np.sin(np.pi*y/2),
lambda x,y:0.0*(x+y),
lambda y,t:0.0*(y+t),
lambda y,t:0.0*(y+t),
lambda x,t:0.0*(x+t),
lambda x,t:0.0*(x+t)]
>>> u,X,Y,t = ode_wave2(0.25,D,2,Fx,40,40,40)
>>> u.shape
(41, 41, 41)

 

若我們想看動態的圖形,可以參考下列範例。

from glophynum import *
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import time

plt.ion()
fig = plt.figure()
ax = axes3d.Axes3D(fig)

D = [0,2,0,2]
Fx = [lambda x,y:0.1*np.sin(np.pi*x)*np.sin(np.pi*y/2),
          lambda x,y:0.0*(x+y),
          lambda y,t:0.0*(y+t),
          lambda y,t:0.0*(y+t),
          lambda x,t:0.0*(x+t),
          lambda x,t:0.0*(x+t)]

u,X,Y,t = ode_wave2(0.25,D,2,Fx,40,40,40)

wframe = None
for i in range(41):
    oldcol = wframe
    wframe = ax.plot_wireframe(u[i,:,:], X, Y, rstride=2, cstride=2)

    if oldcol is not None:
        ax.collections.remove(oldcol)

    ax.view_init(105,20)

    plt.draw()
    time.sleep(0.5)
    
plt.close()

 

為了便利同學把精神放在計算而不是程式,我把關於繪圖的部份寫了兩個靜態繪圖套件EZplot2D() 、EZplot3D(),以及兩個動態繪圖套件movie2D()、movie3D()。透過這幾個套件我們可以很快地看到我們的計算結果。使用範例如下:

>>> D = [0,3,0,2]

>>> Fx = [lambda x,y:0.1*np.sin(np.pi*x)*np.sin(np.pi*y/2), lambda x,y:0.0*(x+y), lambda y,t:0.0*(y+t), lambda y,t:0.0*(y+t), lambda x,t:0.0*(x+t), lambda x,t:0.0*(x+t)]

>>> u,X,Y,t = ode_wave2(0.25,D,2,Fx,40,40,40)

>>> EZplot3D(X,Y,u[0,:,:])

 

我們可以很容易拉動圖形,來觀察細部的變化。若要顯示動態的3D 動畫,我們使用

>>> movie3D(u,X,Y)