大家好！
趁热打铁，在介绍完常微分方程的那一套理论后，我们继续介绍偏微分方程的相关内容。
学过偏微分方程的都知道，偏微分方程一般有三类：椭圆方程，抛物方程和双曲方程。它们每个方程都各有自己的一套理论，要想完全了解相关性质和理论是极为困难的。但是有幸在数值解中，我们可以看到这些方程都可以用一套框架来解决。这也是这篇文章要重点关注的内容。
不过在我们进行偏微分方程的话题之前，我们还会再引入一些之前常微分方程还没有说完的重要的内容和方法。因为最近我还没有更新有限元方法的计划，所以这也是专题的最近的最后一篇文章。
提供之前的笔记：
微分方程数值解法|专题（1）——常微分方程的有限差分方法
我们开始本节的内容。
目录
线性多步法的绝对稳定区
抛物方程的有限差分方法
冯诺依曼分析法
Crank-Nicolson方法
双曲方程的有限差分方法
椭圆方程的有限差分方法
五点格式法
线性多步法的绝对稳定区
我们还是研究方程 
 。无论是最简单的向前Euler法，还是之后的龙格库塔方法，我们的导数项 
 都是没有变的。我们之前也说了它们是显式格式。这是因为我们在计算的时候，
右端的导数项的时间点都在过去。实际情况下，右边的导数项是可以变的，比方说向后Euler法
你会发现，导数项的时间是 
 ，所以你会发现，这个数列通项没有办法再那么容易得到了，因为你在计算时间点 
 的时候，所依赖的导数值与 
 有关。这个时候在进行编程的时候，一般我们采用
迭代法。也就是说，我们把左右两端的 
 都当作未知数 
 ，然后做类似于这样的迭代
然后做出一个关于 
 的收敛数列，不过这不是我们重点关注的内容。
既然我们可以把它改成 
 ，那么我自然可以改成之前某些时间点的线性组合。所以我们假设 
 ，来导出线性多步法的一般形式
设 
 ，我们可以得到我们的一般形式
它对应的特征多项式为 
 ，所以我们只需要考虑让它的根满足一定条件，即可找到它们的绝对稳定区。下面我们来举一些实际的例子，来说明一下如何操作。
Example 1: Forward Euler Scheme
两边写出它的特征多项式，可以得到 
 ，一个一元一次方程（注意这里的 
 不是变量，因为这个方程只是关于 
 的），所以为了要求满足根的条件，我们就要求 
 ，这与我们上一节推的结果完全一致。
Example 2: Backward Euler Scheme
两边化为特征多项式，可以得到 
 ，解 
 可以得到 
 ，那么你也可以看到，只要求 
 就可以，这个式子对于 
 的时候是恒成立的，所以向后Euler其实是无条件稳定的。
Example 3: Trapezoidal Scheme
两边化为特征多项式，有 
 ，解 
 可以得到 
 。学过复变的话，你应该知道这个式子其实是一个保形映射。这样的话，可以得到，如果 
 ，那么 
 ，所以这个的绝对稳定区就是负半平面。
这里只是一些比较简单的例子，但是实际情况下的绝对稳定区的推导其实是相当复杂的。LeVeque的书上有一些篇幅，感兴趣的朋友们可以去查看一下。
抛物方程的有限差分方法
我们研究的就是下面的热方程
学过热方程的都知道，它本质上是一个衡量杆的温度的方程。当然了它本质上其实是一个初边值问题，直观来看是因为，它既需要初值条件，又需要边值条件。
下面我们来看看如何差分吧。在之前我们做过对导数的离散，那么这里 
 相当于对方程针对时间变量的导数做离散，关键在于 
 ，它是一个二阶导数，那么如何去做呢？
事实上方法是完全类似的，在之前我们用向前Euler法来进行近似的时候，其实本质上是下面这个近似。
（注意我们这里因为方程的符号改成了 
 ，所以自然统一的用 
 来表示了）如果将步长标记为单位长度，那么就相当于
注意这里的 
 是
向后差分符号，不是梯度。所以你也能猜到，如果我们要估计 
 ，结果是完全类似的，也就是下面的过程
我们再把步长 
 代回去，就可以得到我们最后的差分格式
其中你也可以看出，这个差分格式因为是对空间变量差分的，所以时间没有变，而空间变量中 
 就对应着上面的 
 。
你也可以看出，从这里开始，因为方程开始有多个变量，所以我们不能够再使用之前的简单的下标标记了。这里的 
 的离散也是容易的，就是下面的结果。
因为我们在空间上设置了离散步长为 
 ，所以我们这里在时间上就用另外一个字母 
 来表示。那么你也容易看出，分别用差分格式代到两边，就能得到我们最终的表达式。而且使用二维泰勒展开，我们也容易得到我们的局部截断误差，分别是 
 和 
 ，其中 
 ，潜台词就是说，它的局部截断误差就是 
 。
当然了，上一节说了那么多稳定性，我们这一节自然也不能错过。我们定义 
 就是在 
 处的
数值解，那么自然的，我们可以将它按照时间的顺序排列，得到下面的一个形式
这里，我们假设 
 ，那么上面的式子就是
你也可以看出来，当我只考虑时间的变化后，固定时间，再考虑空间，其实我们也可以得到一个关于空间变量的数列通项公式。因为我们的变量变成了两个，所以这样的一个通项就不能够轻松的使用特征多项式分析了，但是还好我们还是有招，这就是后面所使用的冯诺依曼分析法。
冯诺依曼 (Von Neumann) 分析法
虽然我们在抛物方程这里引入了它，但实际上它可以应用在各种常见的偏微分方程中。所谓的冯诺依曼分析法，其本质就是用向量，将所有的空间离散点拼起来。那么这样的话，这个矩阵所对应的表达式就仅仅与时间有关了。这样自然分析就会方便很多。
回到上面这个问题，我们针对这个差分格式，就可以化成这样的矩阵。
写的紧凑一点，就是公式 
 。那么为了让这个公式收敛，我们显然需要让 
 ，其中 
 表示它的最大特征值。
这里的 
 你也可以看到，是一个
三对角矩阵，它是数值线性代数里非常重要的一类矩阵。这里我们不加证明的给出一个特殊矩阵的特征值。
Proposition 1:
对于矩阵 
 ，它的特征值为
