偏微分方程数值解
00 分钟
2024-7-28

偏微分方程数值解

1. 模型概览

1.1 学科分属

  • 偏微分方程（Partial differential equation，简称PDE）
偏微分方程是微积分学中来源极广，历史久远的一个重要分支，属于微分方程理论的重要组成部分。
微分方程是伴随着微积分学一起发展起来的。微积分学的奠基人Newton和Leibniz的著作中都处理过与微分方程有关的问题。微分方程的应用十分广泛，可以解决许多与导数有关的问题。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。此外，微分方程在化学、工程学、经济学和人口统计等领域都有应用。
其中，微分方程中，未知函数是一元函数的叫做常微分方程（ODE）；未知函数是多元函数的称为偏微分方程（PDE）。
  • 偏微分方程数值解（Numerical solution of partial differential equations）
偏微分方程数值解是数值分析领域的重要内容。
数值分析(numerical analysis)，为数学的一个分支，是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科。它以计算机求解数学问题的理论和方法为研究对象，为计算数学的主体部分。
偏微分方程数值解是指通过数值计算方法，在计算机上对偏微分方程的近似求解。科学和工程中的大多数实际问题都归结为偏微分方程的定解问题，由于很难求得这些定解问题的解析解（在经典意义下甚至没有解），人们转向求解它们的数值近似解。通常先对问题的求解区域进行网格剖分，然后基于有限元法、有限差分法和有限体积法等数值方法，对原定解问题或其等价形式离散，并归结为一个线性代数方程组，最终在计算机上求得精确解在离散网格点上的近似值。

1.2 历史发展

微积分方程这门学科产生于十八世纪，欧拉在他的著作中最早提出了弦振动的二阶方程，随后不久，法国数学家达朗贝尔也在他的著作《论动力学》中提出了特殊的偏微分方程。不过这些著作当时没有引起多大注意。
1746年，达朗贝尔在他的论文《张紧的弦振动时形成的曲线的研究》中，提议证明无穷多种和正弦曲线不同的曲线是振动的模式。这样就由对弦振动的研究开创了偏微分方程这门学科。
和欧拉同时代的瑞士数学家丹尼尔·贝努利也研究了数学物理方面的问题，提出了解弹性系振动问题的一般方法，对偏微分方程的发展起了比较大的影响。拉格朗日也讨论了一阶偏微分方程，丰富了这门学科的内容。
偏微分方程得到迅速发展是在十九世纪，那时候，数学物理问题的研究繁荣起来了，许多数学家都对数学物理问题的解决做出了贡献。其中法国数学家傅立叶，在从事热流动的研究中，写出了《热的解析理论》，提出了三维空间的热传导方程，也就是一种偏微分方程。他的研究对偏微分方程的发展的影响很大。
数值近似求解的研究由来已久，但只是在二十世纪后期电子计算机产生后，才得到广泛的发展和应用（如有限元理论始于60年代）。目前数值求解的规模也变得更大，例如在航天器设计、湍流模拟、气候预测、油田开发等各种实际问题中，经常过到大规模（网格数至少在百万以上）的运算量问题。偏微分方程的数值求解已渗透到物理、化学、生物等现代科学与工程的各领域，对科技和国民经济的发展有重要作用。

2. 模型介绍

偏微分方程是构建科学、工程学和其他领域的数学模型的主要手段。一般情况下，这些模型都需要用数值方法去求解。借助抛物线型、双曲线型和椭圆型方程常用的有有限差分方法、有限元方法、有限体方法、修正方程分析、辛积分格式、对流扩散问题、多重网络、共轭梯度法。
这里我们主要介绍有限差分方法、有限元方法和有限体方法。

2.1 有限差分法

