目录
一、研究问题
二、C++代码
三、计算结果
一、研究问题
本节我们采用预估校正法(改进欧拉法)求解算例。
预估-校正法的原理及推导请参考:
常微分方程算法之预估-校正法(改进Euler法)_、改进欧拉法-CSDN博客https://blog.csdn.net/L_peanut/article/details/137326300?spm=1001.2014.3001.5501
研究问题依然为
取步长为0.1。
二、C++代码
#include <cmath>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char *argv[])
{
int i,k,N;
double a,b,h,y0,ypredict,err,maxerr;
double *x,*y;
double f(double x, double t);
double exact(double x);
a=0.0; //求解区域左端点
b=1.0; //求解区域右端点
N=10; //总的剖分数
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<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(y[i+1]-exact(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],exact(x[i+1]),err);
if(err>maxerr)
maxerr=err;
}
printf("The max error is %f.\n",maxerr);
return 0;
}
//右端项函数
double f(double x, double y)
{
return y-2*x/y;
}
//精确解
double exact(double x)
{
return sqrt(1.0+2*x);
}
三、计算结果
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.
同样的算例,采用梯形法和欧拉法的计算结果误差为:
+++++++++++++++++++++++++++梯形法+++++++++++++++++++++++++++++
k=3, x[1]=0.1000, y[1]=1.095657, exact=1.095445, err=0.000212.
k=3, x[2]=0.2000, y[2]=1.183596, exact=1.183216, err=0.000380.
k=3, x[3]=0.3000, y[3]=1.265444, exact=1.264911, err=0.000532.
k=3, x[4]=0.4000, y[4]=1.342327, exact=1.341641, err=0.000686.
k=3, x[5]=0.5000, y[5]=1.415064, exact=1.414214, err=0.000850.
k=3, x[6]=0.6000, y[6]=1.484274, exact=1.483240, err=0.001034.
k=3, x[7]=0.7000, y[7]=1.550437, exact=1.549193, err=0.001244.
k=2, x[8]=0.8000, y[8]=1.613948, exact=1.612452, err=0.001496.
k=2, x[9]=0.9000, y[9]=1.675112, exact=1.673320, err=0.001792.
k=2, x[10]=1.0000, y[10]=1.734192, exact=1.732051, err=0.002141.
The max error is 0.002141.
++++++++++++++++++++++++++欧拉法+++++++++++++++++++++++++
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.
可见虽然与梯形法相比,误差较大,但是可以在不增加计算量的前提条件下获得比欧拉法更小的误差。如果取步长为0.05,计算结果为:
x[1]=0.0500, y[1]=1.048869, exact=1.048809, err=0.000060.
x[2]=0.1000, y[2]=1.095561, exact=1.095445, err=0.000116.
x[3]=0.1500, y[3]=1.140345, exact=1.140175, err=0.000169.
x[4]=0.2000, y[4]=1.183437, exact=1.183216, err=0.000221.
x[5]=0.2500, y[5]=1.225017, exact=1.224745, err=0.000273.
x[6]=0.3000, y[6]=1.265236, exact=1.264911, err=0.000325.
x[7]=0.3500, y[7]=1.304219, exact=1.303840, err=0.000378.
x[8]=0.4000, y[8]=1.342075, exact=1.341641, err=0.000434.
x[9]=0.4500, y[9]=1.378897, exact=1.378405, err=0.000492.
x[10]=0.5000, y[10]=1.414767, exact=1.414214, err=0.000553.
x[11]=0.5500, y[11]=1.449756, exact=1.449138, err=0.000618.
x[12]=0.6000, y[12]=1.483927, exact=1.483240, err=0.000687.
x[13]=0.6500, y[13]=1.517337, exact=1.516575, err=0.000762.
x[14]=0.7000, y[14]=1.550035, exact=1.549193, err=0.000842.
x[15]=0.7500, y[15]=1.582067, exact=1.581139, err=0.000928.
x[16]=0.8000, y[16]=1.613472, exact=1.612452, err=0.001021.
x[17]=0.8500, y[17]=1.644289, exact=1.643168, err=0.001122.
x[18]=0.9000, y[18]=1.674551, exact=1.673320, err=0.001231.
x[19]=0.9500, y[19]=1.704288, exact=1.702939, err=0.001350.
x[20]=1.0000, y[20]=1.733530, exact=1.732051, err=0.001479.
The max error is 0.001479.
可见缩小步长,可以取得较理想的计算结果。