求助微分方程Adams算法
那位好心人有微分方程的Adams三步和四步外插法还有四阶预校算法的源程序谢谢了 用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,t;
void Adams_03(float h,float u0,float t0,float T)
{
FILE *fp;
float f0,f1,f2;
u=u0;
t=t0;
fp=fopen("ant.txt","a+");
fprintf(fp,"Adams三步法输出结果为 :\n");
u=u+h*f(t,u);
t=t+h;
u=u+h/2.0*(3.0*f(t,u)-f(t,t));
t=t+h;
for(n=0;n<=(T-t0)/h;n++)
{
t=t+h;
f0=f(t,u);
f1=f(t,u);
f2=f(t,u);
u=u+h/12.0*(23.0*f2-16.0*f1+5.0*f0);
fprintf(fp,"t=%.4f,u[%d]=%f,",t,n,u);
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=u0;
t=t0;
fp=fopen("ant.txt","a+");
fprintf(fp,"Adams四步法输出结果为 :\n");
u=u+h*f(t,u);
u=u+h/2.0*(3.0*f(t,u)-f(t,t));
u=u+h/12.0*(23.0*f(t,u)-16.0*f(t,u)+5.0*f(t,u));
t=t+h;
t=t+h;
t=t+h;
for(n=0;n<=(T-t0)/h;n++)
{
t=t+h;
f0=f(t,u);
f1=f(t,u);
f2=f(t,u);
f3=f(t,u);
u=u+h/24.0*(55.0*f3-59.0*f2+37.0*f1-9.0*f0);
fprintf(fp,"t=%.4f,u[%d]=%f,",t,n,u);
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); }
} 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,double y,double h) //用经典4阶R-K法计算前三个点，以得到较好的初值
{
int i;
double k;
for (i=1;i<4;i++)
{
k=h*f(x,y);
k=h*f(x+h/2,y+k/2);
k=h*f(x+h/2,y+k/2);
k=h*f(x+h,y+k);
y=y+(k+2*k+2*k+k)/6;
x=x+h;
printf("y(%lf)=%.15lf\n",x,y);
}
}
main()
{
int i;
double a,b;
double h;
double c=0;
double p_0=0,p_1;
double x,y,dy;
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=f(x,y);
while(x<b)
{
x=x+h;
p_1=y+h*(55*dy-59*dy+37*dy-9*dy)/24; //预估
c=p_1+251*(c-p_0)/270; //修正
c=f(x,c); //求f
c=y+h*(9*c+19*dy-5*dy+dy)/24; //校正
y=c-19*(c-p_1)/270; //修正
printf("y(%lf)=%.15lf\n",x,y); //输出y
//为利用前一次计算结果并节省存储空间，赋值
dy=dy;
dy=dy;
dy=dy;
dy=f(x,y);
p_0=p_1;
}
end=clock(); //停止计时
t=end-begin;
printf("运行时间为%.15f.\n",t);
}
万分感谢
页:
[1]