常微分方程算法之阿当姆斯法（Adams法）

猿何试Bug个踌

已于 2024-04-05 00:05:47 修改

阅读量2k

点赞数 13

分类专栏：常微分方程文章标签：算法

于 2024-04-05 00:03:34 首次发布

本文链接：https://blog.csdn.net/L_peanut/article/details/137378532

版权

常微分方程专栏收录该内容

9 篇文章 3 订阅

订阅专栏

一、方法背景

欧拉法、梯形法、改进欧拉法以及龙格-库塔法都属于单步法，即计算 $y_{i+1}$ 的时候只用到上一步 $y_{i}$ 的值。如果可以将已经计算出来的结果值 $y_{0},y_{1},y_{2},\cdot \cdot \cdot ,y_{i}$ 用于 $y_{i+1}$ 的计算，就能够增加计算精度。我们尝试使用前面已经计算出来的k个值 $y_{i-k+1},y_{i-k+2},\cdot \cdot \cdot ,y_{i-1},y_{i}$ 来提高求解 $y_{i+1}$ 的精度，即线性多步法的思路。通常线性k步法的公式为：

$y_{i+1}=\sum_{j=0}^{k-1}\alpha_{j}y_{i-j}+h\sum_{j=-1}^{k-1}\beta_{j}f_{i-j} \space\space\space\space (1) \\ \space\space\space\space = \alpha_{0}y_{i}+\alpha_{1}y_{i-1}+\cdot \cdot \cdot +\alpha_{k-1}y_{i-k+1}+h(\beta_{-1}f_{i+1}+\beta_{0}f_{i}+\beta_{1}f_{i-1}+\cdot \cdot \cdot +\beta_{k-1}f_{i-k+1})$

式中 $f_{j}=f(x_{j},y_{i})$ 且 $x_{j}=x_{0}+jh,i=k-1,k,k+1,\cdot \cdot \cdot ,n-1$ ， $\alpha$ 和 $\beta$ 为待定系数。如果 $k=1$ ，就是单步法：

$y_{i+1}=y_{i}+h(\beta_{-1}f_{i+1}+\beta_{0}f_{i})$ ，其中 $\beta_{-1}+\beta_{0}=1$

若 $\beta_{-1}=0,\beta_{0}=1$ ，得到 $y_{i+1}=y_{i}+hf_{i}=y_{i}+hf(x_{i},y_{i})$ ，即为欧拉法；

若 $\beta_{-1}=1,\beta_{0}=0$ ，得到 $y_{i+1}=y_{i}+hf_{i+1}=y_{i}+hf(x_{i+1},y_{i+1})$ ，即为隐式欧拉法；

若 $\beta_{-1}=\beta_{0}=\frac{1}{2}$ ，得到 $y_{i+1}=y_{i}+\frac{1}{2}h(f_{i+1}+f_{i})\\ \space\space\space\space=y_{i}+\frac{1}{2}h(f(x_{i+1},y_{i+1})+f(x_{i},y_{i}))$ ，即为梯形法。

同时，如果 $\beta_{-1}=0$ 时是显格式；如果 $\beta_{-1}\neq 0$ 时是隐格式。

龙格-库塔法的原理是利用函数在某点处的函数值的线性组合来计算导数，而线性多步法的原理则是利用已经算出的函数值的线性组合来计算导数。

二、算法原理及推导过程

我们尝试推导k=3时的隐式线性多步法数值计算公式（即 $\beta_{-1}\neq 0$ ），阶数为四阶：

$y_{i+1}=\alpha_{0}y_{i}+\alpha_{1}y_{i-1}+\alpha_{2}y_{i-2}+h(\beta_{-1}f_{i+1}+\beta_{0}f_{i}+\beta_{1}f_{i-1}+\beta_{2}f_{i-2})$

对于上式，其截断误差为：

$LTE=y(x_{i+1})-[\alpha_{0}y(x_{i})+\alpha_{1}y(x_{i-1})+\alpha_{2}y(x_{i-2})+h(\beta_{-1}f(x_{i+1},y(x_{i+1}))+\beta_{0}f(x_{i},y(x_{i}))+\beta_{1}f(x_{i-1},y(x_{i-1}))+\beta_{2}f(x_{i-2},y(x_{i-2})))]$