事实上，这个矩阵我们可以做一个拆分，得到 
 ，其中 
 ，根据高等代数里有关特征值的结论，容易得到 
 的特征值就是 
 ，对，就是少了一个 
 而已。
回到我们这个例子，你也容易通过变换得到 
 ，所以它的特征值就是 
 。
通过对余弦的估计，你可以知道， 
 ，这样的话根据 
 ，就能得到我们最后的结果 
 。注意
我们一般不取 
 ，这是因为数值算法往往存在误差。
下面我们再介绍一种向后差分方法。我们不再给出推导细节，直接给出它的差分格式
其实变化就是在左边，把它改成了一种向后差分的格式（注意这和之前的向后Euler不是一回事）。那么也容易通过化简得到我们的矩阵方程。
写的紧凑一些，其实就是 
 。那么我们只需要找到保证 
 即可。这是显然的，因为用和上面一模一样的方法，就可以得到特征值为 
 ，所以这个格式是无条件稳定的。
Crank-Nicolson方法
因为之前的两个方法它们的误差是 
 ，这其实会有问题就是主要的误差其实都在时间步长上。我们需要把时间步长 
 取得很小，可能才能够达到我们要的精度。因此Crank-Nicolson方法可以很好地弥补这个缺陷。
这个方法本质上就是用不同的差分方法来近似 
 和 
 。这里我们的 
 使用向后差分，和上面那一个一样。而我们的 
 使用的是混合差分。具体来说，就是用
代替 
 。如何考虑它的稳定性分析呢？如果我们还是一样，按照时间的顺序去排列，可以得到下面的形式
这个方程矩阵化会稍微麻烦一些的，大体上的形式就是 
 ，其中 
 ，且
 。那么现在问题自然就是落到了矩阵 
 的特征值上。
首先，我们把矩阵做一些拆解，容易看出 
 ， 
 。那么这样的话，考虑设 
 为 
 的特征值，并且对应一个特征向量 
 ，那么这就容易得到
然后呢，如果我们有 
 ，那么就一定有 
 。所以根据这个结论，我们容易得到最后的 
 的每一个特征值为
容易验证，这里每一个数其实都是严格小于1的，也就是说格式是无条件稳定的。
但是如果单纯说稳定性，它并没有什么优势。我们通过对它做局部截断误差的误差分析，就可以看出来更多细节。我们这里详细的来走一遍。这里不失一般性，我们假设 
 。
首先还是考虑，使用泰勒展开，首先容易得到的是
（我希望你没有忘记， 
 对应的是 
 的时间点）然后呢，我们再考虑两个中心差分，有
当然只要改一个时间点，就可以得到
如果我们略去最高阶的误差项，那么也能够看出，右边最最关键的部分就是 
 ，而左边的东西就是 
 。我们还是抓住我们的时间和空间都在 
 点展开，就可以得到
剩下的就是拼在一起的事情了，这个证明已经没什么难度了，因为 
这就完成了分析，因为你可以看到， 
 那一项已经被消掉了，所以这个误差是 
 ，这已经比上面的误差要好很多了。
事实上，如果你把左边的时间变量的差分离散改成中心差分，就变成了Richardson格式。这是一个无条件不稳定的格式，它的具体分析就是19年丘成桐杯竞赛的应用数学个人赛的第二题。据北大小伙伴说，如果做对3个题，排名据说就能到银牌……
双曲方程的有限差分方法
这个时候，我们考虑的方程为
怎么差分，我觉得你已经有数了。两个变量都做中心差分即可，也就是格式写成
按照时间顺序排一下，可以得到方程的形式为
这里的 
 。
当你把这个方程写出来之后你可能就发现情况有点不太对了。首先就是初值的情况要单独拉出来，其次就是这个方程即使是按照时间排序，也会与之前的两步时间量有关。那应该怎么办呢？
我们一个一个来解决，首先初值的情况，就是这里的 
 的情况，你知道 
 对应的就是初值的时候，这个值已经由初值条件控制住了。但是 
 你会发现它会与一个 
 的量有关系，这是个什么玩意？也是因为这个，我们需要给出一个估计，一个关于 
 的估计。比较好办的是我们有了关于 
 的导数的表达式，所以我们可以考虑构造一个中心差分格式来近似导数。具体点来说，我们的估计是
代入，合并就可以得到我们的初值方程
事实上，它也可以被写成是一个矩阵形式。而之后的矩阵方程都按照那个格式来写即可，最后我们可以得到下面的矩阵方程组。
如果我们把一般的这个方程拼的紧凑一些，我们就可以得到它的一般情形其实是诸如 
 的，所以核心还是在前面。这是一个多步的数列，它的解决方案类似于高次的常微分方程——
化为方程组。具体点来说，我们应该把方程组写成下面这个形式。
这样我们只需要考察另一个矩阵 
 的特征值就可以了。
用传统的数值线性代数的思路当然没有问题，但是这里我们换一种分析方法。我们考虑设 
 为 
 的特征值，并且设 
 为对应的特征向量对，那么直接走矩阵运算，你能够得到， 
 是 
 的特征值。而 
 的特征值我们用完全一样的思路可以知道它应该是落在 
 这个区间内。那么可以想一想，如果 
 ，那么对应的 
 是完全可以取到一个 
 的值，这就没有办法使得数值解稳定。所以我们让 
 ，因为这个时候，可以发现 
 ，才有可能 
 取到我们的那个区间内（其实也不好说，因为两边习惯上是开区间，但是因为数值误差，写成闭区间倒也无妨）。虽然这个结果不是特别好，但是至少误差不会呈指数增大，所以我们接受这个结果。
通过与上面相同的误差分析，你可以得到，这个格式是 
 的误差。
椭圆方程的有限差分方法
其实你也可以看出，上面的双曲方程和抛物方程的稳定性分析方法都是使用冯诺依曼分析法的，可以说是极为方便的。但是椭圆方程的分析并不是这样，它的方法有些另类。因此我们以它作为我们最后的结束的内容。
首先我们研究的方程一般长这个样子
这里的 
 是拉普拉斯算子，
