Python小白的數學建模課-11.偏微分方程數值解法

偏微分方程可以描述各種自然和工程現象， 是構建科學、工程學和其他領域的數學模型主要手段，

偏微分方程主要有三類：橢圓方程，拋物方程和雙曲方程，

本文采用有限差分法求解偏微分方程，通過案例講解一維平流方程、一維熱傳導方程、二維雙曲方程、二維拋物方程和二維橢圓方程等常見型別的偏微分方程的數值解法，給出了全部例程和運行結果，

歡迎關注『Python小白的數學建模課 @ Youcans』系列，每周持續更新，

1. 偏微分方程基本知識

微分方程是指含有未知函式及其導數的關系式，偏微分方程是包含未知函式的偏導數（偏微分）的微分方程，

偏微分方程可以描述各種自然和工程現象， 是構建科學、工程學和其他領域的數學模型主要手段， 科學和工程中的大多數實際問題都歸結為偏微分方程的定解問題，如：波傳播，流動和擴散，振動，固體力學，電磁學和量子力學，等等，

偏微分方程主要有三類：橢圓方程，拋物方程和雙曲方程，

• 雙曲方程描述變數以一定速度沿某個方向傳播，常用于描述振動與波動問題，
• 橢圓方程描述變數以一定深度沿所有方向傳播，常用于描述靜電場、引力場等穩態問題，
• 拋物方程描述變數沿下游傳播，常用于描述熱傳導和擴散等瞬態問題，

偏微分方程的定解問題通常很難求出決議解，只能通過數值計算方法對偏微分方程的近似求解，常用偏微分方程數值解法有：有限差分方法、有限元方法、有限體方法、共軛梯度法，等等，通常先對問題的求解區域進行網格剖分，然后將定解問題離散為代數方程組，求出在離散網格點上的近似值，

Python 語言求解偏微分方程的功能是比較弱的，主要有 Fipy, FEniCS 等有限元方法的工具包，另外還有機器學習工具如 Tensorflow 也可以進行偏微分方程的仿真模擬，但是，這些工具都不適合 Python 小白學習和使用，因此本篇采用比較簡單的有限差分法對 5種典型的偏微分方程進行編程，通過案例講解偏微分方程的數值解法，

2. 案例一：一維線性平流方程

2.1 一維線性平流方程的數學模型

平流程序是大氣運動中重要的程序，平流方程（Advection equation）描述某一物理量的平流作用而引起局地變化的物理程序，最簡單的形式是一維平流方程，

{ ? u ? t + v ? u ? x = 0 u ( x , 0 ) = F ( x ) 0 ≤ t ≤ t e x a < x < x b \begin{cases} \begin{aligned} &\frac{\partial u}{\partial t} + v \frac{\partial u}{\partial x}= 0\\ &u(x,0)=F(x)\\ &0 \leq t \leq t_e\\ &x_a<x<x_b \end{aligned} \end{cases} ????????????????t?u?+v?x?u?=0u(x,0)=F(x)0≤t≤te?xa?<x<xb???

式中 u 為某物理量，v 為系統速度，x 為水平方向分量，t 為時間，

雖然該方程可以求得決議解：
u ( x , t ) = F ( x ? v ? t ) , 0 ≤ t ≤ t e u(x,t)=F(x-v*t), \ 0 \leq t \leq t_e u(x,t)=F(x?v?t), 0≤t≤te?

考慮一維線性平流偏微分方程的數值解法，采用有限差分法求解，簡單地， 采用一階迎風格式的差分方法（First-order Upwind)，一階導數的差分運算式為：
( ? u ? x ) i , j = u i + 1 , j ? u i , j Δ x + O ( Δ x ) (\frac{\partial u}{\partial x})_{i,j}=\frac{u_{i+1,j}-u_{i,j}}{\Delta x}+O(\Delta x) (?x?u?)i,j?=Δxui+1,j??ui,j??+O(Δx)
于是得到差分方程：
u i , j + 1 = u i , j ? v ? d t / d x ? ( u i , j ? u i ? 1 , j ) u_{i,j+1}=u_{i,j}-v*dt/dx*(u_{i,j}-u_{i-1,j}) ui,j+1?=ui,j??v?dt/dx?(ui,j??ui?1,j?)
即可遞推求得該平流方程的數值解，

2.2 偏微分方程編程步驟

以該題為例一類有限差分法求解一維平流問題偏微分方程的步驟：

1. 匯入numpy、matplotlib 包；
2. 定義初始條件函式 U(x,0)；
3. 輸入模型引數 v, p，定義求解的時間域 (tStart, tEnd) 和空間域 (xMin, xMax)，設定差分步長 dt, nNodes；
4. 初始化；
5. 遞推求解差分方程在區間 [xa,xb] 的數值解，獲得網格節點的處的 u(t) 值，

在例程中，設初值條件為 F ( x ) = s i n ( 2 ? ( x ? p ) 2 ) F(x) = sin(2 * (x-p)^2) F(x)=sin(2?(x?p)2)，取 v = 1.0 , p = 0.25 v= 1.0, p=0.25 v=1.0,p=0.25，空間域 x ∈ ( 0 , π ) x\in(0,\pi) x∈(0,π)，