有限差分方法(Finite Difference Methods，简称FDM）是数值模拟偏微分方程最早采用的方法，至今仍被广泛使用。该方法包括区域剖分和差商代替导数两个步骤。首先将求解区域划分为差分网格，用有限个网格节点代替连续的求解区域。其次，利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替进行离散，从而建立以网格节点上的值为未知量的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且十分成熟的数值方法。
差商代替导数后的格式称为有限差分格式，从格式的精度来考虑，有一阶格式、二阶格式和高阶格式。从差分的空间离散形式来考虑，有中心格式和迎风格式。对于瞬态方程，考虑时间方向的离散，有显格式、隐格式、交替显隐格式等。目前常见的差分格式，主要是以上几种格式的组合，不同的组合构成不同的差分格式。差分方法主要适用于结构网格，网格的大小一般根据问题模型和 Courant稳定条件来决定。
利用有限差分法求偏微分方程的数值解的核心思路在于将微分方程的定解问题转化成离散系统的求解问题，一般步骤如下：
  • 对微分方程及定解条件选择差分近似，列出差分格式;
  • 结合初始条件与边界条件，求解差分格式；
  • 讨论差分格式解对于微分方程解的收敛性及误差估计。
另外，有一些常见概念如下：
  • 有限差分格式的截断误差（Truncation Error）
对于一个偏微分方程Lu = 0​​​​，令方程的真实解为ue​​​​，数值解为u​​​​，离散模型为LΔu = 0​​​​​​。于是截断误差定义为T(xt) = LΔue​​​​，此外可以通过该定义得到局部截断误差（Local Truncation Error）T(xitn) = LΔue(xjtn)​​​​​。
  • 有限差分格式的相容性(Consistency)
即，这意味着差分方程能和微分方程充分接近。
  • 有限差分格式的稳定性（Stability）
对初始误差ϵ0以及范数$||\epsilon^n||_h=\sqrt{\sum\limits_jh(u^n_j)^2}$，则差分格式是稳定的当且仅当||ϵn||h ≤ K||ϵ0||h。
  • 有限差分格式的收敛性（Convergence）
即。直接分析收敛性常常很复杂，一般采用间接的方式。

2.1.1 求解形如 ut + cux = 0 的数值解

确定网格如下图所示：
notion image
notion image
设方程的数值解为 ，则
代入得差分格式
其中 .
当 c > 0 时，仅Backward形式是稳定的；当 c < 0 时，仅Forward形式是稳定的。一般对于一种情景仅选取一种形式，对于边界情况选定所能够使用到的形式单独列出边界条件：
  • 在左边界 x = 0 使用Forward形式，在右边界 x = l 使用Backward形式。即
$$ u_0^{n+1}=u_0^n-\mu(u_1^n-u_0^n)\\ u_l^{n+1}=u_l^n-\mu(u_l^n-u_{l-1}^n) $$
对x与t取适当划分代入即可得到一个有唯一解的线性方程组，出于数值解的稳定性（即不会出现类龙格振荡）一般要求 |μ| ≤ 1 即 $$ \Delta t\le\frac{\Delta x}{|c|} $$

2.1.2 求解形如 ut = cuxx 的数值解

对方程两边使用central形式差分得到 $$ \frac{u_i^{n+1}-u_i^n}{\Delta t}=c\frac{u^n_{i-1}-2u^n_i+u^n_{i+1}}{\Delta x^2} $$ 整理得 uin + 1 = Pui − 1n + (1−2P)uin + Pui + 1n 其中 $P=\displaystyle\frac{c\Delta t}{\Delta x^2}$ 。
为使结果稳定，应满足 $P\le\displaystyle\frac{1}{2}$ ，即 $$ \Delta t\le \frac{\Delta x^2}{2c} $$

2.1.3 求解形如 uxx + uyy = f(xy) 的数值解

设方程的数值解为 ui, j ≈ u(xiyj) ，则有 $$ \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}+\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta y^2}=f(x_i,y_j) $$ 令 h2 = Δx2 = Δy2 整理得五点菱形式，即 $$ f(x_i,y_j)=\frac{1}{h^2}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}) $$

2.2 有限元法

