Gimme梦 Gimme梦
关注数: 16 粉丝数: 78 发帖数: 453 关注贴吧数: 12
非线性常微分方程组的解失真是怎么回事? 如上图,要计算一个非线性常微分方程组,采用的是四阶龙格库塔算法,得到的结果在0-0.6左右计算还是可以的,但到后来的计算结果就完全失真发散了。网上查找了很多资料也没找到什么原因,有人能帮忙看一下吗?附带代码如下: #include <stdio.h> #include <math.h> #define M 3 // 维数,方程组中方程的个数// 算斜率的右端函数 void RightSlopeFun( double *pRightK, //算出来的右端斜率项[out] double iT, //计算点的微分自变量值 double *iY //计算点的函数初始值 ) { pRightK[0] = -(iY[0]-iY[1])-(2*iY[0]*iY[2]+iY[1]*iY[0]+iY[1]*iY[2]-iY[2]*iY[2])+3; pRightK[1] = -(iY[1]+iY[2])-(2*iY[0]*iY[0]/3-iY[0]*iY[1]/3+iY[1]*iY[0]-2*iY[1]*iY[2]/3+2*iY[2]*iY[0]/3+4*iY[2]*iY[1]/3+iY[2]*iY[2]/3)+1; pRightK[2] = -iY[1]-(-iY[0]*iY[0]/3+2*iY[0]*iY[1]/3-iY[0]*iY[2]-iY[1]*iY[0]+iY[1]*iY[2]/3-iY[2]*iY[0]/3-2*iY[2]*iY[1]/3+4*iY[2]*iY[2]/3)-1; } // 龙格库塔法解算程序 double LgktSolution( double *pOY, // 返回函数值[out] double iT, // 积分自变量输入值 double *iY, // 初始值 double h, // 单步步长 int Count=1 // 解算步数,默认形参为1,调用时不给此 // 参数赋值,则只解算一步 ) { double K1[M]; // 第一个斜率 double K2[M]; // 第二个斜率 double K3[M]; // 第三个斜率 double K4[M]; // 第四个斜率 double tempY[M]; // 保存Y的中间值的临时空间 int i; // 计数器 double tempT; // 临时保存x的空间 tempT=iT; // 载入初始的x值 for (i=0;i<M;i++) { tempY[i]=iY[i]; // 载入初始的函数值 pOY[i]=iY[i]; } while (Count>0) { RightSlopeFun(K1,tempT,tempY); // 解算斜率K1 tempT += h/2; for (i=0;i<M;i++) { pOY[i]=tempY[i]+K1[i]*h/2; } RightSlopeFun(K2,tempT,pOY); // 解算斜率K2 for (i=0;i<M;i++) { pOY[i]=tempY[i]+K2[i]*h/2; } RightSlopeFun(K3,tempT,pOY); // 解算斜率K3 tempT += h/2; for (i=0;i<M;i++) { pOY[i]=tempY[i]+K3[i]*h; } RightSlopeFun(K4,tempT,pOY); // 解算斜率K4 for (i=0;i<M;i++) { // 得到函数值 pOY[i]=tempY[i]+h*(K1[i]+2*K2[i]+2*K3[i]+K4[i])/6; // 推算了一步的函数值,此处可以对该数据进行相应的操作 } // 递推数据更新 Count--; if (Count<=0) { break; } for (i=0;i<M;i++) { tempY[i]=pOY[i]; } } return tempT; } void main() { FILE *fp=fopen("g:\\随机振动相关论文及程序代码\\试算\\非平稳解\\测试组\\无耦合项.txt","a"); double iy[M]; double it; int i; iy[0]=1; iy[1]=0.00001; iy[2]=0.00001; it=0; for(i=0;i<1000;i++) { it=LgktSolution(iy,it,iy,0.01); fprintf(fp,"%-8.2f %-8.6f %-8.6f %-8.6f %-8.6f %-8.6f %-8.6f\n",it,iy[0],iy[1],iy[2],iy[3],iy[4],iy[5],iy[6],iy[7],iy[8],iy[9],iy[10],iy[11],iy[12],iy[13],iy[14],iy[15]); } fclose(fp); }
1 下一页