数值计算方法(3) 数值微分方法

+++
date = '2024-12-21T15:45:47+08:00'
draft = true
title = '数值计算方法(3) 数值微分方法'
+++

初次发布于我的个人文档

上一期讲了数值积分方法，这一次自然是要讲数值微分方法的，不然太不完善了。

更何况数值微分方法其实是基于数值积分方法得到的。

我们先从比较简单的估计导数值来开篇。

差商方法

估计函数的导数值意义其实没有前面积分那么大，因为基本初等函数都能很容易地求导，有准确值了近似值的意义就没有那么大了。

不过还是来试试吧。

近似导数的最简单的方法就是用导数的定义。

你看一阶差商$\frac{f(a+h)-f(a)}{h}$，当h趋于0时根据这就是导数的定义啊，所以反过来我们可以让h非常小，这样就是导数的近似值了。

这种差商在这里有一个名字叫向前差商。

有向前自然就有向后$\frac{f(a)-f(a-h)}{h}$。

其实还有中心差商$\frac{f(a+h)-f(a-h)}{2h}$。

中心差商还有一个名字叫中点方法，其实这玩意也能看做是向前差商和向后差商的平均值。

如果是我们手算的话，其实这样就很好了。

我们手算的话h越小误差越小，但这么麻烦的东西自然是不可能手算的。

但是一旦交给计算机算啊，h如果太小，计算机计算的误差就会增加（计算机进行浮点数运算是不准确的，这是由浮点数的存储方式决定的，计算机的浮点数大多都按照IEEE的规范设计的，是存在舍入误差的）。

总之，现在这个h如果太大，从数学上看的误差会变大（我们称为截断误差），h如果太小计算机计算的时候舍入误差会太大。

那么又来了，步长h应该怎么选呢？

这情况和龙贝格算法一模一样啊。

还是用理查德森外推法。

我们可以导出一个一模一样的加速公式。

$G_n(h)=\frac{4^{nG_{n-1}(\frac{h}{2})-G_{n-1}(h)}{4}n-1}$

然后画一模一样的表就可以了。

这里确实没太大区别。

我们先用泰勒公式计算余项：

$G(h)-f'(a)=\alpha_1 h^2+\alpha_2h4+...+\alpha_nh^{2n}+...$

把步长折半

$T(\frac{h}{2})-f'(a)=\alpha_1 \frac{h^{2}{4}+\alpha_2\frac{h}4}{16}+...+\alpha_n\frac{h^{2n}}{4n}+...$

然后消掉h²项得到更好的$f'(a)$的误差估计，每次都这样计算，整理出来的递推公式就是上面的那个式子，由于余项的形式和龙贝格算法是一样的所以最终的递推式也是一样的。（详情还是回到上一节看一下详细的理查德森外推法吧）

插值型求导公式

前面我们说插值是对函数的拟合，所以这里自然可以用插值多项式替代函数f(x)，用插值多项式的导数作为f(x)导数的近似，但是这玩意儿没那么靠谱，因为插值多项式在函数值上比较接近，但是他们的导数值可能差得很远，所以这条路子没那么靠谱，也就没那么多理论了。你就当是半破产了吧。而且操作起来也不难，就是无脑用插值多项式替代f(x)而已，所以不多说了。

接下来才是重点，唠唠怎么求解微分方程。

微分方程数值解

微分方程数值解的求法主要是利用了数值求积方法，在此我们以标准的微分方程初值问题为例介绍几个方法。

在《常微分方程》这门课程中（注意是一门课，不是高等数学的那个章节），你应该听过皮卡逐步逼近定理。

它这里的微分方程近似解是一个多项式，而我们这里说的微分方程数值解不是这样的东西。

针对这样的初值问题，你告诉我你想求微分方程的解这个函数，在哪个点的取值，然后我自己拟定一个步长h，我求出微分方程的解g(x)在x+h，x+2h,...等处的函数值，最后我告诉你你要的那个点的函数值，我们这里说的数值解是这样的东西。

当然，为了保证初值问题有且只有唯一的解，我们仍然要求f(x,y)关于y满足利普希茨条件。（没听过不用管，这是因为你没上过《常微分方程》这门课，总之就是要保证解存在且唯一而已，这个条件后文不会使用）

好接下来开始推导。

对微分方程

$\frac{dy}{dx}=f(x,y)$

我们两边同时积分得到

$\int_{x_n}^{{x_{n+1}}\frac{dy}{dx}=\int_{x_n}}{x_{n+1}}f(x,y)dx$

也就是

$y_{n+1}-y_n=\int_{x_n}^{x_{n+1}}f(x,y)dx$