不是差分符号。考虑一下更特殊的情形，就是
那么如何作差分呢？其实这个地方方法倒是类似的。不过呢我们需要给一些比较有趣的名字。
五点格式法
我们还是一样考虑以 
 来表示 
 处的数值解。这里因为我们有一个额外的函数 
 ，所以我们也需要给它一个记号 
 ，并且同样的我们也假设 
 方向的步长为 
 ， 
 方向的步长为 
 。有了这些假设，我们的差分方法其实就可以写成这个形式
因为两个方向其实都是空间变量，我们有理由相信它们的性质比较类似。基于这个考虑我们不妨假设 
 ，那么组合化简可以得到
下面这张图一定程度上解释了为什么我们称这样的差分方法为“五点格式法”
也许你会想，这个表达式写出来之后，我们也确确实实可以考虑用冯诺依曼分析法。但是问题在于，无论怎么排布，我们都没有办法获得诸如双曲和抛物方程那样好的结构。因此这边提出了一个非常有趣的方案叫做堆积 (stacking)，通俗说就是把所有的，每个点上的数值解排成一个列向量。比方说我们每一个维度上取了 
 个点，那么一共就会有 
 个点，向量自然也就是 
 维的了。另外排布的方式我们自然也有讲究，具体来说就是这样
 ， 
按照我们这样的排列，容易发现我们的矩阵应该长这样
你也可以看出，这实际上是一个极为稀疏的大矩阵。但是还好，它的结构还算完整。而对应的方程其实就是 
 ，注意我们这里的结果因为不再与时间有关，所以实际上它不再是一个迭代的数列，这一点和前面两个有本质上的不同。
下面，我们来做一些稳定性分析吧。首先你也能看出来，因为没有时间变量，我们这里自然不能通过上面的计算谱半径的方法来判断是否稳定。所以这里实际上问题就落在了我们怎么处理这单独的一个矩阵方程上。
最直观的想法就是：考虑构造另外一个方程 
 ，只是说我们这里的 
换成对应时间点的精确解。这当然是可以的，不过这样子得到的结果，后面计算的实际上就是精确解代入得到的一个式子。那么如果我们相减两式，得到 
 ，那么左边的 
 就是我们的每一个点的
整体截断误差。但是右边就需要我们单独计算了。
我们具体写出来右边式子的每一个元素，就是
相信你也知道该怎么做了吧？这就是最经典的局部截断误差的分析，它最后的结果是
你可以看到它的局部截断误差是二阶的。
把这些东西都放在右边，就可以得到我们最后的误差矩阵方程 
 。容易看出，只要 
 有界，那这个值 
 ，它的规模应该是和 
 内的数相同的（想想为什么），所以我们下面就需要检查一个事情：是否 
 有界。而一般来说根据范数的等价性，只要找到一个范数说明一下有界性即可。
我们还是考虑二范数，原因也很简单，这样子更容易分析。
首先这个矩阵是 
 维的，因此也就存在 
 个特征值。我们不加证明的给出它的特征值形式。
注意我们这里的 
 相当于对应的 
 坐标。比方说 
 就相当于是 
 方向都在左（下）边界的对应的特征值。也就是对应 
 这个位置的特征值。
那么你也看出来，这个方程的分析不难的，因为我们有 
  ，这样的话可以看出，如果要证明 
 有界，那么就一定是找到原来的那个矩阵 
 中，
绝对值最小的那个做分析。这里我们自然取 
 。容易得到
这是一个有界的数，所以我们自然可以得到结论说这个格式是稳定的。
最后，我们简单的提一下九点格式法。这个方法写出来是长这样的
这个方法你也可以证明，它可以具有四阶的精度，不过需要要求 
 ，具体的细节我们碍于篇幅，就不在这里说了。
小结
我们用很长的篇幅，给大家完整介绍了偏微分方程的有限差分的数值解法。大家可以看出来，虽然说偏微分方程是一个庞大的理论，但是在数值解领域，是存在这样的方法，可以系统的来解决很多很多问题的。这里主要就体现在贯穿全文的泰勒展开和冯诺依曼分析法中。
事实上，偏微分方程的有限元方法其实是更有活力的一大块。它们的构造其实更有意思，也更加实用。但是我们近期自然是没有机会说这么多了，就希望大家能够在这两篇文章中，详细而系统的窥探到有限差分方法的端倪。
——————————————————————————————————————
本专栏为我的个人专栏，也是我学习笔记的主要生产地。任何笔记都具有著作权，不可随意转载和剽窃。
个人微信公众号：cha-diary，你可以通过它来有效的快速的获得最新文章更新的通知。
专栏目录：笔记专栏|目录
想要更多方面的知识分享吗？欢迎关注专栏：一个大学生的日常笔记。我鼓励和我相似的同志们投稿于此，增加专栏的多元性，让更多相似的求知者受益~

