常微分方程算法之梯形法（隐式Eluer法）

蒲公英的孩子

已于 2024-04-02 16:49:59 修改

阅读量979

点赞数 3

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

于 2024-04-02 16:48:45 首次发布

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

版权

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

8 篇文章 3 订阅

订阅专栏

一、研究对象

二、梯形法原理

三、算例

一、研究对象

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

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

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

二、梯形法原理

对公式（1）两端从 $x_{i}$ 到 $x_{i+1}$ 积分，可得：

$y(x_{i+1})-y(x_{i})=\int_{x_{i}}^{x_{i+1}}f(x,y(x))dx \space\space\space\space (2)$

对于欧拉法而言，是采用左矩形公式进行逼近，即：

$\int_{x_{i}}^{x_{i+1}}f(x,y(x))dx\approx f(x_{i},y(x_{i}))h$

其中，用 $y_{i}$ 作为 $y(x_{i})$ 的近似值，得到的便是欧拉法的差分格式： $y_{i+1}-y_{i}=hf(x_{i},y_{i})$ 。

如果采用右矩形公式进行逼近，可得类似的差分格式： $y_{i+1}-y_{i}=hf(x_{i+1},y_{i+1})$ 。

采用梯形公式计算（2）右端积分，可得：

$\int_{x_{i}}^{x_{i+1}}f(x,y(x))dx \approx \frac{1}{2}h[f(x_{i},y(x_{i})+f(x_{i+1},y_(x_{i+1})]$

用 $y_{i}$ 作为 $y(x_{i})$ 的近似值，可得梯形法的差分格式：

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

梯形法属于隐式格式，需要通过迭代的方式来求解，迭代初值由欧拉法给出。于是：

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

三、算例

求初值问题

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

步长h=0.02，迭代误差限制为 $\varepsilon = 10^{-6}$ 。已知精确解为： $y(x)=(1+2x)^{-0.45}$

代码如下：

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

int main(int argc, char* argv[])
{
        int i,N,k;
        double a,b,h,y0,err,ytemp1,ytemp2,epsilon;
        double *x,*y;
        double f(double x, double y);
        double yexact(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));

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

        y0=1.0;
        y[0]=y0;
        epsilon=1e-6;

        for(i=0;i<N;i++)
        {
                ytemp1=y[i]+h*f(x[i],y[i]);
                k=0;
                do
                {
                        k=k+1;
                        if(k!=1)
                                ytemp1=ytemp2;
                        ytemp2=y[i]+h*(f(x[i],y[i])+f(x[i+1],ytemp1))/2.0;
                }while(fabs(ytemp1-ytemp2)>epsilon);
                y[i+1]=ytemp2;
                err=fabs(y[i+1]-yexact(x[i+1]));
                printf("x[%d]=%.4f,   y[%d]=%f,   exact=%F,   err=%f.\n",i+1,x[i+1],i+1,y[i+1],yexact(x[i+1]),err);
        }

        free(x);
        free(y);

        return 0;
}

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

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

+++++++++++++++++++++++++++梯形法++++++++++++++++++++++++++++++
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.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.

可见，相同求解域内，剖分点数越多，计算误差越小，但所需的计算量越大，耗时越多。

同时与欧拉法对比可知，梯形法计算误差更小，精度更高。