那些做网站的那些软件都叫啥,邯郸哪个公司做网站好,深圳住房建设保障局,百度网站建设的十一个偏微分方程可以描述各种自然和工程现象#xff0c; 是构建科学、工程学和其他领域的数学模型主要手段。
偏微分方程主要有三类#xff1a;椭圆方程#xff0c;抛物方程和双曲方程。
本文采用有限差分法求解偏微分方程#xff0c;通过案例讲解一维平流方程、一维热传导方程…
偏微分方程可以描述各种自然和工程现象 是构建科学、工程学和其他领域的数学模型主要手段。
偏微分方程主要有三类椭圆方程抛物方程和双曲方程。
本文采用有限差分法求解偏微分方程通过案例讲解一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程等常见类型的偏微分方程的数值解法给出了全部例程和运行结果。
欢迎关注『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∂tv∂u∂x0u(x,0)F(x)0≤t≤texaxxb\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_axx_b \end{aligned} \end{cases} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧∂t∂uv∂x∂u0u(x,0)F(x)0≤t≤texaxxb
式中 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,jui1,j−ui,jΔxO(Δx)(\frac{\partial u}{\partial x})_{i,j}\frac{u_{i1,j}-u_{i,j}}{\Delta x}O(\Delta x) (∂x∂u)i,jΔxui1,j−ui,jO(Δx) 于是得到差分方程 ui,j1ui,j−v∗dt/dx∗(ui,j−ui−1,j)u_{i,j1}u_{i,j}-v*dt/dx*(u_{i,j}-u_{i-1,j}) ui,j1ui,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)取 v1.0,p0.25v 1.0, p0.25v1.0,p0.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, xb2*dx, dx)
ux0 funcUx0(x, p)u ux0.copy() # u(j)
ujp ux0.copy() # u(j1)# 时域差分
for i in range(nsteps):plt.clf() # 清除当前 figure 的所有axes, 但是保留当前窗口# 计算 u(j1)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[nNodes2] u[1]# 绘图plt.plot(x, u, b-, labelv1 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,colorgainsboro)plt.pause(0.001)tc dtplt.show()2.4 Python 例程运行结果 3. 案例二一维热传导方程
3.1 一维热传导方程的数学模型
热传导方程描述一个区域内的温度随时间的变化是典型的抛物型偏微分方程也称为扩散方程。
一维热传导方程考虑各向同性的均匀细杆在垂直于 x 轴的截面上的温度相同细杆的圆周与周围环境无热交换杆内无热源则温度 uu(t,x)uu(t,x)uu(t,x) 是时间变量 t 和水平方向空间变量 x 的函数。
{∂u∂tdiv(Uu)λ∂2u∂x2u(x,0)u0(x)u(xa)ua(t),u(xb,t)ub(t)0≤t≤te,xaxxb\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∂udiv(Uu)λ∂x2∂2uu(x,0)u0(x)u(xa)ua(t), u(xb,t)ub(t)0≤t≤te, xaxxb
式中 λ\lambdaλ 为热扩散率取决于材料本身的热传导率、密度和热容。
考虑一维热传导偏微分方程的数值解法采用有限差分法求解。简单地 采用迎风法的三点差分格式二阶导数的差分表达式为 (∂2u∂x2)i,jui1,j−2ui,jui−1,j(Δx)2O(Δx)2(\frac{\partial ^2u}{\partial x^2})_{i,j} \frac{u_{i1,j}-2u_{i,j}u_{i-1,j}}{(\Delta x)^2}O(\Delta x)^2 (∂x2∂2u)i,j(Δx)2ui1,j−2ui,jui−1,jO(Δx)2 于是得到差分方程 ui,j1ui,jλdt/dx2∗(ui1,j−2ui,jui−1,j)u_{i,j1} u_{i,j} \lambda dt/dx^2*(u_{i1,j}-2u_{i,j}u_{i-1,j}) ui,j1ui,jλdt/dx2∗(ui1,j−2ui,jui−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,t0)0u(x,t0) 0u(x,t0)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([nNodes1, nSteps1]) # 建立二维数组# 边界条件
for j in np.arange(0, nSteps1): # 时间序列U[0,j] 7.5 (nSteps-j)/2000 * np.sin(j/1000)/(1np.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,j1] K*U[i1,j] (1-2*K)*U[i,j] K*U[i-1,j]# 绘图
xZone np.arange(0, (nNodes1)*dx, dx) # 建立空间网格
tZone np.arange(0, (nSteps1)*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,nSteps1,round(nSteps/5)):ax1.plot(xZone, U[:,k], labelrx{}.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(locupper right)
rect2 [0.55, 0.2, 0.4, 0.65] # [左, 下, 宽, 高], 0.0~1.0
ax2 plt.axes(rect2)
for k in range(0,nNodes1,round(nNodes/5)): # U[nNodes,k] 0.0ax2.plot(tZone, U[k,:], labelrt{}.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(locupper right)
plt.show()3.4 Python 例程运行结果 4. 案例三二维双曲型方程
4.1 二维波动方程的数学模型
波动方程wave equation是典型的双曲形偏微分方程广泛应用于声学电磁学和流体力学等领域描述自然界中的各种的波动现象包括横波和纵波例如声波、光波和水波。
考虑如下二维波动方程的初边值问题 {∂2u∂t2c2(∂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,0x1,0y1\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∂2uc2(∂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, 0x1, 0y1
式中u 是振幅c 为波的传播速率c 可以是固定常数或位置的函数 c(x,y)也可以是振幅的函数 c(u)。
考虑二维波动偏微分方程的数值解法采用有限差分法求解。简单地 采用迎风法的三点差分格式 将上述的偏微分方程离散为差分方程 化简后得到
即可递推求得二维波动方程的数值解。
为了保证算法的收敛性迎风法的三点差分格式要求步长比小于 1
r4c2Δ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/px2p2u/py2)import numpy as np
import matplotlib.pyplot as plt# 模型参数
c 1.0 # 波的传播速率
tc, te 0.0, 1.0 # 时间范围0tte
xa, xb 0.0, 1.0 # 空间范围xaxxb
ya, yb 0.0, 1.0 # 空间范围yayyb# 初始化
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, (tNodes1)*dt, dt) # 建立空间网格
xZone np.arange(0, (xNodes1)*dx, dx) # 建立空间网格
yZone np.arange(0, (yNodes1)*dy, dy) # 建立空间网格
xx, yy np.meshgrid(xZone, yZone) # 生成网格点的坐标 xx,yy (二维数组)# 步长比检验(r1 则算法不稳定)
r 4 * c2 * dt*dt / (dx*dxdy*dy)
print(dt {:.2f}, dx {:.2f}, dy {:.2f}, r {:.2f}.format(dt,dx,dy,r))
assert r 1.0, Error: r1, unstable step ratio of dt2/(dx2dy2) .
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, projection3d)
# ax2 fig.add_subplot(122, projection3d)# 计算初始值
U np.zeros([tNodes1, xNodes1, yNodes1]) # 建立三维数组
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,:,:], rstride2, cstride2, cmapplt.cm.coolwarm)
# wframe ax2.plot_wireframe(xx, yy, U[0], rstride2, cstride2, linewidth1)# 有限差分法求解
for k in range(2,tNodes1):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,i1,j]) ry*(U[k-1,i,j-1]U[k-1,i,j1])\ 2*(1-rx-ry)*U[k-1,i,j] -U[k-2,i,j]surf ax1.plot_surface(xx, yy, U[k,:,:], rstride2, cstride2, cmaprainbow)# wframe ax2.plot_wireframe(xx, yy, U[k,:,:], rstride2, cstride2, linewidth1, cmaprainbow)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,0x1,0y1\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, 0x1, 0y1
式中 λ\lambdaλ 为热扩散率取决于材料本身的热传导率、密度和热容qvq_vqv 是热源强度可以是恒定值也可以是时间、空间的函数。
考虑二维抛物型偏微分方程的数值解法采用有限差分法求解。将上述的偏微分方程离散为差分方程 化简后得到
即可递推求得二维波动方程的数值解 {Uk1Ukrx∗A∗Ukry∗B∗Ukqv∗dt)A[−21⋯001−21⋯001−2⋯0⋮⋮⋮⋱⋮00⋯1−2](Nx1,Nx1)B[−21⋯001−21⋯001−2⋯0⋮⋮⋮⋱⋮00⋯1−2](Ny1,Ny1)\begin{cases} \begin{aligned} U_{k1} 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} _{(Nx1,Nx1)} 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} _{(Ny1,Ny1)} \end{aligned} \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧Uk1Ukrx∗A∗Ukry∗B∗Ukqv∗dt)A⎣⎢⎢⎢⎢⎢⎡−210⋮01−21⋮0⋯1−2⋮⋯0⋯⋯⋱1000⋮−2⎦⎥⎥⎥⎥⎥⎤(Nx1,Nx1)B⎣⎢⎢⎢⎢⎢⎡−210⋮01−21⋮0⋯1−2⋮⋯0⋯⋯⋱1000⋮−2⎦⎥⎥⎥⎥⎥⎤(Ny1,Ny1) 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)绘制等温云图。
例程求解一个薄板受热的温度分布问题其初始条件为 tIni25tIni 25tIni25边界条件为 tBound25tBound 25tBound25热源为 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/px2p2u/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(zMatzMax)[0][0], np.where(zMatzMax)[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, cmapplt.cm.get_cmap(jet), originlower, levels levels)plt.annotate(showText, xy(x[xMax],y[yMax]), xytext(x[xMax],y[yMax]),fontsize10)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 # 时间范围0tte
xa, xb 0.0, 16.0 # 空间范围xaxxb
ya, yb 0.0, 12.0 # 空间范围yayyb# 初始化
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, xNodes1) # 建立空间网格
yZone np.linspace(ya, yb, yNodes1) # 建立空间网格
xx,yy np.meshgrid(xZone, yZone) # 生成网格点的坐标 xx,yy (二维数组)# 计算 差分系数矩阵 A、B (三对角对称矩阵)差分系数 rx,ry,ft
A (-2) * np.eye(xNodes1, k0) (1) * np.eye(xNodes1, k-1) (1) * np.eye(xNodes1, k1)
B (-2) * np.eye(yNodes1, k0) (1) * np.eye(yNodes1, k-1) (1) * np.eye(yNodes1, k1)
rx, ry, ft c*dt/(dx*dx), c*dt/(dy*dy), qv*dt# 计算 初始值
U uIni * np.ones((yNodes1, xNodes1)) # 初始温度 u0# 前向欧拉法一阶差分求解
for k in range(tNodes1):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 t2.0 else qv# (3) 热源位置随时间变化 Qv(x,y,t)f(x(t),y(t))xt, yt x0vx*t, y0vy*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∂y2f(x,y)\frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2} f(x,y)\\ ∂x2∂2u∂y2∂2uf(x,y)
上式在 f(x,y)0f(x,y)0f(x,y)0 时就是拉普拉斯方程Laplace equation。
考虑二维椭圆型偏微分方程的数值解法采用有限差分法求解。简单地 采用五点差分格式表示二阶导数的差分表达式将上述的偏微分方程离散为差分方程 (ui−1,j−2ui,jui1,j)/Δh2(ui,j−1−2ui,jui,j1)/Δh2fi,j(u_{i-1,j}-2u_{i,j}u_{i1,j})/\Delta h^2 (u_{i,j-1}-2u_{i,j}u_{i,j1})/\Delta h^2 f_{i,j} (ui−1,j−2ui,jui1,j)/Δh2(ui,j−1−2ui,jui,j1)/Δh2fi,j 椭圆型偏微分描述不随时间变化的均衡状态没有初始条件因此不能沿时间步长递推求解。对上式的差分方程可以通过矩阵求逆方法求解但当 h 较小时网格很多矩阵求逆的内存占用和计算量极大。于是可以使用迭代松弛法递推求得二维椭圆方程的数值解 ui,jk1(1−w)∗ui,jkw/4∗(ui,j1kui,j−1kui−1,jkui1,jk−h2∗fi,j)u_{i,j}^{k1} (1-w)*u_{i,j}^k w/4*(u_{i,j1}^k u_{i,j-1}^k u_{i-1,j}^k u_{i1,j}^k -h^2* f_{i,j}) ui,jk1(1−w)∗ui,jkw/4∗(ui,j1kui,j−1kui−1,jkui1,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_xxu_yysimport numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# 求解范围
xa, xb 0.0, 1.0 # 空间范围xaxxb
ya, yb 0.0, 1.0 # 空间范围yayyb# 初始化
h 0.01 # 空间步长, dx dy 0.01
w 0.5 # 松弛因子
nodes round((xb-xa)/h) # x轴 空间网格数# 边值条件
u np.zeros((nodes1, nodes1))
# u[:, 0] 0
# u[:, -1] 0 # -1 表示数组的最后一个值
# u[0, :] 1
# u[-1, :] 1 # -1 表示数组的最后一个值
for i in range(nodes1):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[i1, j] u[i, j-1] u[i, j1]) (1-w) * u[i, j]# 绘图
x np.linspace(0, 1, nodes1)
y np.linspace(0, 1, nodes1)
xx, yy np.meshgrid(x, y)
fig plt.figure(figsize(8,6))
ax fig.add_subplot(111, projection3d)
surf ax.plot_surface(xx, yy, u, cmapplt.get_cmap(rainbow))
fig.colorbar(surf, shrink0.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 中应满足步长比 r1。除非找到专业课程教材或范文中有相关内容可以参考套用否则不建议小白自己摸索这些问题不是调整参数试试就能试出来的。
【本节完】 版权声明
欢迎关注『Python小白的数学建模课 Youcans』 原创作品
原创作品转载必须标注原文链接https://blog.csdn.net/youcans/article/details/119755450
Copyright 2021 Youcans, XUPT
Crated2021-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数模笔记-模拟退火算法
