使用Boost.Odeint求解复数微分方程组的编译问题排查

问题场景 #

尝试用Boost.Odeint库求解2个复数微分方程组，编写的代码如下：

#include <iostream>
#include <boost/array.hpp>
#include <complex>

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

const double sigma = 10.0;
 
const double gamma = 0.1; 
const double J12 = 1; 
const double G = 1; 
const double eta = 0.01; 
const double gamma_s = 0.1; 

const int n = 2; 

typedef boost::array< complex<double> , n > state_type;
//typedef runge_kutta4< complex<double> > stepper_type; 试过这个
//typedef runge_kutta_dopri5< complex<double> > stepper_type; 还有这个，但都报大量错误


void lorenz( const state_type &x , state_type &dxdt , double t )
{
    dxdt[0] =  -0.5*gamma*x[0] 
            - eta * (abs(x[0])*x[0])
            + G * conj(x[0])
            + J12 * gamma_s * (x[0]+x[1]);

    dxdt[1] =  -0.5*gamma*x[1] 
            - eta * (abs(x[1])*x[1])
            + G * conj(x[1])
            + J12 * gamma_s * (x[0]+x[1]);
}

编译时出现大量类型匹配错误，示例错误信息：

include/boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp:142:52:
error: no matching function for call to
    'boost::numeric::odeint::default_operations::scale_sum4<std::complex<double>, double,
     double, double>::scale_sum4(double, std::complex<double>, std::complex<double>,
     std::complex<double>)'
142 |                 typename operations_type::template scale_sum4< value_type ,
      time_type , time_type , time_type >( 1.0 , dt*b41 , dt*b42 , dt*b43 ));

错误原因 #

Boost.Odeint的默认操作类无法处理标量（double）与复数（complex）的混合运算，默认的代数类型（如range_algebra）依赖逐点操作，在类型混合时会出现模板参数不匹配的问题。

解决方法 #

显式指定stepper的完整模板参数，使用vector_space_algebra——该代数类型原生支持向量空间运算，可正确处理标量与复数的乘法、加法操作。

针对Runge-Kutta DOPRI5的修改 #

typedef runge_kutta_dopri5<
    state_type,          // 状态变量类型
    double,              // 时间值类型
    state_type,          // 导数变量类型
    double,              // 时间步长类型
    vector_space_algebra // 适配复数运算的代数类型
> stepper_type;

针对Runge-Kutta4的修改 #

typedef runge_kutta4<
    state_type,
    double,
    state_type,
    double,
    vector_space_algebra
> stepper_type;

补充完整代码示例 #

为确保代码可运行，需添加调用积分器的main函数：

int main()
{
    // 初始化初始状态
    state_type x = {complex<double>(1.0, 0.0), complex<double>(0.0, 1.0)};
    stepper_type stepper;
    
    // 调用积分器：从t=0到t=10，步长0.1
    integrate_const(stepper, lorenz, x, 0.0, 10.0, 0.1);
    
    // 输出最终状态
    cout << "最终状态: " << x[0] << ", " << x[1] << endl;
    return 0;
}

内容的提问来源于stack exchange，提问作者JTT