runge_kutta.cpp

四阶龙格库算法,

float s;
#include <stdio.h>
#include <math.h>
#define M 8		

void rungek(int n,float t,float h,float y[]);
void fct(int n,float x,float y[],float f[]);
void main()
{
	float y[M],t0,h;
	int n,nk,i,k;
	FILE *fp;
	n=M;
	printf("h-步长，T0-自变量初值，NK-计算点数\n");
	scanf("%f %f %d",&h,&t0,&nk);
    printf("输入初值y[i]\n");
	for(i=0;i<n;i++)
	{
		scanf("%f",&y[i]);
	}
    fp=fopen("result.txt","w");
	printf("%6.4f  ",t0);
	fprintf(fp,"%6.4   ",t0);
    for(i=0;i<n;i++)
	{
		printf("%9.6f   ",y[i]);
		fprintf(fp,"%9.6f   ",y[i]);
	}
     printf("\n");
	 fprintf(fp,"\n");
     s=t0;
	 for(k=0;k<nk;k++)
	 {
		 rungek(n,s,h,y);
		 printf("%6.4f   ",s);
		 fprintf(fp,"%6.4f   ",s);
	     for(i=0;i<n;i++)
		  {
		   printf("%9.6f  ",y[i]);
	       fprintf(fp,"%9.6f   ",y[i]);

		  }
	 printf("\n");
	 fprintf(fp,"\n");
	 }	
	 fclose(fp);
} 
void fct(int n,float x,double y[],double f[])
{
     f[0]=y[1];
	 f[1]=18462*y[1]-11*(10^6)*y[0]*y[2];
	 f[2]=y[3];
	 f[3]=18462*y[5]-11*(10^6)*y[0]*y[2];
	 f[4]=y[5];
	 f[5]=18462*y[5]-1.56*(10^8)*(y[6]-y[4])-4.5*(10^10)*y[0]*y[2];
	 f[6]=y[7];
	 f[7]=(1.6*(10^6)*(y[6]-y[4])-0.73*(10^(-8))*y[6]*y[7]*y[6]*y[7])/(0.54-2.4*(10^(-9))*y[6]*y[6]*y[6]);
 
	return;
}
void rungek(int n,float t,float h,float y[])
{
	float u[5],xw;
	double yw[M];
	double yf[M],ym[M];

	int i,j;
	u[0]=0.5*h;
	u[1]=u[0];
	u[2]=h;
	u[3]=h;
	u[4]=u[1];
	xw=t;
    for(i=0;i<n;i++)
	{
		ym[i]=y[i];
		yw[i]=y[i];
	}
	for(j=0;j<4;j++)
	{
		fct(n,t,ym,yf);
		t=xw+u[j];
	 
		for(i=0;i<n;i++)
		{
			ym[i]=yw[i]+u[j]*yf[i];
			y[i]=y[i]+u[j+1]*yf[i]/3.0;
		}
	}
	s=t;
	return;
}