有限元方法(Finite Element Methods，简称FEM）的基础是虚位移原理和分片多项式插值。该方法的构造过程包括以下三个步骤。首先，利用虚位移原理得到偏微分方程的弱形式，将计算区域划分为有限个互不重叠的单元（三角形、四边形、四面体、六面体等)，在每个单元上选择合适的节点作为求解函数的插值点，将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式，得到微分方程的离散形式。利用插值函数的局部支集性质及数值积分可以得到未知量的代数方程组。
有限元方法有较完善的理论基础，具有求解区域灵活（复杂区域)、单元类型灵活（适于结构网格和非结构网格)、程序代码通用（数值模拟软件多数基于有限元方法）等特点。有限元方法最早应用于结构力学，随着计算机的发展已经渗透到计算物理、流体力学与电磁学等各个数值模拟领域。
根据所采用的检验函数（虚位移函数）和插值函数的不同，有限元方法也分为多种计算格式。从检验函数的选择来说，有配置法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多面体网格等，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。
对于有限元方法，其基本思路和解题步骤可归纳为：
  • 建立积分方程，根据虚位移原理或方程余量，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。
  • 区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分采用有限元方法的前处理完成，并给出计算单元和节点编号相互之间的关系、节点的位置坐标，同时还需要列出问题的边界的节点号和相应的边值条件。
  • 确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元的形函数。有限元方法中的形函数是在单元中选取的，由于各单元具有规则的几何形状，在选取形函数时可遵循一定的法则。
  • 单元分析:将各个单元中的求解函数用单元形函数的线性组合表达式进行逼近;再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各节点的函数值)的单元矩阵与荷载。
  • 总体合成:在得出单元矩阵与荷载之后，将区域中所有单元矩阵与荷载按一定法则进行迭加，形成总体有限元方程。
  • 边界条件的处理:一般边界条件有三种形式，分为本质边界条件(Dirichlet边界条件)、自然边界条件(Neumann边界条件)、混合边界条件(Cauchy边界条件)。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则后对总体有限元方程进行修正。
  • 解有限元方程:根据边界条件修正的总体有限元方程组，采用适当的代数方程组求解器，求出各节点的函数值。

2.3 有限体积法

