微分4.c

使用Adams三次外插方法求解微分方程的数值解

 #include<stdio.h>
#include<math.h>
double zhenzhi(double t)
{/*函数真值计算公式*/
 double u;
 u=exp(-8*t)+0.5*t*t-t;
return u;
} 
double f(double u,double t)
{/* f(t,u)的公式*/
 double y;
 y=-8*u+4*t*t-7*t-1;
 return y;
}
double euler(double u,double h,double t)
{/*欧拉法用于计算多步方法的迭代的初始值*/
 u=u+h*f(u,t);
 return u;
}
double wadmas(double u,double h,double f0,double f1,double f2,double f3)
{/*四步外插公式*/
 u=u+(h/24)*((-9)*f0+37*f1-59*f2+55*f3);
 return u;
}
double nadmas(double u,double h,double f0,double f1,double f2 ,double f3)
{/*三步内插公式*/
 u=u+(h/24)*(f0+(-5)*f1+(19)*f2+9*f3);
 return u;
}
 main()
{
double u,h,t,T,a[4],f0,f1,f2,f3,up,uc,temp1,temp2,s;
int i;
printf("please input chizhi u0,and bu chang h\n");
scanf("%lf%lf",&u,&h);
printf("please input (t0,T)\n");
scanf("%lf%lf",&t,&T);
s=zhenzhi(t);
printf("................................................................\n");
printf("近似值（校正算法）：           真值：          相对误差:\n");
printf("u(%lf)=%lf\t u(%lf)=%lf\t%f \n",t,u,t,s,fabs((u-s)/u));
for(i=0;i<4;i++)
{
 a[i]=f(u,t);
 u=euler(u,h,t);
 up=u;
 uc=u;
 t=t+h;
 s=zhenzhi(t); /* s 保留函数的真值*/
 printf("u(%lf)=%lf\t u(%lf)=%lf\t%f \n",t,u,t,s,fabs((u-s)/u));
}
f0=a[0];
f1=a[1];
f2=a[2];
f3=a[3];
printf("\n");
for(i=3;i<=46;i++)
{
 temp1=up;
 temp2=uc;
 up=wadmas(temp1,h,f0,f1,f2,f3);
 up=up+251/270*(temp2-temp1);/*m1*/
 f3=f(up,t);
 uc=nadmas(temp2,h,f0,f1,f2,f3);
 u=uc-19/270*(uc- up);/*m2*/
 f0=f1;
 f1=f2;
 f2=f3;
 t=t+h;
 s=zhenzhi(t);
 printf("u(%lf)=%lf\t u(%lf)=%lf\t%f \n",t,u,t,s,fabs((u-s)/u));
}
}