微分方程数值解法

admin管理员组

文章数量:1794759

微分方程数值解法

大家好！

趁热打铁，在介绍完常微分方程的那一套理论后，我们继续介绍偏微分方程的相关内容。

学过偏微分方程的都知道，偏微分方程一般有三类：椭圆方程，抛物方程和双曲方程。它们每个方程都各有自己的一套理论，要想完全了解相关性质和理论是极为困难的。但是有幸在数值解中，我们可以看到这些方程都可以用一套框架来解决。这也是这篇文章要重点关注的内容。

不过在我们进行偏微分方程的话题之前，我们还会再引入一些之前常微分方程还没有说完的重要的内容和方法。因为最近我还没有更新有限元方法的计划，所以这也是专题的最近的最后一篇文章。

提供之前的笔记：

微分方程数值解法|专题（1）——常微分方程的有限差分方法

我们开始本节的内容。

目录线性多步法的绝对稳定区抛物方程的有限差分方法冯诺依曼分析法Crank-Nicolson方法双曲方程的有限差分方法椭圆方程的有限差分方法五点格式法线性多步法的绝对稳定区

我们还是研究方程 y' = f(t,y) 。无论是最简单的向前Euler法，还是之后的龙格库塔方法，我们的导数项 f(t_i, y_i) 都是没有变的。我们之前也说了它们是显式格式。这是因为我们在计算的时候，右端的导数项的时间点都在过去。实际情况下，右边的导数项是可以变的，比方说向后Euler法

w_{i+1} = w_i + h f(t_{i+1}, w_{i+1})

你会发现，导数项的时间是 t_{i+1} ，所以你会发现，这个数列通项没有办法再那么容易得到了，因为你在计算时间点 t_{i+1} 的时候，所依赖的导数值与 t_{i+1} 有关。这个时候在进行编程的时候，一般我们采用迭代法。也就是说，我们把左右两端的 w_{i+1} 都当作未知数 x ，然后做类似于这样的迭代

x_{n+1} = w_i + hf(t_{i+1}, x_n)

然后做出一个关于 \\{x_n\\} 的收敛数列，不过这不是我们重点关注光环效应的内容。

既然我们可以把它改成 f(t_{i+1},y_{i+1}) ，那么我自然可以改成之前某些时间点的线性组合。所以我们假设 y' = \\lambda y ，来导出线性多步法的一般形式

\\sum_{j =0 } ^ r \\alpha_j w_{n+j} = h \\sum_{j=0}^r \\beta_j \\lamb猫花da w_{n+j}

设 z = h\\lambda ，我们可以得到我们的一般形式

\\sum_{j=0}^{r}\\left(\\alpha_{j}-z \\beta_{j}\\right) w_{n+j}=0

它对应的特征多项式为 R(x) = \\sum_{j=0}^{r}\\left(\\alpha_{j}-z \\beta_{j}\\right) x^{j} ，所以我们只需要考虑让它的根满足一定条件，即可找到它们的绝对稳定区。下面我们来举一些实际的例子，来说明一下如何操作。

Example 1: Forward Euler Scheme w_{i+1} = w_i + h \\lambda w_i

两边写出它的特征多项式，可以得到 x = 1 + z ，一个一元一次方程（注意超兽机神断空我这里的 z 不是变量，因为这个方程只是关于 x 的），所以为了要求满足根的条件，我们就要求 |1 + z | < 1 ，这与我们上一节推的结果完全一致。

Example 2: Backward Euler Scheme w_{i+1} = w_i + h \\lambda w_{i+1}

两边化为特征多项式，可以得到 x = 1 + zx ，解 x 可以得到 x = \\frac1{1-z} ，那么你也可以看到，只要求 |\\frac1{1-z}| < 1 就可以，这个式子对于 z<0 的时候是恒成立的，所以向后Euler其实是无条件稳定的。

Example 3: Trapezoidal Scheme w_{i+1}-w_i= \\frac{h}{2}(\\lambda w_i + \\lambda w_{i+1})

两边化为特征多项式，有 x - 1 = \\frac{z}{2}(1+x) ，解 x 可以得到 x = \\frac{2+z}{2-z} 。学过复变的话，你应该知道这个式子其实是一个保形映射。这样的话，可以得到，如果 |x| < 1 ，那么 \\operatorname{Re}(z) < 0 ，所以这个的绝对稳定区就是负半平面。

这里只是一些比较简单的例子，但是实际情况下的绝对稳定区的推导其实是相当复杂的。LeVeque的盐酸克伦特罗书上有一些篇幅，感兴趣的朋友们可以去查看一下。

抛物方程的有限差分方法

我们研究的就是下面的热方程

