常微分方程算法之编程示例四（龙格-库塔法）

news2024/11/15 23:48:16

一、算例一

1.1 研究问题

1.2 C++代码

1.3 计算结果

二、算例二

2.1 研究问题

2.2 C++代码

2.3 计算结果

一、算例一

本节我们采用龙格-库塔法（Runge-Kutta法）求解算例。

龙格-库塔法的原理及推导请参考：

常微分方程算法之龙格-库塔法（Runge-Kutta法）_龙格库塔法-CSDN博客https://blog.csdn.net/L_peanut/article/details/137336203?spm=1001.2014.3001.5501

1.1 研究问题

研究问题依然为

$\left\{\begin{matrix} \frac{dy}{dx}=y-\frac{2x}{y},0<x\leqslant 1,\\ y(0)=1 \end{matrix}\right.$

取步长为0.1。已知精确解为 $y=\sqrt{1+2x}$ 。

1.2 C++代码

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


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

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

        x0=a;
        y0=1.0;
        i=1;
        do
        {
                x1=x0+h;
                k1=h*f(x0,y0);
                k2=h*f(x0+0.5*h,y0+0.5*k1);
                k3=h*f(x0+0.5*h,y0+0.5*k2);
                k4=h*f(x1,y0+k3);
                y1=y0+(k1+2*k2+2*k3+k4)/6.0;

                yexact=exact(x1);
                err=fabs(yexact-y1);
                printf("x=%.2f, y=%f, exact=%f, error=%f.\n",x1,y1,yexact,err);
                
                i++;
                x0=x1;
                y0=y1;
        }while(i<=N);

        return 0;
}



double f(double x, double y)
{
        return y-2*x/y;
}
double exact(double x)
{
        return sqrt(1+2.0*x);
}

1.3 计算结果

x=0.10, y=1.095446, exact=1.095445, error=0.000000.
x=0.20, y=1.183217, exact=1.183216, error=0.000001.
x=0.30, y=1.264912, exact=1.264911, error=0.000001.
x=0.40, y=1.341642, exact=1.341641, error=0.000002.
x=0.50, y=1.414216, exact=1.414214, error=0.000002.
x=0.60, y=1.483242, exact=1.483240, error=0.000003.
x=0.70, y=1.549196, exact=1.549193, error=0.000003.
x=0.80, y=1.612455, exact=1.612452, error=0.000004.
x=0.90, y=1.673325, exact=1.673320, error=0.000005.
x=1.00, y=1.732056, exact=1.732051, error=0.000006.

与相同条件下的欧拉法、梯形法、预估-校正法计算结果相比：

++++++++++++++++++++++++++欧拉法+++++++++++++++++++++++++
x[1]=0.1000, y[1]=1.100000, exact=1.095445, err=0.004555.
x[2]=0.2000, y[2]=1.191818, exact=1.183216, err=0.008602.
x[3]=0.3000, y[3]=1.277438, exact=1.264911, err=0.012527.
x[4]=0.4000, y[4]=1.358213, exact=1.341641, err=0.016572.
x[5]=0.5000, y[5]=1.435133, exact=1.414214, err=0.020919.
x[6]=0.6000, y[6]=1.508966, exact=1.483240, err=0.025727.
x[7]=0.7000, y[7]=1.580338, exact=1.549193, err=0.031145.
x[8]=0.8000, y[8]=1.649783, exact=1.612452, err=0.037332.
x[9]=0.9000, y[9]=1.717779, exact=1.673320, err=0.044459.
x[10]=1.0000, y[10]=1.784771, exact=1.732051, err=0.052720.

+++++++++++++++++++++++++梯形法++++++++++++++++++++++++++
x[1]=0.1000, y[1]=1.095657, exact=1.095445, err=0.000212.
x[2]=0.2000, y[2]=1.183596, exact=1.183216, err=0.000380.
x[3]=0.3000, y[3]=1.265444, exact=1.264911, err=0.000532.
x[4]=0.4000, y[4]=1.342327, exact=1.341641, err=0.000686.
x[5]=0.5000, y[5]=1.415064, exact=1.414214, err=0.000850.
x[6]=0.6000, y[6]=1.484274, exact=1.483240, err=0.001034.
x[7]=0.7000, y[7]=1.550437, exact=1.549193, err=0.001244.
x[8]=0.8000, y[8]=1.613948, exact=1.612452, err=0.001496.
x[9]=0.9000, y[9]=1.675112, exact=1.673320, err=0.001792.
x[10]=1.0000, y[10]=1.734192, exact=1.732051, err=0.002141.

+++++++++++++++++++++++预估-校正法+++++++++++++++++++++++
x[1]=0.1000, y[1]=1.095909, exact=1.095445, err=0.000464.
x[2]=0.2000, y[2]=1.184097, exact=1.183216, err=0.000881.
x[3]=0.3000, y[3]=1.266201, exact=1.264911, err=0.001290.
x[4]=0.4000, y[4]=1.343360, exact=1.341641, err=0.001719.
x[5]=0.5000, y[5]=1.416402, exact=1.414214, err=0.002188.
x[6]=0.6000, y[6]=1.485956, exact=1.483240, err=0.002716.
x[7]=0.7000, y[7]=1.552514, exact=1.549193, err=0.003321.
x[8]=0.8000, y[8]=1.616475, exact=1.612452, err=0.004023.
x[9]=0.9000, y[9]=1.678166, exact=1.673320, err=0.004846.
x[10]=1.0000, y[10]=1.737867, exact=1.732051, err=0.005817.

相比之下，四阶龙格-库塔法计算结果精度最高，但是计算量较大。

二、算例二

2.1 研究问题

研究问题为

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

取步长为0.25。已知精确解为 $y=x^{2}-2x+2-e^{-x}$ 。

2.2 C++代码

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


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

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

        x0=a;
        y0=1.0;
        i=1;
        do
        {
                x1=x0+h;
                k1=h*f(x0,y0);
                k2=h*f(x0+0.5*h,y0+0.5*k1);
                k3=h*f(x0+0.5*h,y0+0.5*k2);
                k4=h*f(x1,y0+k3);
                y1=y0+(k1+2*k2+2*k3+k4)/6.0;

                yexact=exact(x1);
                err=fabs(yexact-y1);
                printf("x=%.2f, y=%.8f, exact=%f, error=%.4e.\n",x1,y1,yexact,err);
                
                i++;
                x0=x1;
                y0=y1;
        }while(i<=N);

        return 0;
}



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

2.3 计算结果

x=0.25, y=0.78371175, exact=0.783699, error=1.2534e-05.
x=0.50, y=0.64349336, exact=0.643469, error=2.4024e-05.
x=0.75, y=0.59016776, exact=0.590133, error=3.4318e-05.
x=1.00, y=0.63216394, exact=0.632121, error=4.3382e-05.

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/1861045.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！