本文参考的是来自mooc上北京师范大学彭芳麟老师的计算物理基础
基础知识
偏微分方程的三种类型
椭圆型
初始条件：无
抛物型
初始条件：初始温度分布
双曲型
初始条件：初始位移与初始速度
边界条件
Dirchlet边界条件
区域边界的函数值
Neumann边界条件
给出边界上函数的法向导数
混合边界条件
给出边界上函数及其法向导数的线性组合
差分法解热传导方程
热传导方程： 
令 
可以获得显式公式
由此递推公式可以得出下列矩阵
例1
求解程序
x
同时上述的传导方程也有隐式公式
利用 
 得到 
由此递推公式可以得到
引入算符 
则上述矩阵可以化为 
例2 解无量纲化后的薛定谔方程
差分法解弦振动方程
波动方程： 
则可以得到相应的显式差分公式
其中
 解是稳定的
 可能能得到正确数值解
 解是不稳定的
由此递推公式可以得到对应的矩阵形式
初始条件的设定
例3 两端固定的弦振动
求解程序
figure
用差分法解椭圆型方程
差分法：
令 
则可以得到显式差分公式
稳定条件为 
第一类边界条件
中心点用 
 表示，边界点用 
 表示
第二类边界条件
迭代法与松弛法
迭代法解线性方程组
矩阵解法
迭代法
x1
雅可比迭代法
高斯-赛德尔迭代法
松弛法
令 
则有 
为提高运算效率，可以加上松弛权重 
 ，则有 
 时为低松弛， 
 时为超松弛
松弛法迭代公式
启动计算：所有内部点都用边界点的平均值作为启动值
pdetool求解偏微分方程
pdetool中方程的输入格式
边界条件格式
可解问题的分类
解题步骤
设置定解问题
Draw Mode 画求解区域如矩形，椭圆，多边形及其组合
Boundary Mode 定义边界条件
PDE Mode 定义偏微分方程，即给定方程的类型及其系数
解方程
Mesh Mode 将区域分割为三角形网格
Solve Mode 设置初始条件并求解，本征值问题可设搜索本征值范围
将结果可视化输出
Plot Mode
用彩图、高度图、矢量场图、曲面图、网线图等直线图和线头图表现解
对抛物型方程和双曲型方程，可以用动画表现解
用pdetool解椭圆形方程和抛物型方程
用pdetool解波动方程和本征值方程
特殊函数的调用和计算
在matlab中查询特殊函数的方法
help matlabspecfun
legendre
例：计算n阶勒让德函数 
 在x处的值
勒让德函数的值
画勒让德多项式的图像
sym 

	

01
PDE Modeler使用方法介绍
物理学中的偏微分方程(PDE)无处不在，如热传导方程、扩散方程、电磁场方程，甚至量子力学中也能大量遇到偏微分方程——薛定谔方程。偏微分方程的计算十分复杂，而且大部分是没有解析解的。PDE的数值解法有有限差分法等方法，然而这是需要学习相关的知识，通过编程进行求解。而PDE Modeler提供了一个不需要编程就能接偏微分方程的平台。
1.1界面功能介绍
MATLAB提供了很多工具箱，其中PDE Modeler就是一个用来解偏微分方程的工具箱。
图1-1 PDE Modeler的位置
打开后的最初界面是这样的：
图1-2 PDE Modeler的界面
首先映入眼帘的是中间一大块空白，这个是展示方程求解的区域及最后解的分布的。该工具箱可以解不同类型的偏微分方程，也能解一些本征值问题。类型可以自主选择。
图1-3 可求解的问题
从上到下的翻译：
通用标量
通用系统
平面应力
平面应变
静电学
静磁学
交流电源电磁学
导电介质(直流电)
热传导
扩散
界面左上角是功能的核心区。解偏微分方程问题只需要熟悉这些模块的内容。
图1-4 求解功能区
前五个几何图形的图标是用来画求解的边界的。PDE就是确定方程，即输入方程中的一些参数。紧接着两个三角形图标是划分有限元。等号代表求解。后两个依次是画图和放大。
1.2使用步骤
接下来我们通过简单的例子说明一下具体使用步骤。
Step1 选择需要求解的问题
以静电学为例。
图太多了，懒得插图注了
Step2 确定边界
画出一个圆心为(0,0)，半径为1的圆。
Step3 点击偏导符号
点击之后会出现边界。
Step4 输入边界条件
直接点击边界即可输入边界条件，边界条件有两种。
Neumann边界条件
Dirchlet边界条件
本例使用的是Dirchlet边界条件 
Step5 输入PDE
本文以点电荷产生的电势为例。
 ，则  ，其中  。
