Python小白的数学建模课-11.偏微分方程数值解法

偏微分方程可以描述各种自然和工程现象，是构建科学、工程学和其他领域的数学模型主要手段。

偏微分方程主要有三类：椭圆方程，抛物方程和双曲方程。

本文采用有限差分法求解偏微分方程，通过案例讲解一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程等常见类型的偏微分方程的数值解法，给出了全部例程和运行结果。

欢迎关注『Python小白的数学建模课 @ Youcans』系列，每周持续更新。

文章目录

1. 偏微分方程基本知识
2. 案例一：一维线性平流方程
- 2.1 一维线性平流方程的数学模型
- 2.2 偏微分方程编程步骤
- 2.3 Python 例程：偏微分方程（一维平流方程）
- 2.4 Python 例程运行结果
3. 案例二：一维热传导方程
- 3.1 一维热传导方程的数学模型
- 3.2 偏微分方程编程步骤
- 3.3 Python 例程：偏微分方程（一维热传导方程）
- 3.4 Python 例程运行结果
4. 案例三：二维双曲型方程
- 4.1 二维波动方程的数学模型
- 4.2 偏微分方程编程步骤
- 4.3 Python 例程
- 4.4 Python 例程运行结果
5. 案例四：二维抛物型方程
- 5.1 二维热传导方程的数学模型
- 5.2 偏微分方程编程步骤
- 5.3 Python 例程
- 5.4 Python 例程运行结果
6. 案例五：二维椭圆型方程
- 6.1 二维椭圆方程的数学模型
- 6.2 偏微分方程编程步骤
- 6.3 Python 例程
- 6.4 Python 例程运行结果
7. 小结

1. 偏微分方程基本知识

微分方程是指含有未知函数及其导数的关系式，偏微分方程是包含未知函数的偏导数（偏微分）的微分方程。

偏微分方程可以描述各种自然和工程现象，是构建科学、工程学和其他领域的数学模型主要手段。科学和工程中的大多数实际问题都归结为偏微分方程的定解问题，如：波传播，流动和扩散，振动，固体力学，电磁学和量子力学，等等。

偏微分方程主要有三类：椭圆方程，抛物方程和双曲方程。

双曲方程描述变量以一定速度沿某个方向传播，常用于描述振动与波动问题。
椭圆方程描述变量以一定深度沿所有方向传播，常用于描述静电场、引力场等稳态问题。
抛物方程描述变量沿下游传播，常用于描述热传导和扩散等瞬态问题。

偏微分方程的定解问题通常很难求出解析解，只能通过数值计算方法对偏微分方程的近似求解。常用偏微分方程数值解法有：有限差分方法、有限元方法、有限体方法、共轭梯度法，等等。通常先对问题的求解区域进行网格剖分，然后将定解问题离散为代数方程组，求出在离散网格点上的近似值。

Python 语言求解偏微分方程的功能是比较弱的，主要有 Fipy, FEniCS 等有限元方法的工具包，另外还有机器学习工具如 Tensorflow 也可以进行偏微分方程的仿真模拟。但是，这些工具都不适合 Python 小白学习和使用，因此本篇采用比较简单的有限差分法对 5种典型的偏微分方程进行编程，通过案例讲解偏微分方程的数值解法。

2. 案例一：一维线性平流方程

2.1 一维线性平流方程的数学模型

平流过程是大气运动中重要的过程。平流方程（Advection equation）描述某一物理量的平流作用而引起局地变化的物理过程，最简单的形式是一维平流方程。

