odeint v2 - 在 C++ 中求解常微分方程。

下载 odeint-v2.zip - 810.97 KB

引言

本文介绍了 odeint 的第二个版本——一个用于求解常微分方程 (ODE) 的 C++ 框架。它基于模板元编程，不依赖于特定的容器类型，并可与现代图形卡配合使用。在 odeint 的内部，当前版本进行了大量更改。它与第一个版本不直接兼容，但只需进行少量更改即可使版本一和版本二兼容。在此 codeproject 上已简要介绍 odeint 的第一个版本。

当前版本的 odeint-v2 可从 github 获取。完整的文档包含大量教程和示例，也可在 github 上找到。odeint 已得到 Google summer of code 项目的支持。计划在 boost 库中发布 odeint，尽管尚不清楚 odeint 是否能通过 boost 的审查过程。

背景

求解常微分方程在物理、化学、生物甚至社会系统的数学建模中是一项非常重要的任务。分析 ODE 求解方法有着悠久的传统。这些方法总结在数学学科数值分析中。常微分方程总是具有以下形式

d x(t) / d t = f( x(t) , t )

x 被称为 ODE 的状态，可以是向量形式。t 是因变量。我们始终将 t 称为时间，尽管它可能有不同的（物理）含义。函数 f 定义了 ODE。每个 ODE 都关联一个初值问题 (IVP)，该问题由 ODE 和初值 x(t0)=x0 给出。初值问题的解是 x(t) 的时间演化，并附加了 x(t0)=x0 的条件。可以证明，对于大多数与物理相关的 ODE，每个 IVP 都有唯一解。更多细节请参阅第一个版本文章的背景部分，或查阅一些关于微分方程或数值分析的经典教科书。

有什么新功能？

• 代数和运算的引入。 它们允许您更改基本数学运算的执行方式。利用这些概念，可以求解 CUDA 中的常微分方程，或使用 Intel Math Kernel 库或 SIMD 库等复杂的数值库。
• 密集输出。 odeint 现在支持带有密集输出的方法。密集输出意味着解在中间点之间进行插值，通常与数值解具有相同的阶数。已引入一个用于密集输出步进器的专用概念，并且积分函数支持带有密集输出的方法。
• 经典龙格-库塔步进器的元函数和类。 已引入一个用于泛化龙格-库塔求解器的新概念。该概念在很大程度上基于模板元编程。特定的求解器仅由其系数定义。这些方法的性能与手工编写的版本相当。
• 新方法。 已引入新方法，例如 Dormand-Prince 5(4) 方法，求解刚性系统的求解器（如 Rosenbrock-4），以及多步方法（如 Adams-Bashforth-Moulton）。还有一个尝试实现 Taylor 级数方法。
• 对 Fancy 库的支持。 odeint 可以很好地与 Boost.Range 和 Boost.Units 配合使用。前者允许您仅求解完整状态的一部分，而后者在编译时执行数学方程的维度检查。
• 工厂函数。 用于轻松生成受控步进器和密集输出步进器。

示例

在详细介绍 odeint 之前，我们将展示一些示例。

洛伦兹系统的积分

让我们从以前文章中的一个例子开始——著名的洛伦兹系统及其蝴蝶状的翅膀结构。使用 odeint 非常简单，只需包含相应的头文件即可。odeint 完全是头文件形式的。您无需链接其他库。这样做的优点是编译器可以进行强大的优化技术，从而产生非常快速的运行时代码。要使用 odeint，只需包含

#include <boost/numeric/odeint.hpp>

整个库位于命名空间 boost::numeric::odeint 中

using namespace boost::numeric::odeint;

洛伦兹系统的定义要么是一个函数，要么是一个函数对象。但在声明 ODE 之前，我们需要定义状态类型，即 odeint 在其步进器内部用于存储中间值并且在系统函数中使用的类型。

typedef std::vector< double > state_type;

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;

// the system function can be a classical functions
void lorenz( state_type &x , state_type &dxdt , double t )
{                                                         
    dxdt[0] = sigma * ( x[1] - x[0] );
    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
    dxdt[2] = x[0]*x[1] - b * x[2];
}

// the system function can also be a functor
class lorenz_class
{
 public:
    void operator()( state_type &x , state_type &dxdt , double t )
    {
        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = x[0]*x[1] - b * x[2];
    }
};

现在，我们原则上可以开始求解洛伦兹系统了。首先，我们必须选择一个步进器。让我们从四阶经典龙格-库塔步进器开始。我们将系统积分 1000 步，步长为 0.01。

