常微分方程算法之预估-校正法（改进Euler法）

蒲公英的孩子

已于 2024-04-26 11:08:04 修改

阅读量784

点赞数 8

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

于 2024-04-03 09:44:56 首次发布

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

版权

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

8 篇文章 0 订阅

订阅专栏

一、研究对象

二、预估-校正法原理

三、算例

一、研究对象

主要研究对象：一阶常微分方程的初值问题，与前述方法相同。

$\left\{\begin{matrix} \frac{dy}{dx}=f(x,y) & ,a<x\leqslant b, \space\space\space (1)\\ y(a)=\alpha & \end{matrix}\right.$

通过数值解法来近似获得其数值解，并利用计算机编程实现计算过程。

二、预估-校正法原理

该方法是在欧拉法与梯形法的基础上提出的，梯形法精度高，但计算成本大，每次迭代都需要重新计算 $f\left ( x,y \right )$ 的值，然后又需要进行若干次迭代，计算量大。预估-校正法的原理是通过欧拉法求得一个近似值 $\overline{y}_{i+1}$ ，称为预估值。预估值精度较差，再利用梯形法对其进行一次校正，得到 $y_{i+1}$ ，称为校正值。通过这样的方式构建起来的方法即为“预估-校正法”。相应的公式如下：

$\left\{\begin{matrix} \overline{y}_{i+1}=y_{i} +hf(x_{i},y_{i})& ,\space\space\space i=0,1,2,\cdot \cdot \cdot ,m-1, \space\space(2)\\ y_{i+1}=y_{i}+\frac{1}{2}h[f(x_{i},y_{i})+f(x_{i+1},\overline{y}_{i+1})]& \end{matrix}\right.$

公式（2）中第一式称为预估公式，第二式称为校正公式。公式（2）还可以写作：

$\left\{\begin{matrix} y_{i+1}=y_{i}+\frac{1}{2}k_{1}+\frac{1}{2}k_{2}, & \\ k_{1}=hf(x_{i},y_{i}) & ,\space\space\space i=0,1,2,\cdot \cdot \cdot ,m-1\\ k_{2}=hf(x_{i}+h,y_{i}+k_{1}) & \end{matrix}\right.$

三、算例

与欧拉法、梯形法一样，求相同的求初值问题

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

步长h=0.02。已知精确解为： $y(x)=(1+2x)^{-0.45}$

代码如下：

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

int main(int argc, char* argv[])
{
    int i,N;
    double *x,*y;
    double a,b,h,y0,ypredict,err,maxerr;
    double f(double x, double y);
    double exact(double x);

    a=0.0;
    b=0.1;
    N=5;
    h=(b-a)/N;
    
    x=(double*)malloc(sizeof(double)*(N+1));
    y=(double*)malloc(sizeof(double)*(N+1));
    
    y0=1.0;
    y[0]=y0;
    maxerr=0.0;

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

    for(i=0;i<N;i++)
    {
        ypredict=y[i]+h*f(x[i],y[i]);
        y[i+1]=y[i]+h*(f(x[i],y[i])+f(x[i+1],ypredict))/2.0;
        err=fabs(exact(x[i+1])-y[i+1]);
        printf("x[%d]=%.4f,   y[%d]=%f,   exact=%f,   err=%f.\n",i+1,x[i+1],i+1,y[i+1],exact(x[i+1]),err);
        if(err>maxerr)
            maxerr=err;
    }

    printf("The max error is %f.\n",maxerr);

    free(x);
    free(y);

    return 0;
}


double f(double x, double y)
{
    return -0.9*y/(1.0+2*x);
}
double exact(double x)
{
    return pow((1.0+2*x),-0.45);
}

N=5时，运行结果如下（并与欧拉法、梯形法对比）：

+++++++++++++++++++++++++预估-校正法++++++++++++++++++++++++++++
x[1]=0.0200,   y[1]=0.982324,   exact=0.982506,   err=0.000182.
x[2]=0.0400,   y[2]=0.965616,   exact=0.965960,   err=0.000344.
x[3]=0.0600,   y[3]=0.949791,   exact=0.950281,   err=0.000490.
x[4]=0.0800,   y[4]=0.934772,   exact=0.935393,   err=0.000621.
x[5]=0.1000,   y[5]=0.920492,   exact=0.921231,   err=0.000739.
The max error is 0.000739.

+++++++++++++++++++++++++++梯形法++++++++++++++++++++++++++++++
x[1]=0.0200,   y[1]=0.982498,   exact=0.982506,   err=0.000008.
x[2]=0.0400,   y[2]=0.965946,   exact=0.965960,   err=0.000015.
x[3]=0.0600,   y[3]=0.950260,   exact=0.950281,   err=0.000021.
x[4]=0.0800,   y[4]=0.935367,   exact=0.935393,   err=0.000026.
x[5]=0.1000,   y[5]=0.921201,   exact=0.921231,   err=0.000030.

++++++++++++++++++++++++++++欧拉法+++++++++++++++++++++++++++++
x[1]=0.0200,   y[1]=0.982000,   exact=0.982506,   err=0.000506.
x[2]=0.0400,   y[2]=0.965004,   exact=0.965960,   err=0.000957.
x[3]=0.0600,   y[3]=0.948920,   exact=0.950281,   err=0.001360.
x[4]=0.0800,   y[4]=0.933670,   exact=0.935393,   err=0.001723.
x[5]=0.1000,   y[5]=0.919182,   exact=0.921231,   err=0.002049.

N=10时，运行结果如下（并与欧拉法、梯形法对比）：

+++++++++++++++++++++++++预估-校正法+++++++++++++++++++++++++++
x[1]=0.0100,   y[1]=0.991081,   exact=0.991128,   err=0.000047.
x[2]=0.0200,   y[2]=0.982413,   exact=0.982506,   err=0.000092.
x[3]=0.0300,   y[3]=0.973985,   exact=0.974120,   err=0.000135.
x[4]=0.0400,   y[4]=0.965786,   exact=0.965960,   err=0.000175.
x[5]=0.0500,   y[5]=0.957805,   exact=0.958017,   err=0.000213.
x[6]=0.0600,   y[6]=0.950032,   exact=0.950281,   err=0.000248.
x[7]=0.0700,   y[7]=0.942459,   exact=0.942742,   err=0.000283.
x[8]=0.0800,   y[8]=0.935078,   exact=0.935393,   err=0.000315.
x[9]=0.0900,   y[9]=0.927879,   exact=0.928225,   err=0.000346.
x[10]=0.1000,   y[10]=0.920856,   exact=0.921231,   err=0.000375.
The max error is 0.000375.

++++++++++++++++++++++++++++梯形法+++++++++++++++++++++++++++++
x[1]=0.0100,   y[1]=0.991127,   exact=0.991128,   err=0.000001.
x[2]=0.0200,   y[2]=0.982504,   exact=0.982506,   err=0.000002.
x[3]=0.0300,   y[3]=0.974117,   exact=0.974120,   err=0.000003.
x[4]=0.0400,   y[4]=0.965957,   exact=0.965960,   err=0.000004.
x[5]=0.0500,   y[5]=0.958013,   exact=0.958017,   err=0.000004.
x[6]=0.0600,   y[6]=0.950276,   exact=0.950281,   err=0.000005.
x[7]=0.0700,   y[7]=0.942736,   exact=0.942742,   err=0.000006.
x[8]=0.0800,   y[8]=0.935386,   exact=0.935393,   err=0.000006.
x[9]=0.0900,   y[9]=0.928218,   exact=0.928225,   err=0.000007.
x[10]=0.1000,   y[10]=0.921223,   exact=0.921231,   err=0.000008.

++++++++++++++++++++++++++++欧拉法+++++++++++++++++++++++++++++
x[1]=0.0100,   y[1]=0.991000,   exact=0.991128,   err=0.000128.
x[2]=0.0200,   y[2]=0.982256,   exact=0.982506,   err=0.000250.
x[3]=0.0300,   y[3]=0.973756,   exact=0.974120,   err=0.000364.
x[4]=0.0400,   y[4]=0.965488,   exact=0.965960,   err=0.000473.
x[5]=0.0500,   y[5]=0.957442,   exact=0.958017,   err=0.000575.
x[6]=0.0600,   y[6]=0.949609,   exact=0.950281,   err=0.000672.
x[7]=0.0700,   y[7]=0.941978,   exact=0.942742,   err=0.000764.
x[8]=0.0800,   y[8]=0.934541,   exact=0.935393,   err=0.000851.
x[9]=0.0900,   y[9]=0.927290,   exact=0.928225,   err=0.000934.
x[10]=0.1000,   y[10]=0.920218,   exact=0.921231,   err=0.001013.

可见，与欧拉法相比，预估-校正法具有更小的计算误差，与梯形法相比，又不需要进行迭代，计算量小。所以在实际应用中，预估-校正法更加实用。