\\left\\{\\begin{array}{ll}{u_{t}=c u_{x x},} & {a \\leqslant x \\leqslant b \\quad t \\geqslant 0} \\\\ {u(x, 0)=f(x),} & {a \\leqslant x \\leqslant b} \\\\ {u(a, t)=l(t),} & {t \\geqslant 0} \\\\ {u(b, t)=r(t),} & {t \\geqslant 0}\\end{array}\\right.

学过热方程的都知道，它本质上是一个衡量杆的温度的方程。当然了它本质上其实是一个初边值问题，直观来看是因为，它既需要初值条件，又需要边值条件。

下面我们来看看如何差分吧。在之前我们做过对导数的离散，那么这里 u_t 相当于对方程针对时间变量的导数做离散，关键在于 u_{xx} ，它是一个二阶导数，那么如何去做呢？

事实上方法是完全类似的，在之前我们用向前Euler法来进行近似的时候，其实本质上是下面这个近似。

\\frac{u_{n+1} - u_n}{h} \\simeq u'

（注意我们这里因为方程的符号改成了 u ，所以自然统一的用 u 来表示了）如果将步长标记为单位长度，那么就相当于

u' = \\nabla u_{n+1} = u_{n+1} - u_{n}

注意这里的 \\nabla 是向后差分符号，不是梯度。所以你也能猜到，如果我们要估计 u'' ，结果是完全类似的，也就是下面的过程

u'' = \\nabla u_{n+1}- \\nabla u_{n} = u_{n+1} - 2u_{n} + u_{n-1}

我们再把步长 h 代回去，就可以得到我们最后的差分格式

u_{x x}(x, t) \\approx \\frac{1}{h^{2}}(u(x+h, t)-2 u(x, t)+u(x-h, t))

其中你也可以看出，这个差分格式因为是对空间变量差分的，所以时间没有变，而空间变量中 u(x+h,t) 就对应着上面的 u_{n+1} 。

你也可以看出，从这里开始，因为方程开始有多个变量，所以我们不能够再使用之前的简单的下标标记了。这里的 u_t 的离散也是容易的，就是下面的结果。

u_{t}(x, t) \\approx \\frac{1}{k}(u(x, t+k)-u(x, t))

因为我们在空间上设置了离散步长为 h ，所以我们这里在时间上就用另外一个字母 k 来表示。那么你也容易看出，分别用差分格式代到两边，就能得到我们最终的表达式。而且使用二维泰勒展开，我们也容易得到我们的局部截断误差，分别是 h^{2} u_{x x x x}\\left(c_{1}, t\\right) / 12 和 k u_{t t}\\left(x, c_{2}\\right) / 2 ，其中 c_1 \\in (x-h ,x+h), c_2 \\in (t, t+h) ，潜台词就是说，它的局部截断误差就是 O(k) + O(h^2) 。

当然了，上一节说了那么多稳定性，我们这一节自然也不能错过。我们定义 w_{ij} 就是在 (x,t) 处的数值解，那么自然的，我们可以将它按照时间的顺序排列，得到下面的一个形式

w_{i, j+1}=w_{i j}+\\frac{c k}{h^{2}}\\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\\right)

这里，我们假设 \\sigma = \\frac{ck}{h^2} ，那么上面的式子就是

w_{i j}+\\frac{c k}{h^{2}}\\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\\right)=\\sigma w_{i+1, j}+(1-2 \\sigma) w_{i j}+\\sigma w_{i-1, j}

你也可以看出来，当我只考虑时间的变化后，固定时间，再考虑空间，其实我们也可以得到一个关于空间变量的数列通项公式。因为我们的变量变成了两个，所以这样的一个通项就不能够轻松的使用特征多项式分析了，但是还好我们还是有招，这就是后面所使用的冯诺依曼分析法。

冯诺依曼 (Von Neumann) 分析法

虽然我们在抛物方程这里引入了它，但实际上它可以应用在各种常见的偏微分方程中。所谓的冯诺依曼分析法，其本质就是用向量，将所有的空间离散点拼起来。那么这样的话，这个矩阵所对应的表达式就仅仅与时间有关了。这样自然分析就会方便很多。

回到上面这个问题，我们针对这个差分格式，就可以化成这样的矩阵。