state_type x( 3 );
x[0] = x[1] = x[2] = 10.0;
const double dt = 0.01;
integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 10.0 , dt );
// or use the functor:
// integrate_const( runge_kutta4< state_type >() , lorenz_class() , x , 0.0 , 10.0 , dt );

除了 integrate_const 函数，您也可以直接调用步进器。这可能有一些优点，如果您想在连续步骤之间执行非平凡代码，但我们建议您尽可能使用积分函数和观察者。对于受控步进器和密集输出步进器，积分函数实现了适当积分所需的复杂逻辑。尽管如此，上面的代码与以下代码完全等效：

state_type x( 3 );
x[0] = x[1] = x[2] = 10.0;
const double dt = 0.01;
runge_kutta4< state_type > rk4;
double t = 0.0;
for( size_t i=0 ; i<1000 ; ++i,t+=dt )
{
    rk4.do_step( lorenz , x , t , dt );
}

在大多数情况下，您希望在积分过程中对状态执行某些操作，例如将状态写入文件或执行一些花哨的统计方法。您可以创建一个观察者并使用该观察者调用积分函数。

struct streaming_observer
{
     std::ostream &m_out;
     streaming_observer( std::ostream &out ) : m_out( out ) {}

     void operator()( const state_type &x , double t ) const
     {
          m_out << t;
          for( size_t i=0 ; i < x.size() ; ++i )
              m_out << "\t" << x[i];
          m_out << "\n";
     }
};

// ...

integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 10.0 , dt , streaming_observer( std::cout ) );

在之前的示例中，步进器始终是四阶龙格-库塔步进器。这个步进器确实很不错。它众所周知且被广泛使用。但它没有步长控制，更不用说密集输出了。如果您想使用具有自适应积分 ODE 的受控步进器，则需要一个受控步进器。例如，您可以使用 Runge Kutta Cash Karp 方法

typedef controlled_runge_kutta< runge_kutta_cash_karp54< state_type > > stepper_type;
integrate_adaptive( stepper_type() , lorenz , x , 0.0 , 10.0 , dt , streaming_observer( cout ) );

runge_kutta_cash_karp54< state_type > 是一个简单的误差步进器，它估计单步产生的误差。此步进器在 controlled_error_stepper 中使用，该步进器基于 Runge-Kutta Cash Karp 步进器估计的误差实现步长控制。integrate_adaptive 执行自适应积分。这会导致时间演化过程中步长发生变化。

另一个误差步进器是所谓的 Dormand-Prince 5(4) 步进器，它也具有密集输出功能。要使用它，您需要将控制器嵌套在一个通用的密集输出步进器中

typedef controlled_runge_kutta< runge_kutta_dopri5< state_type > > dopri_stepper_type;
typedef dense_output_runge_kutta< dopri_stepper_type > dense_stepper_type;
integrate_const( dense_stepper_type() , lorenz , x , 0.0 , 10.0 , dt , streaming_observer( std::cout ) );

现在，我们再次使用了 integrate_const，它利用了密集输出步进器的所有功能。原则上，也可以在 integrate_const 中使用前一个示例中的受控步进器。但在这种情况下，自适应积分会精确到时间点 0、dt、2dt、3dt、...。对于小的 dt，会损失大量性能。

当然，您也可以调整步长控制的精度，即决定步长过大或过小的度量。对于 Runge-Kutta Cash-Karp 步进器的示例，精度通过以下方式调整：

double abs_error = 1.0e-10 , rel_error = 1.0e-10;
integrate_adaptive( stepper_type( default_error_checker< double >( abs_error , rel_error ) ) , lorenz , x , 0.0 , 10.0 , dt , streaming_observer( std::cout ) );

abs_error 确定绝对误差的容差，而 rel_error 是相对误差的容差。

工厂函数

为了简化受控步进器或密集输出步进器的创建，引入了工厂函数。例如，可以使用以下方法轻松创建 Dormand-Prince 5(4) 步进器的密集输出步进器：

integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) , system , x , t1 , t2 , dt );

make_dense_output 的两个参数是绝对和相对误差容差。对于受控步进器，工厂函数称为 make_controlled，并且可以使用与 make_dense_output 完全相同的方式使用。

在 GPU 上求解常微分方程