题主想问的是常微分方程(ODE)和偏微分方程(PDE)的数值方法区别呢还是微分方程这个领域和微分方程数值解领域的区别呢？按照前面@赵永峰 的回答，我也按照前者理解吧。毕竟后者的一些区别是显而易见的。
先说一点共性。微分方程的数值方法，无论是ODE还是PDE，都是将连续的、无限未知数的问题近似为离散的、有限未知数的问题求解。从经典数值分析的角度，通常会关心下面一些问题：相容性、稳定性、收敛性、收敛阶、计算量等等。相容性是指格式在局部是不是做出了正确的近似；稳定性是说局部的近似误差会不会随着计算而积累放大；收敛性是说当离散尺度无穷小的时候数值解是否会趋向于真实解；收敛阶则刻画了收敛的速度，高阶的格式可以用较大的离散尺度获得较好的数值结果，但是代价通常是单步下稍多的计算量。因此数值方法的最终表现需要在误差和计算量之间找到一个平衡。
先说说ODE。在这个领域里，无论是初值问题还是边值问题，有限差分方法都是最常用的方法，比如说著名的Runge-Kutta方法。最常用的RK4方法就有稳定性条件比较宽泛、收敛阶很高(4阶)、计算量较小的优点。ODE数值方法中，差分方法是绝对的主流。尽管有限元方法、谱方法等等也可以用于解ODE，但是差分法依然更受欢迎。即便是边值问题，基于差分法的打靶法也比有限元更受欢迎。
由于ODE的解行为通常比较好，只要右端项满足一定的Lipschitz连续性，解就存在唯一，对初值参数连续依赖。所以ODE数值方法的特点是有限差分法是一种适用面非常广泛的方法。也就是说，如果你是一个工程师，对数值方法并不熟悉。你在实际工作用需要求解一个(规模不太大的)ODE，那么你闭着眼睛把这个方程扔给一个RK4标准程序，效果一般不会太差……
实际应用中ODE数值方法面临的最主要问题是刚性。简单说，如果把方程组理解为一组粒子的运动，那么这些粒子的运动存在时间尺度的分离，而你的数值方法应该要抓住最小的时间尺度，这就意味着超大的计算量。这种问题在分子动力学模拟(MD)中特别常见。本来MD就要计算10^6量级的粒子，再有很强的刚性就会使得模拟几乎无法进行。实际中，无论是从理论上做渐近分析或是平均化(averaging)抑或是数值上构造稳定性条件更加宽松的数值格式都是非常有挑战性的工作。
ODE数值解面对的另一个困难时长时间模拟。再好的数值格式也会有误差，误差总会随着时间积累，时间充分长之后总会让数值解变得不可信。尤其是如果方程的解包含周期结构的时候数值误差很容易在长时间上破坏解的周期性(一个典型的例子是用Euler法求解地球轨道方程，数值解最终会远离太阳而去)。因此一个很有挑战性的问题就是如何在长时间的计算中保持数值解的某种结构，比如说能量守恒。如何构造这种满足特殊要求的数值格式同时还能尽量保持高精度是需要仔细设计的。实际中如果面对超大规模方程的长时间模拟，计算量的限制使得高阶格式都难以应用的时候，其结果的可信度基本属于玄学……
除此之外，ODE数值解还有一些具体的问题。比如说不适定问题的求解、方程在临近分岔时的精确求解等等。总的来说，ODE数值解的领域相对成熟，理论比较完善，有一些可以作为标准方法的解法。实际应用中，可以根据实际问题的特点在这些标准方法上做出改进。
说到PDE数值解，那简直就是天坑……这个领域太大了，即便你说PDE数值解就是全部的计算数学，错的也不算离谱。教授们如果不注意维护自己的个人主页，很容易发现一所高校计算数学系教授的研究兴趣都是偏微分方程数值解……
还是简单说几句好了。从方法构造上，前面@赵永峰 的答案中提到的有限差分法、有限元方法和谱方法确实是最主要的几种方法。有限差分法依然是最基础的。差分法有直观清楚、构造简单、易于编程的优点，对于没有受过专门数值方法训练的工程师来说，差分法依然是最好的选择。精心构造的差分方法可以非常高效。比如在求解流体力学方程的时候，守恒型差分格式有非常成熟的理论和方法。有限差分法的缺点主要是只能用于比较规则的区域，对于复杂区域边界的处理不但困难，而且很容易损失精度，进而影响数值解在全局的精度。
一种改进的方式是有限体积法(Finite Volume Method)。有限体积法的做法是将微分方程写成积分方程，在每一个小区域中用数值积分来近似精确积分，进而求解方程组。因为数值积分的方法比较灵活，有限体积法对于区域的要求宽松许多，并且可以选择合适的积分法来保持方程的物理性质。缺点则是如果使用较高阶的数值积分方法，那么计算量将非常大，甚至需要求解非线性方程组；而如果使用较低阶的数值积分法，又不如差分法简洁。
差分法的思想是在局部用差商代替微商，这是一个局部的近似。从全局看，差分法相当于用分片常数近似导数，也就是用分片线性函数近似精确解。而分片线性函数在全局其实是不可导的，所以我们通常在连续函数的最大值范数下来考察收敛性。而有限元方法(Finite Element Method)则是用分片多项式来近似精确解，我们不但可以在整体上考虑函数值的收敛性，还可以考虑导数的收敛性。有限元方法的优点在于可以用于不规则的一般区域，原则上可以构造出非常高阶的格式，收敛性和收敛阶有比较成熟的理论，缺点则是有限元的构造比较困难，也不容易写程序。在一些汉译文献中经常混淆有限体积法和有限元方法两个术语，需要特别注意。(一个特别有名的例子，LeVeque的名著“Finite Volume Methods for Hyperbolic Problems”就被翻译成了有限元方法……)
谱方法则是一种无网格方法。它不像差分法和有限元那样需要首先将区域做剖分，而是将解按照一组正交基做展开(也就是广义的Fourier展开)，截取有限项作为近似，需要求解的是对应的Fourier系数。谱方法的好处是高精度，以及搭配一些快速算法(比如快速Fourier变换)计算速度很快，缺点则是一般只适用于非常规则的区域，并且对边界条件有比较苛刻的限制。此外，将谱方法和有限元方法结合起来的谱元法也是当下比较热门的领域。
可以看到，和ODE不同，PDE数值解没有一种占绝对优势地位的方法乃至于框架。一般来说，我们需要针对不同的方程设计不同的数值方法。所以PDE方程的数值求解是一件技术含量比较高的事情。如果你是一个对数值方法不熟悉的工程师，在实际应用中需要求解一个PDE，那么最好还是找一本书简单学习一下。即便是最简单的方程、最简单的差分法，也需要一些知识来设计合适的格式(举两个学生作业中常见的例子，对流方程的差分格式需要满足CFL条件，对流占优的对流扩散方程也需要仔细设计格式来避免数值耗散对解的污染，即便这些方程都是常系数的)。
PDE数值解的困难主要在于PDE的解表现出的行为太丰富了。很多时候，我们对要求解的方程性质都缺少基本的认识，更说不上根据方程的特点设计有效的算法。实际中我们只能针对一类方程来设计一类格式，这一类格式对另一类方程很可能根本就不灵。我们都知道，
 和
 一个符号之差就是两种完全不一样的方程。适用于前者的格式根本就解不了后者。
