用adams三步四步法求解微分方程.c

用Adams三步四步法求解微分方程

#include "stdio.h"
#include "math.h"
#define f(t,u) -8.0*(u)+4.0*(t)*(t)-7.0*(t)-1.0
#define N 1000
int i,n;
float u[N],t[N];

void Adams_03(float h,float u0,float t0,float T)
{
    FILE *fp;
    float f0,f1,f2;
    u[0]=u0;
    t[0]=t0;
    fp=fopen("ant.txt","a+");
     fprintf(fp,"Adams三步法输出结果为 :\n");
    u[1]=u[0]+h*f(t[0],u[0]);
    t[1]=t[0]+h;
    u[2]=u[1]+h/2.0*(3.0*f(t[1],u[1])-f(t[0],t[0]));
    t[2]=t[1]+h;
    for(n=0;n<=(T-t0)/h;n++)
     {
       t[n+3]=t[n+2]+h;
       f0=f(t[n],u[n]);
       f1=f(t[n+1],u[n+1]);
       f2=f(t[n+2],u[n+2]);
       u[n+3]=u[n+2]+h/12.0*(23.0*f2-16.0*f1+5.0*f0);
       fprintf(fp,"t=%.4f,u[%d]=%f,",t[n],n,u[n]);
       if(n%2==0)
       fprintf(fp,"\n");
}
    fprintf(fp,"\n");  
    fclose(fp);
}

void Adams_04(float h,float u0,float t0,float T)
{
    FILE *fp;
    float f0,f1,f2,f3;
    u[0]=u0;
    t[0]=t0;
    fp=fopen("ant.txt","a+");
     fprintf(fp,"Adams四步法输出结果为 :\n");
    u[1]=u[0]+h*f(t[0],u[0]);
    u[2]=u[1]+h/2.0*(3.0*f(t[1],u[1])-f(t[0],t[0]));
    u[3]=u[2]+h/12.0*(23.0*f(t[2],u[2])-16.0*f(t[1],u[1])+5.0*f(t[0],u[0]));
    t[1]=t[0]+h;
    t[2]=t[1]+h;
    t[3]=t[2]+h;
    for(n=0;n<=(T-t0)/h;n++)
     {
       t[n+4]=t[n+3]+h;
       f0=f(t[n],u[n]);
       f1=f(t[n+1],u[n+1]);
       f2=f(t[n+2],u[n+2]);
       f3=f(t[n+3],u[n+3]);
       u[n+4]=u[n+3]+h/24.0*(55.0*f3-59.0*f2+37.0*f1-9.0*f0);
       fprintf(fp,"t=%.4f,u[%d]=%f,",t[n],n,u[n]);
       if(n%2==0)
       fprintf(fp,"\n");
}
    fprintf(fp,"\n");  
    fclose(fp);
}
main()
{
    FILE *fp;
   int k;
    float h,u0,t0,T;
    fp=fopen("ant.txt","a+");
    printf("h,u0,t0,T=");
    scanf("%f,%f,%f,%f",&h,&u0,&t0,&T);
    printf("使用Adams K步法:");
    scanf("%d",&k);
    fprintf(fp,"用Adams%d步法求解的结果为:\n",k);
    fprintf(fp,"初值：u(0)=%f\nt的取值范围:%f<t<=%f\n",u0,t0,T);
    fprintf(fp,"h=%f\n",h);
    fclose(fp);
    if(k==3)
       {  Adams_03(h,u0,t0,T);   }
    if(k==4)
       { Adams_04(h,u0,t0,T);    }
    }