第二个示例展示了如何使用 odeint 在 GPU 上求解常微分方程。当然，并非所有问题都适合移植到 GPU。例如，上一节的洛伦兹系统不是一个好的选择，它太小了。使用 GPU 仅对非常大的系统有意义。在 GPU 上的一次调用速度很慢，但这次调用可以并行执行大量操作。适合 GPU 的问题示例包括大型网格，如偏微分方程 (PDE) 的离散化或 ODE 的集合。

在示例中，我们将展示如何积分整个未耦合的洛伦兹系统集合，其中一个参数被改变。此示例可用于参数研究。在官方文档的教程中，对同一示例进行了更详细的解释，这里我们只展示如何使用 GPU。

目前，odeint 支持通过 Thrust 的 CUDA。原则上，也可以使用 OpenCL，但此方法尚未实现。Thrust 是用于原生 CUDA API 的类似 STL 的接口。它遵循函数式编程范例。其标准用法是调用一个通用算法，如 thrust::for_each，并使用一个执行基本操作的特定函数对象。例如：

struct add_two
{
    template< class T >
    __host__ __device__
    void operator()( T &t ) const
    {
        t += T( 2 );
    }
};

// ...

thrust::device_vector< int > x( N );
thrust::for_each( x.begin() , x.end() , add_two() );

该函数通用地将 2 加到容器 x 的每个元素。Thrust 的另一个非常好的特性是代码无需修改即可与 OpenMP 一起使用，这意味着您可以将计算分布在 CPU 的内核上。要使用此功能，您只需使用 -Xcompiler -fopenmp -DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP 编译代码。

.

要将 Thrust 与 odeint 结合使用，您需要：

• 将 thrust::device_vector< double > 用作状态类型，
• 使用专门的代数和专门的运算，即 thrust_algebra、thrust_operations（稍等片刻，看看它是如何工作的），
• 使您的系统函数为 Thrust 就绪。

前两点可以通过定义正确的状态类型和正确的步进器轻松解决

// change to float if your GPU does not support doubles
typedef double value_type;
typedef thrust::device_vector< double > state_type;
typedef runge_kutta4< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper_type;

第三点是最困难的，尽管 Thrust 使这项任务相对容易。我们想积分 N 个洛伦兹系统的一个集合，其中每个子系统都有不同的参数 R。单个洛伦兹系统的维度为三，因此状态类型的维度为 3*N。我们将使用一种内存布局，其中状态向量的前 N 个条目是 x 分量，第二个 N 个条目是 y 分量，第三个 N 个条目是 z 分量。此外，我们想改变洛伦兹系统中的参数 beta。因此，我们还需要一个具有 N 个条目的参数向量。我们的洛伦兹类的布局是：

lorenz_system
{
    // ...

    lorenz_system( size_t N , const state_type &beta )
    : m_N( N ) , m_beta( beta ) { }

    void operator()( const state_type &x , state_type &dxdt , double t )  const
    {
         // ...
    }

    size_t m_N;
    const state_type &m_beta;
};

现在，我们已准备好计算 dxdt。为此，我们使用 Thrust 的一个非常好的特性——zip 迭代器。Zip 迭代器允许我们一次迭代一堆向量，并对这堆向量执行一个操作。

lorenz_system
{
    // ...

    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {
        thrust::for_each(
                thrust::make_zip_iterator( thrust::make_tuple(
                        x.begin() , x.begin() + m_N , x.begin() + 2 * m_N ,
                        m_beta.begin() ,
                        dxdt.begin() , dxdt.begin() + m_N , dxdt.begin() + 2 * m_N ) ) ,
                thrust::make_zip_iterator( thrust::make_tuple(
                        x.begin() + m_N , x.begin() + 2 * m_N , x.begin() + 3 * m_N ,
                        m_beta.begin() ,
                        dxdt.begin() + m_N , dxdt.begin() + 2 * m_N , dxdt.begin() + 3 * m_N  ) ) ,
                lorenz_functor() );
    }

    // ...
};

因此，我们对 x、y、z、R、dxdt、dydt、dzdt 进行迭代。值被放入一个 thrust::tuple 中。在迭代过程中，使用元组 x、y、z、R、dxdt、dydt、dzdt 调用 lorenz_functor()。我们所需要了解的关于单个系统的一切都包含在这个元组中，我们可以计算洛伦兹系统的右侧（r.h.s.）：

