目录
一、研究问题
二、C++代码
三、计算结果
一、研究问题
本节我们采用阿当姆斯法(Adams法)求解算例。
阿当姆斯法的原理及推导请参考:
常微分方程算法之阿当姆斯法(Adams法)_四步四阶adams显格式;三步四阶adams隐格式;四阶adams预估-校正格式-CSDN博客https://blog.csdn.net/L_peanut/article/details/137378532 研究问题依然为
取步长为0.1。已知精确解为。
二、C++代码
这里我们使用四阶阿当姆斯显格式和四阶阿当姆斯预估-校正格式进行求解,如下:
#include <cmath>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char *argv[])
{
int i,N;
double a,b,h,k1,k2,k3,k4,upredict;
double *x,*y,*u;
double f(double x, double y);
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));
u=(double *)malloc(sizeof(double)*(N+1));
for(i=0;i<=N;i++)
x[i]=a+i*h;
y[0]=1.0;
u[0]=1.0;
for(i=0;i<=2;i++) //用龙格-库塔法求解y1,y2,y3作为初始值
{
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;
u[i+1]=y[i+1];
}
for(i=3;i<N;i++)
{
//四阶阿当姆斯显格式
y[i+1]=y[i]+(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]))*h/24.0;
//四阶阿当姆斯预估-校正格式
upredict=u[i]+(55*f(x[i],u[i])-59*f(x[i-1],u[i-1])+37*f(x[i-2],u[i-2])-9*f(x[i-3],u[i-3]))*h/24.0;
u[i+1]=u[i]+(9*f(x[i+1],upredict)+19*f(x[i],u[i])-5*f(x[i-1],u[i-1])+f(x[i-2],u[i-2]))*h/24.0;
}
for(i=0;i<=N;i++)
printf("x=%.2f, 显示解=%f, 预估-校正解=%f, exact=%f.\n",x[i],y[i],u[i],exact(x[i]));
free(x);free(y);free(u);
return 0;
}
double f(double x, double y)
{
return y-2*x/y;
}
double exact(double x)
{
return sqrt(1+2.0*x);
}
三、计算结果
x=0.00, 显示解=1.000000, 预估-校正解=1.000000, exact=1.000000.
x=0.10, 显示解=1.095446, 预估-校正解=1.095446, exact=1.095445.
x=0.20, 显示解=1.183217, 预估-校正解=1.183217, exact=1.183216.
x=0.30, 显示解=1.264912, 预估-校正解=1.264912, exact=1.264911.
x=0.40, 显示解=1.341552, 预估-校正解=1.341641, exact=1.341641.
x=0.50, 显示解=1.414046, 预估-校正解=1.414214, exact=1.414214.
x=0.60, 显示解=1.483019, 预估-校正解=1.483240, exact=1.483240.
x=0.70, 显示解=1.548919, 预估-校正解=1.549193, exact=1.549193.
x=0.80, 显示解=1.612116, 预估-校正解=1.612452, exact=1.612452.
x=0.90, 显示解=1.672917, 预估-校正解=1.673320, exact=1.673320.
x=1.00, 显示解=1.731570, 预估-校正解=1.732051, exact=1.732051.
从计算结果可知,预估-校正法求解结果有效数字更多,结果更为精确,优于显示法。