c++ - 嵌入式 Runge-Kutta 方法的步长控制器的正确实现

Question

我一直在尝试用 C++ 为嵌入式 Runge-Kutta 方法（目前是显式和 Rosenbrock）编写代码。这个想法是保持代码简单和通用，以便人们可以通过他们的 Butcher 画面（任何顺序）并运行它。

我已经验证了代码通常可以正常工作，但是在某些情况下（当我有一个非常复杂的 4 个微分方程系统时）步长控制无法适应（我得到恒定的步长或通常只是错误的）。

我使用的步长控制（我在本文中发现它是：

//beta is some safety parameter
fac=beta;


if(Delta<=1) {  //if current step is accepeted 
        if(delta_rej<=1){fac*=h/h_old; } //if previous step was rejected

        fac*=std::pow(Delta, -0.65/( (LD) method.p + 1.) );   
        fac*=std::pow( delta_acc/Delta, 0.3/ ( (LD) method.p + 1. ) );   
        h_stop=true ; //this is used exit a loop inside which the stepsize control is called

    }else{ //if current step is rejected
        fac*=std::pow( Delta , -1./((LD) method.p +1. ) );
    }
    
    //don't allow h to increase or decrease very much  
    if(fac> fac_max){fac = fac_max;}
    if(fac< fac_min){fac = fac_min;}
    
    //update h
    h=h*fac;

这里，h_old是之前接受的步长，当前试行步的步长是h。此外，Delta[1] 是当前尝试的相对（局部）误差估计（控制器尝试做出 ~1），delta_rej是前一次尝试delta_acc的 Delta，是上一个接受步骤的 Delta，并且method.p是方法（LD 是一个可以是 double 或 long double 的宏）。

我尝试过使用这个的简单版本（即 just fac*=std:: pow( Delta, -1./((LD) method.p +1. ) );），但似乎前一个版本更稳定一些。

例如，这些是我得到的直方图，显示了我的代码与 scipy 所采取的步骤数： Explicit RKF， Rosenbrock。如您所见，它们很接近，并且差异可能是由实现细节的差异引起的。

话虽如此，我仍然不确定，我真正想知道的是我是否正确使用了控制器。

谢谢

[1]：这是Delta的定义

标签： c++differential-equationsrunge-kutta