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