2.3 Python 例程：偏微分方程（一維平流方程）

# mathmodel13_v1.py
# Demo10 of mathematical modeling algorithm
# Solving partial differential equations
# 偏微分方程數值解法

import numpy as np
import matplotlib.pyplot as plt

# 1. 一維平流方程 (advection equation)
# U_t + v*U_x = 0

# 初始條件函式 U(x,0)
def funcUx0(x, p): 
    u0 = np.sin(2 * (x-p)**2)
    return u0

# 輸入引數
v1 = 1.0  # 平流方程引數，系統速度
p = 0.25  # 初始條件函式 u(x,0) 中的引數

tc = 0  # 開始時間
te = 1.0  # 終止時間: (0, te)
xa = 0.0  # 空間范圍: (xa, xb)
xb = np.pi
dt = 0.02  # 時間差分步長
nNodes = 100  # 空間網格數

# 初始化
nsteps = round(te/dt)
dx = (xb - xa) / nNodes
x = np.arange(xa-dx, xb+2*dx, dx)
ux0 = funcUx0(x, p)

u = ux0.copy()  # u(j)
ujp = ux0.copy()  # u(j+1)

# 時域差分
for i in range(nsteps):
    plt.clf()  # 清除當前 figure 的所有axes, 但是保留當前視窗

    # 計算 u(j+1)
    for j in range(nNodes + 2):
        ujp[j] = u[j] - (v1 * dt/dx) * (u[j] - u[j-1])

    # 更新邊界條件
    u = ujp.copy()
    u[0] = u[nNodes + 1]
    u[nNodes+2] = u[1]

    # 繪圖
    plt.plot(x, u, 'b-', label="v1= 1.0")
    plt.axis((xa-0.1, xb + 0.1, -1.1, 1.1))
    plt.xlabel("x")
    plt.ylabel("U(x)")
    plt.legend(loc=(0.05,0.05))
    plt.title("Advection equation with finite difference method, t = %1.f" % (tc + dt))
    plt.text(0.05,0.9,"youcans-xupt",color='gainsboro')
    plt.pause(0.001)
    tc += dt

plt.show()

2.4 Python 例程運行結果

3. 案例二：一維熱傳導方程

3.1 一維熱傳導方程的數學模型

熱傳導方程描述一個區域內的溫度隨時間的變化，是典型的拋物型偏微分方程，也稱為擴散方程，

一維熱傳導方程考慮各向同性的均勻細桿，在垂直于 x 軸的截面上的溫度相同，細桿的圓周與周圍環境無熱交換，桿內無熱源，則溫度 u = u ( t , x ) u=u(t,x) u=u(t,x) 是時間變數 t 和水平方向空間變數 x 的函式，

{ ? u ? t = d i v ( U u ) = λ ? 2 u ? x 2 u ( x , 0 ) = u 0 ( x ) u ( x a ) = u a ( t ) , u ( x b , t ) = u b ( t ) 0 ≤ t ≤ t e , x a < x < x b \begin{cases} \begin{aligned} &\frac{\partial u}{\partial t} = div(U_u) = \lambda \frac{\partial ^2u}{\partial x^2}\\ &u(x,0)=u_0(x)\\ &u(x_a)=u_a(t),\ u(x_b,t)=u_b(t)\\ &0\leq t \leq t_e, \ x_a < x < x_b \end{aligned} \end{cases} ????????????????t?u?=div(Uu?)=λ?x2?2u?u(x,0)=u0?(x)u(xa?)=ua?(t), u(xb?,t)=ub?(t)0≤t≤te?, xa?<x<xb???

式中 λ \lambda λ 為熱擴散率，取決于材料本身的熱傳導率、密度和熱容，

考慮一維熱傳導偏微分方程的數值解法，采用有限差分法求解，簡單地， 采用迎風法的三點差分格式，二階導數的差分運算式為：
( ? 2 u ? x 2 ) i , j = u i + 1 , j ? 2 u i , j + u i ? 1 , j ( Δ x ) 2 + O ( Δ x ) 2 (\frac{\partial ^2u}{\partial x^2})_{i,j} = \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}+O(\Delta x)^2 (?x2?2u?)i,j?=(Δx)2ui+1,j??2ui,j?+ui?1,j??+O(Δx)2
于是得到差分方程：
u i , j + 1 = u i , j + λ d t / d x 2 ? ( u i + 1 , j ? 2 u i , j + u i ? 1 , j ) u_{i,j+1} = u_{i,j} + \lambda dt/dx^2*(u_{i+1,j}-2u_{i,j}+u_{i-1,j}) ui,j+1?=ui,j?+λdt/dx2?(ui+1,j??2ui,j?+ui?1,j?)
即可遞推求得一維熱傳導方程的數值解，

3.2 偏微分方程編程步驟

以該題為例一類有限差分法求解一維平流問題偏微分方程的步驟：

1. 匯入numpy、matplotlib 包；
2. 輸入模型引數 L, lam，定義求解的時間域 (0, te) 和空間域 (0, L)，設定差分步長 dt, dx；
3. 初始化；
4. 計算每一時點的邊界條件 U[0,j]，U[L,j]、每一位置的初始條件U[i,0]；
5. 遞推求解差分方程在空間域 [0, L] 的數值解，獲得網格節點的處的 U(x,t) ，