ODE中我们提到的困难对于PDE都存在，比如刚性，比如长时间行为。但是这都不是PDE数值解的主要问题。因为PDE的数值解还远到不了讨论这么精细问题的程度，当务之急还是在有限的计算时间内解出来。对ODE数值解要求4、5阶的精度不算过分，但是PDE数值解能有时空2阶精度就非常令人满意了。
和ODE相比，PDE的数值解更加强调对方程物理性质的保持。因为PDE问题通常都来自物理背景。计算流体力学中要求保持物理量的守恒性，还要能够准确的捕捉激波。既要利用数值粘性来避免数值振荡，还要尽量减小数值粘性来保持解的守恒性。这些使得某一种PDE的数值求解都变成一门需要深入研究的学问。泛泛的谈PDE的数值解通常是谈不出什么来的。
PDE数值解的另一个巨大困难就是维数灾难(curse of dimensionality)。一般的说，PDE需要求解的未知数数量是随着问题维数指数增加的。这就意味着合理的计算量根本处理不了高维的问题。现今，无论是差分法、有限元还是谱方法，一般都只能处理三维以下的问题。超过三维，如果没有可以利用的对称性，基本可以宣告放弃了。然而高维的PDE求解在统计物理中随处可见。即便要求解Boltzmann方程，也是7维的，远远超出了传统方法的能力范围。
对于一类特殊的PDE，我们可以将它视作是某个随机变量的期望，然后利用Monte Carlo方法来计算这个期望。众所周知，Monte Carlo方法的优点就是计算量对维数的增加不敏感，可以针对少量特殊点求解方程而不必在全局解出整个解，可并行化程度高，是求解高维PDE的一种很有吸引力的方法。当然，Monte Carlo方法的缺点也很多。比如说收敛慢(通常只有半阶)、精度低、随机误差不可避免、对问题形式要求严苛等等。
总的来说，PDE的求解通常是根据具体问题设计具体方法的，泛泛地说PDE的数值方法很难深入下去。PDE求解的问题和困难非常之多，如果说解ODE的时候闭着眼睛上RK4是个不算糟糕的方案，那么解PDE就一定要对待求解的方程和数值方法理论本身都有基本的认识。

大家好！
最近运气可能比较好发现自己在丘成桐大学生数学竞赛（简称“丘赛”）的应用数学与计算数学（Applied and Computational Math）这个track下拿到了第12名，也就是压线入围的意思。这可是把我这个半吊子吓坏了哦……也是因为这个，学院安排了一些额外的学习任务（大雾），而这一系列的文章也是抓大放小，主要关注方法论和实用性。诸如微分方程数值解法，其实是有套路可循的，而掌握这一种套路往往就能够用来解决很多很多问题。所以这篇专题我们将关注数值解的相关方法，并给出相关的应用举例。
微分方程数值解法这一个学科需要的前置知识，除了数分高代以外，还需要数值分析和常微分方程，偏微分方程的一套理论，所以确实对于初学者来说已经完全不能算是“零门槛”了。所以这一系列专题，我们不会特别的关注理论（因为我不可能为了这个专门搬出一本数值分析）。相反，我相信直截了当的介绍方法，反而会让那些不是特别钻理论的应用者们有所受益。
我们follow的课本是Timothy, Sauer的Numerical Analysis（好吧，我还是搬出来了一本……），李荣华的《微分方程数值解法》，LeVeque的Numerical Methods of Partial Differential Equations 以及我们学院许传炬老师的课程材料。因此在任何单独一本教材上，你都无法找到这份笔记上的所有内容。但是也因为是专题的缘故，参考的书本较多，尽管我在写作前已经事先组织过，但确实可能逻辑线上不如查只查阅一本教材那么缜密。还请大家谅解。
废话就说到这里，我们开始正题吧。
目录
一阶常微分方程初值问题的有限差分方法与误差分析
向前Euler法及误差分析
梯形法
龙格库塔方法简介
稳定性
绝对稳定区
刚性常微分方程的有限差分方法
龙格库塔——切比雪夫方法
一阶常微分方程初值问题的有限差分方法与误差分析
在这里，我们主要分析的是下面这个方程
你可以看到，它的导数是一个函数 
 ，我们称它为“初值问题”的原因是，这个方程本质上描述的是一个物体在给定导数（也就是速度）后，随着时间的变化趋势。也就是说它主要是与时间变量有关。
所谓数值解法，本质上就是一种离散。我通过对时间步长作分割，通过一种特殊的数列迭代的公式，使得在我给定数列的前几项之后，通过数列的迭代公式得到我每一个分点上的数值解。因为很多时候，我们需要用到方程的解在某些点上的值，而传统的偏微分方程，很多时候都没有办法找出显式解，这就是数值解所起作用的地方了。
你也可以看到，数值解必然不可能和精确解完全相同的，会产生误差。我们下面给出与它有关的误差定义。首先很明显的就是整体截断误差，它一般定义为
Definition 1: Global Truncation Error
其中 
 就是我们的
精确解。而 
 就是给定计算格式（也就是数列的迭代公式）之后，你所计算出来的数值解。
第二个误差是局部截断误差，我们之后来解释为什么我们需要这个概念。它的定义是
Definition 2: Local Truncation Error
首先我们需要解释， 
 是下面这个初值问题的解。
你可以看到，这个初值的时间是从 
 开始的，它相当于是一小段时间内的 
 的表现。值得关注的是，它的初值设成了 
 ，也就是说它本质上的意思就是，当
我的初值精确解就是我的数值解的时候，我在下一个时间点的值应该满足的方程。所以我们这个时候假设的误差，其实是相当于假设之前的计算都是精确的情况下，下一步计算产生的误差。这也是“局部”的含义。一般情况下，“局部”带来的误差是容易估计的，这样的话，整体的误差就可以通过局部的误差来累加得到，所以引入这个局部截断误差，多是为了之后的误差分析方便。
有了这些误差的定义，我们就可以来分析相关的计算格式了。
向前Euler法及误差分析
我们先介绍向前Euler方法。在此之前，我们设时间是等步长的，并且步长为 
 。
Definition 3: Forward Euler Scheme
这样子直接看可能会糊，我们其实稍稍修改一下，可以得到下面的这个结果。
如果你把这些数值解都用精确解去代替的话，式子就会变成这样
你可以看到，当我令 
 的时候，右边其实就是导数的定义表达式。所以数值格式的本质，就是使用我的解的线性组合，去
估计导数。而如何估计导数其实方法各有不同，而这也自然产生了不同的格式。
它的误差好分析吗？其实都是一个套路，这个套路叫做泰勒展开。我们用它举例。
注意到 
 ，所以通过泰勒展开，我们可以得到
其中 
 。你可以看到，这里的 
 与 
 是不一定相等的。所以如果直接推导总体的误差，其实是会有些麻烦的（不过麻烦的点并不是这两个，而是之后的 
 与 
 ）。所以我们退而求其次，先求局部截断误差，再考虑之后的分析。
