求助微分方程Adams算法

hayrhen

那位好心人有微分方程的Adams三步和四步外插法还有四阶预校算法的源程序
谢谢了

风花雪月

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

1. #include "stdio.h"
   
2. #include "math.h"
   
3. #define f(t,u) -8.0*(u)+4.0*(t)*(t)-7.0*(t)-1.0
   
4. #define N 1000
   
5. int i,n;
   
6. float u[N],t[N];
   
7. 
   
8. void Adams_03(float h,float u0,float t0,float T)
   
9. {
   
10.     FILE *fp;
    
11.     float f0,f1,f2;
    
12.     u[0]=u0;
    
13.     t[0]=t0;
    
14.     fp=fopen("ant.txt","a+");
    
15.      fprintf(fp,"Adams三步法输出结果为 :\n");
    
16.     u[1]=u[0]+h*f(t[0],u[0]);
    
17.     t[1]=t[0]+h;
    
18.     u[2]=u[1]+h/2.0*(3.0*f(t[1],u[1])-f(t[0],t[0]));
    
19.     t[2]=t[1]+h;
    
20.     for(n=0;n<=(T-t0)/h;n++)
    
21.      {
    
22.        t[n+3]=t[n+2]+h;
    
23.        f0=f(t[n],u[n]);
    
24.        f1=f(t[n+1],u[n+1]);
    
25.        f2=f(t[n+2],u[n+2]);
    
26.        u[n+3]=u[n+2]+h/12.0*(23.0*f2-16.0*f1+5.0*f0);
    
27.        fprintf(fp,"t=%.4f,u[%d]=%f,",t[n],n,u[n]);
    
28.        if(n%2==0)
    
29.        fprintf(fp,"\n");
    
30. }
    
31.     fprintf(fp,"\n");  
    
32.     fclose(fp);
    
33. }
    
34. 
    
35. void Adams_04(float h,float u0,float t0,float T)
    
36. {
    
37.     FILE *fp;
    
38.     float f0,f1,f2,f3;
    
39.     u[0]=u0;
    
40.     t[0]=t0;
    
41.     fp=fopen("ant.txt","a+");
    
42.      fprintf(fp,"Adams四步法输出结果为 :\n");
    
43.     u[1]=u[0]+h*f(t[0],u[0]);
    
44.     u[2]=u[1]+h/2.0*(3.0*f(t[1],u[1])-f(t[0],t[0]));
    
45.     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]));
    
46.     t[1]=t[0]+h;
    
47.     t[2]=t[1]+h;
    
48.     t[3]=t[2]+h;
    
49.     for(n=0;n<=(T-t0)/h;n++)
    
50.      {
    
51.        t[n+4]=t[n+3]+h;
    
52.        f0=f(t[n],u[n]);
    
53.        f1=f(t[n+1],u[n+1]);
    
54.        f2=f(t[n+2],u[n+2]);
    
55.        f3=f(t[n+3],u[n+3]);
    
56.        u[n+4]=u[n+3]+h/24.0*(55.0*f3-59.0*f2+37.0*f1-9.0*f0);
    
57.        fprintf(fp,"t=%.4f,u[%d]=%f,",t[n],n,u[n]);
    
58.        if(n%2==0)
    
59.        fprintf(fp,"\n");
    
60. }
    
61.     fprintf(fp,"\n");  
    
62.     fclose(fp);
    
63. }
    
64. main()
    
65. {
    
66.     FILE *fp;
    
67.    int k;
    
68.     float h,u0,t0,T;
    
69.     fp=fopen("ant.txt","a+");
    
70.     printf("h,u0,t0,T=");
    
71.     scanf("%f,%f,%f,%f",&h,&u0,&t0,&T);
    
72.     printf("使用Adams K步法:");
    
73.     scanf("%d",&k);
    
74.     fprintf(fp,"用Adams%d步法求解的结果为:\n",k);
    
75.     fprintf(fp,"初值：u(0)=%f\nt的取值范围:%f<t<=%f\n",u0,t0,T);
    
76.     fprintf(fp,"h=%f\n",h);
    
77.     fclose(fp);
    
78.     if(k==3)
    
79.        {  Adams_03(h,u0,t0,T);   }
    
80.     if(k==4)
    
81.        { Adams_04(h,u0,t0,T);    }
    
82.     }

复制代码

风花雪月

4阶Adams预估-校正算法

1. #include<stdio.h>
   
2. #include<math.h>
   
3. #include<time.h>
   
4. double f(double x,double y)
   
5. {
   
6. double z;
   
7. z=-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x);
   
8. return(z);
   
9. }
   
10. void RK(double x[4],double y[4],double h) //用经典4阶R-K法计算前三个点，以得到较好的初值
    
11. {
    
12. int i;
    
13. double k[4];
    
14. for (i=1;i<4;i++)
    
15. {
    
16.   k[0]=h*f(x[i-1],y[i-1]);
    
17.   k[1]=h*f(x[i-1]+h/2,y[i-1]+k[0]/2);
    
18.   k[2]=h*f(x[i-1]+h/2,y[i-1]+k[1]/2);
    
19.   k[3]=h*f(x[i-1]+h,y[i-1]+k[2]);
    
20.   y[i]=y[i-1]+(k[0]+2*k[1]+2*k[2]+k[3])/6;
    
21.   x[i]=x[i-1]+h;
    
22.   printf("y(%lf)=%.15lf\n",x[i],y[i]);
    
23. }
    
24. }
    
25. main()
    
26. {
    
27. int i;
    
28. double a,b;
    
29. double h;
    
30. double c=0;
    
31. double p_0=0,p_1;
    
32. double x[4],y[4],dy[100];
    
33. clock_t begin,end;
    
34. double t;
    
35. printf("请输入求解区间（a,b）:\n");
    
36. scanf("%lf,%lf",&a,&b);
    
37. printf("请输入步长h:\n");
    
38. scanf("%lf",&h);
    
39. printf("请输入初值（x_0,y_0）:\n");
    
40. scanf("%lf,%lf",&x,&y);
    
41. begin=clock();                         //计时开始
    
42. RK(x,y,h);                            //使用R-K法开始，计算前三个点
    
43. for (i=0;i<4;i++)
    
44.     dy[i]=f(x[i],y[i]);
    
45. while(x[3]<b)
    
46. {
    
47.   x[3]=x[3]+h;
    
48.   p_1=y[3]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24;        //预估
    
49.   c=p_1+251*(c-p_0)/270;                                     //修正
    
50.   c=f(x[3],c);                                               //求f
    
51.   c=y[3]+h*(9*c+19*dy[3]-5*dy[2]+dy[1])/24;                  //校正
    
52.   y[3]=c-19*(c-p_1)/270;                                     //修正
    
53.   printf("y(%lf)=%.15lf\n",x[3],y[3]);                       //输出y
    
54.   //为利用前一次计算结果并节省存储空间，赋值
    
55.   dy[0]=dy[1];
    
56.   dy[1]=dy[2];
    
57.   dy[2]=dy[3];
    
58.   dy[3]=f(x[3],y[3]);
    
59.   p_0=p_1;
    
60. }
    
61. end=clock();                                           //停止计时
    
62. t=end-begin;
    
63. printf("运行时间为%.15f.\n",t);
    
64. }

复制代码

hayrhen