在例程中，設初值條件為 u ( x , t = 0 ) = 0 u(x,t=0) = 0 u(x,t=0)=0，邊界條件為 u ( x a , t ) = F ( t ) , u ( x b , t ) = 0 u(x_a,t)=F(t), \ u(x_b,t)=0 u(xa?,t)=F(t), u(xb?,t)=0，在時間域 t ∈ ( 0 , 2.0 ) t\in(0,2.0) t∈(0,2.0)、空間域 x ∈ ( 0 , 1.0 ) x\in(0,1.0) x∈(0,1.0) 求數值解即溫度分布，
（1）例程中的初始條件 U[i,0] 為常數，如果初始條件是 x 的函式 u(x,0)，將該函式在初始條件陳述句賦值即可（參加例程中注釋的陳述句），
（2）例程中的邊界條件，一端是 t 的函式 u(0,t)，另一端是常數 u(L,t) =0，這些條件也可以根據具體問題設定為相應的常數或函式，

3.3 Python 例程：偏微分方程（一維熱傳導方程）

# mathmodel13_v1.py
# Demo10 of mathematical modeling algorithm
# Solving partial differential equations
# 偏微分方程數值解法

import numpy as np
import matplotlib.pyplot as plt

# 2. 一維熱傳導方程（拋物型偏微分方程）
# pu/pt = l*p2u/px2

# 模型引數
L = 1.0  # 細桿長度
lam = 1.0  # 熱擴散率
tc = 0  # 開始時間
te = 10.0  # 終止時間: (0, te)

# 初始化
dx = 0.05  # 空間步長
dt = 0.001  # 時間步長
nNodes = round(L/dx)  # 空間網格數
nSteps = round(te/dt)  # 時序網格數
K = lam * dt/(dx**2)  # lambda * dt/dx^2
U = np.zeros([nNodes+1, nSteps+1])  # 建立二維陣列

# 邊界條件
for j in np.arange(0, nSteps+1):  # 時間序列
    U[0,j] = 7.5 + (nSteps-j)/2000 * np.sin(j/1000)/(1+np.exp(-j))
    U[nNodes,j] = 0.0  # 每一時點的邊界條件

# 初始條件
for i in np.arange(0, nNodes):  # 空間序列
    # U[i,0]= 0.2*i*h*(L-i*h)  # 初始條件是 x 的函式
    U[i,0]= 0  # 每一位置的初始條件

# 時域差分法求解
for j in np.arange(0, nSteps):  # 時間步長
    for i in np.arange(1, nNodes):  # 空間步長
        U[i,j+1] = K*U[i+1,j] + (1-2*K)*U[i,j] + K*U[i-1,j]

# 繪圖
xZone = np.arange(0, (nNodes+1)*dx, dx)  # 建立空間網格
tZone = np.arange(0, (nSteps+1)*dt, dt)  # 建立空間網格
fig = plt.figure(figsize=(10, 6))
rect1 = [0.05, 0.2, 0.4, 0.65]  # [左, 下, 寬, 高], 0.0~1.0
ax1 = plt.axes(rect1)
for k in range(0,nSteps+1,round(nSteps/5)):
    ax1.plot(xZone, U[:,k], label=r"x={}".format(k/nSteps))
ax1.set_ylabel('u(x,t)')
ax1.set_xlabel('x')
ax1.set_xlim(0,L)
ax1.set_ylim(-1,12)
ax1.set_title("Temperature distribution along t-axis")
ax1.legend(loc='upper right')
rect2 = [0.55, 0.2, 0.4, 0.65]  # [左, 下, 寬, 高], 0.0~1.0
ax2 = plt.axes(rect2)
for k in range(0,nNodes+1,round(nNodes/5)):  # U[nNodes,k] = 0.0
    ax2.plot(tZone, U[k,:], label=r"t={}".format(k/nNodes))
ax2.set_ylabel('u(x,t)')
ax2.set_xlabel('t')
ax2.set_xlim(0,te)
ax2.set_ylim(-1,12)
ax2.set_title("Temperature distribution along x-axis")
ax2.legend(loc='upper right')
plt.show()

3.4 Python 例程運行結果

4. 案例三：二維雙曲型方程

4.1 二維波動方程的數學模型

波動方程（wave equation）是典型的雙曲形偏微分方程，廣泛應用于聲學，電磁學，和流體力學等領域，描述自然界中的各種的波動現象，包括橫波和縱波，例如聲波、光波和水波，