struct lorenz_system
{
    struct lorenz_functor
    {
        template< class T >
        __host__ __device__
        void operator()( T t ) const
        {
            // unpack the parameter we want to vary and the Lorenz variables
            value_type R = thrust::get< 3 >( t );
            value_type x = thrust::get< 0 >( t );
            value_type y = thrust::get< 1 >( t );
            value_type z = thrust::get< 2 >( t );
            thrust::get< 4 >( t ) = sigma * ( y - x );
            thrust::get< 5 >( t ) = R * x - y - x * z;
            thrust::get< 6 >( t ) = -b * z + x * y ;

        }
    };

    // ...
};

就是这样。现在，我们初始化状态和 beta 并执行积分：

// initialize beta
vector< value_type > beta_host( N );
const value_type beta_min = 0.0 , beta_max = 56.0;
for( size_t i=0 ; i < N ; ++i )
    beta_host[i] = beta_min + value_type( i ) * ( beta_max - beta_min ) / value_type( N - 1 );
state_type beta = beta_host;

state_type x( 3 * N );

// initialize x,y,z
thrust::fill( x.begin() , x.begin() + 3 * N , 10.0 );

lorenz_system lorenz( N , beta );

integrate_const( stepper_type() , lorenz , x , 0.0 , 10.0 , dt );

当然，对于实际应用来说，这个示例并没有多大意义，因为这个洛伦兹系统集合没有任何事情发生。在 odeint 的教程部分，使用相同的示例计算每个系统的最大李雅普诺夫指数。使用李雅普诺夫指数可以识别参数空间中的不同区域。对于洛伦兹系统，可以识别出混沌、规则（振荡）和消失的动力学区域。使用 Thrust 和 CUDA 进行参数研究的方法非常高效，因为许多系统是并行求解的。与 CPU 相比，性能可以轻松提高 20 倍甚至更多。

设计和结构

在上一节中，我们看到了 odeint 如何用于求解常微分方程的两个示例。本节将介绍 odeint 的技术细节。此外，它还概述了 odeint 中现有的方法和步进器。

odeint 原则上由四部分组成：

• 集成函数
• 步进器
• 代数和运算
• 实用程序

集成函数负责执行 ODE 的迭代。它们进行步长控制，并可能利用密集输出。集成函数附带一个观察者概念，该概念允许您在迭代过程中观察您的解。步进器是 odeint 的主要构建块。此处实现了各种方法。存在多种用于不同目的的步进器，例如：

• 经典龙格-库塔步进器，
• 用于哈密顿系统的辛步进器，
• 隐式方法，以及
• 刚性求解器。

代数和运算构建了一个抽象层，允许您更改 odeint 执行基本数值操作（如加法、乘法等）的方式。在实用工具部分，您会找到用于调整和复制状态类型等功能。接下来，我们将更详细地介绍各个部分。

步进器概念

odeint 中的步进器执行单步。odeint 中存在多种不同的步进器类型。最简单的一种由 Stepper 概念描述。它只有一个 do_step 方法，用于执行步进。例如，四阶经典龙格-库塔步进器。

runge_kutta4< state_type > rk4;
rk4.do_step( system , x , t , dt );

大多数步进器都有一组模板参数来调整其行为。在上面的龙格-库塔 4 示例中，您必须明确指定状态类型。此类型也用于存储中间值。对于同一个步进器，还有许多其他参数，它们带有某些默认值。例如，您可以明确指定值类型（在计算过程中使用的数值）、导数类型、时间类型、代数和运算。此外，还有一个用于调整大小的策略参数。在大多数情况下，默认参数都是不错的选择。对于使用 Boost.Units 或奇异代数等特殊应用，您需要显式指定它们。odeint 的文档提供了一些示例。

do_step 方法的第一个参数指定 ODE。期望系统是函数或函数对象，并且必须是具有以下签名的可调用对象：

system( x , dxdt , t )

另一个概念是 ErrorStepper。它的基本用法与 Stepper 概念非常相似：

runge_kutta54_ck< state_type > rk54;
rk54.do_step( system , x , t , dt , xerr );

主要区别在于它估计单步产生的误差并将此误差存储在 xerr 中。另一个概念是 ControlledStepper 概念，它尝试执行单步。ControlledStepper 概念的一个例子是 controlled_runge_kutta，它有一个模板参数，指定底层的龙格-库塔步进器。例如，可以使用 Runge-Kutta-Cask-Karp 步进器：

controlled_runge_kutta< runge_kutta54_ck< state_type > > stepper;
controlled_step_result res = stepper.try_step( system , x , t , dt );

实现 ControlledStepper 概念的步进器通常具有一些必须满足的误差界限。如果调用 try_step，步进器将以给定的步长 dt 执行一步，并检查该步产生的误差。如果误差在所需的容差范围内，它则接受该步，并可能为下一次评估增加步长。如果达到误差界限，则拒绝该步并减小步长。

