如上图，要计算一个非线性常微分方程组，采用的是四阶龙格库塔算法，得到的结果在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);
}