考慮如下二維波動方程的初邊值問題：
{ ? 2 u ? t 2 = c 2 ( ? 2 u ? x 2 + ? 2 u ? y 2 ) ? ? t u ( 0 , x , y ) = 0 u ( x , y , 0 ) = u 0 ( x , y ) u ( 0 , y , t ) = u a ( t ) , u ( 1 , y , t ) = u b ( t ) u ( x , 0 , t ) = u c ( t ) , u ( x , 1 , t ) = u d ( t ) 0 ≤ t ≤ t e , 0 < x < 1 , 0 < y < 1 \begin{cases} \begin{aligned} &\frac{\partial^2 u}{\partial t^2} = c^2 (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2})\\ &\frac{\partial }{\partial t} u(0,x,y)= 0\\ &u(x,y,0)=u_0(x,y)\\ &u(0,y,t)=u_a(t),\ u(1,y,t)=u_b(t)\\ &u(x,0,t)=u_c(t),\ u(x,1,t)=u_d(t)\\ &0\leq t \leq t_e, \ 0 < x < 1, \ 0< y < 1 \end{aligned} \end{cases} ??????????????????????????????t2?2u?=c2(?x2?2u?+?y2?2u?)?t??u(0,x,y)=0u(x,y,0)=u0?(x,y)u(0,y,t)=ua?(t), u(1,y,t)=ub?(t)u(x,0,t)=uc?(t), u(x,1,t)=ud?(t)0≤t≤te?, 0<x<1, 0<y<1??

式中：u 是振幅；c 為波的傳播速率，c 可以是固定常數，或位置的函式 c(x,y)，也可以是振幅的函式 c(u)，

考慮二維波動偏微分方程的數值解法，采用有限差分法求解，簡單地， 采用迎風法的三點差分格式， 將上述的偏微分方程離散為差分方程 ：

化簡后得到：

即可遞推求得二維波動方程的數值解，

為了保證演算法的收斂性，迎風法的三點差分格式要求步長比小于 1：

r = 4 c 2 Δ t 2 Δ x 2 + Δ y 2 ≤ 1 r = \frac{4c^2\Delta t^2}{\Delta x^2 + \Delta y^2} \leq 1 r=Δx2+Δy24c2Δt2?≤1

4.2 偏微分方程編程步驟

以該題為例一類有限差分法求解二維波動問題偏微分方程的步驟：

1. 匯入numpy、matplotlib 包；
2. 輸入模型引數 c，定義求解的時間域 (0, te) 和空間域 (0,1)；
3. 初始化，設定差分步長 dt, dx, dy，檢驗步長比以保證演算法的穩定性；
4. 計算初值條件U[0]、U[1]；
5. 遞推求解差分方程在空間域 [0, 1]*[0, 1] 的數值解，獲得網格節點的處的 U(t,x,y) ，

在例程中，取初始條件為 u ( x , y , 0 ) = s i n ( 6 π x ) + c o s ( 4 π y ) u(x,y,0)=sin(6\pi x)+cos(4\pi y) u(x,y,0)=sin(6πx)+cos(4πy)，邊界條件為 u ( 0 , y , t ) = u ( 1 , y , t ) = 0 u(0,y,t)=u(1,y,t)=0 u(0,y,t)=u(1,y,t)=0， u ( x , 0 , t ) = u ( x , 1 , t ) = 0 u(x,0,t)=u(x,1,t)=0 u(x,0,t)=u(x,1,t)=0，在時間域 t ∈ ( 0 , 1.0 ) t\in(0,1.0) t∈(0,1.0)、空間域 x ∈ ( 0 , 1.0 ) , y ∈ ( 0 , 1.0 ) x\in(0,1.0),\ y\in(0,1.0) x∈(0,1.0), y∈(0,1.0) 求數值解即振動狀態，

4.3 Python 例程

# mathmodel13_v1.py
# Demo10 of mathematical modeling algorithm
# Solving partial differential equations
# 偏微分方程數值解法

# 4. 二維波動方程（雙曲型二階偏微分方程）
# p2u/pt2 = c^2*(p2u/px2+p2u/py2)

import numpy as np
import matplotlib.pyplot as plt

# 模型引數
c = 1.0  # 波的傳播速率
tc, te = 0.0, 1.0  # 時間范圍，0<t<te
xa, xb = 0.0, 1.0  # 空間范圍，xa<x<xb
ya, yb = 0.0, 1.0  # 空間范圍，ya<y<yb

# 初始化
c2 = c*c  # 方程引數
dt = 0.01  # 時間步長
dx = dy = 0.02  # 空間步長
tNodes = round(te/dt)  # t軸 時序網格數
xNodes = round((xb-xa)/dx)  # x軸 空間網格數
yNodes = round((yb-ya)/dy)  # y軸 空間網格數
tZone = np.arange(0, (tNodes+1)*dt, dt)  # 建立空間網格
xZone = np.arange(0, (xNodes+1)*dx, dx)  # 建立空間網格
yZone = np.arange(0, (yNodes+1)*dy, dy)  # 建立空間網格
xx, yy = np.meshgrid(xZone, yZone)  # 生成網格點的坐標 xx,yy (二維陣列)

# 步長比檢驗(r>1 則演算法不穩定)
r = 4 * c2 * dt*dt / (dx*dx+dy*dy)
print("dt = {:.2f}, dx = {:.2f}, dy = {:.2f}, r = {:.2f}".format(dt,dx,dy,r))
assert r < 1.0, "Error: r>1, unstable step ratio of dt2/(dx2+dy2) ."
rx = c*c * dt**2/dx**2
ry = c*c * dt**2/dy**2