{∂u∂t+v∂u∂x=0u(x,0)=F(x)0≤t≤texa<x<xb\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≤texa<x<xb

式中 u 为某物理量，v 为系统速度，x 为水平方向分量，t 为时间。

虽然该方程可以求得解析解：
u(x,t)=F(x−v∗t),0≤t≤teu(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=ui+1,j−ui,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)
于是得到差分方程：
ui,j+1=ui,j−v∗dt/dx∗(ui,j−ui−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 偏微分方程编程步骤

以该题为例一类有限差分法求解一维平流问题偏微分方程的步骤：

导入numpy、matplotlib 包；
定义初始条件函数 U(x,0)；
输入模型参数 v, p，定义求解的时间域 (tStart, tEnd) 和空间域 (xMin, xMax)，设置差分步长 dt, nNodes；
初始化；
递推求解差分方程在区间 [xa,xb] 的数值解，获得网格节点的处的 u(t) 值。

在例程中，设初值条件为 F(x)=sin(2∗(x−p)2)F(x) = sin(2 * (x-p)^2)F(x)=sin(2∗(x−p)2)，取 v=1.0,p=0.25v= 1.0, p=0.25v=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 += dtplt.show()

2.4 Python 例程运行结果

3. 案例二：一维热传导方程

3.1 一维热传导方程的数学模型

热传导方程描述一个区域内的温度随时间的变化，是典型的抛物型偏微分方程，也称为扩散方程。

一维热传导方程考虑各向同性的均匀细杆，在垂直于 x 轴的截面上的温度相同，细杆的圆周与周围环境无热交换，杆内无热源，则温度 u=u(t,x)u=u(t,x)u=u(t,x) 是时间变量 t 和水平方向空间变量 x 的函数。

{∂u∂t=div(Uu)=λ∂2u∂x2u(x,0)=u0(x)u(xa)=ua(t),u(xb,t)=ub(t)0≤t≤te,xa<x<xb\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∂2uu(x,0)=u0(x)u(xa)=ua(t), u(xb,t)=ub(t)0≤t≤te, xa<x<xb

式中 λ\lambdaλ 为热扩散率，取决于材料本身的热传导率、密度和热容。

考虑一维热传导偏微分方程的数值解法，采用有限差分法求解。简单地，采用迎风法的三点差分格式，二阶导数的差分表达式为：
(∂2u∂x2)i,j=ui+1,j−2ui,j+ui−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
于是得到差分方程：
ui,j+1=ui,j+λdt/dx2∗(ui+1,j−2ui,j+ui−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 偏微分方程编程步骤

以该题为例一类有限差分法求解一维平流问题偏微分方程的步骤：

导入numpy、matplotlib 包；
输入模型参数 L, lam，定义求解的时间域 (0, te) 和空间域 (0, L)，设置差分步长 dt, dx；
初始化；
计算每一时点的边界条件 U[0,j]，U[L,j]、每一位置的初始条件U[i,0]；
递推求解差分方程在空间域 [0, L] 的数值解，获得网格节点的处的 U(x,t) 。

在例程中，设初值条件为 u(x,t=0)=0u(x,t=0) = 0u(x,t=0)=0，边界条件为 u(xa,t)=F(t),u(xb,t)=0u(x_a,t)=F(t), \ u(x_b,t)=0u(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.0ax2.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）是典型的双曲形偏微分方程，广泛应用于声学，电磁学，和流体力学等领域，描述自然界中的各种的波动现象，包括横波和纵波，例如声波、光波和水波。

考虑如下二维波动方程的初边值问题：
{∂2u∂t2=c2(∂2u∂x2+∂2u∂y2)∂∂tu(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\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=4c2Δt2Δx2+Δy2≤1r = \frac{4c^2\Delta t^2}{\Delta x^2 + \Delta y^2} \leq 1 r=Δx2+Δy24c2Δt2≤1

4.2 偏微分方程编程步骤

以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤：

导入numpy、matplotlib 包；
输入模型参数 c，定义求解的时间域 (0, te) 和空间域 (0,1)；
初始化，设置差分步长 dt, dx, dy，检验步长比以保证算法的稳定性；
计算初值条件U[0]、U[1]；
递推求解差分方程在空间域 [0, 1]*[0, 1] 的数值解，获得网格节点的处的 U(t,x,y) 。

在例程中，取初始条件为 u(x,y,0)=sin(6πx)+cos(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)=0u(0,y,t)=u(1,y,t)=0u(0,y,t)=u(1,y,t)=0， u(x,0,t)=u(x,1,t)=0u(x,0,t)=u(x,1,t)=0u(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=λ(∂2u∂x2+∂2u∂y2)+qvu(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\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)+qvu(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λ 为热扩散率，取决于材料本身的热传导率、密度和热容；qvq_vqv 是热源强度，可以是恒定值，也可以是时间、空间的函数。

考虑二维抛物型偏微分方程的数值解法，采用有限差分法求解。将上述的偏微分方程离散为差分方程：

化简后得到：

即可递推求得二维波动方程的数值解：
{Uk+1=Uk+rx∗A∗Uk+ry∗B∗Uk+qv∗dt)A=[−21⋯001−21⋯001−2⋯0⋮⋮⋮⋱⋮00⋯1−2](Nx+1,Nx+1)B=[−21⋯001−21⋯001−2⋯0⋮⋮⋮⋱⋮00⋯1−2](Ny+1,Ny+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⋮01−21⋮0⋯1−2⋮⋯0⋯⋯⋱1000⋮−2⎦⎥⎥⎥⎥⎥⎤(Nx+1,Nx+1)B=⎣⎢⎢⎢⎢⎢⎡−210⋮01−21⋮0⋯1−2⋮⋯0⋯⋯⋱1000⋮−2⎦⎥⎥⎥⎥⎥⎤(Ny+1,Ny+1)

5.2 偏微分方程编程步骤

以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤：

导入numpy、matplotlib 包；
输入热传导参数、热源参数等模型参数，定义求解的时间域 (0, te) 和空间域；
初始化，设置差分步长 dt, dx, dy，计算三对角系数矩阵 A、B，差分系数 rx, ry, ft；
计算初始条件 U0U_0U0；
递推求解差分方程在空间域的数值解，更新边界条件，获得网格节点的处的 U(x,y)U(x,y)U(x,y)；
绘制等温云图。

例程求解一个薄板受热的温度分布问题，其初始条件为 tIni=25tIni =25tIni=25，边界条件为 tBound=25tBound = 25tBound=25，热源为 QvQvQv，在空间域 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 pltdef 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] = uBoundU[0,:] = U[-1,:] = uBound# 差分求解U = U + rx * np.dot(U,A) + ry * np.dot(B,U) + Qv*dtif 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 二维椭圆方程的数学模型

椭圆型偏微分方程是一类重要的偏微分方程，主要用来描述物理的平衡稳定状态，如定常状态下的电磁场、引力场和反应扩散现象等，广泛应用于流体力学、弹性力学、电磁学、几何学和变分法中。
考虑如下二维泊松方程：
∂2u∂x2+∂2u∂y2=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)=0f(x,y)=0f(x,y)=0 时就是拉普拉斯方程（Laplace equation）。

考虑二维椭圆型偏微分方程的数值解法，采用有限差分法求解。简单地，采用五点差分格式表示二阶导数的差分表达式，将上述的偏微分方程离散为差分方程：
(ui−1,j−2ui,j+ui+1,j)/Δh2+(ui,j−1−2ui,j+ui,j+1)/Δh2=fi,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 较小时网格很多，矩阵求逆的内存占用和计算量极大。于是，可以使用迭代松弛法递推求得二维椭圆方程的数值解：
ui,jk+1=(1−w)∗ui,jk+w/4∗(ui,j+1k+ui,j−1k+ui−1,jk+ui+1,jk−h2∗fi,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 偏微分方程编程步骤

以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤：

导入numpy、matplotlib 包；
输入模型参数，定义空间域 (0,1)；
初始化，设置差分步长 h、松弛因子 w；
计算边值条件 u[0,:], u[1,:], u[:,0], u[:,1]；
迭代松弛法递推求解差分方程在空间域 [0, 1]*[0, 1] 的数值解。

在例程中，取边界条件为 u(x,0)=u(x,1)=0,u(0,y)=u(1,y)=1u(x,0)=u(x,1)=0, \; u(0,y)=u(1,y)=1u(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=simport 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. 小结

偏微分方程是数学中的重要学科分支，其数值解法的研究和方法可谓博大精深，远非本文所能涵盖。
本文针对 Python小白，采用有限差分法求解偏微分方程，通过案例介绍了一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程的数值解法，给出了例程和运行结果，供 Python 小白参考，以便进行数模学习。
需要特别说明的是：一阶差分格式的精度相对较低，往往需要较多网格计算。随着精度的不断提高，可以推导出各种差分表达式。高阶精度的差分公式可以给出质量更高的解，但需要使用更多的网格点进行计算，程序更复杂，计算量很大。类似地，显式方法的建立和编程简单，但要求步长较小才能保持稳定性，隐式方法的建立和编程更复杂，运算量大，但稳定性好。
需要特别说明的是：关于偏微分方程数值解法的稳定性、收敛性和误差分析，是非常专业的问题。例如步长选择不恰当可能导致算法不稳定，如 4.2 中应满足步长比 r>1。除非找到专业课程教材或范文中有相关内容可以参考套用，否则不建议小白自己摸索，这些问题不是调整参数试试就能试出来的。

【本节完】

欢迎关注『Python小白的数学建模课 @ Youcans』原创作品

原创作品，转载必须标注原文链接：https://blog.csdn.net/youcans/article/details/119755450

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数模笔记-模拟退火算法