Step6 差分
点击三角形图标，第一个是粗略划分，第二个是精细划分。之后区域内会划分为很多小三角形。
Step7 解方程
点击等号得出方程的解，图中会显示解情况。到这步为止，方程已经解出。该偏微分方程的解为点电荷的电势分布。
Step8 画图功能
该工具箱还能画不同类型的图，如等高线图、梯度图、三维图、动态图等，点击菜单栏的Plot——Parameters可以进入，具体功能大家可以自行探索，也就是每个勾都选一选试试有什么功能。
02
以数模2018A题为例进行演示
接下来我以数学建模-2018A题的第一小题为例进行演示，顺便补充一些细节。需要注意的是，该方法可能并不适用当年的比赛，因为该方法无法进行遍历，因此难以求出最优解。
本文为解出数值解提供了一种方法，但目的是为了让大家学习pde工具箱的用法，而不是完美地解决18年的数学建模问题。因为这种方法一次只能解出一个，并不能写循环、遍历(当然，可以自己一个个输入，手动循环)。
第一小题叙述如下：
“在高温环境下工作时，人们需要穿着专用服装以避免灼伤。专用服装通常由三层织物材料构成，记为I、II、III层，其中I层与外界环境接触，III层与皮肤之间还存在空隙，将此空隙记为IV层。为设计专用服装，将体内温度控制在37ºC的假人放置在实验室的高温环境中，测量假人皮肤外侧的温度。第一题：专用服装材料的某些参数值由附件1给出，对环境温度为75ºC、II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟的情形开展实验，测量得到假人皮肤外侧的温度(见附件2)。建立数学模型，计算温度分布，并生成温度分布的Excel文件(文件名为problem1.xlsx)。2.1建立偏微分方程组
这题给我的感觉是：建模极其简单，计算怀疑人生。
我相信任何一个学过偏微分方程的人(尤其是物理专业)，都能很快地列出这组偏微分方程。然而解析解是不可能解出来的了，但数值解并不好算。
偏微分方程：
衔接条件：
初始条件：
边界条件：
2.2利用PDE Modeler求解
画出防热服的各层
输入边界条件
由于原题是一维问题，而这个工具箱只能解决二维问题，所以我们画成条状的图形。
同时，两侧需要将边界条件设置为0，才能达到一维解的效果。
工作服外侧与内侧边界条件如下图：
设置时间及初始条件
点击菜单栏的Solve——Parameters进入，可输入时间和初始条件。
解出方程
2.3误差分析
点击菜单栏的Solve——Export可导出解的数值。
我们画出四个分界面的温度变化：
将解出的皮肤温度变化与题目已知的皮肤温度变画做比较，这里计算了差的绝对值。结果如下图，最高差2℃，效果还行。
如果想保存结果，可以点击File——Save As保存为m格式文件。我解的方程的m文件代码在附录中给出。
内容来源：知乎同名文章，作者：拉格朗日L2附录 MATLAB代码：欢迎加入QQ交流群大学生数学竞赛与考研群：645654483数学建模交流群：875120648ACM、CCPC交流群：740862242Python数据分析交流群：8181495302021保研交流群：1140271992

文章目录
前言
调用示例
例题
求解
命令介绍
具体实现
步骤1：化标准式
步骤 2：编写偏微分方程的系数向量函数
步骤3：编写初始条件函数
步骤 4：编写边界条件函数
步骤 5： 取点
主程序

前言
在python3安装fipy失败之后，懒得下载python2的我还是选择了matlab。
调用示例
例题

求解

命令介绍



具体实现
步骤1：化标准式



步骤 2：编写偏微分方程的系数向量函数
function[c,f,s]=pdefun(x,t,u,ux)  %建立偏微分方程函数
c=[1;1];
y=u(1)-u(2);
F=exp(5.73*y)-exp(-11.46*y);
s=F*[-1;1];
f=[0.024*ux(1);0.017*ux(2)];

