數值偏微分方程

 

Elliptic PDE

這一節我們只討論一個特別形式的 elliptic equation 稱為Helmholtz's equation,它的形式如下 \(\frac{\partial^2 u(x,y)}{\partial x^2} + \frac{\partial^2 u(x,y)}{\partial y^2} + g(x,y)u(x,y) = f(x,y)\) 在Domain \(D=\{ (x,y) | x_0 \leq x \leq x_f, y_0 \leq y \leq y_f \}\) 並且有下列邊界條件 \(u(x,y_0) = b_{y_0}(x)\), \(u(x,y_f) = b_{y_f}(x)\),\(u(x_0,y) = b_{x_0}(y)\) 以及 \(u(x_f,y) = b_{x_f}(y)\) 若 \(g(x,y) = f(x,y) = 0\),我們稱這個方程式為 Laplace's equation。 為了用差分的方法,我們把 \(x\) 方向的定義域切割成 \(M_x\) 等分,每一段的間距是 \(\Delta x = \frac{x_f - x_0}{M_x}\),同樣地,我們等分 \(y\) 方向的定義域為 \(M_y\) 等分,每一段的間距是 \(\Delta y = \frac{y_f - y_0}{M_y}\)。

我們用下列的差分式來代替微分值

\( \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1})}{\Delta x^2} \simeq \frac{\partial^2 u(x,y)}{\partial x^2}|_{xj,yi}\) \( \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta y^2} \simeq \frac{\partial^2 u(x,y)}{\partial y^2}|_{xj,yi}\)

這裡的 \(x_j = x_0 + j \Delta x\), \(y_i = y{0} + i \Delta y\),\(1\leq i\leq M_y-1\) 以及 \(1\leq j\leq M_x-1\)。 我們得到下列 finite difference equation

\(\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta x^2} + \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta y^2} + g_{i,j} u_{i,j} = f_{i,j}\)

其中 \(u_{i,j} = u(x_j,y_i)\),\(f_{i,j} = f(x_j,y_i)\) 以及 \(g_{i,j} = g(x_j,y_i)\) 這是一個 \((M_x-1)(M_y-1)\) 個未知數和方程式構成的線性方程組,可想而知,當 \(M_x\) 與 \(M_y\) 稍微放大一些,我們就糊了。這和先前三對角矩陣的情況不同,我們可以想像,不論我們如何把所有未知數排序,對應的矩陣不是一個三對角或幾對角類型的 sparse 矩陣。如果電腦硬體夠強,我們可以試著把矩陣寫開,然後真的使用軟體或是 linear solver套件去解它。如果不是呢?我們就要借用疊代法的技巧,不真的把矩陣寫開來計算。

要寫出一個疊代式我們把上式二次差分中的 \(u_{i,j}\) 保留在等號的左邊,其餘項移到等號右邊,我們得到

\(u_{i,j} = r_y (u_{i,j+1}+u_{i,j-1}) + r_x(u_{i+1,j}+u_{i-1,j}) + r_{xy}(g_{i,j}u_{i,j} - f_{i,j})\)

其中

\(r_y = \frac{\Delta y^2}{2(\Delta x^2+\Delta y^2)}\),\(r_x = \frac{\Delta x^2}{2(\Delta x^2+\Delta y^2)}\), \(r_{xy} = \frac{\Delta x^2 \Delta y^2}{2(\Delta x^2+\Delta y^2)}\), \(u_{i,0} = b_{x_0}(y_i)\), \(u_{i,Mx} = b_{x_f}(y_i)\), \(u_{0,j} = b_{y_0}(x_j)\) 以及 \(u_{My,j} = b_{y_f}(x_j)\)。

和解ODE的遞迴式一樣,他們的寫法可以有很多種變形,重點是我們會把它看做是

\(x^{n+1} = A x^n + b = A^n x^0 + (A^{n-1}+A^{n-2}+...+A)b\)

的遞迴形式。這樣的形式如何保證會收斂呢?答案當然是對應的 \(A\) 矩陣每一項 eigenvalue 都小於1 。如果我們知道所寫出來的遞迴式是一個會收斂的式子,再來的問題是起始值怎麼決定?因為我們有僅有的真解是邊界條件,把每一個 \(u_{i,j}\) 設成邊界條件的平均值似乎是一個合理的選擇。