\\left[ \\begin{array}{c}{w_{1, j+1}} \\\\ {\\vdots} \\\\ {w_{m, j+1}}\\end{array}\\right]=\\left[ \\begin{array}{ccccc}{1-2 \\sigma} & {\\sigma} & {0} & {\\cdots} & {0} \\\\ {\\sigma} & {1-2 \\sigma} & {\\sigma} & {\\ddots} & {\\vdots} \\\\ {0} & {\\sigma} & {1-2 \\sigma} & {\\ddots} & {0} \\\\ {\\vdots} & {\\ddots} & {\\ddots} & {\\ddots} & {\\sigma} \\\\ {0} & {\\cdots} & {0} & {\\sigma} & {1-2 \\sigma}\\end{array}\\right] \\left[ \\begin{array}{c}{w_{1 j}} \\\\ {\\vdots} \\\\ {w_{m j}}\\end{array}\\right] + \\sigma \\left[ \\begin{array}{c}{w_{0, j}} \\\\ {0} \\\\ {\\vdots} \\\\ {0} \\\\ {w_{m+1, j}}\\end{array}\\right]

写的紧凑一点，就是公式 \\mathbf{w}_{j+1} = A \\mathbf{w}_j + \\mathbf{s}_{j} 。那么为了1942票房让这个公式收敛，我们显然需要让 \\rho(A)<1 ，其中 \\rho(A) 表示它的最大特征值。

这里的 A 你动物实验室也可以看到，是一个三对角矩阵，它是数值线性代数里非常重要的一类矩阵。这里我们不加证明的给出一个特殊矩阵的特征值。

Proposition 1:对于矩阵 \\boldsymbol{T}=\\left[ \\begin{array}{ccccc}{1} & {-1} & {0} & {\\cdots} & {0} \\\\ {-1} & {1} & {-1} & {\\ddots} & {\\vdots} \\\\ {0} & {-1} & {1} & {\\ddots} & {0} \\\\ {\\vdots} & {\\ddots} & {\\ddots} & {\\ddots} & {-1} \\\\ {0} & {\\cdots} & {0} & {-1} & {1}\\end{array}\\right] ，它的特征值为\\lambda_{j}=1-2 \\cos \\pi j /(m+1), j=1, \\cdots, m

事实上，这个矩阵我们可以做一个拆分，得到 \\mathbf{T} = \\mathbf{I + S} ，其中 \\mathbf{S} = \\left[ \\begin{array}{ccccc}{0} & {-1} & {0} & {\\cdots} & {0} \\\\ {-1} & {0} & {-1} & {\\ddots} & {\\vdots} \\\\ {0} & {-1} & {0} & {\\ddots} & {0} \\\\ {\\vdots} & {\\ddots} & {\\ddots} & {\\ddots} & {-1} \\\\ {0} & {\\cdots} & {0} & {-1} & {0}\\end{array}\\right] 变频电机与普通电机的区别，根据高等代数里有关特征值的结论，容易得到 \\mathbf{S} 的特征值就是 \\lambda_{j}=-2 \\cos (\\pi j /(m+1)), j=1, \\cdots, m ，对，就是少了一个 1 而已。

回到我们这个例子，你也容易通过变换得到 A = (1-2\\sigma)\\mathbf{I} + \\sigma \\mathbf{S} ，所以它的特征值就是

\\lambda_j氦质谱检漏仪 = 2 \\sigma(\\cos ( \\pi j /(m+1) )-1)+1, j=1, \\cdots, m 。

通过对余弦的估计，你可以知道， \\rho(A) \\in (-4\\sigma+1, 1) ，这样的话根据 |\\rho(A)| < 1 ，就能得到我们最后的结果 \\sigma < \\frac12 。注意我们一般不取 \\sigma = \\frac12 ，这是因为数值算法往往存在误差。

下面我们再介绍一种向后差分方法。我们不再给出推导细节，直接给出它的差分格式

\\frac{1}{k}\\left(w_{i j}-w_{i, j-1}\\right)=\\frac{c}{h^{2}}\\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\\right)

其实变化就是在左边，把它改成了一种向后差分的格式（注意这和之前的向后Euler不是一回事）。那么也容易通过化简得到我们的矩阵方程。

\\left[ \\begin{array}{ccccc}{1+2 \\sigma} & {-\\sigma} & {0} & {\\cdots} & {0} \\\\ {-\\sigma} & {1+2 \\sigma} & {-\\sigma} & {\\ddots} & {\\vdots} \\\\ {0} & {-\\sigma} & {1+2 \\sigma} & {\\ddots} & {0} \\\\ {\\vdots} & {\\vdots} & {\\vdots} & {\\vdots} & {-\\sigma} \\\\ {护理论文0} & {\\cdots} & {0} & {-\\sigma} & {1+2 \\sigma}\\end{array}\\right] \\left[ \\begin{array}{c}{w_{1 j}} \\\\ {\\vdots} \\\\ {w_{m j}}\\end{array}\\right]=\\left[ \\begin{array}{c}{w_{1, j-1}} \\\\ {\\vdots} \\\\ {w_{m j}-1}\\end{array}\\right]+\\sigma \\left[ \\begin{array}{c}{w_{0 j}} \\\\ {0} \\\\ {\\vdots} \\\\ {0} \\\\ {w_{m+1, j}}\\end{array}\\right]

