4阶adams预估-校正算法.c

4阶Adams预估-校正算法

#include<stdio.h>
#include<math.h>
#include<time.h>
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];
 clock_t begin,end;
 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);
 begin=clock();                         //计时开始
 RK(x,y,h);                            //使用R-K法开始，计算前三个点
 for (i=0;i<4;i++)
	dy[i]=f(x[i],y[i]);
 while(x[3]<b)
 {
  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;                                     //修正
  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;
 }
 end=clock();                                           //停止计时
 t=end-begin;
 printf("运行时间为%.15f.\n",t);
}