# 繪圖
fig = plt.figure(figsize=(8,6))
ax1 = fig.add_subplot(111, projection='3d')
# ax2 = fig.add_subplot(122, projection='3d')

# 計算初始值
U = np.zeros([tNodes+1, xNodes+1, yNodes+1])  # 建立三維陣列
U[0] = np.sin(6*np.pi*xx)+np.cos(4*np.pi*yy)  # U[0,:,:]
U[1] = np.sin(6*np.pi*xx)+np.cos(4*np.pi*yy)  # U[1,:,:]
surf = ax1.plot_surface(xx, yy, U[0,:,:], rstride=2, cstride=2, cmap=plt.cm.coolwarm)
# wframe = ax2.plot_wireframe(xx, yy, U[0], rstride=2, cstride=2, linewidth=1)

# 有限差分法求解
for k in range(2,tNodes+1):
    if surf:
        ax1.collections.remove(surf)  # 更新三維影片視窗

    for i in range(1,xNodes):
        for j in range(1,yNodes):
            U[k,i,j] = rx*(U[k-1,i-1,j]+U[k-1,i+1,j]) + ry*(U[k-1,i,j-1]+U[k-1,i,j+1])\
                     + 2*(1-rx-ry)*U[k-1,i,j] -U[k-2,i,j]

    surf = ax1.plot_surface(xx, yy, U[k,:,:], rstride=2, cstride=2, cmap='rainbow')
    # wframe = ax2.plot_wireframe(xx, yy, U[k,:,:], rstride=2, cstride=2, linewidth=1, cmap='rainbow')
    ax1.set_xlim3d(0, 1.0)
    ax1.set_ylim3d(0, 1.0)
    ax1.set_zlim3d(-2, 2)
    ax1.set_title("2D wave equationt (t= %.2f)" % (k*dt))
    ax1.set_xlabel("x-youcans")
    ax1.set_ylabel("y-XUPT")
    plt.pause(0.01)

plt.show()

4.4 Python 例程運行結果

5. 案例四：二維拋物型方程

5.1 二維熱傳導方程的數學模型

熱傳導方程（heat equation）是典型的拋物形偏微分方程，也成為擴散方程，廣泛應用于聲學，電磁學，和流體力學等領域，描述自然界中的各種的波動現象，包括橫波和縱波，例如聲波、光波和水波，

考慮如下二維熱傳導方程的初邊值問題：
{ ? u ? t = λ ( ? 2 u ? x 2 + ? 2 u ? y 2 ) + q v u ( x , y , 0 ) = u 0 ( x , y ) u ( 0 , y , t ) = u a ( t ) , u ( 1 , y , t ) = u b ( t ) u ( x , 0 , t ) = u c ( t ) , u ( x , 1 , t ) = u d ( t ) 0 ≤ t ≤ t e , 0 < x < 1 , 0 < y < 1 \begin{cases} \begin{aligned} &\frac{\partial u}{\partial t} = \lambda (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) + q_v\\ &u(x,y,0)=u_0(x,y)\\ &u(0,y,t)=u_a(t),\ u(1,y,t)=u_b(t)\\ &u(x,0,t)=u_c(t),\ u(x,1,t)=u_d(t)\\ &0\leq t \leq t_e, \ 0 < x < 1, \ 0< y < 1 \end{aligned} \end{cases} ??????????????????????t?u?=λ(?x2?2u?+?y2?2u?)+qv?u(x,y,0)=u0?(x,y)u(0,y,t)=ua?(t), u(1,y,t)=ub?(t)u(x,0,t)=uc?(t), u(x,1,t)=ud?(t)0≤t≤te?, 0<x<1, 0<y<1??

式中 λ \lambda λ 為熱擴散率，取決于材料本身的熱傳導率、密度和熱容； q v q_v qv? 是熱源強度，可以是恒定值，也可以是時間、空間的函式，

考慮二維拋物型偏微分方程的數值解法，采用有限差分法求解，將上述的偏微分方程離散為差分方程 ：

化簡后得到：

即可遞推求得二維波動方程的數值解：
{ U k + 1 = U k + r x ? A ? U k + r y ? B ? U k + q v ? d t ) A = [ ? 2 1 ? 0 0 1 ? 2 1 ? 0 0 1 ? 2 ? 0 ? ? ? ? ? 0 0 ? 1 ? 2 ] ( N x + 1 , N x + 1 ) B = [ ? 2 1 ? 0 0 1 ? 2 1 ? 0 0 1 ? 2 ? 0 ? ? ? ? ? 0 0 ? 1 ? 2 ] ( N y + 1 , N y + 1 ) \begin{cases} \begin{aligned} &U_{k+1} = U_{k} + r_x*A*U_k + r_y*B*U_k + q_v*dt)\\ &A= \begin{bmatrix} -2 & 1 & \cdots & 0 & 0\\ 1 & -2 & 1 & \cdots & 0\\ 0 & 1 & -2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots& \vdots\\ 0 & 0 & \cdots & 1 & -2\\ \end{bmatrix} _{(Nx+1,Nx+1)} B= \begin{bmatrix} -2 & 1 & \cdots & 0 & 0\\ 1 & -2 & 1 & \cdots & 0\\ 0 & 1 & -2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots& \vdots\\ 0 & 0 & \cdots & 1 & -2\\ \end{bmatrix} _{(Ny+1,Ny+1)} \end{aligned} \end{cases} ?????????????????????Uk+1?=Uk?+rx??A?Uk?+ry??B?Uk?+qv??dt)A=?????????210?0?1?21?0??1?2???0???1?000??2?????????(Nx+1,Nx+1)?B=?????????210?0?1?21?0??1?2???0???1?000??2?????????(Ny+1,Ny+1)???