步骤3：编写初始条件函数
function[u0]=pdeic(x)   %建立偏微分方程的初始条件函数
u0=[1;0];

步骤 4：编写边界条件函数
function[pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)   %建立偏微分方程的边界条件函数
pa=[0;ua(2)];qa=[1;0];
pb=[ub(1)-1;0];qb=[0;1];

步骤 5： 取点
由于此问题的端点均受边界条件的限制， 且时间t 很小时状态的变动很大(由多次求

解后的经验得知)，故在两端点处的点可稍微密集些。同时对于t 小处亦可取密一些。例

如，
x=[0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t=[0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

主程序
以上几个主要步骤编写完成后，事实上就可直接完成主程序来求解。此问题的参考

程序如下：
x=0:0.05:1;
t=0:0.05:2;
m=0;

sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
u1=sol(:,:,1);
u2=sol(:,:,2);

figure;
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure;
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')

之前介绍的 Physics Informed Deep Learning 解偏微分方程时都是使用深度神经网络得到一个参数化的解，然后通过最小化解对初始条件，边界条件以及偏微分方程的破坏程度，训练学习得到近似解的参数。缺点是对每一个初始条件都要重新训练一次神经网络。
今天介绍一个使用神经网络解偏微分方程的新思路。这种方法训练一次神经网络之后，可能适用于不同的初始条件。具体方法是对给定的偏微分方程，学习最好的差分近似离散化方案。此离散化方案依赖于具体的偏微分方程，时空格点以及场量的局域变化曲率。
Learning data driven discretizations for partial differential equations​arxiv.org
结果
专业术语的解释放在方法部分。
文章图（a) 在格子16倍粗粒化时，1阶迎风(1st order)算法已经发散；放大32倍时，WENO算法在第20个时间步也开始发散，而神经网络提供的差分离散化方案一直到第40个时间步都工作正常。
文章图 (b) 对比了7个区域的精确解与神经网络解在40个时间步长上的演化，结果还算比较精确。横轴是时间。
图(c), 定量对比了不同算法的平均绝对误差，发现在32倍粗粒化时，传统方法的误差是神经网络解的 10 倍。
方法
好的离散化方案应该类似于重整化群，在更粗粒度的离散化方案中得到细粒度下物理的有效近似。在流体力学里面，这相当于将高频信号和局部涨落积分，得到系统的长波近似。在光学中，当系统尺度远大于光的波长，几何光学也是对麦克斯韦方程组的很好近似。所以，如果神经网络能够在粗粒度下学到对底层微观动力学的有效近似，可以解决很多问题。
多体系统的涌现（集体）现象需要微观机制，但如果从第一原理即微观机制出发，模拟整个物理系统，则会受到很多非相关细节的干扰。如何构造粗粒度下的有效理论近似微观机制的集体效应，可能推动整个物理学领域的发展。重整化群即是这样一种思想。本文立意很高，试图使用神经网络学到粗粒化的差分离散化方案，有效近似微观动力学带来的宏观效应。
差分近似入门
如果我们要数值求解下面这个偏微分方程（假设v是常数），该怎么做？
这个方程告诉我们，流体密度 
 是时间 t 和空间坐标 x 的函数。
如果将 
 时刻流体密度的空间位置分布 
 定义为初始条件，那么偏微分方程的解就是 
 的任意时刻流体密度的空间位置分布 
 。
如果将 x 离散化为 
 ，空间步长为 
 ，将时间离散化，时间步长取为 
 , 则连续的偏微分方程可以化为求解这些离散时空点上的流体密度。此时可以定义，
相应的就有
根据泰勒展开，
可以知道 
对 x 的偏微分用同样的方法化为差分，得到这一章最开头的偏微分方程最简单的离散化形式为，
简化一下，将 n+1 时刻需要计算的量放左边，将 n 时刻所有已知的量移到右边，得到
初始条件 
 知道，从 n 到 n+1 时刻的迭代方程知道，此偏微分方程的数值解也就完全确定了。
差分近似进阶
可能你会觉得差分近似太简单了，一看就会。但现实情况是一用就错。
比如上面这个差分近似，如果常数 
 , 那么很快会出现数值发散！如果 
 ，又可正常迭代。对于 
 , 必须选用迎风格式 (upwind), 即
注意上节末尾为了计算 i 点的空间微分，使用了 (i+1, i) 两个格点的信息；而迎风格式使用了 (i-1, i) 两个格点的信息。就这一点点不同，导致迎风格式收敛，而另一个发散。最直观的理解是，当 v>0 时，信息从左往右传播，计算 i 点的微分使用 i+1 点的流体密度违背了因果律。
为何  (i-1, i) 两点也能用来计算空间的偏微分呢？还是泰勒展开，
同样可以得到 
这里的迎风格式还是一阶精度，即近似到 
 ，
如何得到更高精度的差分近似呢？泰勒展开到更高阶，联立求解便是！
都展开到二阶的时候，两个相减则得到，
这种差分方案称作中心差分，具有二阶精度。
将两者相加，消掉 f 的一阶偏导，则可以得到二阶微分的差分形式，
同理，对 
 ,  
 ,  
 ,  
 泰勒展开到更高阶，并联立求解，可以得到更高阶精度的一阶微分，或更高阶微分的差分离散化方案。
微分的差分近似很像卷积核，比如这个二阶微分，化为卷积操作则系数为，
将这个卷积核从左到右点乘到 
 上，就能得到每个空间格点的二阶微分。
有限体积法
有限体积法使用每个格子内部场量的积分平均值来代替格子边界上的值，更好的保证物理量的守恒。本质上与有限差分法没有太大区别。
本文做法
上文所示差分离散化方案一般是人为指定的，在数值求解偏微分方程时不随方程和时空的不同区域发生变化。
这样就有一些潜在的问题，比如低阶差分近似一般会带来大的耗散，即 
 对全空间的积分不守恒，演化时间越长，因为耗散丢失的质量越多。
高阶差分耗散小，但又会引入大的色散。一个小的数值涨落，随时间的演化不是消失，而是逐渐增大，使得演化发散。
近代计算流体力学发展了各种各样的有限差分和有限体积法，比如 FCT, Godunov, TVD, NND，ENO, WENO, COMPACT 算法等等，都是为了解决低阶与高阶，耗散与振荡的矛盾。
ENO 和 WENO 算法会根据局部的曲率，从一系列预先定义好的差分形式中选择最适合的离散核 
。本文使用的方法类似 ENO 与 WENO 的推广，使用多层神经网络来输出n阶导的差分近似系数 
 。这些系数依赖于所处的时间和空间位置，以及局部场量的数值（比如曲率）。
对空间的不同区域使用不同的差分格式非常合理，比如说在一个激波的内部和外部，应该使用不同精度的差分近似及差分系数。
机器学习和数据驱动的方法如下：首先生成高分辨率的训练数据，然后从数据中学习微分的离散近似。产生高分辨率的训练数据可能花费很高，但是如果在小的系统生成解在流形空间的局域近似，然后应用于大的系统求解，可能会显著减少计算时间。
在有限积分格式下，Burgers 方程完全受每个cell边界上流入，流出的场量决定。所以唯一的挑战是根据cell的平均值精确的估计边界上的流（Flux）。这一步由神经网络给出。
计算分成3个步骤：
使用神经网络构造Cell边界上的空间微分，（利用边界两边Cell的场量均值）。
使用计算得到的空间微分去计算 Flux J。
使用旧的 Cell 均值减去左右边界上向外的 Flux J 来计算 Cell 均值的时间微分。并将时间微分（而非空间微分）的预测误差作为loss function进行优化。
作者尝试过使用传统的单调算法（此处单调指的是没有色散引起的振荡），比如 Godunov Flux，没有带来比神经网络模型更好的提升。
这篇文章提出本方法的缺陷是：（1）速度比普通差分近似慢很多。普通差分的卷积核一般用到左右5个格子，而本方法神经网络参数约几千个。（2）本文方法是在规则格点上的差分方案，如何推广到不规则格点呢？图神经网络与点云神经网络显然是很好的候选。这个需要进一步的研究。
展望
如果神经网络能够根据偏微分方程，边界条件和初始条件生成不规则的离散化格子，在曲率大的区域构造更精细的格子，在平坦区域构造粗粒度格子，并在不规则格子上构建最优的差分离散化方案，可能会带来计算流体力学的新突破。
代码：
旧版 https://github.com/google/data-driven-discretization-1d 
新版 https://github.com/google-research/data-driven-pdes  
新版开源代码中有 Tutorial 和 Example