用微信扫码二维码

分享至好友和朋友圈

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 t1。除非找到专业课程教材或范文中有相关内容可以参考套用，否则不建议小白自己摸索，这些问题不是调整参数试试就能试出来的。

特别声明：以上内容(如有图片或视频亦包括在内)为自媒体平台“网易号”用户上传并发布，本平台仅提供信息存储服务。

Notice: The content above (including the pictures and videos if any) is uploaded and posted by a user of NetEase Hao, which is a social media platform and only provides information storage services.