写的紧凑一些，其实就是 \\boldsymbol{w}_{j}=A^{-1} \\boldsymbol{w}_{j-1}+\\boldsymbol{b} 。那么我们只需要找到保证 \\rho(A)>1 即可。这是显然的，因为用和上面一模一样的方法，就可以得到特征值为 1+2 \\sigma-2 \\sigma \\cos \\frac{\\pi j}{m+1} >1 ，所以这个格式是无条件稳定的。

Crank-Nicolson方法

因为之前的两个方法它们的误差是 O(k+h^2) ，这其实会有问题就是主要的误差其实都在时间步长上。我们需要把时间步长 k 取得很小，可能才能够达到我们要的精度。因此Crank-Nicolson方法可以很好地弥补这个缺陷。

这个方法本质上就是用不同的差分方法来近似 u_{xx} 和 u_t 。这里我们的 u_t 使用向后差分，和上面那一个一样。而我们的 u_{xx} 使用的是混合差分。具体来说，就是用

\\frac{1}{h^{2}}\\left(\\frac{1}{2}\\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\\right)+\\frac{1}{2}\\left(w_{i+1, j-1}-2 w_{i, j-1}+w_{i-1, j-1}\\right)\\right)

代替 u_{xx} 。如何考虑它的稳定性分析呢？如果我们还是一样，按照时间的顺序去排列，可以得到下面的形式

-\\sigma w_{i-1, j}+(2+2 \\sigma) w_{i j}-\\sigma w_{i+1, j}=\\sigma w_{i-1, j-1}+(2-2 \\sigma) w_{i, j-1}+\\sigma w_{i+1, j-1}

这个方程矩阵化各板块龙头股会稍微麻烦一些的，大体上的形式就是 A \\mathbf{w}_j = B\\mathbf{w}_{j-1} + \\sigma (\\mathbf{s}_{j-1} + \\mathbf{s}_j) ，其中 A=\\left( \\begin{array}{ccccc}{2+2 \\sigma} & {-\\sigma} & {0} & {\\cdots} & {0} \\\\ {-\\sigma} & {2+2 \\sigma} & {-\\sigma} & {\\ddots} & {\\vdots} \\\\ {0} & {-\\sigma} & {2+2 \\sigma} & {\\ddots} & {0} \\\\ {\\vdots} & {\\ddots} & {\\ddots} & {\\ddots} & {-\\sigma} \\\\ {0} & {\\cdots} & {0} & {-\\sigma} & {2+2 \\sigma}\\end{array}\\right) ，且B=\\left( \\begin{array}{ccccc}{2-2 \\sigma} & {\\sigma} & {0} & {\\cdots} & {0} \\\\ {\\sigma} & {2-2 \\sigma} & {\\sigma} & {\\ddots} & {\\vdots} \\\\ {0} & {\\sigma} & {2-2 \\sigma} & {\\ddots} & {0} \\\\ {\\vdots} & {\\ddots} & {\\ddots} & {\\ddots} & {\\sigma} \\\\ {0} & {\\cdots} & {0} & {\\sigma} & {2-2 \\sigma}\\end{array}\\right) 。那么现在问题自然就是落到了矩阵 A^{-1}B 的特征值上。

首先，我们把矩阵做一些拆解，容易看出 A = (2+2\\sigma)\\mathbf{I} - \\sigma \\mathbf{S} ， B = (2-2\\sigma)\\mathbf{I} + \\sigma\\mathbf{S} 。那么这样的话，考虑设 \\lambda_j 为 \\mathbf{S} 的特征值，并且对应一个特征向量 v_j ，那么这就容易得到

Bv_j = (2-2\\sigma)v_j + \\sigma\\mathbf{S}v_j = (2-2\\sigma + \\sigma \\lambda_j)v_j

然后呢，如果我们有 Av_j = \\lambda_j v_j ，那么就一定有 A^{-1}v_j = \\frac1\\lambda v_j 。所以根据这个结论，我们容易得到最后的 A^{-1}B 的每一个特征值为

\\mu_j = \\fr中科院计算所ac{2-2\\sigma + \\sigma \\lambda_j}{2+2\\sigma - \\sigma\\lambda_j}, \\lambda_{j}=-2 \\cos \\pi j /(m+1)

容易验证，这里每一个数其实都是严格小于法国全称1的，也就是说格式是无条件稳定的。

