#include
void main()
{
void rkt1(double t,double y[],int n,double h,int k,double z[3][15001]);
int i,j;
double t,h,y
#include
void main()
{
void rkt1(double t,double y[],int n,double h,int k,double z[3][15001]);
int i,j;
double t,h,y[3],z[3][15001]={{0},{0},{0}};
y[0]=-1.0;
y[1]=0.0;
y[2]=1.0;
t=0.0;
h=0.01;
rkt1(t,y,3,h,15001,z);
printf("
");
for(i=0;i<=15000;i )
{
t=i*h;
printf("t=%5.2f
",t);
for(j=0;j<=2;j )
printf("y(%d)=%e ",j,z[j][i]);
printf("
");
}
}
void rkt1f(double t,double y[],int n,double d[])
{
double time,number;
time=t;
number=n;
d[0]=y[1];
d[1]=-y[0];
d[2]=-y[2];
}
void rkt1(double t,double y[],int n,double h,int k,double z[3][15001])
{
void rkt1f(double t,double y[],int n,double d[]);
int i,j,l;
double a[4],tt,b[32767],d[32767];
a[0]=h/2.0;
a[1]=a[0];
a[2]=h;
a[3]=h;
for(i=0;i<=n-1;i ) z[i][0]=y[i];
for(l=1;l<=k-1;l )
{
rkt1f(t,y,n,d);
for(i=0;i<=n-1;i ) b[i]=y[i];
for(j=0;j<=2;j )
{
for(i=0;i<=n-1;i )
{
y[i]=z[i][l-1] a[j]*d[i];
b[i]=b[i] a[j 1]*d[i]/3.0;
}
tt=t a[j];
rkt1f(tt,y,n,d);
}
for(i=0;i<=n-1;i )
y[i]=b[i] h*d[i]/6.0;
for(i=0;i<=n-1;i )
z[i][l]=y[i];
t=t h;
}
}
展开
全部