从而

$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)dx$

而如果我们可以用数值方法求出这个积分，不就解出了这个微分方程？

例如，我们直接用左矩形公式求积分，得到

$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)=y_n+(x_{n+1}-x_n)f(x_n,y_n)$

这就是欧拉格式。(注意，数值微分里我们称为格式而不是公式，这里没有打错字)

你可能会发现我这里给的欧拉格式和网上一般的不一样，其实你只要设$h=x_{n+1}-x_n$（当我设了h就默认等步长了）就能得到一般的欧拉格式了。

$y_{n+1}=y_n+hf(x_n,y_n)$

同样地，我们可以用右矩形公式得到

$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)=y_n+(x_{n+1}-x_n)f(x_{n+1},y_{n+1})$

类似地，我们设$h=x_{n+1}-x_n$就有，

$y_{n+1}=y_n+hf(x_{n+1},y_{n+1})$

这就是隐式欧拉格式了。

那如果我用中矩形公式，还能得到两步欧拉格式：

对微分方程

$\frac{dy}{dx}=f(x,y)$

我们两边同时积分得到

$\int_{x_n-1}^{{x_{n+1}}\frac{dy}{dx}=\int_{x_n-1}}{x_{n+1}}f(x,y)dx$

注意，这次是从$x_{n-1}$积到$x_n$，并且我们要求步长为h，则$x_n$是积分区间中点。

有$y_{n+1}-y_{n-1}=2hf(x_n,y_n)$

从而，$y_{n+1}=y_{n-1}+2hf(x_n,y_n)$

这就是两步欧拉格式。

还可以用梯形公式来求积分得到梯形格式。

$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})]$

注意力好的人应该已经发现了，梯形格式恰好是显式欧拉格式和隐式欧拉格式的算术平均。

预报-矫正系统与改进欧拉

前面扯了这么久感觉非常平淡啊，就是套用一下求积公式，下面来点有意思的。

你看梯形格式里是不是说需要$f(x_{n+1},y_{n+1})$？

你有没有想过这玩意儿哪里来？

是不是前面看得特爽但是这个最基本的东西没发现。

$f(x_{n+1},y_{n+1})$明明是下一次的东西，你怎么能让我现在就给出来呢？

其实是这样的，我们先用前面的一个显式格式预报一个$f(x_{n+1},y_{n+1})$出来，然后再用隐式格式得到一个矫正后的$y_{n+1}$，这样精度会更高。

例如，使用显式欧拉格式预报$y_{n+1}=y_n+hf(x_n,y_n )$

然后再把这个预报的$y_{n+1}$代入隐式欧拉，这就是改进欧拉格式了。

这样的一套系统我们称为预报-矫正系统。

你其实是可以任意组合和重复的。

你先通过一个显式格式预报第一次$y_{n+1}$。

然后利用一个隐式格式预报$y_{n+1}$，

再用任何一个隐式格式预报，哪怕是刚刚用过的也可以。

这样一直预报下去，我们得到的$y_{n+1}$的精度就会不断提高。

但是这样的方法仍然没有脱离求积公式的桎梏，而下面的方法就是重新开了一个方向了。

龙格-库塔方法

这个龙格-库塔方法从另一个角度看问题，他是这么说的。

函数f的一阶差商$\frac{f(x_{n+1})-f(x_{n})}{x_{n+1}-x_n}$，因为我们微分方程数值解指的是你给我一个步长h，我求出微分方程的解g(x)在x+h，x+2h,...等处的函数值，所以分母就是步长h，分子则可以运用拉格朗日中值定理，从而$\frac{f(x_{n+1})-f(x_{n})}{h}=f'(\alpha)$

于是，$f(x_{n+1})=f(x_n)+hf'(\alpha)$

你们呐，折腾这么多其实就是在估计$f'(\alpha)$！

这个$f'(\alpha)$被称为函数f在区间内的平均斜率。

你可以试试看，显式欧拉就是用函数在左端点的斜率值估计平均斜率，改进欧拉是用两个端点处的斜率的平均值估计斜率。

这些估计方式对平均斜率的估计显然精度没有那么高，所以这些方法在计算数值解时的精度自然也不高。

那我如果多用一些点来估计函数的平均斜率，精度不就提高了？

龙格-库塔方法就是这样的方法。

我们在$x_n,x_{n+1}$之间取p个不同点的斜率，计算他们的加权平均，把这个当做平均斜率的估计。

那取哪些点呢？这个就是龙格-库塔方法精度高的精髓了。

我们以二阶龙格-库塔方法为例进行推导。

