Adams预测-校正系统

2012-2013（1）专业课程实践论文

Adams 预测-校正系统

高强强，0818180212，R数学08-2班

王岳，0818180108，R数学08-1班

一、算法理论

显示格式和隐式格式都具有四阶精度，两种格式可匹配成下列Adams 预测-校正

系统：

Adams 预测 - 校正

预测 ()3-2-1-n 19-3759f -5524

n n n n n f f f h y y ++=+, ()111,+++=n n n y x f f ;

校正 ()

2-1-115-19924n n n n n n f f f f h y y +++=++ ; ()111,+++=n n n y x f f .

Adams 内插外插公式联合使用称作 Adams 预测 - 校正系统，用外插公式计算

预测，用内插公式进行校正。计算时需要注意以下两点：

1 、外插公式为显式，内插公式为隐式。故用内插公式计算需要进行迭代。

2 、这种预测 - 校正法是四步法，计算1+n y 时，不但用到前一步的信息n x ，n y ，

而且要用到更前三步的信息 1-n f ,2-n f ,3-n f ，因此它不是自开始的，实际计算时必

须借助于某种单步法，用Runge-Kutta 格式为它提供开始值 1y ,2y ,3y 。

依据局部截断误差公式得

()()()n n n x y h y x y 5511720

251-≈++, ()()()n n n x y h y x y 5511720

19-≈

++. 进一步将Adams 预测-校正系统加工成下列方案：

预测 ()3-2-1-19-3759-5524

n n n n n n f f f f h y p ++=+, 改进 ()n n n n p c p m -72025111+=++, 计算 ()111,+++=n n n m x f m ,

校正 ()2-1-115-19924n n n n n n f f f m h y c ++'+

=++, 改进 ()1111

-72019-++++=n n n n p c c y , 计算 ()111,+++='n n n

y x f y ; 运用上述方案计算1+n y 时，要用到先一步的信息n y ，n

y '，n n c p -和更前一步的

信息1-n y ，因此在计算机之前必须给出初值1y 和11-c p ，1y 可用其他单步法来计

算，11-c p 则一般令它等于零。

开始

输入区间（a,b),步长h,初始值

x[0],y[0]

调用RK函数求出y[1],y[2,y[3]]

计算dy[1],dy[2],dy[3]

开始迭代

p_0=0

x[3]

y

预估值p_1,修正值c,并打印预估值，调用函数f求出dy[i]校正值c，求出修正值y[i],输出y[i],p_0=p_1

x[3]=x[3]+h

j结束

#include

#include

double f(double x,double y)

{

double z;

z=-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x);

return(z);

}

void RK(double x[4],double y[4],double h) //用经典4阶R-K法计算前三个点，以得到较好的初值

{

int i;

double k[4];

for (i=1;i<4;i++)

{

k[0]=h*f(x[i-1],y[i-1]);

k[1]=h*f(x[i-1]+h/2,y[i-1]+k[0]/2);

k[2]=h*f(x[i-1]+h/2,y[i-1]+k[1]/2);

k[3]=h*f(x[i-1]+h,y[i-1]+k[2]);

y[i]=y[i-1]+(k[0]+2*k[1]+2*k[2]+k[3])/6;

x[i]=x[i-1]+h;

printf("y(%lf)=%.15lf\n",x[i],y[i]);

}

}

main()

{

int i;

double a,b;

double h;

double c=0;

double p_0=0,p_1;

double x[4],y[4],dy[100];

double t;

printf("请输入区间（a,b）:\n");

scanf("%lf,%lf",&a,&b);

printf("请输入步长h:\n");

scanf("%lf",&h);

printf("请输入初值（x_0,y_0）:\n");

scanf("%lf,%lf",&x,&y); //计时开始

RK(x,y,h); //使用R-K法开始，计算前三个点

for (i=0;i<4;i++)

dy[i]=f(x[i],y[i]);

while(x[3]

{

x[3]=x[3]+h;

p_1=y[3]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24; //预估

c=p_1+251*(c-p_0)/270; //修正printf("预测值y(%lf)=%.15lf\n",x[3],c);

c=f(x[3],c); //求f

c=y[3]+h*(9*c+19*dy[3]-5*dy[2]+dy[1])/24; //校正y[3]=c-19*(c-p_1)/270; //修正printf("校正值y(%lf)=%.15lf\n",x[3],y[3]); //输出y //为利用前一次计算结果并节省存储空间，赋值

dy[0]=dy[1];

dy[1]=dy[2];

dy[2]=dy[3];

dy[3]=f(x[3],y[3]);

p_0=p_1;

}

}