我们在前面提到说，局部截断误差的本质就是假设之前的数值解都是精确的。因此这里相当于 
 。这样的话，我们可以得到
这个时候，这里的 
 其实就对应了上面的 
 。
我们自然还可以把数值格式写出来，得到下面的式子
我们把两式相减就可以得到 
 。
我们有了局部的截断误差，自然需要考虑整体截断误差。刚开始显然是没问题的，然后我们做了一步之后，得到的整体截断误差其实就和局部截断误差没差，关键是第二步之后。因为第二步之后，除了局部截断误差，你在第一步得到的整体截断误差也会变大，因此这实际上是好几个因素的综合。所以我们需要拆分来看。
我们定义 
 是假设 
 时间点精确的情况下得到的方程的精确解。那么显然如果我们需要估计 
 ，最直接的方式就是补一个 
 进去。因为 
 我们是知道的，就是局部截断误差。所以不知道的是 
 。我们先把这些推导写在下面。
为了估计出这个式子，我们回想一下常微分方程的一些要求。在考虑常微分方程的连续依赖性的时候，我们有给出过方程所需要满足的Lipschitz条件。我们写在下面。
Definition 3: Lipschitz Condition
设存在 
 为常数，使得 
 对 
 中任意的 
 都有 
 ，则称 
 在 
 中关于变量 
 满足Lipschitz条件。
我们引入这个条件的原因相信你也能看出来：通过导数的差距，衡量出在不同的初值下，方程行为的差异程度。其实如果你的眼睛够尖，你就能看出来，虽然 
 比 
 多了一个假设，但是它们的微分方程的形式是完全相同的，只是初值不同。有了这个思路，我们给出下面的性质
Proposition 1:
假设 
 是对应微分方程 
 的不同初值，那么有 
简单说明一下这个结论。设 
 ，那么它显然要不是一直正，要不是一直负的。这样的话，我们可以去掉两边的绝对值。再考虑 
 ，这样的话就会有下面的结果
也就是说 
 （左边去掉绝对值就好），再根据中值定理可以得到
化简即可得到我们最后的结果。
使用这个估计，我们再回头来看 
 。你可以看到，其实 
 就是比 
 变了个初值而已。那么为什么这个定理有用呢？你看上面的定理，其实就相当于说，如果我们把时间从 
 “退回”到 
 ，那么就会使得误差降低一个膨胀因子 
 。所以我们有下面的结果
首先，我们考虑的是第二步，那么往回退一步（这里也就是退了一个步长 
 ）， 
 就退回到了 
 ，而 
 就退回到了 
 。所以你可以看到，这个时候两个值的差就恰好是
第一步的整体截断误差。而第一步的整体截断误差恰好就等于第一步的局部截断误差。所以当我们假设了之前的数值解都精确的时候，我们其实只能够给出第一个式子 
 的误差估计，也就是局部截断误差，而其它的误差（这里指的是 
 ）其实是从之前步的整体截断误差再扩大所得到的。
有了这个思路，相信你也不难做后面的推导，我们直接写出一般的结果
所以通过这个方法，我们就将整体的误差，用一系列局部截断误差估计出来了。
如果我们进一步假设 
 （我们之后说假设这个干啥），就可以得到下面的一些不等式估计。我就直接抄书了……
 （注意 
 ）
所以你可以看出来，对于这种方程，如果局部截断误差是 
 阶的，就可以得到整体截断误差是 
 阶的。
梯形法
考虑到我自己在开始分析梯形法的格式时也是处处碰壁，这里还是决定提一下这个格式。
Definition 4: Trapezoidal Scheme
你看着感觉挺复杂的，但其实这第二个方程只是取了两个时间的中值而已。
为了笔记的多样性，我们用另外一个方法来解释这个格式的来源。其实你只需要移个项，得到
（为了方便我们省去了第二个参数）然后我们还是，用精确解去代替这里的数值格式，得到
这个估计为什么合理呢？除了用导数的估计来看，还可以通过数值积分来看。首先注意到 
 ，那么这个时候，你可以看到，右边的式子其实就是一个数值积分格式（梯形格式），因此这个格式也被称为梯形格式。
我们再走一遍流程，做一遍误差分析。你可以看到上面其实我们就是重点关注了两个式子：精确解的泰勒展开，数值解的迭代格式。所以我们这里也是一样的思路，首先精确解的泰勒展开我们不难得到，可以写成
当你写到这个式子的时候其实你就应该会感觉分析上会有些麻烦了。数值格式上既有 
 又有 
 ，展开的话究竟在哪一个点展开？又怎么样与精确解的泰勒展开式子相消？所以我们这里最关键的，自然就是对 
 做处理，这就需要我们下一个式子了。
有些复杂，但它就是链式法则而已。
既然你发现了我们使用的本质上是二维的泰勒展开，那么相同的思路自然也可以用到我们的数值格式上。显然我们推导局部误差是方便的，所以只需要假设 
 ，那么我们会得到下面的结果
显然第二个式子可以用相同的方法作展开，也就是说
这样的话，其实可以看到，再两项作相减之后，都只会有 
 阶项存在，所以可以认为局部截断误差 
 。那么整体方法就是一个 
 阶方法。
关于隐式格式，我们不再做详细的介绍，虽然我们之后还是会用到它们。
龙格库塔（Runge-Kutta）方法简介
这个方法本质上是用导数在各个点的取值来构造数值积分，估计函数的端点值之差。也就是估计
右边的这个积分。龙格库塔法的意思就是：取定 
 ，设 
 ，那么想法就是找到一组系数 
 ，使得 
 ，也就是尽可能的去逼近导数。
一个显然的问题是：如何给出 
 ？这个确实是没有办法的，因为它可能会和 
 有关，而这个你一般是不知道的。那么这个时候，我们给出了这么一个逼近：先给出
（左边界上的点的值我们知道），然后对于 
 时间，用向前Euler估计，得到下面的估计
（注意他们一步一步都是约等于过来的）问题就落在了：如何取这些系数？
那么我们推广上面的思想，可以得到下面的这个推导公式。
其中 
 。
