我需要使用自适应RK4方法集成一个ODES系统，并通过步长加倍技术进行步长控制。
问题是程序不断地将步长缩小到机器精度，而不提前时间。
的想法是一步一步的解决方案，也通过两个连续的半步，比较结果的差异，并将其存储在eps。所以eps是误差的度量。现在我想根据eps是否大于指定精度的eps0(如"数值食谱"一书中所述)来确定下一步的步长
RK4Step(double t, double* Y, double *Yout, void (*RHSFunc)(double, double *, double *),double h)将解向量Y按h步进，并使用RHSFunc函数将结果放入Yout。
#define NEQ 4 //problem dimension
int main(int argc, char* argv[])
{
ofstream frames("./frames.dat");
ofstream graphs("./graphs.dat");
double Y[4] = {2.0, 2.0, 1.0, 0.0}; //initial conditions for solution vector
double finaltime = 100; //end of integration
double eps0 = 10e-5; //error to compare with eps
double t = 0.0;
double step = 0.01;
while(t < finaltime)
{
double eps = 0.0;
double Y1[4], Y2[4]; //Y1 will store half step solution
//Y2 will store double step solution
double dt = step; //cache current stepsize
for(;;)
{
//make a step starting from state stored in Y and
//put solution into Y1. Then from Y1 make another half step
//and store into Y1.
RK4Step(t, Y, Y1, RHS, step); //two half steps
RK4Step(t+step,Y1, Y1, RHS, step);
RK4Step(t, Y, Y2, RHS, 2*step); //one long step
//compute eps as maximum of differences between Y1 and Y2
//(an alternative would be quadrature sums)
for(int i=0; i<NEQ; i++)
eps=max(eps, fabs( (Y1[i]-Y2[i])/15.0 ) );
//if error is within tolerance we grow stepsize
//and advance time
if(eps < eps0)
{
//stepsize is accepted, grow stepsize
//save solution from Y1 into Y,
//advance time by the previous (cached) stepsize
Y[0] = Y1[0]; Y[1] = Y1[1];
Y[2] = Y1[2]; Y[3] = Y1[3];
step = 0.9*step*pow(eps0/eps, 0.20); //(0.9 is the safety factor)
t+=dt;
break;
}
//if the error is too big we shrink stepsize
step = 0.9*step*pow(eps0/eps, 0.25);
}
}
frames.close();
graphs.close();
return 0;
}
您从未在内循环中重置eps。这可能是你的问题的直接原因。虽然实际误差随着步长的减小而减小，但eps中的最大值保持不变，并且高于eps0。这将导致步长更新中的常数减小因子，而不会有任何机会打破循环。
又一个"错误的事情";误差估计和公差是不相容的。容错eps0是一个误差密度或单位步长误差。要将错误估计eps转换成这种格式，您需要将eps除以step。或者换句话说，目前你正在强迫实际的步骤误差接近0.5*eps0，所以全局误差是0.5*eps0乘以所采取的步骤数，而步骤数与eps0^0.2大致成正比。在使用单位步长误差的版本中，局部误差被强制为"动态的"。接近0.5*eps0*step，使得全局误差约为5*eps0乘以积分间隔的长度。我想说，第二种变体更符合对预期行为的直觉。
这不是一个临界错误，但可能导致次优步长和实际的全局误差，偏离期望的误差容忍度。
你也有一个编码不一致，因为在状态的传播和状态向量的声明中，你在状态向量中硬编码了4个组件，而在错误计算中，你在方程和组件的变量数NEQ上有一个循环。在使用c++时，可以使用状态向量类在内部处理所有依赖维度的循环。(如果做得太过，频繁地分配生命周期短的实例可能会导致效率问题。)