假设 $y_{i}=y(x_{i})$ ，则 $f_{i}=f(x_{i},y_{i})=f(x_{i},y(x_{i}))=y^{'}(x_{i})$ ，在 $x_{i}$ 处对 $y(x_{i+1}),y(x_{i-1})$ 分别使用泰勒公式，将得到的展开式带入到截断误差计算公式中可得：

$LTE=(1-\alpha_{0}-\alpha_{1}-\alpha_{2})y(x_{i})+\\(1+\alpha_{1}+2\alpha_{2}-\beta_{-1}-\beta_{0}-\beta_{1}-\beta_{2})hy^{'}(x_{i})+\\ (\frac{1}{2}-\frac{1}{2}\alpha_{1}-2\alpha_{2}-\beta_{-1}+\beta_{1}+2\beta_{2})h^{2}y^{''}(x_{i})+\\(\frac{1}{6}+\frac{1}{6}\alpha_{1}+\frac{4}{3}\alpha_{2}-\frac{1}{2}\beta_{-1}-\frac{1}{2}\beta_{1}-2\beta_{2})h^{3}y^{'''}(x_{i})+\\ (\frac{1}{24}-\frac{1}{24}\alpha_{1}-\frac{2}{3}\alpha_{2}-\frac{1}{6}\beta_{-1}+\frac{1}{6}\beta_{1}+\frac{4}{3}\beta_{2})h^{4}y^{(4)}(x_{i})+\\ (\frac{1}{120}+\frac{1}{120}\alpha_{1}+\frac{4}{15}\alpha_{2}-\frac{1}{24}\beta_{-1}-\frac{1}{24}\beta_{1}-\frac{2}{3}\beta_{2})h^{5}y^{5}(x_{i})+O(h^{6})$

为了得到四阶格式，需要满足：

$\left\{\begin{matrix} \alpha_{0}+\alpha_{1}+\alpha_{2}=1,\\ \alpha_{1}+2\alpha_{2}-\beta_{-1}-\beta_{0}-\beta_{1}-\beta_{2}=-1,\\ \frac{1}{2}\alpha_{1}+2\alpha_{2}+\beta_{-1}-\beta_{1}-2\beta_{2}=\frac{1}{2},\\ \frac{1}{6}\alpha_{1}+\frac{4}{3}\alpha_{2}-\frac{1}{2}\beta_{-1}-\frac{1}{2}\beta_{1}-2\beta_{2}=-\frac{1}{6},\\ \frac{1}{24}\alpha_{1}+\frac{2}{3}\alpha_{2}+\frac{1}{6}\beta_{-1}-\frac{1}{6}\beta_{1}-\frac{4}{3}\beta_{2}=\frac{1}{24},\\ \frac{1}{240}\alpha_{1}+\frac{4}{15}\alpha_{2}-\frac{1}{24}\beta_{-1}-\frac{1}{24}\beta_{1}-\frac{2}{3}\beta_{2}\neq -\frac{1}{120},\\ \beta_{-1}\neq0 \end{matrix}\right.$

上式有无穷多解，只要满足上式条件的系数 $\alpha$ 、 $\beta$ 所组成的隐式线性多步格式都是四阶的。

其中，满足 $\alpha$ 和 $\beta$ 系数条件的个别情况下的著名方法有：

1、辛普森（Simpson）方法

取 $\alpha_{0}=0,\alpha_{1}=1,\alpha_{2}=0,\beta_{-1}=\frac{1}{3},\beta_{0}=\frac{4}{3},\beta_{1}=\frac{1}{3},\beta_{2}=0$ ，可得：

$y_{i+1}=y_{i-1}+\frac{1}{3}h(f_{i+1}+4f_{i}+f_{i-1})$

2、汉明（Hamming）方法

取 $\alpha_{0}=\frac{9}{8},\alpha_{1}=0,\alpha_{2}=-\frac{1}{8},\beta_{-1}=\frac{3}{8},\beta_{0}=\frac{3}{4},\beta_{1}=-\frac{3}{8},\beta_{2}=0$ ，可得：

$y_{i+1}=\frac{1}{8}(9y_{i}-y_{i-2})+\frac{3}{8}h(f_{i+1}+2f_{i}-f_{i-1})$

3、隐式四阶阿当姆斯（Adams）方法

取 $\alpha_{0}=1,\alpha_{1}=0,\alpha_{2}=0,\beta_{-1}=\frac{9}{24},\beta_{0}=\frac{19}{24},\beta_{1}=-\frac{5}{24},\beta_{2}=\frac{1}{24}$ ，可得：

$y_{i+1}=y_{i}+\frac{1}{24}h(9f_{i+1}+19f_{i}-5f_{i-1}+f_{i-2})$

在线性多步法中，形如：

$y_{i+1}=y_{i}+h\sum^{k-1}_{j=-1}\beta_{j}f_{i-j}=y_{i}+h(\beta_{-1}f_{i+1}+\beta_{0}f_{i}+\beta_{1}f_{i-1}+\cdot \cdot \cdot +\beta_{k-1}f_{i-k+1})$

的k步法称为阿当姆斯法。当 $\beta_{-1}=0$ 时为阿当姆斯显格式，称为Adams-Bashforth公式；当 $\beta_{-1}\neq0$ 时为阿当姆斯隐格式，称为Adams-Monlton公式。

4、预估-校正阿当姆斯法

与改进欧拉法类似，可以通过显示法预测初值，然后利用隐式法进行迭代。阿当姆斯法也可以使用预估-校正的方式进行计算：

显式阿当姆斯法为：

$y_{i+1}=y_{i}+\frac{1}{24}h(55f_{i}-59f_{i-1}+37f_{i-2}-9f_{i-3})$

隐式阿当姆斯法为：

$y_{i+1}=y_{i}+\frac{1}{24}h(9f_{i+1}+19f_{i}-5f_{i-1}+f_{i-2})$

于是，预估-校正阿当姆斯法为：

①预估: $\bar{y}_{i+1}=y_{i}+\frac{1}{24}h(55f_{i}-59f_{i-1}+37f_{i-2}-9f_{i-3})$

②校正： $y_{i+1}=y_{i}+\frac{1}{24}h(9f(x_{i+1},\bar{y}_{i+1})+19f_{i}-5f_{i-1}+f_{i-2})$

三、算例实现

求解初值问题：

$\left\{\begin{matrix} y^{'}=-y+x+1 &, 0<x\leqslant 1,\\ y(0)=1 & \end{matrix}\right.$

步长h=0.2。已知精确解为：

代码如下：

#include <cmath>
#include <stdlib.h>
#include <stdio.h>


int main(int argc, char* argv[])
{
        int i,N;
        double a,b,h,y0,k1,k2,k3,k4,err,y_predict;
        double *x,*y;
        double exact(double x);
        double f(double x, double y);

        a=0.0;
        b=1.0;
        N=5;
        h=(b-a)/N;

        x=(double*)malloc(sizeof(double)*(N+1));
        y=(double*)malloc(sizeof(double)*(N+1));

        for(i=0;i<=N;i++)
                x[i]=a+i*h;

        y0=1.0;
        y[0]=y0;

        for(i=0;i<3;i++)
        {
                k1=h*f(x[i],y[i]);
                k2=h*f(x[i]+0.5*h,y[i]+0.5*k1);
                k3=h*f(x[i]+0.5*h,y[i]+0.5*k2);
                k4=h*f(x[i]+h,y[i]+k3);
                y[i+1]=y[i]+(k1+2*k2+2*k3+k4)/6.0;
        }

        for(i=3;i<N;i++)
        {
                //显示预测初值
                y_predict=y[i]+h*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])-9*f(x[i-3],y[i-3]))/24.0;
                //隐式求解
                y[i+1]=y[i]+h*(9*f(x[i+1],y_predict)+19*f(x[i],y[i])-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]))/24.0;
        }

        for(i=0;i<=N;i++)
        {
                printf("x[%d]=%.2f,   y=%f,   exact=%f,   err=%.4e.\n",i,x[i],y[i],exact(x[i]),fabs(y[i]-exact(x[i])));

        }