若我們要用硬算的方式,可以用下列的Python 程式:

import numpy as np
import numpy.linalg as lg

def pde_Elliptic(f,g,BFx,Domain,Mx,My):
    """
    The matrix generator of the Elliptic PDE

    u_xx + u_yy + g(x,y)u = f(x,y) over the region
    D = [x0,xf]x[y0,yf] with the boundary conditions:
    u(x0,y) = bx0(y), u(xf,y) = bxf(y)
    u(x,y0) = by0(x), u(x,yf) = byf(x)
    Mx : the number of subintervals along x axis
    My : the number of subintervals along y axis

    BFx = [bx0,bxf,by0,byf]
    Domain = [x0,xf,y0,yf]
    
    The output A, b the matrix and vector of the corresponding linear system
    """

    n = (Mx-1)*(My-1) #The number of variables
    A = np.zeros([n,n])
    b = np.zeros([n,1])
    
    bx0, bxf, by0, byf = BFx[0], BFx[1], BFx[2], BFx[3]
    x0, xf, y0, yf = Domain[0], Domain[1], Domain[2], Domain[3]

    hx = float(xf-x0)/Mx
    hy = float(yf-y0)/My

    x = [x0+j*hx for j in range(1,Mx)]
    y = [y0+i*hy for i in range(1,My)]

    h2x = 1.0/(hx**2)
    h2y = 1.0/(hy**2)
    
    def ij2k(i,j): #Transfer (i,j) 2D-index to 1D index
        return i*(Mx-1)+j
    l = -1
    for i in range(My-1):
        for j in range(Mx-1):
            l += 1 
            k = ij2k(i,j)
            b[k] += f(x[j],y[i])
            #set the center coefficient
            A[l][k] = g(x[j],y[i])-2*(h2x+h2y)
            #set the up coefficient
            if i==0:
                b[k] -= h2y*by0(x[j])
            else:
                ku = ij2k(i-1,j)
                A[l][ku] += h2y
            #set the low coefficient
            if i==(My-2):
                b[k] -= h2y*byf(x[j])
            else:
                kd = ij2k(i+1,j)
                A[l][kd] += h2y
            #set the left coefficient
            if j==0:
                b[k] -= h2x*bx0(y[i])
            else:
                kl = ij2k(i,j-1)
                A[l][kl] += h2x
            #set the right coefficient
            if j==(Mx-2):
                b[k] -= h2x*bxf(y[i])
            else:
                kr = ij2k(i,j+1)
                A[l][kr] += h2x
    return A, b

 

若我們要解的PDE 為Laplace's equation,並且邊界條件為

\(b_{x_0}(y) = e^y-cos(y)\), \(b_{x_f}(y) = e^y cos(4)-exp(4)cos(y)\)

\(b_{y_0}(x) = cos(x) -e^x\), \(b_{y_f}(x) = e^4 cos(x)-e^x cos(4)\)

且Domain 為 \([0,4]\times [0,4]\)

則我們用下列Python 程式求解。

>>> BFx = [lambda y:np.exp(y)-np.cos(y), lambda y:np.exp(y)*np.cos(4)-np.exp(4)*np.cos(y),
.... lambda x:np.cos(x)-np.exp(x), lambda x: np.exp(4)*np.cos(x) - np.exp(x)*np.cos(4)]
>>> Domain = [0,4,0,4]
>>> A,b = pde_Elliptic(lambda x,y:0.0, lambda x,y:0.0, BFx,Domain,20,20)
>>> A
array([[-100., 25., 0., ..., 0., 0., 0.],
[ 25., -100., 25., ..., 0., 0., 0.],
[ 0., 25., -100., ..., 0., 0., 0.],
...,
[ 0., 0., 0., ..., -100., 25., 0.],
[ 0., 0., 0., ..., 25., -100., 25.],
[ 0., 0., 0., ..., 0., 25., -100.]])
>>> A.shape
(361, 361)
>>> len(b)
361
>>> x = lg.solve(A,b)
>>> X = x.reshape([19,19])
>>> X.shape
(19, 19)
>>> A,b = pde_Elliptic(lambda x,y:0.0, lambda x,y:0.0, BFx,Domain,100,100)
>>> A.shape
(9801, 9801)
>>> x = lg.solve(A,b)

 