在这里，龙格-库塔方法的阶数指的就是用的点的数量。

在这里我们设要取的点的斜率为$k_1,k_2$，要计算他们的加权平均，所以还需要假设权。

所以我们再设平均斜率$k^*=(1-\lambda)k_1+\lambda k_2$

我们先取其中一个端点为左端点$x_n$试试水，那么就有$k_1=f(x_n,y_n)$了。

接下来我们要找第二个点$x_{n+p}=x_n+ph,p \in(0,1]$然后计算他的斜率为$k_2$。

那怎么找呢？似乎没有头绪了，其实我也没有。不过好在数学家想到了，可以继续用预报-矫正系统啊。

我们先用显式欧拉预报一个$y_{n+p}=y_n+phk_1$，那么$k_2=f(x_{n+p},y_{n+p})=f(x_{n+p},y_n+phk_1)$了。

呕吼，搞了一圈合着还是俩未知量，一个没消啊。

这个方法里真正消元的方法是泰勒展开+比较系数。

我们先对$k_2=f(x_{n+p},y_n+phk_1)$进行泰勒展开，由于$x_{n+p}=x_n+ph$，所以$k_2$实际上是一个二元函数$f(x,y)$的泰勒展开。

这个二元函数的泰勒展开，好像只有数学分析介绍了，如果你只学过高等数学的话你应该不知道，所以我这边给一下结论。

如果函数f在点$P(x_0,y_0)$的邻域U(P)内有n阶连续的偏导数，则对U(P)内的任意一点$(x_0+h,y_0+k)$，存在$\theta \in (0,1)$，使得

$f(x_0+h,y_0+k)=\sum_{p=0}^n\frac{1}{p!}(h\frac{\partial}{\partial x}+k\frac{\partial}{\partial y})^{pf(x_0,y_0)+o(\rho),\rho}n=\sqrt{h^2+k2}$

这是二元函数的n阶泰勒公式。

现在我们来展开一下但是这里只需要展开两项，

$k_2=f(x_n+ph,y_n+phk_1)=f(x_n,y_n)+(ph\frac{\partial}{\partial x}+phk_1\frac{\partial}{\partial y})f(x_n,y_n)+o(\rho ^2)$

现在只是单纯的代入了公式，我们来进行计算。

先关注偏导数部分。

$(ph\frac{\partial}{\partial x}+phk_1\frac{\partial}{\partial y})f(x_n,y_n)$

这里前面一个括号写的是偏微分算子，其实我们直接把函数$f(x,y)$乘进去就可以了。

也就得到$ph\frac{\partial}{\partial x}f(x_n,y_n)+phk_1\frac{\partial}{\partial y}f(x_n,y_n)$

由于没有给出f(x,y)具体的表达式，所以这里就不能计算f(x,y)的偏导了。

把这个结果代入就是最后的展开式了。

$k_2=f(x_n+ph,y_n+phk_1)=f(x_n,y_n)+[ph\frac{\partial f(x_n,y_n)}{\partial x}+phk_1\frac{\partial f(x_n,y_n)}{\partial y}]+o(\rho ^2)$

等等，$k_1$是啥？

回到前面看看，$k_1=f(x_n,y_n)$

而我们解的微分方程是啥来着？$\frac{dy}{dx}=f(x,y)$

$x_n,y_n$既然是微分方程解上的点自然也就满足这个微分方程了喽，所以$k_1=f(x_n,y_n)=\frac{dy}{dx}|_{x=x_n}$

这么写有点难看，我们把式子弄得漂亮一点。微分方程最后的解是一个关于x的函数，我们设为$y(x)$，那么他的导数就是$y'(x)$，他在$x_n$处的导数自然就是$y'(x_n)$了。

所以$k_1=y'(x_n)$

也因此

$k_2=f(x_n+ph,y_n+phk_1)=f(x_n,y_n)+[ph\frac{\partial f(x_n,y_n)}{\partial x}+phy'(x_n)\frac{\partial f(x_n,y_n)}{\partial y}]+o(\rho ^2)$

再等等，刚刚是不是还说了$k_1=f(x_n,y_n)=y'(x_n)$来着？

这下不是又能代入了？

$k_2=f(x_n+ph,y_n+phk_1)=y'(x_n)+[ph\frac{\partial f(x_n,y_n)}{\partial x}+phy'(x_n)\frac{\partial f(x_n,y_n)}{\partial y}]+o(\rho ^2)$

这时候有些注意力比较集中的人就发现了。

