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

四阶龙格库塔法解一阶二元微分方程//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#includestdio.h#includeconio.h#includemath.h#includestdlib.h#defineN1000//定义运算步数；#defineh0.01//定义步长；floata,b,c;//定义全局变量常数a，b，c//定义微分方程：doublefx(doublex[],doubledx,doubley[],doubledy,doublez[],inti,doublek,doublexavg){intj;doublexi,yi;xi=x[i]+dx;yi=y[i]+dy;returnc*(xi-pow(xi,3)/3+yi)+k*(xavg-xi)+c*z[i];}doublefy(doublex[],doubledx,doubley[],doubledy,inti){doublexi,yi;xi=x[i]+dx;yi=y[i]+dy;return(xi-b*yi+a)/c;}voidmain(){doubleKx[3][4],Ky[3][4],x[3]={1,2,3},y[3]={2,3,4},xavg,k=0;//定义x,y的初值；doublez[3]={0};inti,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(Failedtoopenfile.\n);getch();return;}printf(Inputthevalueofconsta,b,c(seperatedby,eg0.1,0.2,0.3):);scanf(%f,%f,%f,&a,&b,&c);printf(Inputthethreevaluesofz(seperatedbyspacekey,eg0.10.20.3):);for(m=0;m3;++m){scanf(%lf,&z[m]);}printf(InputthevalueofStepstogetdifferentvaluesofxt,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;jN;++j){printf(%.3lf,j*h);fprintf(fp,%.3lf,j*h);xavg=0;for(i=0;i3;++i){xavg+=x[i];}xavg=xavg/3.0;//四阶龙格库塔法：for(i=0;i3;++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;i3;++i){printf(\t%.3lf,x[i]);fprintf(fp,\t%.3lf,x[i]);}for(i=0;i3;++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;n3;++n){fprintf(fp1,\t%.3lf,x[n]);}for(n=0;n3;++n){fprintf(fp1,\t%.3lf,y[n]);}break;}continue;}k+=0.1;}fclose(fp1);fclose(fp);printf(\nFinished.\nDatashavebeensavedto\sjy.txt,sjykxy.txt\.\n);getch();}