但是如果单纯说稳定性，它并没有什么优势。我们通过对它做局部截断误差的误差分析，就可以看出来更多细节。我们这里详细的来走一遍。这里不失一般性，我们假设 c = 1 。

首先还是考虑，使用泰勒展开，首先容易得到的是

\\frac{u(x,t) - u(x, t-k)}{k} = u_t(x, t) - \\frac k2u_{tt}(x,t) + \\frac{k^2}{6} u_{ttt}(x,\\xi_1)

（我希望你没有忘记， (x,t) 对应的是 w_{ij} 的时间点）然后呢，我们再考虑两个中心差分，有

\\frac{1}{h^2}\\left(u(x+h,t)-2 u(x,t)+u(x-h,t)\\right) = u_{xx}(x,t) + \\frac{h^2}{24}(u_{xxxx}(\\xi_2,t) + u_{xxxx}(\\xi_3,t))

当然只要改一个时间点，就可以得到

\\frac{1}{h^2}\\left(u(x+h,t-k)-2 u(x,t-k)+u(x-h,t-k)\\right) = u_{xx}(x,t-k) + \\frac{h^2}{24}(u_{xxxx}(\\xi_4,t-k) + u_{xxxx}(\\xi_5,t-k))

如果我们略去最高阶的误差项，那么也能够看出，右边最最关键的部分就是 \\frac{u_{xx}(x,t) + u_{xx}(x,t-k) }{2} ，而左边的东西就是 u_t(x, t) - \\frac k2u_{tt}(x,t) 。我们还是抓住我们的时间和空间都在 (x,t) 点展开，就可以得到

\\frac{u_{xx}(x,t-k)}{2} = \\frac{u_{xx}(x,t) - ku_{xxt}(x,t)}2 + O(k^2)

剩下的就是拼在一起的事情了，这个证明已经没什么难度了，因为

u_{xxt}(x,t) =\\frac{\\partial{u_{xx}}(x,t)}{ \\partial t} = \\frac{\\partial{u_{t}}(x,t)}{ \\partial t} =u_{tt}(x,t)

这就完成了分析，因为你可以看到， O(k) 那一项已经被消掉了，所以这个误差是 O(k^2+h^2) ，这已经比上面的误差要好很多了。

事实上，如果你把左边的时间变量的差分离散改成中心差分，就变成了Richardson格式。这是一个无条件不稳定的格式，它的具体分析就是19年丘成桐杯竞赛的应用数学个人赛的第二题。据北大小伙伴说，如果做对3个题，排名据说就能到银牌……

双曲方程的有限差分方法

这个时候，我们考虑的方程为

