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

Posted 2021-09-06 youcans

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了Python小白的数学建模课-11.偏微分方程数值解法相关的知识，希望对你有一定的参考价值。

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

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

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

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

文章目录

1. 偏微分方程基本知识

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

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

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

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

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

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

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

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

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

$\\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}$

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

虽然该方程可以求得解析解：
$\\ 0 \\leq t \\leq t_e$

考虑一维线性平流偏微分方程的数值解法，采用有限差分法求解。简单地，采用一阶迎风格式的差分方法（First-order Upwind)，一阶导数的差分表达式为：
$(\\frac{\\partial u}{\\partial x})_{i,j}=\\frac{u_{i+1,j}-u_{i,j}}{\\Delta x}+O(\\Delta x)$
于是得到差分方程：
$u_{i,j+1}=u_{i,j}-v*dt/dx*(u_{i,j}-u_{i-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)$ ，取 $v = 1.0, p = 0.25$ ，空间域 $x\\in(0,\\pi)$ 。

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)$ 是时间变量 t 和水平方向空间变量 x 的函数。

$\\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}$

Python小白的数学建模课-10 微分方程边值问题

Python小白的数学建模课-A1.国赛赛题类型分析

Python小白的数学建模课-B5. 新冠疫情 SEIR模型

Python小白的数学建模课-B6. 新冠疫情 SEIR 改进模型