boost odeint求解结果与Python3 scipy odeint差异过大问题及修复

我尝试使用boost odeint求解一个非常简单的常微分方程（ODE），部分场景下计算结果与Python的scipy odeint函数输出一致或误差极小，但更换初始条件后两者结果差异极大。
求解的方程为：d(uhat) / dt = - alpha^2 * kappa^2 * uhat，其中alpha取值为1.0，kappa为随测试场景变化的常量（具体取值见下文）。
我测试了boost提供的多种ODE求解器，均未能解决该问题。
更新：下方代码现已可以正常运行。
测试用例中，第一个场景输出几乎完全一致，第二个为验证性简单用例输出正常，第三个场景原先的C++版本会输出错误结果。

问题原因 #

原先手动调用runge_kutta_dopri5的do_step方法时，使用的是无自适应步长控制的定步长求解模式，而scipy的odeint默认使用带自适应步长控制的LSODA算法，修正后的代码用make_dense_output包装求解器后也开启了自适应步长控制。
第三个测试用例中kappa值很大，对应的时间衰减常数约为0.0012，远小于设定的定步长0.1，定步长求解会出现数值不稳定问题，导致结果完全偏离理论值，开启自适应步长后求解器会自动调整步长满足精度要求，结果就和scipy输出一致。

C++ 实现代码 #

#include <boost/numeric/odeint.hpp>
#include <cstdlib>
#include <iostream>
#include <vector>

typedef boost::numeric::odeint::runge_kutta_dopri5<double> Stepper_Type;

struct ResultsObserver
{
    std::ostream& m_out;
    ResultsObserver( std::ostream &out ) : m_out( out ) { }
    void operator()(const double& x , double t ) const
    {
      m_out << t << "  :  " << x << std::endl;
    }
};

// 方程右侧定义:  d_uhat_dt = - alpha^2 * kappa^2 * uhat
class Eq {
public:
  Eq(double alpha, double kappa)
    : m_constant(-1.0 * alpha * alpha * kappa * kappa) {}
  
  void operator()(double uhat, double& d_uhat_dt, const double t) const
  {
    d_uhat_dt = m_constant * uhat;
  }

private:
  double m_constant;
};

void integrate(double kappa, double initValue)
{
  const double dt = 0.1;
  const double alpha = 1.0;
  const double startTime = 0.0;
  const double endTime = 10.0;

  double uhat = initValue;      // 初始条件
  Eq rhs(alpha, kappa);         // ODE右侧表达式

  // 原错误写法：无自适应步长的定步长求解
  // boost::numeric::odeint::runge_kutta_dopri5<double> stepper;
  // std::vector<double> uhats;
  // for(unsigned step = 0; step < 100; ++step) {
  //   uhats.push_back(uhat);
  //   stepper.do_step(rhs, uhat, step*dt, dt);
  // }

  // 修正后写法：开启自适应步长控制
  integrate_const(
     boost::numeric::odeint::make_dense_output<Stepper_Type>( 1E-12, 1E-6 ),
     rhs, uhat, startTime, endTime, dt, ResultsObserver(std::cout) 
  );

  std::cout << "kappa = " << kappa << ", initial value = " << initValue << std::endl;
  std::cout << "---" << std::endl << std::endl;
}

int main() {
  const double kappa1 = 0.062831853071796;
  const double initValue1 = -187.097241230045967;
  integrate(kappa1, initValue1);

  const double kappa2 = 28.274333882308138;
  const double initValue2 = 0.000000000000;
  integrate(kappa2, initValue2);

  const double kappa3 = 28.337165735379934;
  const double initValue3 = -0.091204068895190;
  integrate(kappa3, initValue3);
  
  return EXIT_SUCCESS;
}

Python3 实现代码 #

#!/usr/bin/env python3
import numpy as np
from scipy.integrate import odeint

def Eq(uhat, t, kappa, a):
    d_uhat = -a**2 * kappa**2 * uhat
    return d_uhat

def integrate(kappa, initValue):
    dt = 0.1
    t = np.arange(0,10,dt)
    a = 1.0
    print("kappa =", kappa)
    print("initValue =", initValue)
    uhats = odeint(Eq, initValue, t, args=(kappa,a))
    print(uhats)
    print("---\n")

kappa1 = 0.062831853071796
initValue1 = -187.097241230045967
integrate(kappa1, initValue1)

kappa2 = 28.274333882308138
initValue2 = 0.000000000000
integrate(kappa2, initValue2)

kappa3 = 28.337165735379934
initValue3 = -0.091204068895190
integrate(kappa3, initValue3)

内容的提问来源于stack exchange，提问作者Madeleine P. Vincent