为什么说这是我们上面的那个思路的一种推广呢？这是因为我们这里其实每一步都运用了前面的信息，并且要估计未知的 
 以便我们计算 
 。那么你也可以看到，这个时候呢，我们自然需要的是：当 
 中的时间参数 
 在某一个点 
 时，它的函数值参数 
也必须要是 
 的一个线性估计
。这就是 
 的含义（不然时间就不同步了）。
然而，这只是最最一般的情况，我们很显然不能就到此为止了。因为我们要逼近导数，所以我们考虑把真值的表达式写成下面的这种形式
而这里的 
 是一个表达式，因为这里的 
 都是精确解，所以这个表达式是可以通过
泰勒展开得到的。很明显，你想要几阶的精度，就把它展开到几阶即可。但是要注意的是每多展开一阶，计算的复杂度就会多一个量级。比方说我们把它展开到三阶，那么表达式就可以写成下面这样
其中
你就会发现其实这个系数就已经相当复杂了。
那么精确解是这样的一个形式，对应的数值格式呢？其实就是把公式改成了下面这个形式
这里的 
 就是时间为 
 时刻的数值解， 
 啥意思自然不必多说。你可以看到其实本质上的差异就是 
 和 
 。而这个 
 其实就是我们之前假设的那个 
 。这个计算事实上也不是特别容易的，如果你要求 
 的数量比较大的话。我们这里随便计算两个表达式如下
其中 
 参考上面的表达式。所以比方说你要求 
 ，那你就需要把 
 计算出来，但是你也能看出来这个计算量已经是非常大了。
最后，我们所需要做的事情，就是在确定了我们要的误差阶数之后，根据不同的 
 ，去匹配不同的系数，
使得 
 在系数上匹配
。这便是龙格库塔方法的全部思想，而它的具体格式实现，大家可以通过wikipedia寻找更多内容。
最后，提一下，我们对于龙格库塔方法其实是要控制两个参数 
 ，而这对应得到的格式，我们就称他们为 
 级 
 阶龙格库塔方法。
稳定性
我们用刚性的常微分方程来看我们为什么需要稳定性。一个比较简单的方程形式是
其中 
 是一个绝对值很大的负数（比如说 
 ）。你自然会疑惑，这也算是上面的方程的一个例子，我们如果使用向前Euler方法，你可以知道它具有一阶的精度（误差为 
 ），也就是说你肯定不怕它不收敛。既然收敛了，为什么还要担心其它的东西？
这就牵涉到一个概念叫做稳定性。也就是说，这个格式确实收敛没有错，但是很有可能因为你的其它参数（比如步长）选的不好，导致它在给定的步数下没有办法计算到给定的精度。所以我们先介绍稳定性的相关概念，再在之后观察为什么刚性问题会有其它需要注意的地方。
绝对稳定区
我们还是分析之前分析过的常微分方程，但是稍微写的具体一点，像这样
如果我们使用向前Euler法进行计算，我们就可以得到下面的计算格式
同样的，我们对精确解的式子做泰勒展开，可以得到
其中 
 就是式子最后的那个误差项。相减之后就可以得到我们下面的这个误差式子
其中这里的 
 （我们没有再加绝对值）。所以你可以看出来，尽管我们存在一个局部的误差，但是如果这里的膨胀因子 
 过大，你就自然会发现，一步一步的迭代，其实误差是会发散的，并不是我们之前想的“收敛”。这可能会与之前的推导矛盾，你会觉得，明明我之前都找到阶数了，结果反而发散了，这不是胡扯吗？但是要注意的是，我们给出阶数，其实
只是相当于给出了误差与步长的关系。只有你的步长趋于0的时候，我们才能得到收敛性。那么如果你的 
 过大，一般情况下都是 
 取得太大了，那自然不存在收敛这一说。
运用几乎相同的分析方法，你可以观察一下，如果回到上面一般的方程 
 ，那么Euler格式一般就是有诸如
的形式，这个形式的分析方法和上面类似的，但是这里的膨胀因子变成了 
 ，这里的 
 就是函数 
 的Lipschitz常数，也就是取使得
的最小的常数 
 。
那么根据这些特殊的例子，你也能看出来如何定义绝对稳定区了。
Definition 5: Absolute Stability Area
定义特征多项式 
 ，那么满足 
 的区域称为绝对稳定区。
这个定义是很一般的，不同的数值格式会对应有不同的特征多项式。比如说上面的向前Euler法，它的特征多项式就是 
 。也就是 
 之前那个我们说的“膨胀因子”。所以你也能看出来，这个格式的绝对稳定区其实就是 
 ，它在复平面上就是一个圆，当然了如果限制在实数上，对应的 
 的取值就是 
 。
在这部分最后我们简单提一下数值格式的各种可能的稳定性。
Definition 6: Stability
若对于任意的步长 
 ，数值格式计算得到的解都会收敛，则称为无条件稳定。若需要步长满足一定的条件才可稳定，则称为条件稳定。否则称为无条件不稳定。
一些具体的例子会牵涉到多步法的概念，但是多步法的解决方法和框架等其实应该是偏微分方程数值解里更为常见的，所以我们这里不会介绍这样的方法，但是我们还是用一个具体的多步法例子，来介绍一般的线性差分方程，如何判断它的稳定性。
Example 1:
这是一种多步法，在进行之前，我们先看看这个格式为什么是存在的。老方法，用精确解去替代数值解，我们可以得到
也就是说
我们对两个不是在时间 
 的式子做一下泰勒展开，看一下我们的结果
所以代入之后，你会发现 
 的项都被消掉了，而 
 项之前的系数就是1，和右边的正好匹配上，剩下的就是一个小量了。所以这个式子是一个合理的导数的估计。
好的，那我们如何判断呢？这里给大家介绍一个数列的特征方程法。如果你是竞赛党出身，你绝对不可能对这个名词感到陌生。
针对这个数值例子，我们忽略所有的与右端导数项有关的内容，只考察差分的数列项。那么针对这个例子，它对应的特征方程其实就是 
 ，可以得到它的两个根为 
 。这样的话，理论证明了它的根的形式就是 
 。你看这个 
 就可以说明了这个格式其实是无条件不稳定的，我们不应该使用它。
