1.偏微分方程基本知识
微分方程是指含有未知函数及其导数的关系式，偏微分方程是包含未知函数的偏导数（偏微分）的微分方程。
偏微分方程可以描述各种自然和工程现象， 是构建科学、工程学和其他领域的数学模型主要手段。科学和工程中的大多数实际问题都归结为偏微分方程的定解问题，如：波传播，流动和扩散，振动，固体力学，电磁学和量子力学，等等。
偏微分方程主要有三类：椭圆方程，抛物方程和双曲方程。
双曲方程描述变量以一定速度沿某个方向传播，常用于描述振动与波动问题。
椭圆方程描述变量以一定深度沿所有方向传播，常用于描述静电场、引力场等稳态问题。
抛物方程描述变量沿下游传播，常用于描述热传导和扩散等瞬态问题。
偏微分方程的定解问题通常很难求出解析解，只能通过数值计算方法对偏微分方程的近似求解。常用偏微分方程数值解法有：有限差分方法、有限元方法、有限体方法、共轭梯度法，等等。通常先对问题的求解区域进行网格剖分，然后将定解问题离散为代数方程组，求出在离散网格点上的近似值。
Python 语言求解偏微分方程的功能是比较弱的，主要有 Fipy, FEniCS 等有限元方法的工具包，另外还有机器学习工具如 Tensorflow 也可以进行偏微分方程的仿真模拟。但是，这些工具都不适合 Python 小白学习和使用，因此本篇采用比较简单的有限差分法对 5种典型的偏微分方程进行编程，通过案例讲解偏微分方程的数值解法。
2. 案例一：一维线性平流方程
2.1 一维线性平流方程的数学模型
平流过程是大气运动中重要的过程。平流方程（Advection equation）描述某一物理量的平流作用而引起局地变化的物理过程，最简单的形式是一维平流方程。
式中 u 为某物理量，v 为系统速度，x 为水平方向分量，t 为时间。
虽然该方程可以求得解析解：
考虑一维线性平流偏微分方程的数值解法，采用有限差分法求解。简单地， 采用一阶迎风格式的差分方法（First-order Upwind)，一阶导数的差分表达式为：
于是得到差分方程：
即可递推求得该平流方程的数值解。
2.2 偏微分方程编程步骤
以该题为例一类有限差分法求解一维平流问题偏微分方程的步骤：
1.导入numpy、matplotlib 包；
2.定义初始条件函数 U(x,0)；
3.输入模型参数 v, p，定义求解的时间域 (tStart, tEnd) 和空间域 (xMin, xMax)，设置差分步长 dt, nNodes；
4.初始化；
5.递推求解差分方程在区间 [xa,xb] 的数值解，获得网格节点的处的 u(t) 值。
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 的函数。
式中 λ \lambdaλ 为热扩散率，取决于材料本身的热传导率、密度和热容。
考虑一维热传导偏微分方程的数值解法，采用有限差分法求解。简单地， 采用迎风法的三点差分格式，二阶导数的差分表达式为：
于是得到差分方程：
即可递推求得一维热传导方程的数值解。
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) = 0u(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)=0u(x a ,t)=F(t), u(x b,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）是典型的双曲形偏微分方程，广泛应用于声学，电磁学，和流体力学等领域，描述自然界中的各种的波动现象，包括横波和纵波，例如声波、光波和水波。
考虑如下二维波动方程的初边值问题：
式中：u 是振幅；c 为波的传播速率，c 可以是固定常数，或位置的函数 c(x,y)，也可以是振幅的函数 c(u)。
考虑二维波动偏微分方程的数值解法，采用有限差分法求解。简单地， 采用迎风法的三点差分格式， 将上述的偏微分方程离散为差分方程 ：
化简后得到：
即可递推求得二维波动方程的数值解。
为了保证算法的收敛性，迎风法的三点差分格式要求步长比小于 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 ) = 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)=0u(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)=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
xa, xb = 0.0, 1.0 # 空间范围，xa
ya, yb = 0.0, 1.0 # 空间范围，ya
# 初始化
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）是典型的抛物形偏微分方程，也成为扩散方程。广泛应用于声学，电磁学，和流体力学等领域，描述自然界中的各种的波动现象，包括横波和纵波，例如声波、光波和水波。
考虑如下二维热传导方程的初边值问题：
式中 λ \lambdaλ 为热扩散率，取决于材料本身的热传导率、密度和热容；q v 是热源强度，可以是恒定值，也可以是时间、空间的函数。
考虑二维抛物型偏微分方程的数值解法，采用有限差分法求解。将上述的偏微分方程离散为差分方程 ：
化简后得到：
即可递推求得二维波动方程的数值解：
5.2 偏微分方程编程步骤
以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤：
导入numpy、matplotlib 包；
输入热传导参数、热源参数等模型参数，定义求解的时间域 (0, te) 和空间域；
初始化，设置差分步长 dt, dx, dy，计算三对角系数矩阵 A、B，差分系数 rx, ry, ft；
计算初始条件 U 0 U_0U
0
递推求解差分方程在空间域的数值解，更新边界条件，获得网格节点的处的 U ( x , y ) U(x,y)U(x,y)；
绘制等温云图。
例程求解一个薄板受热的温度分布问题，其初始条件为 t I n i = 25 tIni =25tIni=25，边界条件为 t B o u n d = 25 tBound = 25tBound=25，热源为 Q v QvQv，在空间域 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
xa, xb = 0.0, 16.0 # 空间范围，xa
ya, yb = 0.0, 12.0 # 空间范围，ya
# 初始化
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 二维椭圆方程的数学模型
椭圆型偏微分方程是一类重要的偏微分方程，主要用来描述物理的平衡稳定状态，如定常状态下的电磁场、引力场和反应扩散现象等，广泛应用于流体力学、弹性力学、电磁学、几何学和变分法中。
考虑如下二维泊松方程：
上式在 f ( x , y ) = 0 f(x,y)=0f(x,y)=0 时就是拉普拉斯方程（Laplace equation）。
考虑二维椭圆型偏微分方程的数值解法，采用有限差分法求解。简单地， 采用五点差分格式表示二阶导数的差分表达式，将上述的偏微分方程离散为差分方程 ：
椭圆型偏微分描述不随时间变化的均衡状态，没有初始条件，因此不能沿时间步长递推求解。对上式的差分方程，可以通过矩阵求逆方法求解，但当 h 较小时网格很多，矩阵求逆的内存占用和计算量极大。于是，可以使用迭代松弛法递推求得二维椭圆方程的数值解：
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)=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=s
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 求解范围
xa, xb = 0.0, 1.0 # 空间范围，xa
ya, yb = 0.0, 1.0 # 空间范围，ya
# 初始化
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。除非找到专业课程教材或范文中有相关内容可以参考套用，否则不建议小白自己摸索，这些问题不是调整参数试试就能试出来的。
热门跟贴