龙格库塔法 C语言程序

输入文件runge_kutta_input.txt中内容：
0
1
1
10

以下为代码部分：

#include
#include "stdlib.h"
#include
#include

void size(int *n) ;
void input(float *a,float *b,float *y0,int *n);
void RK(int n,float y0,float a,float b,float x[],float y[]);
float f(float x,float y);
float af(float x);



void main()
{
int n;
float a,b,y0,*x,*y;
size(&n) ;
x=(float*)malloc((n+1)*sizeof(float));
y=(float*)malloc((n+1)*sizeof(float));
input(&a,&b,&y0,&n);
RK(n,y0,a,b,x,y);
system("pause");
free(x);free(y);
}


void size(int *n) /* 读取维数n */
{
FILE *fp;
if((fp=fopen("runge_kutta_input.txt","r"))==NULL)
{
printf("can not open file.\n");

exit(0);
}
fscanf(fp,"%d",n);
fclose(fp);
}



void input(float *a,float *b,float *y0,int *n)
{ FILE *fp;

if((fp=fopen("runge_kutta_input.txt","r"))==NULL)
{
printf("cannot open file runge_kutta.txt!\n");
}
printf("\n二阶Runge_Kutta法求初值：\n");
printf("\n区间：");
fscanf(fp,"%f",a);
fscanf(fp,"%f",b);
printf("[a,b]=[%f,%f]",*a,*b) ;
printf("\n") ;
fscanf(fp,"%f",y0);
printf("\n初值：");
printf("y0=%f\n",*y0) ;
fscanf(fp,"%d",n);
printf("\n计算步数：");
printf("n=%d\n",*n) ;
fclose(fp);
}



void RK(int n,float y0,float a,float b,float x[],float y[])
{

int i;
float x1,y1,h;
float K1,K2;
float *yy,*dy;
FILE *fp;
char *filename="runge_kutta_output.txt";
yy=(float*)malloc((n+1)*sizeof(float));
dy=(float*)malloc((n+1)*sizeof(float));
if((fp=fopen(filename,"w"))==NULL)
{
printf("can not open file.\n");
getchar();
exit(0);
}
h=0.1;
x1=a;
y1=y0;
x[0]=x1;
y[0]=y1;
yy[0]=y0;
dy[0]=0;
printf("\n常微分方程y-2*x/y=y'的数值解为：\n");
fprintf(fp,"常微分方程y-2*x/y=y'的数值解为：\n");
printf("\n\txn\t\t");
fprintf(fp,"\n\txn\t\t");
printf("yn\t\t");
fprintf(fp,"yn\t\t");
printf("y(xn)\t\t");
fprintf(fp,"y(xn)\t\t");
printf("|y(xn)-yn|\n\n");
fprintf(fp,"|y(xn)-yn|\n\n");
printf("x[0]=%f\t",x[0]) ;
fprintf(fp,"x[0]=%f\t",x[0]) ;
printf("y[0]=%f\t",y[0]) ;
fprintf(fp,"y[0]=%f\t",y[0]) ;
printf("y(x[0])=%f\t",yy[0]) ;
fprintf(fp,"y(x[0])=%f\t",yy[0]) ;
printf("dy(x[0])=%f\n",dy[0]) ;
fprintf(fp,"dy(x[0])=%f\n",dy[0]) ;
for(i=1;i<=n;i++)
{
/*K1=f(x1,y1);
K2=f((x1+h),(y1+h*K1));
y1=y1+h*((K1/2)+(K2/2)); */
K1=f(x1,y1);
K2=f((x1+3*h/4),(y1+3*h*K1/4));
y1=y1+h*((K1/3)+(2*K2/3));
x1=a+i*h;
x[i]=x1;
printf("x[%d]=%f\t",i,x[i]);
fprintf(fp,"x[%d]=%f\t",i,x[i]);

y[i]=y1;
printf("y[%d]=%f\t",i,y[i]);
fprintf(fp,"y[%d]=%f\t",i,y[i]);

yy[i]=af(x[i]);
printf("y(x[%d])=%f\t",i,yy[i]);
fprintf(fp,"y(x[%d])=%f\t",i,yy[i]);

dy[i]=fabs(y[i]-yy[i]) ;
printf("dy(x[%d])=%f\n",i,dy[i]);
fprintf(fp,"dy(x[%d])=%f\n",i,dy[i]);
}
free(yy);free(dy) ;
fclose(fp);
}

float f(f

loat x,float y)
{
float fun;
fun=y-2*x/y;
return fun;
}

float af(float x)
{
float fun;
fun=sqrt(1+2*x);
return fun;
}