四阶龙格库塔法解一阶二元微分方程

四阶龙格库塔法解一阶二元微分方程
//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>
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;
#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;
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);
}
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");
FILE *fp1,*fp;
fp=fopen("sjy.Leabharlann Baiduxt","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");
}
break;
}
continue;
}k+=0.1;
}
fclose(fp1);
fclose(fp);
printf("\nFinished.\nDatas have been saved to \"sjy.txt,sjykxy.txt\".\n");
getch();
}
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);
{
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);
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]);
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;
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):");
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)