rk.c

四阶R-K方法求解微分方程.

//dxi/dt=c*(xi-xi^3/3+yi)+K*(X-xi)+c*zi
//dyi/dt=(xi-b*yi+a)/c
//i=1,2,3
//X=sum(xi)/N

#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>

#define N 1000 //定义运算步数；
#define h 0.01 //定义步长；

float a,b,c;//定义全局变量常数a，b，c

//定义微分方程：
double fx(double x[],double dx,double y[],double dy,double z[],int i,double k,double xavg){
int j;
double xi,yi;
xi=x[i]+dx;
yi=y[i]+dy;
return c*(xi-pow(xi,3)/3+yi)+k*(xavg-xi)+c*z[i];
}

double fy(double x[],double dx,double y[],double dy,int i){
double xi,yi;
xi=x[i]+dx;
yi=y[i]+dy;
return (xi-b*yi+a)/c;
}

void main(){
double Kx[3][4],Ky[3][4],x[3]={1,2,3},y[3]={2,3,4},xavg,k=0;//定义x,y的初值；
double z[3]={0};
int i,j,m,n,S;
FILE *fp1,*fp;
fp=fopen("sjy.txt","w");
fp1=fopen("sjykxy.txt","w");
fprintf(fp1,"k\tx1\tx2\tx3\ty1\ty2\ty3\n");
if(fp==NULL||fp1==NULL){
   printf("Failed to open file.\n");
   getch();
   return;
}
   printf("Input the value of const a,b,c(seperated by ,eg 0.1,0.2,0.3):");
   scanf("%f,%f,%f",&a,&b,&c);
   printf("Input the three values of z(seperated by spacekey,eg 0.1 0.2 0.3):");
   for (m=0;m<3;++m)
   {
    scanf("%lf",&z[m]);
   }
   printf("Input the value of Steps to get different values of xt,yt(S):");
   scanf("%d",&S);
while(k<=1)
{ 
fprintf(fp,"k=%.3f\n",k);
    fprintf(fp,"t\tx1\tx2\tx3\ty1\ty2\ty3\n");
fprintf(fp1,"\n%.3f",k);
    for(j=1;j<N;++j)
    {
   printf("%.3lf",j*h);
   fprintf(fp,"%.3lf",j*h);
   xavg=0;
   for(i=0;i<3;++i)
     {
    xavg+=x[i];
   }
   xavg=xavg/3.0;
   //四阶龙格库塔法：
   for(i=0;i<3;++i)
   {
    Kx[i][0]=fx(x,0,y,0,z,i,k,xavg);
    Ky[i][0]=fy(x,0,y,0,i);
    Kx[i][1]=fx(x,h/2*Kx[i][0],y,h/2*Ky[i][0],z,i,k,xavg);
    Ky[i][1]=fy(x,h/2*Kx[i][0],y,h/2*Ky[i][0],i);
    Kx[i][2]=fx(x,h/2*Kx[i][1],y,h/2*Ky[i][1],z,i,k,xavg);
    Ky[i][2]=fy(x,h/2*Kx[i][1],y,h/2*Ky[i][1],i);
    Kx[i][3]=fx(x,h*Kx[i][2],y,h*Ky[i][2],z,i,k,xavg);
    Ky[i][3]=fy(x,h*Kx[i][2],y,h*Ky[i][2],i);
    x[i]=x[i]+(Kx[i][0]+2*Kx[i][1]+2*Kx[i][2]+Kx[i][3])/6*h;
    y[i]=y[i]+(Ky[i][0]+2*Ky[i][1]+2*Ky[i][2]+Ky[i][3])/6*h;
   }
   for(i=0;i<3;++i){
    printf("\t%.3lf",x[i]);
    fprintf(fp,"\t%.3lf",x[i]);
   }
   for(i=0;i<3;++i){
    printf("\t%.3lf",y[i]);
    fprintf(fp,"\t%.3lf",y[i]);
   }
   printf("\n");
   fprintf(fp,"\n");
        //取第S步，即时间为S*h的x1，x2，x3，y1，y2，y3随k值的变化；
   while(j==S)
   {
    for(n=0;n<3;++n)
    {
     fprintf(fp1,"\t%.3lf",x[n]);
    }
    for(n=0;n<3;++n)
    {
     fprintf(fp1,"\t%.3lf",y[n]);
    }
    break;
   }
   continue;
}k+=0.1;
}
fclose(fp1);
fclose(fp);
printf("\nFinished.\nDatas have been saved to \"sjy.txt,sjykxy.txt\".\n");
getch();
}