k2的式子里面啊，中括号的内容是$ph\frac{\partial f(x_n,y_n)}{\partial x}+phy'(x_n)\frac{\partial f(x_n,y_n)}{\partial y}$

我们把ph提出来，就是$ph[\frac{\partial f(x_n,y_n)}{\partial x}+y'(x_n)\frac{\partial f(x_n,y_n)}{\partial y}]$

你看中括号里面的东西。

额，也许你看不出来。我给个提示。

观察二元函数f(x,y)，这里我们已经假设了y是关于x的函数y(x)，所以这时候我们可以求f(x,y)关于x的全导数，利用复合函数求导法则就可以了。

或者画个路径图。

f(x,y)首先是关于x,y的二元函数，y又是关于x的一元函数，

所以f(x,y)关于x的全导数$\frac{df(x,y)}{dx}=\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}\frac{dy}{dx}$

我把式子写全

$\frac{df(x,y)}{dx}=\frac{\partial f(x,y)}{\partial x}+\frac{\partial f(x,y)}{\partial y}\frac{dy(x)}{dx}$

我取f(x,y)在点$(x_n,y_n)$的全导数，

这特么，不就是括号里的$\frac{\partial f(x_n,y_n)}{\partial x}+y'(x_n)\frac{\partial f(x_n,y_n)}{\partial y}$吗？

所以又能化简了。

$k_2=f(x_n+ph,y_n+phk_1)=y'(x_n)+ph\frac{df(x_n,y_n)}{dx}+o(\rho ^2)$

我再问一遍，我们要求解的微分方程是什么？

$\frac{dy}{dx}=f(x,y)$

所以，$f(x_n,y_n)=\frac{dy(x_n)}{dx}$

那$f(x_n,y_n)$再对x求一阶导是啥？

不就是$y(x)$在$x_n$处的二阶导$y''(x_n)$吗？

所以又能化简了。

$k_2=f(x_n+ph,y_n+phk_1)=y'(x_n)+phy''(x_n)+o(\rho ^2)$

最后还有个小尾巴，$o(\rho^2)$是啥？

回去看看泰勒公式，$\rho^2=(ph)2+(phk_1)^{2$，它和$h}2$是同阶无穷小，所以$o(\rho^2)=o(h2)$

所以到了最后，

$k_2=f(x_n+ph,y_n+phk_1)=y'(x_n)+phy''(x_n)+o(h^2)$

$k_1$我们说了是$f(x_n,y_n)=f'(x_n)$

代入我们的迭代格式$y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2]$

得到

$y_{n+1}=y_n+hy'(x_n)+\lambda ph^{2y''(x_n)+o(h}3)$

而$y_{n+1}$实际上是y(x)在$x_{n+1}$的取值$y(x_{n+1})=y(x_n+h)$

这时候我们还可以把一元函数$y(x)$泰勒展开，展开三项得到

$y_{n+1}=y(x_n+h)=y(x_n)+hy'(x_n)+\frac{h^{2}{2}y''(x_n)+o(h}3)$

由于泰勒展开是唯一的，所以我们前后得到的两个泰勒展开应该是一样的才对，也就是$\lambda p=\frac{1}{2}$

总之，最终结论是$2\lambda p=1$，并且二阶龙格-库塔方法的阶段误差是$h^3$这个数量级的。

也就是说，只要满足$2\lambda p=1$的迭代格式都是二阶龙格-库塔格式，误差都是$h^3$这个数量级的。

我们取$p=1,\lambda = \frac{1}{2}$得到的迭代格式就是改进欧拉格式，取$p=\frac{1}{2},\lambda=1$得到的迭代格式我们称为变形的欧拉格式，或者也称为中点格式。

$y_{n+1}=y_n+hk_2$

$k_1=f(x_n,y_n)$

$k2=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}k_1)$

如果我们令$p=\frac{2}{3},\lambda=\frac{3}{4}$，得到的就是休恩公式。

如果进一步增加点就可以得到三阶、四阶乃至更高阶的龙格-库塔迭代格式不过这个迭代式子啊比较难看了而且也没有什么规律，想了解的话可以自己上网搜寻。他的推导思路倒是和二阶的没什么区别。

龙格-库塔方法的问题呢已经很小了。不过我们还是可以吹毛求疵出来，在计算$y_{n+1}$的时候我们已经算出了$y_1,y_2,...,y_n$而龙格-库塔方法没有充分地利用这些已经算出来的数据。

亚当姆斯方法

如果要利用这些信息你就得到了亚当姆斯方法。

亚当姆斯方法可以当做龙格-库塔方法的推广继续用泰勒展开得到，但是你懂的，这非常麻烦。