第四个步进器概念是 DenseOutputStepper。密集输出步进器有两个基本方法：do_step 和 calc_state。第一个方法执行一步。ODE 的状态在内部管理，并且通常通过受控步进器来执行自适应步进。第二个函数计算解的中间值。通常，中间值的精度与解的阶数相同。例如，Dormand-Prince 步进器，可以使用如下方式：

typedef controlled_runge_kutta< runge_kutta_dopri5< state_type > > stepper_type;
dense_output_runge_kutta< stepper_type > stepper;
stepper.initialize( x0 , t0 , dt0 );
std::pair< double , double > time_interval = stepper.do_step( system );
stepper.calc_state( ( time_interval.first + time_interval.second ) / 2.0 , x );

集成函数

集成函数负责使用指定的步进器集成 ODE。因此，它们负责调用 do_step 或 try_step 方法。如果您使用受控步进器或密集输出步进器，这尤其有用，因为在这种情况下，您需要一些调用步进方法的逻辑。当然，您也可以始终手动调用步进方法。如果您使用的是没有步长控制和密集输出功能的经典步进器，这可能比调用集成方法更简单。

odeint 定义了四种不同的集成方法：

• integrate_const
• integrate_adaptive
• integrate_times
• integrate_n_steps

integrate_const 以恒定步长执行积分。在每一步都可以调用一个观察者。例如：

struct streaming_observer
{
    std::ostream &m_out;
    streaming_observer( std::ostream &out ) : m_out( out ) { };

    void operator()( const state_type &x , double t ) const
    {
         m_out << t;
         for( size_t i=0 ; i < x.size() ; ++i )
             m_out << "\t" << x[i];
    }
};

integrate_const( stepper , system , x , t_start , t_end , dt , streaming_observer( cout ) );

integrate_const 实现为使用步进器的所有功能。如果您使用密集输出步进器调用它，它将执行步进并借助密集输出功能计算状态。对于受控步进器也是如此。在这里，步长在每一步都会被自适应调整。

其他集成函数与 integrate_const 函数类似。例如，integrate_adaptive 执行 ODE 的自适应积分。integrate_times 期望一个时间范围，您希望在这些时间点获得解，而 integrate_n_steps 则精确地积分 n 步。这些方法仅适用于 Stepper 概念。

代数和运算

代数和运算引入了一个自由度，允许您更改和调整 odeint 执行基本数值计算的方式。一般来说，代数定义了您如何迭代状态类型集，而运算则对状态类型中的条目执行加法、减法、乘法等数值运算。当然，可用的代数和运算很大程度上取决于状态类型。

并非所有步进器都提供更改代数和运算的可能性。例如，Rosenbrock 步进器和隐式欧拉步进器则不提供。它们需要执行矩阵计算，因此依赖于提供基本向量和矩阵算法的 boost::numeric::ublas::vector 和 boost::numeric::ublas::matrix。但是，龙格-库塔步进器以及辛步进器和多步方法都使用代数和运算。

每个代数都由一组成员函数组成：for_each1( c1 , op )、for_each2( c1 , c2 , op )、for_each3( c1 , c2 , c3 , op )... 和 accumulate( c1 , op )。变量 c1、c2、c3... 代表 ODE 的状态或其一阶导数。它们通常与步进器的 state_type 相同，但不必如此。op 是在状态元素上执行的函数或函数对象。在大多数情况下，此函数对象取自运算。运算类应具有适当的成员结构 scale_sum1、scale_sum2、scale_sum3...。每个结构都拥有一个执行运算的 operator()。有关代数和运算的更多详细信息，请参阅 odeint 的文档或阅读源代码。

代数和运算作为额外的模板参数进入步进器。例如，所有龙格-库塔步进器都具有以下参数：

runge_kutta_stepper< state_type , value_type , deriv_type , time_type , algebra , operations >

在此，state_type 表示 ODE 的状态，value_type 是基本数值类型，如 double 或 complex<double>。deriv_type 表示状态类型的导数，即主方程的左侧，而 time_type 是因变量的类型。在大多数情况下，deriv_type 与 state_type 相同，time_type 与 value_type 相同。接下来的两个参数是代数和运算。除 state_type 外，所有类型都有默认值，选择这些默认值是为了适合大多数情况。