5.2 偏微分方程編程步驟

以該題為例一類有限差分法求解二維波動問題偏微分方程的步驟：

1. 匯入numpy、matplotlib 包；
2. 輸入熱傳導引數、熱源引數等模型引數，定義求解的時間域 (0, te) 和空間域；
3. 初始化，設定差分步長 dt, dx, dy，計算三對角系數矩陣 A、B，差分系數 rx, ry, ft；
4. 計算初始條件 U 0 U_0 U0?；
5. 遞推求解差分方程在空間域的數值解，更新邊界條件，獲得網格節點的處的 U ( x , y ) U(x,y) U(x,y)；
6. 繪制等溫云圖，

例程求解一個薄板受熱的溫度分布問題，其初始條件為 t I n i = 25 tIni =25 tIni=25，邊界條件為 t B o u n d = 25 tBound = 25 tBound=25，熱源為 Q v Qv Qv，在空間域 x ∈ ( 0 , 16 ) , y ∈ ( 0 , 12 ) x\in(0,16),\ y\in(0,12) x∈(0,16), y∈(0,12) ，時間域 t ∈ ( 0 , 5 ) t\in(0,5) t∈(0,5) 求數值解即溫度分布，

對于外加熱源，例程中給出了三種情況：（1）恒定熱源，（2）熱源功率（或開關）隨時間變化，（3）熱源位置隨時間變化（從 (x0,y0) 以速度 (xv,yv) 移動，以模擬焊接加熱的情況），

5.3 Python 例程

# mathmodel13_v1.py
# Demo10 of mathematical modeling algorithm
# Solving partial differential equations
# 偏微分方程數值解法

# 5. 二維熱傳導方程（拋物型二階偏微分方程）
# pu/p2 = c*(p2u/px2+p2u/py2)

import numpy as np
import matplotlib.pyplot as plt

def showcontourf(zMat, xyRange, tNow):  # 繪制等溫云圖
    x = np.linspace(xyRange[0], xyRange[1], zMat.shape[1])
    y = np.linspace(xyRange[2], xyRange[3], zMat.shape[0])
    xx,yy = np.meshgrid(x,y)
    zMax = np.max(zMat)
    yMax, xMax = np.where(zMat==zMax)[0][0], np.where(zMat==zMax)[1][0]
    levels = np.arange(0,100,1)
    showText = "time = {:.1f} s\nhotpoint = {:.1f} C".format(tNow, zMax)

    plt.clf()  # 清除當前圖形及其所有軸，但保持視窗打開
    plt.plot(x[xMax],y[yMax],'ro')  # 繪制最高溫度點
    plt.contourf(xx, yy, zMat, 100, cmap=plt.cm.get_cmap('jet'), origin='lower', levels = levels)
    plt.annotate(showText, xy=(x[xMax],y[yMax]), xytext=(x[xMax],y[yMax]),fontsize=10)
    plt.colorbar()
    plt.xlabel('Xupt')
    plt.ylabel('Youcans')
    plt.title('Temperature distribution of the plate')
    plt.draw()

# 模型引數
uIni = 25  # 初始溫度值
uBound = 25.0  # 邊界條件
c = 1.0  # 熱傳導引數
qv = 50.0  # 熱源功率
x0, y0 = 0.0, 3.0  # 熱源初始位置
vx, vy = 2.0, 1.0  # 熱源移動速度
# 求解范圍
tc, te = 0.0, 5.0  # 時間范圍，0<t<te
xa, xb = 0.0, 16.0  # 空間范圍，xa<x<xb
ya, yb = 0.0, 12.0  # 空間范圍，ya<y<yb

# 初始化
dt = 0.002  # 時間步長
dx = dy = 0.1  # 空間步長
tNodes = round(te/dt)  # t軸 時序網格數
xNodes = round((xb-xa)/dx)  # x軸 空間網格數
yNodes = round((yb-ya)/dy)  # y軸 空間網格數
xyRange = np.array([xa, xb, ya, yb])
xZone = np.linspace(xa, xb, xNodes+1)  # 建立空間網格
yZone = np.linspace(ya, yb, yNodes+1)  # 建立空間網格
xx,yy = np.meshgrid(xZone, yZone)  # 生成網格點的坐標 xx,yy (二維陣列)

# 計算 差分系數矩陣 A、B (三對角對稱矩陣)，差分系數 rx,ry,ft
A = (-2) * np.eye(xNodes+1, k=0) + (1) * np.eye(xNodes+1, k=-1) + (1) * np.eye(xNodes+1, k=1)
B = (-2) * np.eye(yNodes+1, k=0) + (1) * np.eye(yNodes+1, k=-1) + (1) * np.eye(yNodes+1, k=1)
rx, ry, ft = c*dt/(dx*dx), c*dt/(dy*dy), qv*dt