\\left\\{\\begin{array}{ll}u_{tt} = c^2u_{xx} \\\\{u(x, 0)=f(x),} & {a \\leqslant x \\leqslant b} \\\\ {u_{t}(x, 0)=g(x),} & {a \\leqslant x \\leqslant b} \\\\ {u(a, t)=l(t),} & {t \\geqslantbp公司 0} \\\\ {u(b, t)=r(t),} & {t \\geqslant 0}\\end{arrvx毒气ay}\\right.

怎么差分，我觉得你已经有数了。两个变量都做中心差分即可，也就是格式写成

\\frac{w_{i, j+1}-2 w_{i j}+w_{i, j-1}}{k^{2}}声环境-c^{2} \\frac{w_{i-1, j}-2 w_{i j}+w_{i+1, j}}{h^{2}}=0

按照时间顺序排一下，可以得到方程的形式为

w_{i, j+1}=\\left(2-2 \\sigma^{2}\\right) w_{i j}+\\sigma^{2} w_{i-1, j}+\\sigma^{2} w_{i+1, j}-w_{i, j-1}

这里的 \\sigma = \\frac{ck}{h} 。

当你把这个方程写出来之后你可能就发现情况有点不太对了。首先就是初值的情况要单独拉出来，其次就是这个方程即使是按照时间排序，也会与之前的两步时间量有关。那应该怎么办呢？

我们一个一个来解决，首先初值的情况，就是这里的 j = 1 的情况，你知道 j=0 对应的就是初值的时候，这个值已经由初值条件控制住了。但是 j=1 你会发现它会与一个 j= -1 的量有关系，这是个什么玩意？也是因为这个，我们需要给出一个估计，一个关于 w_{i, -1} 的估计。比较好办的是我们有了关于 u 的导数的表达式，所以我们可以考虑构造一个中心差分格式来近似导数。具体点来说，我们的估计是

w_{i,-1} \\approx w_{i 1}-2 k g\\left(x_{i}\\right)

代入，合并就可以得到我们的初值方程

w_{i 1}=\\left(1-\\sigma^{2}\\right) w_{i 0}+k g\\left(x_{i}\\right)+\\frac{\\sigma^{2}}{2}\\left(w_{i-1,0}+w_{i+1,0}\\right)

事实上，它也可以被写成是一个矩阵形式。而之后的矩阵方程都按照那个格式来写即可，最后我们可以得到下面的矩阵方程组。

\\left( \\begin{array}{c}{w_{11}} \\\\ {\\vdots} \\\\ {w_{m 1}}\\end{array}\\right)=\\frac{1}{2} \\boldsymbol{A} \\left( \\begin{array}{c}{f\\left(x_{1}\\right)} \\\\ {\\vdots} \\\\ {f\\left(x_{m}\\right)}\\end{array}\\right)+k \\left( \\begin{array}{c}{g\\left(x_{1}\\right)} \\\\ {\\vdots} \\\\ {g\\left(x_{m}\\right)}\\end{array}\\right)+\\frac{1}{2} \\sigma^{2} \\left( \\begin{array}{c}{l\\left(t_{0}\\right)} \\\\ {0} \\\\ {\\vdots} \\\\ {0} \\\\ {r\\left(t_{0}\\right)}\\end{array}\\right) \\left( \\begin{array}{c}{w_{1, j+1}} \\\\ {\\vdots} \\\\ {w_{m, j+1}}\\end{array}\\right)=\\boldsymbol{A} \\left( \\begin{array}{c}{w_{1 j}} \\\\ {\\vdots} \\\\ {w_{m j}}\\end{array}\\right)-\\left( \\begin{array}{c}{w_{1, j-1}} \\\\ {\\vdots} \\\\ {w_{m, j-1}}\\end{array}\\right)+\\sigma^{2} \\left( \\begin{array}{c}{l\\left(t_{j}\\right)} \\\\ {0} \\\\ {\\vdots} \\\\ {0} \\\\ {r\\left(t_{j}\\right)}\\end{array}\\right) \\boldsymbol{A}=\\left( \\begin{array}{ccccc}{2-2 \\sigma^{2}} & {\\sigma^{2}} & {0} & {\\cdots} & {0} \\\\ {\\sigma^{2}} &缺口 {2-2 \\sigma^{2}} & {\\sigma^{2}} & {\\ddots} & {\\vdots} \\\\ {0} & {\\sigma^{2}} & {2-2 \\sigma^{2}} & {\\ddots} & {0} \\\\ {\\vdots} & {\\ddots} & {\\ddots} & {\\ddots} & {\\sigma^{2}} \\\\ {0} & {\\cdots} & {0} & {\\sigma^{2}} & {2-2 \\sigma^{2}}\\end{array}\\right)

如果我们把一般的这个方程拼的紧凑一些，我们就可以得到它的一般情形其实是诸如 \\mathbf{w}_{j+1} = A \\mathbf{w}_j - \\mathbf{w}_{j-1} + \\sigma^2 \\mathbf{s}_j 的，所以核心还是在前面。这是一个多步的数列，它的解决方案类似于高次的常微分方程——化为方程组。具体点来说，我们应该把方程组写成下面这个形式。

\\left( \\begin{array}{c}{\\mathbf{w}_{j+1}} \\\\ {\\mathbf{w}_{j}}\\end{array}\\right)=\\left( \\begin{array}{cc}{A} & {-\\mathbf{I}} \\\\ {\\mathbf{I}} & {\\mathbf{0}}\\end{array}\\right) \\left( \\begin{array}{c}{\\mathbf{w}_{j}} \\\\ {\\mathbf{w}_{j-1}}\\end{array}\\right)+\\sigma^{2} \\left( \\begin{array}{c}{\\mathbf{s}_{j}} \\\\ {\\mathbf{0}}\\end{array}\\right)

这样我们只需要考察另一个矩阵 A' = \\left( \\begin{array}{cc}{A} & {-\\mathbf{I}} \\\\ {\\mathbf{I}} & {\\mathbf{0}}\\end{array}\\right) 的特征值就可以了。

用传统的数值线性代数的思路当然没有问题，但是这里我们换一种分析方法。我们考虑设 \\lambda 为 A' 的特征值，并且设 (v_1, v_2)' 为对应的特征向量对，那么直接走矩阵运算，你能够得到， \\mu = \\lambda + \\frac1\\lambda 是 A 的特征值。而 A 的特征值我们用完全一样的思路可以知道它应该是落在 (2 - 4 \\sigma^2 , 2) 这个区间内。那么可以想一想，如果 2-4\\sigma^2 < -2 ，那么对应的 \\lambda 是完全可以取到一个 <-1 的值，这就没有办法使得数值解稳定。所以我们让 \\sigma^2 \\le 1 ，因为这个时候，可以发现 |\\lambda| = 1 ，才有可能 \\mu 取到我们的那个区间内（其实也不好说，因为两边习惯上是开区间，但是因为数值误差，写成闭区间倒也无妨）。虽然这个结果不是特别好，但是至少误差不会呈指数增大，所以我们接受这个结果。

通过与上面相同的误差分析，你可以得到，这个格式是 O(k^2+h^2) 的误差。

椭圆方程的有限差分方法

其实你也可以看出，上面的双曲方程和抛物方程的稳定性分析方法都是使用冯诺依曼分析法的，可以说是极为方便的。但是椭圆方程的分析并不是这样，它的方法有些另类。因此我们以它作为我们最后的结束的内容。

首先杭州人才市场我们研究的方程一般长这个样子

\\Delta u = f

这里的 \\Delta 是拉普拉斯算子，不是差分符号。考虑一下更特殊的情形，就是

u_{xx} + u_{yy} = f

那么如何作差分呢？其实这个地方方法倒是类似的。不过呢我们需要给一些比较有趣的名字。

五点格式法

我们还是一样考虑以 w_{ij} 来表示 (x,y) 处的数值解。这里因为我们有一个额外的函数 f ，所以我们也需要给它一个记号 f_{ij} ，并且同样的我们也假设 x 方向的步长为 h ， y 方向的步长为 k 。有了这些假设，我们的差分方法其实就可以写成这个形式

\\frac1{h^2}(w_{i-1, j} - 2 w _{i ,j} + w_{i+1, j}) + \\frac1{k^2}(w_{i, j-1} - 2w_{i ,j} + w_{i, j+1}) = f_{ij}

因为两个方向其实都是空间变量，我们有理由相信它们的性质比较类似。基于这个考虑我们不妨假设 k = h ，那么组合化简可以得到

\\frac{1}{h^{2}}\\left(w_{i-1, j}+w_{i+1, j}+w_{i, j-1}+w_{i, j+1}-4 w_{i j}\\right)=f_{i j}

下面这张图一定程度上解释了为什么我们称这样的差分方法为“五点格式法”

也许你会想，这个表达式写出来之后，我们也确确实实可以考虑用冯诺依曼分析法。但是问题在于，无论怎么排布，我们都没有办法获得诸如双曲和抛物方程那样好的结构。因此这边提出了一个非常有趣的方案叫做堆积 (stacking)，通俗说就是把所有的，每个点上的数值解排成一个列向量。比方说我们每一个维度上取了 m 个点，那么一共就会有 m^2 个点，向量自然也就是 m^2 维的了。另外排布的方式我们自然也有讲究，具体来说就是这样

w=\\left[ \\begin{array}{c}{w^{[1]}} \\\\ {w^{[2]}} \\\\ {\\vdots} \\\\ {w^{[m]}}\\end{array}\\right] ， w^{[j]}=\\left[ \\begin{array}{c}{w_{1 j}} \\\\ {w_{2 j}} \\\\ {\\vdots} \\\\ {w_{m j}}\\end{array}\\right]

按照我们这样的排列，容易发现我们的矩阵应该长这样

A=\\frac{1}{h^{2}} \\left普通话证书怎么考[ \\begin{array}{ccccc}{T} & {I} & { } & { } & { } \\\\ {I} & {T} & {I} & { } & { } \\\\ { } & {I} & {T} & {I} & { } \\\\ { } & { } & {\\ddots} & {\\ddots} & {\\ddots} \\\\ { } & { } & { } & {I} & {T}\\end{array}\\right] T=\\left[ \\begin{array}{ccccc}{-4} & {1} & { } & { } & { } \\\\ {1} & {-4} & {1} & { } & { } \\\\ { } & {1} & {-4} & {1} & {} \\\\ { } & { } & {\\ddots} & {\\ddots} & {\\ddots} \\\\ { } & { } & { } & {1} & {-4}\\end{array}\\right]

你也可以看出，这实际上是一个极为稀疏的大矩阵。但是还好，它的结构还算完整。而对应的方程其实就是 A\\mathbf{W} = \\mathbf{F} ，注意我们这里的结果因为不再与时间有关，所以实际上它不再是一个迭代的数列，这一点和前面两个有本质上的不同。

下面，我们来做一些稳定性分析吧。首先你也能看出来，因为没有时间变量，我们这里自然不能通过上面的计算谱半径的方法来判断是否稳定。所以这里实际上问题就落在了我们怎么处理这单独的一个矩阵方程上。

最直观的想法就是：考虑构造另外一个方程 A\\mathbf{W_2} = \\mathbf{F_2} ，只是说我们这里的 \\mathbf{W_2} 换成对应时间点的精确解。这当然是可以的，不过这样子得到的结果，后面计算的实际上就是精确解代入得到的一个式子。那么如果我们相减两式，得到 A(\\mathbf{W_2-W} )= \\m春丽athbf{F_2-F} ，那么左边的 \\mathbf{W_2-W} 就是我们的每一个点的整体截断误差。但是右边就需要我们单独计算了。

我们具体写出来右边式子的每一个元素，就是

\\tau_{i j}=\\frac{1}{h^{2}}\\left(u\\left(x_{i-1}, y_{j}\\right)+u\\left(x_{i+1}, y_{j}\\right)+u\\left(x_{i}, y_{j-1}\\right)+u\\left(x_{i}, y_{j+1}\\right)-4 u\\left(x_{i}, y_{j}\\right)\\right)-f\\left(x_{i}, y_{j}\\right)

相信你也知道该怎么做了吧？这就是最经典的局部截断误差的分析，它最后的结果是

\\tau_{i j}=\\frac{1}{12} h^{2}\\left(u_{x x x x}+u_{y y y y}\\right)+O\\left(h^{4}\\right)

你可以看到它的局部截断误差是二阶的。

把这些东西都放在右边，就可以得到我们最后的误差矩阵方程 A\\mathbf{E} = -\\tau 。容易看出，只要 \\|A^{-1}\\| 有界，那这个值 \\mathbf{E} ，它的规模应该是和 - \\tau 内的数相同的（想想为什么），所以我们下面就需要检查一个事情：是否 \\|A^{-1}\\| 有界。而一般来说根据范数的等价性，只要找到一个范数说明一下有界性即可。

我们还是考虑二范数，原因也很简单，这样子更容易分析。

首先这个矩阵是 m^2 \\times m^2 维的，因此也就存在 m^2 个特征值。我们不加证明的给出它的特征值形式。

\\lambda_{p, q}=\\frac{2}{h^{2}}((\\cos (p \\pi h)-1)+(\\cos (q \\pi h)-1))

注意我们这里的 (p, q) 相当于对应的 x,y 坐标。比方说 天鹅绒丝袜\\lambda_{1,1} 就相当于是 x,y 方向都在左（下）边界的对应的特征值。也就是对应 w_{11} 这个位置的特征值。

那么你也看出来，这个方程的分析不难的，因为我们有 h = \\frac1m ，这样的话可以看出，如果要证明 ||A^{-1}|| 有界，那么就一定是找到原来的那个矩阵 A 中，绝对值最小的那个做分析。这里我们自然取 \\lambda_{1,1} 。容易得到

\\lambda_{1,1}=-2 \\pi^{2}+O\\left(h^{2}\\right)

这是一个有界的数，所以我们自然可以得到结论说这个格式是稳定的。

最后，我们简单的提一下九点格式法。这个方法写出来是长这样的

\\begin{aligned} w_{i j}=\\frac{1}{6 h^{2}} &\\left[4 w_{i-1, j}+ w_{i+1, j}+4 w_{i, j-1}+4 w_{i, j+1}\\right.\\\\ &+w_{i-1, j-1}+w_{i-1, j+1}+w_{i+1, j-1}+w_{i+1, j+1}-20 w_{i j} ] \\end{aligned}

这个方法你也可以证明，它可以具有四阶的精度，不过需要要求 f = 0 ，具体的细节我们碍于篇幅，就不在这里说了。

小结

我们用很长的篇幅，给大家完整介绍了偏微分方程的有限差分的数值解法。大家可以看出来，虽然说偏微分方程是一个庞大的理论，但是在数值解领域，是存在这样的方法，可以系统的来解决很多很多问题的。这里主要就体现在贯穿全文的泰勒展开和冯诺依曼分析法中。

事实上，偏微分方程的有限元方法其实是更有活力的一大块。它们的构造其实更有意思，也更加实用。但是我们近期自然是没有机会说这么多了，就希望大家能够在这两篇文章中，详细而系统的窥探到有限差分方法的端倪。

——————————————————————————————————————

本专栏为我的个人专栏，也是我学习笔记的主要生产地。任何笔记都具有著作权，不可随意转载和剽窃。

个人微信公众号：cha-diary，你可以通过它来有效的快速的获得最新文章更新的通知。

专栏目录：笔记专栏|目录

想要更多方面的知识分享吗？欢迎关注专栏：一个大学生的日常笔记。我鼓励和我相似的同志们投稿于此，增加专栏的多元性，让更多相似的求知者受益~

本文标签： 微分方程解法数值