我们换一个简单的方法——利用插值多项式得到。

针对微分方程$\frac{dy}{dx}=f(x,y)$，我们两边积分（熟悉的开局）

但是这次积分区间不一样，看你想用几个历史信息点了，比如说我想用r+1个历史信息，那就从$x_{n-r+1}$积到$x_{n+1}$，那么左边是函数值的差，右边是一个积分。

这个积分呢，我用一个r次插值多项式来近似。

那插值节点选哪r个呢？答案是任意。如果我选择包含$x_{n+1}$的前r+1个插值节点得到的就是r+1阶隐式亚当姆斯格式，我选择$x_n$开始的r+1个插值节点得到的就是r+1阶显式亚当姆斯格式。

例如，我想得到2阶隐式亚当姆斯方法，我就从$x_n$积到$x_{n+1}$

得到$y_{n+1}-y_n=\int_{x_n}^{x_{n+1}}f(x,y)dx$

再用过$x_{n+1},x_n$的线性插值多项式$\frac{x-x_{n+1}}{x_n-x_{n+1}}f(x_n,y_n)+\frac{x-x_n}{x_{n+1}-x_{n}}f(x_{n+1},y_{n+1})$近似替代f(x,y)然后算出积分就可以了。

结果是$y_{n+1}=y_n+\frac{h}{2}(f_n+f_{n+1})$

那如果我取插值节点是$x_n,x_{n-1}$就能得到二阶显式亚当姆斯格式，结果是$y_{n+1}=y_n+\frac{h}{2}(3f_n-f_{n-1})$

亚当姆斯预报-矫正系统

很经典的场景出现了。

我们同时有显式和隐式迭代格式，而隐式格式里需要提前知道$y_{n+1}$，那$y_{n+1}$哪里来呢？

可以由同阶显式迭代格式预报得来。

这就是亚当姆斯预报-矫正系统，没啥花头。

前面说了这么多都是关于迭代格式，也就是算法的讨论，那这部分有没有理论呢？

迭代格式的收敛性和数值稳定性

其实有的。

一个最严重的问题是收敛性问题，也就是我们叭了叭叭bb了老半天，你迭代的$y_n$是否真的趋于准确值。

这些方法我都拿出来讲了你应该明白，肯定是收敛的。

这玩意儿的证明又涉及到了利普希茨条件，又要搞泰勒展开的，比较麻烦，我就给个文章吧。

知乎的@123456写的关于一阶显式和隐式方法收敛性的理论分析，地址https://zhuanlan.zhihu.com/p/515389776，当然你去网上找其他资料也可以。

我这里只给个结论。

显式欧拉格式的全局截断误差$|y(x_n)-y_n|\le e^{{L(b-a)}|y(x_0)-y_0|+\frac{Mh}{2L}(e})$，其中$M=\sup|y''(x)|$

这个全局截断误差有两项，我们说后一项是局部截断误差的累计，而在前一项里，$|y(x_0)-y_0|$是初值条件的误差称为条件误差，e^{L(b-a)}$是条件误差放大倍数，他只于初值问题本身的利普希茨系数L有关。当这个L很大的时候，注意，哪怕是微小的条件误差也会导致显式欧拉格式有很大的全局截断误差，这种情况我们称原问题数值稳定性差。

当然，这只是定性的说法，如果要定量的话是这样说的。

我们给一个扰动$\delta$，如果迭代格式产生的各个节点的误差均在$\delta$以内，我们就说这个迭代方法是稳定的。

但是这个稳定性问题是非常复杂的，这个问题在其他学科还有一个名字叫混沌。

2021年的诺贝尔物理学奖就颁给了混沌相关的研究，咱们就不蹚这深水了。

我可以告诉你的是，显式欧拉格式是条件稳定的，只有当步长h比较小的时候才稳定。

关于方程组和边值问题

我们扯了那么多，都是在说微分方程的初值问题。那方程组呢？

方程组的话其实迭代格式不用改变的，方程组的初值问题不就是$\frac{d\vec{y} }{dx}=f(x,\vec{y})$嘛，你看写成向量不就熟悉了？

迭代格式不用发生任何变化。

总之，这样就能求解常微分方程组了。

而如果是高阶常微分方程，我们知道可以通过换元的方法转化成方程组的求解所以也解决了。

而微分方程的边值问题的数值解，上面的方法都用不了，而且求解起来比较复杂。

甚至哪怕是上了《常微分方程》课程的人也未必知道微分方程的边值问题是啥吧。

所以我就只给个名词了，可以用差商法、有限元法、打靶法等方法求解。