一个所有这些类型都不同的示例是 Boost.Units，它可用于实现方程维度的编译时验证，请参阅 odeint 文档中的教程。代数和运算与默认值不同的另一个示例是上面 Thrust 的示例，它需要特殊的 thrust_algebra 以及特殊的 thrust_operations。还有其他代数，例如用于 Intel math kernel library 的代数，或者您可以轻松实现自己的代数来处理自己的容器。

龙格-库塔元算法

对于显式龙格-库塔步进器的实现，使用了一种新的模板元编程方法。使用此方法，龙格-库塔算法仅由其系数指定。当然，也可以使用经典的编程技术来完成此任务，但这无疑会牺牲性能。模板元编程算法避免了这种缺点。这里只使用编译时容器和算法，以便编译器能够执行强大的优化。这些算法的性能与其手工编写的版本相当，请参阅下一节。完整的算法在此处（仅德语） [5] 和此处 [7] 有详细描述。

性能

odeint 的性能已与 GNU Scientific library (gsl) 或 Numerical Recipes (NR) 的方法等现有解决方案进行了仔细比较。结果表明，odeint 中引入的抽象不会影响其性能，前提是使用了带有优化的现代 C++ 编译器。下面图表中显示了洛伦兹系统和 boost::array 作为状态类型的性能结果。百分比值与 odeint 的第一个版本有关。odeint 代表 odeint 的第一个版本，generic 代表 odeint v2 中的当前实现，NR 代表 Numerical Recipes 的方法，gsl 代表 GNU Scientific library，rt_gen 代表所考虑方法的通用运行时实现。

总结和展望

本文介绍了 odeint v2——一个用于求解常微分方程的高级 C++ 库。与第一个版本相比，它的主要增强功能是使用了代数和运算，这两个概念允许您更改 odeint 中基本运算的执行方式。利用这些技术，您可以在现代图形卡上求解 ODE，或者您可以使用编译时容器和高级维度分析类，如 Boost.Fusion 和 Boost.Units。此外，还实现了更多的数值方法并引入了密集输出功能。

下表显示了 odeint 中当前实现的所有步进器：

方法	类名	注释
显式欧拉	euler	非常简单，仅用于演示目的
龙格-库塔 4	runge_kutta4	经典的龙格-库塔方案，良好的通用方案，无误差控制
Cash-Karp	runge_kutta_cash_karp54	良好的通用方案，具有误差估计
Dormand-Prince 5	runge_kutta_dopri5	标准的带误差控制和密集输出的方法
Fehlberg 78	runge_kutta_fehlberg78	良好的高阶方法，具有误差估计
Adams-Bashforth-Moulton	adams_bashforth_moulton	高性能多步方法
受控误差步进器	controlled_runge_kutta	龙格-库塔步进器的误差控制
密集输出步进器	dense_output_runge_kutta	龙格-库塔步进器的密集输出
Bulirsch-Stoer	bulirsch_stoer	具有步长、阶数控制和密集输出的步进器。对于需要高精度的情况非常有用。
隐式欧拉	implicit_euler	基本隐式例程
Rosenbrock 4	rosenbrock4	带误差控制和密集输出的刚性系统求解器
辛欧拉	symplectic_euler	用于可分离哈密顿系统的基本辛求解器
辛 RKN McLachlan	symplectic_rkn_sb3a_mclachlan	用于六阶可分离哈密顿系统的辛求解器

odeint 目前正在大力开发中。因此，某些类名和头文件可能仍会更改，并且可能会实现新功能和方法。本文提供的源代码是当前开发的一个快照。

如果您想为该库做贡献，或有任何有趣的观点、批评或建议，我们诚挚地邀请您。

参考文献

1. [1] - Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., 2009.
2. [2] - Ernst Hairer, and Gerhard Wanner, (Author) Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., 2010
3. [3] - W. H. Press, S. T. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C: the Art of Scientific Computing. 1992.
4. [4] - odeint.aspx
5. [5] - 数值分析和模板元编程（仅德语）
6. [6] - Karsten Ahnert, Mario Mulansky, Odeint - Solving ordinary differential equations in C++, IP Conf. Proc., Volume 1389, pp. 1586-1589, 2011, Arxiv 链接
7. [7] - Mario Mulansky, Karsten Ahnert, Metaprogramming Applied to Numerical Problems, IP Conf. Proc., Volume 1389, pp. 1582-1585, 2011, Arxiv 链接
8. [8] - A. Alexandrescu. Modern C++ Design. 2001.

历史

• 2011年10月19日 - 初始文档。