文档介绍：四阶龙格库塔算法//状态方程/X=AX+BU化为微分方程求解#include<>#include<>#include<>#include<>//定义运算步数;//定义步长;#defineN1000000#[9]={1,2,3,4,5,6,7,8,9};doubleu[4]={0};//定义微分方程:doublefx(doublex[],doubledx,inti){//矩阵A和BdoubleA[9][9]={{0,0,0,0,0,0,1,,},{0,0,0,0,0,0,0,,-},{0,0,0,0,0,0,0,,},{,-,0,-,,,,,},{,-,0,-,-,,-,,},{-,-,0,,,-,,,-},{0,0,0,-,-,,-,,},{0,0,0,,-,-,-,-,-},{0,0,0,,,-,,-,-}};doubleB[9][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{,-,,-},{-,-,,-},{,-,-,0},{,-,,-},{-,,-,-},{0,-,-,}};doublea[9]={};doubleb[9]={};for(intj=0;j<9;j++){a[i]+=A[i][j]*(x[j]+dx);}for(intk=0;k<4;k++){b[i]+=B[k][0]*u[k];}returna[i]+b[i];}//主函数voidmain(){doubleinteger[3]={0};doubleKx[9][4];intcount=0;intn,S;FILE*fp1,*fp,*fp_integer;fp=fopen("","w");fp1=fopen("","w");fp_integer=fopen("","w");fprintf(fp1,"x1\tx2\tx3\tx4\tx5\tx6\tx7\tx8\tx9\n");fprintf(fp_integer,"S1\tS2\tS3\n");if(fp==NULL||fp1==NULL){printf("Failedtoopenfile.\n");getch();return;}printf("Inputthevalueofconstu1,u2,u3,u4(seperatedby,,,,):");scanf("%f,%f,%f,%f",&u[0],&u[1],&u[2],&u[3]);printf("InputthevalueofStepstogetdifferentvaluesofxt,yt(S):");scanf("%d",&S);for(intj=1;j<N;++j){printf("%.3lf",j*h);fprintf(fp,"%.3lf",j*h);//四阶龙格库塔法:for(inti=0;i<9;i++){Kx[i][0]=fx(x,0,i);Kx[i][1]=fx(x,h/2*Kx[i][0],i);Kx[i][2]=fx(x,h/2*Kx[i][1],i);Kx[i][3]=fx(x,h*Kx[i][2],i);x[i]=x[i]+(Kx[i][0]+2*Kx[i][1]+2*Kx[i][2]+Kx[i][3])/6*h;}for(i=0;i<9;++i){//三位有效数字printf("\t%.3lf",x[i]);fprintf(fp,"\t%.3lf",x[i]);}printf("\n");fprintf(fp,"\n");count++;if(count<=S){integer[0]+=x[3]*h;int