龙格库塔 C语言程序
下面的程序运行800多部的时候出现错误,请帮忙确认一下
#include "stdio.h"
#include "math.h"
#include<graphics.h>
#include<conio.h>
void fun(double y[4],double d[4])
{
d[0]=(0.01-y[0]-2*((y[2]>0)?y[2]:0)+5-2.5*y[1])/2;
d[1]=(0-y[1]+((y[0]>0)?y[0]:0))/12;
d[2]=(0-y[2]-2*((y[0]>0)?y[0]:0)+5-2.5*y[3])/2;
d[3]=(0-y[3]+((y[2]>0)?y[2]:0))/12;
}
void rkg(double h,double y[4],double d[4])
{
double a21=h/2;
double a31=(sqrt(2.0)-1)/2*h;
double a32=(2-sqrt(2.0))/2*h;
double a42=-sqrt(2.0)/2*h;
double a43=(2+sqrt(2.0))/2*h;
double b1=1*h/6;
double b2=(2-sqrt(2.0))*h/6;
double b3=(2+sqrt(2.0))*h/6;
double b4=1*h/6;
double k1[4],k2[4],k3[4],k4[4];
double temp_y[4];
int i=0;
for(i=0;i<4;i++)
temp_y[i]=y[i];
fun(temp_y,k1);
for(i=0;i<4;i++)
temp_y[i]=y[i]+a21*k1[i];
fun(temp_y,k2);
for(i=0;i<4;i++)
temp_y[i]=y[i]+a31*k1[i]+a32*k2[i];
fun(temp_y,k3);
for(i=0;i<4;i++)
temp_y[i]=y[i]+a42*k2[i]+a43*k3[i];
fun(temp_y,k4);
for(i=0;i<4;i++)
d[i]=y[i]+b1*k1[i]+b2*k2[i]+b3*k3[i]+b4*k4[i];
}
void main()
{ FILE *fp;
int i;
double t,h,z[1000][5];
z[0][1]=0.0; z[0][2]=0.0; z[0][3]=0.0;z[0][0]=0.0;
t=0.0;h=0.05;
for(i=0;i<1000;i++)
{
rkg(h,z[i],z[i+1]);
fprintf(fp,"%f,%f,%f,%f,%f\n",t+h*i,z[i][0],z[i][1],z[i][2],z[i][3]);
}
for(i=0;i<1000;i++)
{printf("%f,%f,%f,%f,%f\n",t+h*i,z[i][0],z[i][1],z[i][2],z[i][3]);
/*getch();*/}
}
void graph()
{
int gdrive,gmode=VGAHI;
gdrive=VGA;
initgraph(&gdrive,&gmode,"d:\tc20");
closegraph();
}