# 計算 初始值
U = uIni * np.ones((yNodes+1, xNodes+1))  # 初始溫度 u0

# 前向歐拉法一階差分求解
for k in range(tNodes+1):
    t = k * dt  # 當前時間

    # 熱源條件
    # (1) 恒定熱源：Qv(x0,y0,t) = qv, 功率、位置 恒定
    # Qv = qv
    # (2) 熱源功率隨時間變化 Qv(x0,y0,t)=f(t)
    # Qv = qv*np.sin(t*np.pi) if t<2.0 else qv
    # (3) 熱源位置隨時間變化 Qv(x,y,t)=f(x(t),y(t))
    xt, yt = x0+vx*t, y0+vy*t  # 熱源位置變化
    Qv = qv * np.exp(-((xx-xt)**2+(yy-yt)**2))  # 熱源方程

    # 邊界條件
    U[:,0] = U[:,-1] = uBound
    U[0,:] = U[-1,:] = uBound

    # 差分求解
    U = U + rx * np.dot(U,A) + ry * np.dot(B,U) + Qv*dt

    if k % 100 == 0:
        showcontourf(U, xyRange, k*dt)  # 繪制等溫云圖
        print('t={:.2f}s\tTmax={:.1f}  Tmin={:.1f}'.format(t, np.max(U), np.min(U)))

5.4 Python 例程運行結果

6. 案例五：二維橢圓型方程

6.1 二維橢圓方程的數學模型

橢圓型偏微分方程是一類重要的偏微分方程，主要用來描述物理的平衡穩定狀態，如定常狀態下的電磁場、引力場和反應擴散現象等，廣泛應用于流體力學、彈性力學、電磁學、幾何學和變分法中，
考慮如下二維泊松方程：
? 2 u ? x 2 + ? 2 u ? y 2 = f ( x , y ) \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)\\ ?x2?2u?+?y2?2u?=f(x,y)

上式在 f ( x , y ) = 0 f(x,y)=0 f(x,y)=0 時就是拉普拉斯方程（Laplace equation），

考慮二維橢圓型偏微分方程的數值解法，采用有限差分法求解，簡單地， 采用五點差分格式表示二階導數的差分運算式，將上述的偏微分方程離散為差分方程 ：
( u i ? 1 , j ? 2 u i , j + u i + 1 , j ) / Δ h 2 + ( u i , j ? 1 ? 2 u i , j + u i , j + 1 ) / Δ h 2 = f i , j (u_{i-1,j}-2u_{i,j}+u_{i+1,j})/\Delta h^2 + (u_{i,j-1}-2u_{i,j}+u_{i,j+1})/\Delta h^2 = f_{i,j} (ui?1,j??2ui,j?+ui+1,j?)/Δh2+(ui,j?1??2ui,j?+ui,j+1?)/Δh2=fi,j?
橢圓型偏微分描述不隨時間變化的均衡狀態，沒有初始條件，因此不能沿時間步長遞推求解，對上式的差分方程，可以通過矩陣求逆方法求解，但當 h 較小時網格很多，矩陣求逆的記憶體占用和計算量極大，于是，可以使用迭代松弛法遞推求得二維橢圓方程的數值解：
u i , j k + 1 = ( 1 ? w ) ? u i , j k + w / 4 ? ( u i , j + 1 k + u i , j ? 1 k + u i ? 1 , j k + u i + 1 , j k ? h 2 ? f i , j ) u_{i,j}^{k+1} = (1-w)*u_{i,j}^k + w/4*(u_{i,j+1}^k + u_{i,j-1}^k + u_{i-1,j}^k + u_{i+1,j}^k -h^2* f_{i,j}) ui,jk+1?=(1?w)?ui,jk?+w/4?(ui,j+1k?+ui,j?1k?+ui?1,jk?+ui+1,jk??h2?fi,j?)

6.2 偏微分方程編程步驟

以該題為例一類有限差分法求解二維波動問題偏微分方程的步驟：

1. 匯入numpy、matplotlib 包；
2. 輸入模型引數，定義空間域 (0,1)；
3. 初始化，設定差分步長 h、松弛因子 w；
4. 計算邊值條件 u[0,:], u[1,:], u[:,0], u[:,1]；
5. 迭代松弛法遞推求解差分方程在空間域 [0, 1]*[0, 1] 的數值解 ，

在例程中，取邊界條件為 u ( x , 0 ) = u ( x , 1 ) = 0 , ?? u ( 0 , y ) = u ( 1 , y ) = 1 u(x,0)=u(x,1)=0, \; u(0,y)=u(1,y)=1 u(x,0)=u(x,1)=0,u(0,y)=u(1,y)=1，為了讓圖形更有趣，也可以參考例程中選擇不同的邊界條件，

6.3 Python 例程

# mathmodel13_v1.py
# Demo10 of mathematical modeling algorithm
# Solving partial differential equations
# 偏微分方程數值解法