一般的，我们有下面的性质：
Proposition 2:
对于特征方程，如果它的根绝对值均小于1，则格式无条件稳定。若根均绝对值大于1，或者存在一个以上的根绝对值为1，则格式无条件不稳定。如果存在一个根绝对值为1，其余根绝对值小于1，则格式条件稳定。
具体为什么，我们在之后再说。
刚性常微分方程的有限差分方法
有了稳定性的概念，我们再回头看方程 
 。如果我们用向前Euler法做离散，可以得到它的对应的值 
 ，那么你容易发现，若要求 
 ，其实对于 
 的要求是相当高的（如果 
 ，那么 
 ）。这个时候迭代的话，在机器上受到误差的影响，其实很多时候会严重影响数值解的有效性。究竟如何修改格式和解，其实是一个大的主题。我们这里只介绍部分方法和背后的思想，作为这一节的最后一部分。
我们之前提到了绝对稳定区的概念。那么你也能看到，如果它的绝对稳定区涵盖整个复平面的左半部分，那么无论我的 
 怎么取，在这个例子中它都是稳定的。这个格式是存在的，比如说向后Euler格式，它的稳定性要求是 
 ，所以至少 
 的时候它绝对是可以满足的。而 
 其实是一个非常松的要求，因为 
 （如果 
 ，那么你也能看到，这个方程本身就没有办法得到一个收敛的数值解，这不是数值解法喜欢关注的内容）。
所以你也可以看到，如果它的稳定区范围很广，那么自然刚性常微分方程就能解。但是这种方法多半是隐式方法，而隐式方法因为通过迭代法解初值，所以计算复杂度要比显式方法大得多。所以有没有显式的格式，能够用于解刚性方程？答案是肯定的。这就是我们下面要介绍的内容。
龙格库塔——切比雪夫方法
首先我们回忆一下龙格库塔方法（上面虽然写了，但是可能比较抽象？）。一个具体的例子如下：
Example 2:
我们可以看出它就是取了一个 
 的中点。然后我直接用了中间点的导数值去近似估计了 
 点的导数值。
如果我们设 
 ，代入我们可以得到， 
 ，而 
因此归纳法容易得到，如果我们的方法是 
 步的，就可以将我们的数值格式写成这样的形式
其中 
 ，而 
 。
如果我们需要它具有一阶精度，这里我们显然需要 
 。那么你可以看到，这个要求等价于要求多项式满足 
 。学过数值逼近的同学肯定已经开始有感觉了。这可以通过正交多项式来拟合。因此我们的这个方法带上了切比雪夫的名字。
我们不加证明的给出我们考虑的多项式
容易发现，如果要求格式稳定，那么必须要有 
 。而切比雪夫多项式的性质告诉我们，它的满足的稳定区间范围为 
 ，你也可以看到，通过改变 
 ，那么即使是显式的方法，我们也能够达到我们想要的稳定性结果。
当然，如果存在扰动，那就又会有很多其它的问题，书上提到了一个位移切比雪夫多项式法 (Shifted Chebyshev Polynomial)，碍于篇幅我们不再介绍，感兴趣的可以寻找相关资料，或直接参考LeVeque的书。
小结
和《笔记整理》系列文章不同，《专题》系列文章会显得有点“左一块右一块”的意思。而且因为我们更加注重于直接性，用于铺垫的语句也做了大量的删减。但还好，我们算是用了不短的篇幅，介绍完了大部分常微分方程的有限差分方法。包括显式格式，龙格库塔方法，稳定性分析和扩展的刚性常微分方程的解决方案。
关于隐式格式和与稳定性有关的更多概念，虽然在常微分方程数值解法中就已经出现，但是它们的应用更多的其实是在偏微分方程的差分方法中。因此我们会在之后再提到这些概念。
——————————————————————————————————————
本专栏为我的个人专栏，也是我学习笔记的主要生产地。任何笔记都具有著作权，不可随意转载和剽窃。
个人微信公众号：cha-diary，你可以通过它来有效的快速的获得最新文章更新的通知。
专栏目录：笔记专栏|目录
想要更多方面的知识分享吗？欢迎关注专栏：一个大学生的日常笔记。我鼓励和我相似的同志们投稿于此，增加专栏的多元性，让更多相似的求知者受益~

参考《常微分方程》第三版（王高雄）
主要内容：一阶线性偏微分方程的解可看成积分曲线，而首次积分表示为特征曲线，由特征曲线构成积分曲线。因此 ，可以用常微分方程方法解一阶线性偏微分方程。
7.1基本概念
对于自变量 
 与未知函数 
 的一阶偏导数 
一阶偏微分方程： 
      （7.1）
一阶线性偏微分方程： 
 （7.2）a
一阶齐次线性偏微分方程： 
             （7.2）b
一阶拟线性偏微分方程： 
 （7.3）
解 
 可以想象为在空间 
 中的一张n维曲面，通常称为偏微分方程（7.1）的
积分曲面。
称 
 （7.5）为（7.2）b的特征方程。（7.3）的求解问题则可化为（7.2）b的方程类型来处理。
7.2 一阶线性偏微分方程与常微分方程组的关系
考虑初值问题（p344）
首次积分：
n个首次积分称为彼此独立的：
常微分方程组与一阶线性偏微分方程之间的关系：
7.3 利用首次积分求解常微分方程组
考虑方程组
通积分：方程组（7.11）的n个彼此独立的首次积分的全体（7.12）称为（7.11）的通积分。
如何求首次积分：（p350）
（1）如果能得到常微分方程的含任意常数的通解 
 或 
 ，则当 
 或 
 时，可解出 
 ,
即为首次积分。
（2）还可以通过构造全微分 
 得到首次积分 
 。
7.4 一阶线性偏微分方程的解法
这一节将讨论一阶齐次线性和拟线性偏微分方程的通解的结构.
对于（7.2）b的通解结构，有
关于一阶拟线性偏微分方程的通解结构，有如下定理：
7.5 柯西问题
所谓柯西问题，用几何的语言说，就是求一阶偏微分方程（7.29）或（7.30）的通过某一给定曲线的积分曲面.这里必须指出，对于某些曲线（譬如特征曲线）柯西问题是不确定的，因为对一条特征曲线而言，可以有无穷多个特征曲面经过它；而对于另外一些曲线，柯西问题甚至没有解存在.但是我们有下面的一般结果.
关于拟线性一阶偏微分方程的柯西问题类似地有：
2020.11.26