Boost C++ 库

首先，您必须指定表示系统 x 状态的数据类型。 在数学上，这通常是一个 n 维向量，其中实数或复数作为标量对象。 对于 odeint 来说，最自然的方式是使用 vector< double > 或 vector< complex< double > > 来表示系统状态。 但是，odeint 也可以处理其他容器类型，例如 std::array< double , N >，只要它满足下面定义的一些要求。

为了数值积分微分方程，还必须定义方程 x' = f(x) 的 rhs。 在 odeint 中，您以对象的术语提供此函数，该对象使用特定的参数结构实现 () 运算符。 因此，直接的方法是只定义一个函数，例如

/* The type of container used to hold the state vector */
typedef std::vector< double > state_type;

const double gam = 0.15;

/* The rhs of x' = f(x) */
void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
{
    dxdt[0] = x[1];
    dxdt[1] = -x[0] - gam*x[1];
}

该函数的参数必须遵循上面的示例，其中 x 是当前状态，这里是一个包含振荡器的位置 q 和动量 p 的两分量向量，dxdt 是导数 x'，应由函数填充为 f(x)，t 是当前时间。 请注意，在本例中，计算 f 不需要 t，但是 odeint 期望函数签名正好有三个参数（也有例外，稍后讨论）。

更复杂的方法是将系统实现为一个类，其中 rhs 函数被定义为该类的 () 运算符，其参数结构与上述相同

/* The rhs of x' = f(x) defined as a class */
class harm_osc {

    double m_gam;

public:
    harm_osc( double gam ) : m_gam(gam) { }

    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
    {
        dxdt[0] = x[1];
        dxdt[1] = -x[0] - m_gam*x[1];
    }
};

odeint 可以处理此类类的实例，而不是纯函数，这使得代码更简洁。

数值积分以迭代方式工作，这意味着您从状态 x(t) 开始，并执行长度为 dt 的时间步长，以获得近似状态 x(t+dt)。 存在许多不同的方法来执行这样的时间步长，每种方法都有一定的阶数 q。 如果一个方法的阶数为 q，则它精确到项 ~dt^q，这意味着此类步骤产生的 x 中的误差为 ~dt^q+1。 odeint 提供了几种不同阶数的步进器，请参阅 步进器概述。

上面表格中的一些步进器是特殊的：有些需要 ODE 的雅可比矩阵，另一些是为特殊的 ODE 系统（如哈密顿系统）构建的。 我们将在本教程中展示典型的示例和用例，以及应应用哪种步进器。

基本步进器仅执行一个时间步长，并且不给您任何关于所产生误差的信息（除了您知道它是 q+1 阶的）。 此类步进器与恒定步长一起使用，应选择足够小的步长以获得合理的误差。 但是，您应该对结果应用某种有效性检查（例如观察守恒量），因为您无法控制其他误差。 以下示例定义了一个基于经典 4 阶龙格-库塔方案的基本步进器。 步进器的声明需要状态类型作为模板参数。 现在可以使用 odeint 中的 integrate_const( Stepper, System, state, start_time, end_time, step_size ) 函数完成积分

runge_kutta4< state_type > stepper;
integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

此调用使用 RK4 方法从 t=0 到 10 积分由 harmonic_oscillator 定义的系统，步长为 dt=0.01，初始条件在 x 中给出。 结果 x(t=10) 存储在 x 中（就地）。 每个步进器都定义了一个 do_step 方法，也可以直接使用。 因此，您可以将上面的示例写成

const double dt = 0.01;
for( double t=0.0 ; t<10.0 ; t+= dt )
    stepper.do_step( harmonic_oscillator , x , t , dt );

	提示
	如果您有启用 C++11 的编译器，您可以轻松地使用 lambda 表达式来创建系统函数 { runge_kutta4< state_type > stepper; integrate_const( stepper , []( const state_type &x , state_type &dxdt , double t ) { dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; } , x , 0.0 , 10.0 , 0.01 ); }

为了改进数值结果并进一步最小化计算工作量，建议应用步长控制。 步长控制通过步进器算法实现，该算法还提供应用步长的误差估计。 odeint 提供了许多这样的 ErrorSteppers，我们将在 explicit_error_rk54_ck 的示例中展示它们的用法 - 一种 5 阶龙格-库塔方法，具有 4 阶误差估计和由 Cash 和 Karp 引入的系数。

typedef runge_kutta_cash_karp54< state_type > error_stepper_type;

给定误差步进器，仍然需要一个实例来检查误差并相应地调整步长。 在 odeint 中，这是通过 ControlledSteppers 完成的。 对于 runge_kutta_cash_karp54 步进器，存在一个 controlled_runge_kutta 步进器，可以通过以下方式使用

typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;
integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

与上面一样，这积分了由 harmonic_oscillator 定义的系统，但现在使用基于龙格-库塔 Cash-Karp 54 方案的自适应步长方法，从 t=0 到 10，初始步长为 dt=0.01（将调整），初始条件在 x 中给出。 结果 x(t=10) 也将存储在 x 中（就地）。

在上面的示例中，误差步进器嵌套在受控步进器中。 这是一种很好的技术； 然而，一个缺点是始终需要定义两个步进器。 也可以将受控步进器的实例化写入 integrate 函数的调用中，但仍然需要完全了解底层步进器类型。 另一点是，步长控制的误差容限不容易包含在受控步进器中。 这两个问题都可以通过使用 make_controlled 来解决

integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) ,
                    harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

make_controlled 可以与 odeint 的许多步进器一起使用。 第一个参数是绝对误差容限 eps_abs，第二个参数是相对误差容限 eps_rel，在积分期间使用。 模板参数确定应从哪个误差步进器实例化受控步进器。 make_controlled 的另一种语法是

integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) ,
                    harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

对于龙格-库塔控制器，一个步骤中产生的误差与 eps_abs + eps_rel * ( a_x * |x| + a_dxdt * dt * |dxdt| ) 进行比较。 如果误差小于此值，则接受当前步骤，否则拒绝该步骤并减小步长。 请注意，如果误差相对于上述关系的 rhs 太小，步长也会增加。 因此，controlled_runge_kutta 的完整实例化及其所有参数是

double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
controlled_stepper_type controlled_stepper(
    default_error_checker< double , range_algebra , default_operations >( abs_err , rel_err , a_x , a_dxdt ) );
integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

当使用 make_controlled 时，参数 a_x 和 a_dxdt 使用其标准值 1。

在下表中，可以找到所有与 make_controlled 和 make_dense_output 一起使用的步进器，后者是密集输出步进器的类似物。

当使用 make_controlled 或 make_dense_output 时，应注意使用的确切类型以及步长控制的工作方式。

odeint 支持使用迭代器来求解 ODE。 也就是说，您实例化一对迭代器，而不是使用带有适当观察器的 integrate 例程，而是将迭代器放入 C++ 标准库或 Boost.Range 中的算法之一。 一个例子是

std::for_each( make_const_step_time_iterator_begin( stepper , harmonic_oscillator, x , 0.0 , 0.1 , 10.0 ) ,
               make_const_step_time_iterator_end( stepper , harmonic_oscillator, x ) ,
               []( std::pair< const state_type & , const double & > x ) {
                   cout << x.second << " " << x.first[0] << " " << x.first[1] << "\n"; } );