有限体积法(Finite Volume Method，简称FVM）又称为控制体积法。其基本思路是:将计算区域划分为一系列互不重叠的控制体，并使每个网格点周围有一个控制体;将待求解的微分方程对每一个控制体积积分，便得出一组离散方程。该方法的未知量为网格点上的函数值。为了求出控制体积的积分，须假定函数值在网格点控制体边界上的变化规律。从积分区域的选取方法来看，有限体积法属于有限元方法中检验函数取分片常数插值的子区域法;从未知量的近似方法看来，有限体积法属于采用局部近似多项式插值逼近。
有限体积法的基本思路易于理解，能够保持物理量在控制体上的守恒性质，也即离散方程保持了微分方程物理量在控制体满足某种守恒原理的物理意义这是有限体积法吸引人的优点。此外，在有限体积法中，插值函数只用于计算控制体积的积分，因此可以对微分方程中不同的项采取不同的插值函数。
一般步骤如下：
  • 将计算区域划分为一系列不重复的控制体积,每一个控制体积都有一个节点作代表,将待求的守恒型微分方程在任一控制体积及一定时间间隔内对空间与时间作积分。
  • 对待求函数及其导数对时间及空间的变化型线或插值方式作出假设。
  • 对步骤1中各项按选定的型线作出积分并整理成一组关于节点上未知量的离散方程。

2.4 三种方法的优缺点

有限差分方法直观，经验丰富，格式众多。但是不规则区域处理繁琐，虽然网格生成可以使FDM应用于不规则区域，但是对区域的形状有较大的限制，并且使用不方便，FDM是三种方法中计算量最少的一种，并且易于编程。
有限元方法适合处理复杂区域和各种边值条件，但程序复杂，编程量和计算量是三种方法之首。
有限体积法主要用于流体与传热传质的计算，可以简单应用于非结构网格，能处理复杂区域，但是精度较低，对边值处理较繁琐，不如有限元灵活，程序量与计算量皆居中。

3. 模型应用

3.1 常见应用场景

偏微分方程数值解常用于热传导在三维的等方向均匀介质里的传播，其可用以下方程表达：
$$ \displaystyle\frac{\partial}{\partial x}\Big(k\frac{\partial T}{\partial x}\Big)+\frac{\partial}{\partial y}\Big(k\frac{\partial T}{\partial y}\Big)+\frac{\partial}{\partial z}\Big(k\frac{\partial T}{\partial z}\Big)+q=\rho c_p \frac{\partial T}{\partial t} $$ 其中 t 为时间； x, y, z 为坐标轴； ρ 为密度； cp 为定压比热容；热扩散方程表明：在介质中任意一点处，由传导进入单位体积的净导热速率加上单位体积的热能产生速率必定等于单位体积内所贮存的能量变化速率。
除了热量的扩散，对于物质的扩散问题，在石油开采、环境污染、疾病流行、化学反应、新闻传播、煤矿瓦斯爆炸、水里工程、房屋基建、神经传导、药物在人体内分布等情况大都能由线性或非线性的偏微分方程来定量或者定性的解决。

3.2 数学建模竞赛应用

3.2.1 CUMCM2020A 炉温曲线

为了方便求解，将炉内温度简化为一维热稳态问题。 在一维情况下， 3.1.1 中提到的热传导方程转化为 $$ \displaystyle\frac{\partial^2 T(x,t)}{\partial x^2}=\frac{c\rho}{\lambda}\cdot\frac{\partial T}{\partial t}\\ $$ 对于热稳态问题有 $\displaystyle\frac{\partial T}{\partial t}=0$ ，从而有 $\displaystyle\frac{\partial^2 T(x,t)}{\partial x^2}=0$ ，炉内温度的空间变化率 $K_x=\displaystyle\frac{dT(x)}{dx}=\frac{t_1-t_2}{D}$ （设厚度为 D 的均匀介质两侧分别保持恒温 t1 与 t2 ）
由于电路板在炉内匀速直线运动，则 x = vt，dx = vdt ，从而温度的时间变化率 $K_t=\displaystyle\frac{dT}{dt}=\frac{dT}{dx/v}=vK_x$
又 $$ \displaystyle\frac{dQ}{dt}=hA(T_0+K_t(t-t_0))=cm\frac{dT}{dt}=\rho cdA\frac{dT}{dt} $$ 解得 Tt = (T−T0)e−α(t−t0) + T0
之后经参数分段拟合差值，从第一个温区开始，作出理论值与实验值之差与时间的关系 ，对参数开始进行二分法计算可得具体的温度曲线，通过6至7次迭代可以得到符合精度要求的参数。如图所示：
notion image
notion image

3.2.2 CUMCM2018A 高温作业专用服装设计

针对“环境-防护服-人体”系统，提出高温下织物以及空气层的热传递数学模型。
notion image
notion image
  • 建立热力学传导方程：
$$ \begin{cases} \displaystyle\frac{\partial T_1}{\partial t}=D_1\frac{\partial^2 T_1}{\partial x^2}, & 0\le x<x_1\\ \displaystyle\frac{\partial T_2}{\partial t}=D_2\frac{\partial^2 T_2}{\partial x^2}, & x_1\le x<x_2\\ \displaystyle\frac{\partial T_3}{\partial t}=D_3\frac{\partial^2 T_3}{\partial x^2}, & x_2\le x<x_3\\ \displaystyle\frac{\partial T_4}{\partial t}=D_4\frac{\partial^2 T_4}{\partial x^2}, & x_3\le x<x_4\\ \end{cases} $$
  • 初始条件：对于含有时间变量 t 的数理方程来说，其未知函数将随t 不同而有不同的值，故必然要反映某一时刻物理量与相邻时刻的同一物理量之间的关系，所以求解问题过程中必须追溯到早先某个所谓“初始”时刻的状况，即建立初始条件。t=0时四个介质层的初始温度都为 Th = 37 ℃ ，可得四个区域的热传导方程初始条件： T1(x,0) = T2(x,0) = T3(x,0) = T4(x,0) = Th
  • 衔接条件：
$$ \left\{ \begin{aligned} T_1(x_1-0,t)=T_2(x_1+0,t)\\ T_2(x_2-0,t)=T_3(x_2+0,t)\\T_3(x_3-0,t)=T_4(x_3+0,t)\\ k_1\frac{\partial T_1}{\partial x}\Bigg|_{x=x_1-0}=k_2\frac{\partial T_2}{\partial x}\Bigg|_{x=x_1+0}\\ k_2\frac{\partial T_2}{\partial x}\Bigg|_{x=x_2-0}=k_3\frac{\partial T_3}{\partial x}\Bigg|_{x=x_2+0}\\ k_3\frac{\partial T_3}{\partial x}\Bigg|_{x=x_3-0}=k_4\frac{\partial T_4}{\partial x}\Bigg|_{x=x_3+0}\\ \end{aligned} \right. $$
  • 边界条件：
$$ \left\{ \begin{aligned} T_1(0,t)=T_s\\ \frac{\partial T_4}{\partial x}\Bigg|_{x=x_4}=\delta[T_4(x_4,0)-T_h] \end{aligned} \right. $$
在本问题中采用有限差分法处理求解出偏微分方程的数值解。得到对应条件的差分方程：
  • 热传导差分方程：
$$ \displaystyle\frac{T(k,j)-T(k,j-1)}{\Delta t}=D\frac{T(k+1,j)-2T(k,j)+T(k-1,j)}{(\Delta x)^2} $$ 引入变量 r = DΔt/(Δx)2 : T(kj) + T(kj−1) = rT(k+1,j) + 2rT(kj) + rT(k−1,j) 。此方程为三对角问题，应用三对角矩阵算法（追赶法）即可得到T T(kj) 而不需要对矩阵直接求逆。
  • 初始条件： T(xk,0) = Ts
  • 衔接条件：
$$ k_1\displaystyle\frac{T_1(x_1-0,j)-T_1(x_1+0,j)}{\Delta t}=k_2\frac{T_1(x_1-0,j)-T_2(x_1+0,j)}{\Delta t} $$
  • 边界条件：
$$ \displaystyle\frac{T_4(x_4,j)-T_4(x_4,j)}{\Delta t}=\delta(t)[T_4(x_4,0)-T_h] $$
依据附件中数据即可得到所求温度分布。
notion image
notion image

3.2.3 MCM1990A The Brain-Drug Problem & MCM1994B The Concrete Slab Problem

这两个问题都属于衰减类的扩散问题，均可以用线性抛物型方程解决。
设 u(xyzt) 是 t 时刻点(xyz)处一种物质的浓度。任取一个闭曲面 S ，它所围的区域是 Ω ，由于扩散，从 t 到 t + Δt 时刻这段时间内，通过 S 流入 Ω 的质量为
$$ M_1=\displaystyle\int_t^{t+\Delta t}\iint_S\Big(a^2\frac{\partial u}{\partial x}cos\alpha+b^2\frac{\partial u}{\partial y}cos\beta+c^2\frac{\partial u}{\partial z}cos\gamma\Big)dSdt $$
由高斯公式， $M_1=\displaystyle\int_t^{t+\Delta t}\iiint_{\Omega}\Big(a^2\frac{\partial u}{\partial x}+b^2\frac{\partial u}{\partial y}+c^2\frac{\partial u}{\partial z}\Big)dxdydzdt$
其中， a2, b2, c2 分别是沿 x, y, z 方向的扩散系数。
由于衰减（例如吸收、代谢等）， Ω 内的质量减少为
M2 = ∫tt + Δt∭Ωk2udxdydzdt
其中 k2 为衰减系数.
从而积存在 Ω 内的质量为 M1 − M2 ，故有
$$ \displaystyle\frac{\partial u}{\partial t}=a^2\frac{\partial^2u}{\partial x^2}+b^2\frac{\partial^2u}{\partial y^2}+c^2\frac{\partial^2u}{\partial z^2}-k^2u $$
此方程为衰减的扩散过程的数学模型，对于具体问题与相应的初始条件与边界条件匹配求得确定情况下的数值解。

4.软件/程序介绍

4.1 模拟介绍

MATLAB有大量成熟的偏微分方程函数与工具包，功能简洁强力；同时MATLAB资料库成熟，可及性强，故以下着重介绍MATLAB相关程序。

4.1.1 差分求解

使用差分解法进行求解。
先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。
以Poisson方程为例：
设x，y都在0到1之间，x，y方向的步长均取为 $\frac{1}{3}$ 代码示例如下。
常见的差分格式有古典显式格式、古典隐式格式、杜福特—弗兰克尔（DoFort—Frankel）格式等。
方法同样可解抛物线型、双曲型方程等，但差分的实现代码需要一定熟练度。以下介绍另一种方便的解法，运用MATLAB自带工具包与函数进行运算。

4.1.2 使用函数进行求解

4.1.2.1 一维状态空间的偏微分方程的 MATLAB 解法

MATLAB 有专门的指令pdepe解一维的抛物型或椭圆型方程组的初边值问题，使用指令pdepe要写三个辅助函数文件，即微分方程函数文件、初始条件函数文件、边界条件函数文件。
sol = pdepe(mpdefunicfunbcfunxmeshtspanoptions)
  • 其中xmesh，tspan分别为空间变量x和时间变量t的网格点位置向量。
  • xmesh取值对精确度影响很大，应观察解的特性，在变化较快处更加密集地取点。
  • pdepe求解的算法主要是将原来的椭圆型和拋物型偏微分方程转化为一组常微分方程。此转换的过程是基于使用者所指定的 mesh 点，以二阶空间离散化(spatial discretization)技术为之(Keel and Berzins,1990)，然后以 ode15s 的指令求解。
步骤如下：
  1. 将各方程化为标准形式。 $$ c(x,t,u,\partial u/\partial x)\frac{\partial u}{\partial t}=f(x,t,u,\partial u/\partial x)\frac{\partial u}{\partial x}+s(x,t,u,\partial u/\partial x) $$
pa(xtu) + qa(xt)f(xtu,∂u/∂x) = 0, x = a, t0 ≤ t ≤ tf
pb(xtu) + qb(xt)f(xtu,∂u/∂x) = 0, x = b, t0 ≤ t ≤ tf
  1. 取点
  1. 利用pdepe求解
  1. 显示结果
上述方法可直接推广到方程组，只需把对应系数改为向量即可。注意其中 c(xtu,∂u/∂x) 为偏微分方程的系数对角矩阵。若某一对角线元素为0，则表示该偏微分方程为椭圆型偏微分方程，若正值(不为0)，则为拋物型偏微分方程。c的对角线元素一定不全为0。

4.1.2.2 二维状态空间的偏微分方程的 MATLAB 解法

MATLAB 工具箱可以解决下列类型的偏微分方程
  • 椭圆型偏微分方程
  • 抛物型偏微分方程
  • 双曲型偏微分方程
  • 特征值问题
  • 非线性椭圆偏微分方程
  • 方程组
根据不同方程类型可选择不同函数，如assempde（解决二次椭圆问题）、pdenonlin（非线性问题）、parabolic（解决抛物线问题）、hyperbolic（解决双曲型问题）等。
一个典型求解过程如下：

4.1.2.3 偏微分方程的 PDE Modeler 解法

对于一般的区域，任意边界条件的偏微分方程，我们可以利用 MATLAB 中 PDE Modeler 提供的偏微分方程用户图形界面解法。
图形界面解法步骤大致上为：
  1. 定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。
  1. 产生离散化点，并将原 PDE 方程离散化。
  1. 利用有限元素法(finite element method；FEM)求解并显示答案。
具体步骤为：
  1. 利用绘图(draw)模式，定义欲解问题的空间范围(domain)Ω 。
  1. 利用 boundary 模式，指定边界条件。
  1. 利用 PDE 模式，指定 PDE 系数，即输入 c，a，f 和 d 等 PDE 模式中的系数。
  1. 在 mesh 模式下，产生 mesh 点，以便将原问题离散化。
  1. 在 solve 模式下，求解。
  1. 最后，在 Plot 模式下，显示答案。
解示例：
notion image
notion image
PDE Modeler适用于二维问题的求解。但一维问题同样可以通过画出条状图转化求解。

4.1.3 利用Python求解

Python求解思路与MATLAB中差分解法思路一致。
用Python数值求解偏微分方程 - 知乎 (zhihu.com)
之后利用matplotlib进行可视化即可。

4.2 其他特殊软件

暂无了解，查阅资料得知lingo对于偏微分方程这种大型问题能力有限，故主要考虑使用功能成熟的MATLAB进行计算。

修改记录

  • 2021-08-20，杨伟程校对语法
  • 2021-08-17，陈冠宇、许琳怡、汪蕊完成写作
  • 2021-08-10，郑鸿晓更新模板
  • 2021-08-04，张嘉乐建立模板
上一篇
海明校验的公式推导及讨论
下一篇
数学建模のwiki
评论
Loading...