龙格-库塔法

概要

微分方程：含参数、未知函数、未知函数的导数（或者微分）的方程
数值求解：用若干离散点计算 近似值 来代替准确值
分类：单步法、多步法；隐式法、显示法
欧拉法 (欧拉折线法)，也是一阶龙格-库塔法：以矩阵面积代替曲面梯形面积
改进欧拉法 (预估——校正法)，也是二阶龙格-库塔法：以梯形面积代替曲面梯形面积
龙格-库塔法（重点）：yn+1 的值用 f(x,y) 在某些点上函数值的线性组合成来计算，计算 f(x,y) 的次数是龙格库塔的阶。从几何意义看：用多个斜率加平均计算叠加，从而逼近准确值
一般常用四阶龙格-库塔法。
阶数越高、步长越小，则计算精度越高，但是计算量越大
补充：对于2阶微分方程，则需要拆分成2个一阶微分方程，此时的 f(x,y) 变成了 f(x,y,y') 需要分别用R-K法计算 y 和 y'
局部截断误差：准确值-近似值
局部误差的阶（也是龙格-库塔法的阶）
亚当姆斯法（重点）：是线性多步法，利用已求出的多个值进行计算
拉格朗日插值：构造函数多项式 Ln(x) 来近似代替 y=f(x)
亚当姆斯法用拉格朗日插值的函数多项式 Ln(x) 近似替代 y′(x)
亚当姆斯公式有：显式公式、隐式公式
用高阶的亚当姆斯法（一般四阶）计算，得到的计算结果一般更精确，且计算量比龙格-库塔法少，唯一缺点是最前面几个值需要用其他方法求出

一、微分方程

微分方程的定义：凡是含有参数、未知函数、未知函数的导数（或者微分）的方程。

定义式： F(x,y,y′,y″,...,y(n))

常微分方程：未知函数是一元函数
偏微分方程：未知函数是多元函数

微分方程的阶：未知函数最高阶导数的阶数（注意和后面龙格-库塔法的阶做区分）。

例01：

$\frac{dy}{dx}=2x,y|_{x=1}=2$ 就是一个一阶微分方程。

例02：

$\frac{d^{2}S}{d^{2}t}=-0.4,S|_{t=0}=0,(\frac{dS}{dt})|_{t=0}=20$ 就是一个二阶微分方程。

二、常微分方程的数值求解

求解微分方程的困难:

求解微分方程往往计算量很大，甚至无法求解
在实际应用中，可以解出 近似值 即可——数值求解

2.1 数值求解

概念：利用给定常微分方程及边界条件，求解函数 y(x) 在若干 离散点 的 近似值（再拟合成折现或曲线）的方法

即：在区间 [a,b] 上有若干离散点： a= $x_{0}$ < $x_{1}$ <...< $x_{n}$ =b ，求出函数 y(x) 在离散点 $x_{k}$ 处的近似值 $y_{k}$ ，作为精确值 y( $x_{k}$ ) 的近似。

（这里强调一下： $x_{k}$ 是第k个离散点； $y_{k}$ 是近似值； y( $x_{k}$ ) 是精确值。）

2.2 数值求解法分类

1、分类方式一

单步法：计算 $y_{i+1}$ 时，只用到了前面的一个近似值 $y_{i}$ 。
多步法：计算 $y_{i+1}$ 时，用到了前面的多个近似值 $y_{i}$ 、 $y_{i-1}$ 、 $y_{i-2}$ ...。
一般单步法的运算量少，但是精度不高。

2、分类方式二

隐式法：方程两端均含有 $y_{i+1}$ ,（很难将其合并入公式的一端），即需要求解 隐函数方程 才能计算出 $y_{i+1}$ 。
显示法：方程可以直接计算出 $y_{i+1}$ 。

三、常用数值求解的具体方法

3.1 欧拉法（欧拉折线法）

令常微分方程：

${y}'(x)=f(x,y),y|_{x=0}=y_{0},(a=x_{0}\leq x_{1}\leq ...\leq x_{n}=b)$

起始点 $y_{0}$ 已知， h 等距步长 = $x_{k+1}$ − $x_{k}$ ，步长取的越小，运算精度越高，但运算量增大。

欧拉公式： $y_{k+1}=y_{k}+h*f(x_{k},y_{k})$

1、几何意义一：折线

用步长和前面一个点的斜率计算后面一个点。

（注意下图的纵坐标是 y ）

2、几何意义二：用矩形面积代替曲边梯形面积

注意：图中的纵坐标是 f(x,y) 即上一张图的斜率， y 是面积。

例03（欧拉法数值求解微分方程）

3.2 改进欧拉法（预估——校正法）

改进欧拉公式：

先用欧拉公式计算出 $\tilde{y_{k+1}}$ ——预估
再用 $\tilde{y_{k+1}}$ 求出真正的近似值 ${y_{k+1}}$ ——校正

$\tilde{y_{k+1}}=y_{k}+h*f(x_{k},y_{k}),y_{k+1}=y_{k}+\frac{h}{2}[f(x_{k},y_{k})+f(x_{k+1},\tilde{y_{k+1}})]$

为什么需要多加一步校验：因为改进欧拉法需要用到 ${y_{k+1}}$ ，是隐式法，需要解隐函数，计算量大。

几何意义：用 梯形面积 代替 曲边梯形 面积

改进欧拉法 的运算精度一般比 欧拉法 高。

3.3 龙格——库塔法（Runger—Kutta法，简称R-K法）（重点）

将函数 y(x) 在 xn 处泰勒展开：

$y(x)=y(x_{n})+{y(x_{n})}'(x-x_{n})+\frac{{y(x_{n})}''}{2!}(x-x_{n})^{2}+\frac{y^{(3)}(x_{n})}{(x-x_{n})^{3}}+...$

n阶龙格-库塔法：取泰勒展开的前n阶导数项

将 $x={x_{n+1}}$ 带入 $y(x)$ 计算出近似值 $y_{n+1}$ 。

1、欧拉法和改进欧拉法是龙格-库塔法的一种：

一阶龙格-库塔法：欧拉法
只取前2项，即线性部分 $y(x)\approx y(x_{n})+{y(x_{n})}'(x-x_{n})$
二阶龙格-库塔法：改进欧拉法
取前3项

具体证明过程如下：（数学不好的人可以忽略）------------------------------------------------

--------------------------------------------------------------------------------------------

2、用龙格-库塔法的一般形式表示欧拉法、改进欧拉法

（1）改写欧拉法： $y_{n+1}=y_{n}+h*f(x_{n},y_{n})$
改写成 $y_{n+1}=y_{n}+hK_{1},K_{1}=f(x_{n},y_{n})$

（2）改写改进欧拉法：

$y_{n+1}=y_{n}+hf(x_{n},y_{n})+\frac{h}{2}f(x_{n+1},y_{n+1})$
改写成 $y_{n+1}=y_{n}+\frac{h}{2}(K_{1}+K_{2}),K_{1}=f(x_{n},y_{n}),K_{2}=f(x_{n}+h,y_{n}+hK_{1})$

龙格-库塔法计算顺序： $K_{1},K_{2},...,K_{n}$ 最后计算 $y_{n+1}$ 。

3、误差分析——局部截断误差和阶

假设第n步的计算是准确的，即 $y_{n+1}=y_{n}$ 。

局部截断误差： $o(h^{p+1})=y(x_{n+1})-y_{n+1}$ ，即：精确值-近似值。

阶：上边公式中的 $p$ 。

阶越大，局部误差就越小，近似值就越接近精确值，但同时运算量更大。

欧拉法和改进欧拉法的局部截断误差和阶

（1）欧拉法：没有泰勒展开到二阶导数的，所以精确值-近似值就会出现二阶导数的部分，二阶导数的部分为 $\frac{{y(x_{n})}''}{2!}(x-x_{n})^{2}$ ，有 $(x-x_{n})^{2}$ ，即 $h^{2}$ 。

局部截断误差： $o(h^{2})$
阶：1阶

（2）改进欧拉法：

局部截断误差： $o(h^{3})$
阶：2阶

3.3.4 龙格-库塔法的一般形式

p阶龙格-库塔法： $y_{n+1}$ 的值用 $f(x,y)$ 在某些点上函数值的线性组合成来计算，计算 $f(x,y)$ 的次数是龙格-库塔法的阶，也是局部截断误差的阶。

从几何意义上看：龙格-库塔法可以看成用 $p$ 个斜率 $K_{r}=f(x_{r},y_{r})$ 的一种 加权平均逼近 。如：欧拉法只用了一个斜率 $f(x_{n},y_{n})$ ，即 $K_{1}$ ；而改进欧拉法是 $\frac{h}{2}[f(x_{k},y_{k})+f(x_{k+1},\tilde{y_{k+1}})]$ ，即 $\frac{h}{2}K_{1}+\frac{h}{2}K_{2}$ 。

四阶龙格-库塔法（常用）：

$y_{n+1}=y_{n}+\frac{h}{6}(K_{1}+2K_{2}+2K_{3}+K_{4}),K_{1}=f(x_{n},y_{n}),K_{2}=f(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{1}),K_{3}=f(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{1}),K_{4}=f(x_{n}+h,y_{n}+hK_{3})$

（再次强掉： $K_{1},K_{2},K_{3},K_{4}$ 前面的系数是 加权平均 得到，计算顺序是从 $K_{1}$ 到 $K_{4}$ ，根据这个规律特征，可以推出后续更高阶的龙格-库塔法）。

问：为什么需要推出，不能直接写出其 $p$ 阶的一般表达式吗？
回答：和”斐波那契数列“是一个道理，因为一般表达式很复杂，还不如推的简单。其一般形式如下：

例04（4阶龙格-库塔法数值求解微分方程）

在用4阶龙格-库塔法时，可以适当增大步长h，减少计算容量，一般精度也不会损失很少。

实际选择方法时，可以根据具体情况选择适当的阶数和 步长h。

--------------------------------------------------------------------------------------------

补充：二阶微分方程的 4阶R-K法

对于2阶微分方程，则需要拆分成2个一阶微分方程，此时的 $f(x,y)$ 变成了 $f(x,y,{y}')$ 需要分别用R-K法计算 $y$ 和 ${y}'$ 。

已知初值 $y(0)$ 和 ${y}'(0)$ ,以及 $y{}''$ 关于 $x,y,y'$ 的函数，即 ${y}''=f(x,y,{y}')$ 。

例题：用4阶R-K法求下面的二阶微分方程的数值解，步长 h 取0.1。

$2\frac{d^{2}y}{dt}+(y^{2}-1)\frac{dy}{dt}+y=0,y(0)=0.25,{y}'(0)=0$

即

${y}''=-\frac{y+{y}'(y^{2}-1)}{2},y(0)=0.25,{y}'(0)=0$

解：将二阶微分方程拆成2个一阶微分方程，R-K法计算 $y$ 和 ${y}'$ 。

其中： $K_{mn}$ 表示第 $n$ 次迭代， $m$ 次导数作为 $f(x,y)$ ；

其中： $f_{1}={y}',f_{2}={y}''=-\frac{y+{y}'(y^{2}-1)}{2}$

然后根据上边的步骤，一步步计算出 $y_{n+1}$ 和 ${y}'_{n+1}$ 。

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

龙格-库塔法的缺点：其是单步法，计算 $y_{n+1}$ 时只用到了 $y_{n}$ ,因此需要迭代多次才能得到较精确的结果。

四、亚当姆斯法（Adams法）

4.1 线性多步法

利用已求出的 $y_{n}$ 、 $y_{n-1}$ 、 $...$ 、 $y_{n-q}$ （通过单步法，如：龙格-库塔法）共 $q+1$ 个点处的函数值构造 $y_{n+1}$ 的线性计算式。

4.2 拉格朗日插值

在实际项目中，函数 $y=f(x)$ 可能比较复杂，因此用一个简单的函数 $L(x)$ 近似代替 $f(x)$ 。

我们知道任何一个解析函数（即函数在给定复平面内处处可微）都可以用一个唯一的多项式表示：

$L_{n}(x)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n}x^{n}$

我们就是用这个函数多项式 $L_{n}(x)$ 近似替代 $y=f(x)$ ，即对于每一个 $x_{i}$ ，有 $L_{n}(x_{i})=y_{i},(i=0,1,2,...,n)$

可以发现， $n$ 越大，则 $L_{n}(x)$ 的表达式越复杂，但是更能逼近原函数。

怎么求 $L_{n}(x)$ ？

$L_{n}(x)$ 的具体表达式是通过已知点（前面求过的点）的坐标表示的。

如：一次插值( n=1 )，已知 $(x_{k},y_{k})$ 、 $(x_{k+1},y_{k+1})$

则 $L_{1}(x)=\frac{x-x_{k+1}}{x_{k}-x_{k+1}}y_{k}+\frac{x-x_{k}}{x_{k+1}-x_{k}}y_{k+1}$ 是一条直线。

二次插值( n=2 )，已知 $(x_{k-1},y_{k-1})$ 、 $(x_{k},y_{k})$ 、 $(x_{k+1},y_{k+1})$

$L_{2}(x)=\frac{(x-x_{k})(x-x_{k+1})}{(x_{k-1}-x_{k})(x_{k-1}-x_{k+1})}y_{k-1}+\frac{(x-x_{k-1})(x-x_{k+1})}{(x_{k}-x_{k-1})(x_{k}-x_{k+1})}y_{k}+\frac{(x-x_{k-1})(x-x_{k})}{(x_{k+1}-x_{k-1})(x_{k+1}-x_{k})}y_{k+1}$ 是一条抛物线。

拉格朗日插值多项式一般形式：

4.3 亚当姆斯法：线性多步法之一（插值方法）

还是认为 $f(x,y)={y}'(x)$

我们用拉格朗日插值多项式 $L_{n}(x)$ 表示 $f(x,y)$ 。

注意：这里的 $f(x,y)$ 是原函数的导数。

则：根据面积 $y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(x,y)dx$

用 $L_{n}(x)$ 替代上边公式的 $f(x,y)$ 。

4.4 亚当姆斯法有2种公式：显式Adams公式、隐式Adams公式

如：要求 $y_{n+1}$

显式公式：拉格朗日多项式不使用 $y_{n+1}$ ，根据 $y_{n}$ 、 $y_{n-1}$ 、 $...$ 求。一般用于求未求过的 $y_{n+1}$
隐式公式：拉格朗日多项式使用 $y_{n+1}$ ，根据 $y_{n+1}$ 、 $y_{n}$ 、 $...$ 求。一般用于预测-校验之前已经求过的 $y_{n+1}$

例：多项式都取一阶插值，即公式是二阶的

1.二阶隐式Adams公式——梯形公式

取一次插值 $L_{1}(x)=\frac{x-x_{n+1}}{x_{n}-x_{n+1}}{y}'(x_{n})+\frac{x-x_{n}}{x_{n+1}-x_{n}}{y}'(x_{n+1})$
(还是提醒一下： $f(x,y)={y}'(x)$ ）
则： $\int_{x_{n}}^{x_{n+1}}L_{1}(x)dx=\frac{h}{2}[{y}'(x_{n})+{y}'(x_{n+1})]$
$y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n+1})]$
可以发现，二阶隐式Adams公式，也是梯形公式

图形解释：用梯形面积 (阴影) 近似替代梯形曲面

注意下图的纵坐标是原函数的斜率 $f(x,y)$ ,即图中的面积才是 $y(x)$ 。

用阴影面积近似替代梯形曲面

2.二阶显式Adams公式

取一次插值 $L_{1}(x)=\frac{x-x_{n-1}}{x_{n}-x_{n-1}}{y}'(x_{n})+\frac{x-x_{n}}{x_{n-1}-x_{n}}{y}'(x_{n-1})$
则： $\int_{x_{n}}^{x_{n+1}}L_{1}(x)dx=\frac{h}{2}[3{y}'(x_{n})-{y}'(x_{n-1})]$
$y_{n+1}=y_{n}+\frac{h}{2}[3f(x_{n},y_{n})-f(x_{n-1},y_{n-1})]$