# 7. 二維橢圓方程(Laplace equation)
# u_xx+u_yy=s

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 求解范圍
xa, xb = 0.0, 1.0  # 空間范圍，xa<x<xb
ya, yb = 0.0, 1.0  # 空間范圍，ya<y<yb

# 初始化
h = 0.01  # 空間步長, dx = dy = 0.01
w = 0.5  # 松弛因子
nodes = round((xb-xa)/h)  # x軸 空間網格數

# 邊值條件
u = np.zeros((nodes+1, nodes+1))
# u[:, 0] = 0
# u[:, -1] = 0  # -1 表示陣列的最后一個值
# u[0, :] = 1
# u[-1, :] = 1  # -1 表示陣列的最后一個值
for i in range(nodes+1):
    u[i, 0] = 1.0 + np.sin(0.5*(i-50)/np.pi)
    u[i, -1] = -1.0 + 0.5*np.sin((i-50)/np.pi)
    u[0, i] = -1.0 + 0.5*np.sin((i-50)/np.pi)
    u[-1, i] = 1.0 + np.sin(0.5*(50-i)/np.pi)

# 迭代松弛法求解
for iter in range(100):
    for i in range(1, nodes):
        for j in range(1, nodes):
            u[i, j] = w/4 * (u[i-1, j] + u[i+1, j] + u[i, j-1] + u[i, j+1]) + (1-w) * u[i, j]

# 繪圖
x = np.linspace(0, 1, nodes+1)
y = np.linspace(0, 1, nodes+1)
xx, yy = np.meshgrid(x, y)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(xx, yy, u, cmap=plt.get_cmap('rainbow'))
fig.colorbar(surf, shrink=0.5)
ax.set_xlim3d(0, 1.0)
ax.set_ylim3d(0, 1.0)
ax.set_zlim3d(-2, 2.5)
ax.set_title("2D elliptic partial differential equation")
ax.set_xlabel("Xupt")
ax.set_ylabel("Youcans")
plt.show()

6.4 Python 例程運行結果

7. 小結

1. 偏微分方程是數學中的重要學科分支，其數值解法的研究和方法可謂博大精深，遠非本文所能涵蓋，
2. 本文針對 Python小白，采用有限差分法求解偏微分方程，通過案例介紹了一維平流方程、一維熱傳導方程、二維雙曲方程、二維拋物方程和二維橢圓方程的數值解法，給出了例程和運行結果，供 Python 小白參考，以便進行數模學習，
3. 需要特別說明的是：一階差分格式的精度相對較低，往往需要較多網格計算，隨著精度的不斷提高，可以推匯出各種差分運算式，高階精度的差分公式可以給出質量更高的解，但需要使用更多的網格點進行計算，程式更復雜，計算量很大，類似地，顯式方法的建立和編程簡單，但要求步長較小才能保持穩定性，隱式方法的建立和編程更復雜，運算量大，但穩定性好，
4. 需要特別說明的是：關于偏微分方程數值解法的穩定性、收斂性和誤差分析，是非常專業的問題，例如步長選擇不恰當可能導致演算法不穩定，如 4.2 中應滿足步長比 r>1，除非找到專業課程教材或范文中有相關內容可以參考套用，否則不建議小白自己摸索，這些問題不是調整引數試試就能試出來的，

【本節完】

著作權宣告：

歡迎關注『Python小白的數學建模課 @ Youcans』 原創作品

原創作品，轉載必須標注原文鏈接：https://blog.csdn.net/youcans/article/details/119755450

Copyright 2021 Youcans, XUPT

Crated：2021-08-16

歡迎關注 『Python小白的數學建模課 @ Youcans』 系列，持續更新
Python小白的數學建模課-01.新手必讀
Python小白的數學建模課-02.資料匯入
Python小白的數學建模課-03.線性規劃
Python小白的數學建模課-04.整數規劃
Python小白的數學建模課-05.0-1規劃
Python小白的數學建模課-06.固定費用問題
Python小白的數學建模課-07.選址問題
Python小白的數學建模課-09.微分方程模型
Python小白的數學建模課-10.微分方程邊值問題
Python小白的數學建模課-11.偏微分方程數值解法
Python小白的數學建模課-12.非線性規劃
Python小白的數學建模課-15.圖論的基本概念
Python小白的數學建模課-16.最短路徑演算法
Python小白的數學建模課-17.條件最短路徑演算法
Python小白的數學建模課-18.最小生成樹問題
Python小白的數學建模課-19.網路流優化問題
Python小白的數學建模課-20.網路流優化案例
Python小白的數學建模課-21.關鍵路徑法
Python小白的數學建模課-22.插值方法
Python小白的數學建模課-23.資料擬合全集
Python小白的數學建模課-A1.國賽賽題型別分析
Python小白的數學建模課-A2.2021年數維杯C題探討
Python小白的數學建模課-A3.12個新冠疫情數模競賽賽題及短評
Python小白的數學建模課-B2. 新冠疫情 SI模型
Python小白的數學建模課-B3. 新冠疫情 SIS模型
Python小白的數學建模課-B4. 新冠疫情 SIR模型
Python小白的數學建模課-B5. 新冠疫情 SEIR模型
Python小白的數學建模課-B6. 新冠疫情 SEIR改進模型
Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計回歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火演算法