這裡的 \(x\) 是解出來的向量,我們用 reshape() 函式把他排成矩陣的形式,方便我們繪圖。若我們把 \(M_x\) 和 \(M_y\) 增加到100,我們看到矩陣的大小變成9801x9801。最後一步解線性方程,以我手上的電腦大約花了一分鐘的時間(雖然無法忍受,但相較於Matlab 需要六分鐘的時間,也證明Python 比Matlab 快多了),我們可以想像,若網格點變得更細緻時,一般地電腦基本上是無法負荷的,因此我們需要用到遞迴解的方式來代替展開矩陣。

以下是遞迴解的Python 程式:

def pde_poisson(f,g,BFx,Domain,Mx,My,tol=10**(-8),MaxIter=100):
    """
    The PDE solver of u_xx + u_yy + g(x,y)u = f(x,y) over the region
    D = [x0,xf]x[y0,yf] with the boundary conditions:
    u(x0,y) = bx0(y), u(xf,y) = bxf(y)
    u(x,y0) = by0(x), u(x,yf) = byf(x)
    Mx : the number of subintervals along x axis
    My : the number of subintervals along y axis

    tol : the tolerance
    MaxIter : the Maximal number of the iteration
    
    The output
    u : u(x_j,y_i)
    x : the uniform grids of x-axis
    y : the uniform grids of y-axis 
    """

    bx0, bxf, by0, byf = BFx[0], BFx[1], BFx[2], BFx[3]
    x0, xf, y0, yf = Domain[0], Domain[1], Domain[2], Domain[3]
    
    hx = float(xf-x0)/Mx
    hy = float(yf-y0)/My

    x = [x0+j*hx for j in range(1,Mx)]
    y = [y0+i*hy for i in range(1,My)]

    Mx, My = Mx+1, My+1

    u = np.zeros([My,Mx])
    F = np.zeros([My,Mx])
    G = np.zeros([My,Mx])
    u0 = np.zeros([My,Mx])

    #set boundary condition
    j = 1
    for xj in x:
        u[0,j], u[My-1,j] = by0(xj), byf(xj)
        j+=1
    i = 1
    for yi in y:
        u[i,0], u[i,Mx-1] = bx0(yi), bxf(yi)
        i+=1

    #initialize as the average of boundary values
    sum_of_bv = sum(u[0,:])+sum(u[My-1,:])+sum(u[1:My-1,0])+sum(u[1:My-1,Mx-1])

    u[1:My-1,1:Mx-1] = float(sum_of_bv)/(2*(Mx+My-2))

    #set the f(xj,yi) & g(xj,yi)
    for i in range(1,My-1):
        for j in range(1,Mx-1):
            F[i,j], G[i,j] = f(x[j-1],y[i-1]), g(x[j-1],y[i-1])

    dx2, dy2 = hx**2, hy**2
    dxy2 = 2*(dx2+dy2)
    rx, ry = dx2/dxy2, dy2/dxy2
    rxy = rx*dy2
    for itr in range(MaxIter):
        for i in range(1,My-1):
            for j in range(1,Mx-1):
                u[i,j] = ry*(u[i,j+1]+u[i,j-1])+rx*(u[i+1,j]+u[i-1,j])+rxy*(G[i,j]*u[i,j]-F[i,j])
        Err = abs(u-u0)
        
        if (itr>1) & (Err.max()<tol):
            break
        u0 = u

    u = u[1:My,1:Mx]
    return u, x, y

 

我們用下列程式求先前的Laplace's equation

>>> BFx = [lambda y:np.exp(y)-np.cos(y), lambda y:np.exp(y)*np.cos(4)-np.exp(4)*np.cos(y),
.... lambda x:np.cos(x)-np.exp(x), lambda x: np.exp(4)*np.cos(x) - np.exp(x)*np.cos(4)]
>>> Domain = [0,4,0,4]
>>> u,x,y = pde_poisson(lambda x,y:0.0, lambda x,y:0.0, BFx,Domain,20,20,10**(-8),50)
>>> u.shape
(20, 20)

 

