求高手找一下错!!!
这是一个四阶RK法的程序,本人做了一些改变,可不知哪里出了错误,计算的结果有问题!#include"stdlib.h"
#include"math.h"
void rkt3(double t,double h,double y[],int n,void(*f)())
{int i,j,k;
double tt,a[4],*g,*b,*c,*d,*e;
g=malloc(n*sizeof(double));
b=malloc(n*sizeof(double));
c=malloc(n*sizeof(double));
d=malloc(n*sizeof(double));
e=malloc(n*sizeof(double));
a[0]=h/2.0;a[1]=a[0];
a[2]=h;a[3]=h;
for(i=0;i<=n-1;i++)
{g[i]=y[i];y[i]=c[i];}
{(*f)(t,y,n,d);
for(i=0;i<=n-1;i++)
{b[i]=y[i];e[i]=y[i];}
for(k=0;k<=2;k++)
{for(i=0;i<=n-1;i++)
{y[i]=e[i]+a[k]*d[i];
b[i]=b[i]+a[k+1]*d[i]/3.0;
}
tt=t+a[k];
(*f)(tt,y,n,d);
}
for(i=0;i<=n-1;i++)
y[i]=b[i]+h*d[i]/6.0;
t=t+h;
}
free(g);free(b);free(c);free(d);free(e);
return;
}
主程序为:
#include <stdio.h>
#include<rkt3.c>
#include<math.h>
main()
{int i;
double t,h,y[3];
void rkt3f(double,double[],int,double[]);
y[0]=-1.0;y[1]=0.0;y[2]=1.0;
t=0.0;h=0.1;
printf("\n");
for(i=0;i<=2500;i++)
{rkt3(t,h,y,3,rkt3f);
t=t+h;
printf("%7.3f %e %e %e",t,y[0],y[1],y[2]);
printf("\n");
}
}
void rkt3f(double t,double y[],int n,double d[])
{t=t;n=n;
d[0]=y[1];
d[1]=(0.3+9.0*y[2]*y[2])*sin(y[0])*cos(y[0])-(0.3+9.8/1.5)*sin(y[0])-0.4*y[0];
d[2]=(1.1*cos(y[0])-0.3)/1.2-0.8*sin(1.0*t);
return;
}