free(x);
free(y);

        return 0;
}

double f(double x, double y)
{
        return -y+x+1;
}
double exact(double x)
{
        return x+exp(-x);
}

N=5时，运行结果如下：

x[0]=0.00,   y=1.000000,   exact=1.000000,   err=0.0000e+00.
x[1]=0.20,   y=1.018733,   exact=1.018731,   err=2.5803e-06.
x[2]=0.40,   y=1.070324,   exact=1.070320,   err=4.2251e-06.
x[3]=0.60,   y=1.148817,   exact=1.148812,   err=5.1888e-06.
x[4]=0.80,   y=1.249323,   exact=1.249329,   err=6.4190e-06.
x[5]=1.00,   y=1.367866,   exact=1.367879,   err=1.3775e-05.

N=10时，运行结果如下：

x[0]=0.00,   y=1.000000,   exact=1.000000,   err=0.0000e+00.
x[1]=0.10,   y=1.004838,   exact=1.004837,   err=8.1964e-08.
x[2]=0.20,   y=1.018731,   exact=1.018731,   err=1.4833e-07.
x[3]=0.30,   y=1.040818,   exact=1.040818,   err=2.0132e-07.
x[4]=0.40,   y=1.070320,   exact=1.070320,   err=1.2779e-07.
x[5]=0.50,   y=1.106530,   exact=1.106531,   err=3.9130e-07.
x[6]=0.60,   y=1.148811,   exact=1.148812,   err=6.0354e-07.
x[7]=0.70,   y=1.196585,   exact=1.196585,   err=7.7242e-07.
x[8]=0.80,   y=1.249328,   exact=1.249329,   err=9.0367e-07.
x[9]=0.90,   y=1.306569,   exact=1.306570,   err=1.0029e-06.
x[10]=1.00,   y=1.367878,   exact=1.367879,   err=1.0751e-06.