有了這兩種算法,有興趣的讀者可以進一步比較這兩種算法所計算的結果有哪些差異。先前我們談到,使用遞迴解的方式,要保證收斂所對應的線性系統的矩陣的特徵值的絕對值要比1小。接下來我們要驗證剛才的遞迴式所對應的矩陣特徵值是否真的如此。

產生這個對應矩陣的Python 程式如下:

def pde_poisson_MG(f,g,BFx,Domain,Mx,My):
    """
    The Matrix generator of the iteration method to solve u_xx + u_yy + g(x,y)u = f(x,y) over the region
    D = [x0,xf]x[y0,yf] with the boundary conditions:
    u(x0,y) = bx0(y), u(xf,y) = bxf(y)
    u(x,y0) = by0(x), u(x,yf) = byf(x)
    Mx : the number of subintervals along x axis
    My : the number of subintervals along y axis
    """

    bx0, bxf, by0, byf = BFx[0], BFx[1], BFx[2], BFx[3]
    x0, xf, y0, yf = Domain[0], Domain[1], Domain[2], Domain[3]
    
    hx = float(xf-x0)/Mx
    hy = float(yf-y0)/My

    x = [x0+j*hx for j in range(1,Mx)]
    y = [y0+i*hy for i in range(1,My)]

    Mx, My = Mx+1, My+1

    u = np.zeros([My,Mx])
    F = np.zeros([My,Mx])
    G = np.zeros([My,Mx])
    u0 = np.zeros([My,Mx])

    #set boundary condition
    j = 1
    for xj in x:
        u[0,j], u[My-1,j] = by0(xj), byf(xj)
        j+=1
    i = 1
    for yi in y:
        u[i,0], u[i,Mx-1] = bx0(yi), bxf(yi)
        i+=1

    #initialize as the average of boundary values
    sum_of_bv = sum(u[0,:])+sum(u[My-1,:])+sum(u[1:My-1,0])+sum(u[1:My-1,Mx-1])

    u[1:My-1,1:Mx-1] = float(sum_of_bv)/(2*(Mx+My-2))

    #set the f(xj,yi) & g(xj,yi)
    for i in range(1,My-1):
        for j in range(1,Mx-1):
            F[i,j], G[i,j] = f(x[j-1],y[i-1]), g(x[j-1],y[i-1])

    dx2, dy2 = hx**2, hy**2
    dxy2 = 2*(dx2+dy2)
    rx, ry = dx2/dxy2, dy2/dxy2
    rxy = rx*dy2
    
    def ij2k(i,j): #Transfer (i,j) 2D-index to 1D index
        return (i-1)*(Mx-2)+(j-1)

    n = (Mx-2)*(My-2) #The number of variables
    A = np.zeros([n,n])

    l=-1
    for i in range(1,My-2):
        for j in range(1,Mx-2):
            l += 1 
            k = ij2k(i,j)
            #set the center coefficient
            A[l][k] = rxy*G[i,j]
            #set the up coefficient
            if i>0:
                ku = ij2k(i-1,j)
                A[l][ku] += rx
            #set the low coefficient
            if i<(My--1):
                kd = ij2k(i+1,j)
                A[l][kd] += rx
            #set the left coefficient
            if j>0:
                kl = ij2k(i,j-1)
                A[l][kl] += ry
            #set the right coefficient
            if j<(Mx-1):
                kr = ij2k(i,j+1)
                A[l][kr] += ry
    return A

 

我們用同樣的Laplace's equation 條件來計算這個矩陣。

>>> BFx = [lambda y:np.exp(y)-np.cos(y), lambda y:np.exp(y)*np.cos(4)-np.exp(4)*np.cos(y),
.... lambda x:np.cos(x)-np.exp(x), lambda x: np.exp(4)*np.cos(x) - np.exp(x)*np.cos(4)]
>>> Domain = [0,4,0,4]
>>> A = pde_poisson_MG(lambda x,y:0.0, lambda x,y:0.0, BFx,Domain,20,20)
>>> D,V = lg.eig(A)
>>> max(D)
(0.91231809995413671+0j)

瞧!最大的eigenvalue 是0.9 多,所以會收斂,只可惜這個特徵值是0.9,若能更靠近0 收斂速度會更快。不知道大家通了沒有,我們可以套用解ODE 的Relaxation Method 自己用調整 \(w\) 的方式設計出收斂速度更快的方法。