在 C++ 中求解常微分方程。
4.82/5 (32投票s)
本文介绍了一个基于模板元编程的常微分方程求解框架。
注意： 存在一个新版本的 odeint，此处有详细介绍。
引言
本文介绍了一个用于求解常微分方程 (ODEs) 的 C++ 框架 odeint，该框架基于模板元编程。我认为这个框架比现有的 ODE 代码具有一些不错的优势，并且它以一种非常优雅的方式使用了模板。此外，odeint 具有类似 STL 的语法，并且独立于特定容器。整个项目位于 boost sandbox 中；本文提供的代码或多或少是当前开发的一个快照。
背景
我将只非常简要地介绍常微分方程。如果大家对更详细的 ODE 解释感兴趣，我可以在文章的后续版本中扩展这部分内容。
常微分方程
常微分方程用一个量的导数来描述该量的演化。它通常采用以下形式：
d x(t) / d t = f( x(t) , t )
函数 f 定义了 ODE，而 x 和 f 可以是向量。与每个 ODE 相关联的是一个初值问题 (IVP)，即 ODE 本身，以及一个初值 x(t0)=x0。初值问题的解是 x(t) 的时间演化，并附加条件 x(t0)=x0。可以证明，每个 IVP 都有唯一解。当然，常微分方程不限于时间问题，因此变量 t 可以被替换为其他量，例如空间坐标。ODE 及其相关的 PDE（偏微分方程）在几乎所有科学领域都非常重要。物理学中的所有主要方程都属于此类，例如经典物理学的牛顿定律、电磁学的麦克斯韦方程组、量子世界的薛定谔方程及其相对论推广，或者广义相对论的爱因斯坦方程。非常频繁地，PDE 的空间轴会被离散化，从而产生格点上的 ODE。
在下图所示的例子中，展示了一个来自混沌理论的 ODE：著名的洛伦兹吸引子。它是确定性混沌系统的第一个例子，并引发了大量科学研究。洛伦兹吸引子由一组耦合的常微分方程描述：
d x1 / dt = sigma * ( x2 - x1 )
d x2 / dt = R * x1 - x2 - x1 * x3
d x3 / dt = x1 * x2 - b * x3
用
sigma = 10
R = 28
b = 8 / 3
该方程没有解析解，因此只能通过数值方法求解。
ODE 的数值解
为了数值求解 ODE，存在各种方法；所有这些方法都会离散化时间。最简单的方法无疑是显式欧拉法，它将导数表示为差商：
d x(t) / d t = x(t+dt) - x(t) / dt
然后，一些基本的代数运算得到：
x(t+dt) = x(t) + dt*f (x(t) , t )
这样，就可以在知道时间 t 处的 x 的情况下获得 x(t+dt)。再次注意，x 不必是标量，可以是向量。其余的非常简单。选择一个初值 x(t=0)=x0，然后迭代地将欧拉步应用于 x(t) 以获得 x(t) 的演化。当然，存在更高级的求解器，其中最常用的是四阶龙格-库塔方法。
在下面的源代码中，定义 ODE 的函数 f(x(t)) 被同义地称为系统或动力系统。
主要思想
现在，我们来讨论库的设计方面。显然需要两个主要构建块：
- 一个步进函数，它从 x(t) 的知识计算出 x(t+d t)。
- 一个积分函数，它执行迭代，从而多次应用步进函数。
整个库位于命名空间 boost::odeint::numeric 中。.
步进方法
我们试图避免使用 抽象 类，因为这会产生一些额外的性能开销。相反，我们使用泛型编程和概念来描述类型或类应该如何表现，以及其方法应该可用。对于步进函数，我们需要以下要素：
template<
class Container ,
class Time = double ,
class Traits = container_traits< Container >
>
class ode_step
{
// provide basic typedefs
//
public:
typedef unsigned short order_type;
typedef Time time_type;
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
// public interface
public:
// return the order of the stepper
order_type order_step() const;
// standard constructor
ode_step( void );
// adjust the size of the internal arrays
void adjust_size( const container_type &x );
// performs one step
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
const container_type &dxdt ,
time_type t ,
time_type dt );
// performs one step
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
time_type t ,
time_type dt );
};
容器类型定义了状态的存储方式。它必须满足一些基本要求，例如值类型，并且应该是默认可构造的。container_traits 类型抽象了 序列概念 的 resize 功能。此外，它知道如何获取适当的 begin() 和 end() 迭代器。其用法将在下面展示。adust_size 方法负责调整用于存储中间计算结果的任何内部容器的大小。
do_step 函数从 x(t) 的知识计算 x(t+dt)，并就地转换 x(t)，这意味着 x(t) 在 do_step 中被 x(t+dt) 替换。参数 system 定义了 ODE。这可以是一个带有三个参数的简单函数：system( state , dxdt , t )，或者一个带有适当括号运算符的类 operator()(state,dxdt,t)。这样的类也称为函数对象。存在两个版本的 do_step：一个在显式知道 f(x(t)) 的情况下计算 x(t+dt)，另一个在内部计算该值。原则上，第二个版本只调用第一个版本。这两个版本都存在，因为在许多情况下，用户需要知道 f(x(t)) 并会调用带有 x(t) 的 system 方法。为了避免重复调用 system，引入了带有 f(x(t)) 显式声明的 do_step 版本。
注意，上面定义的类并不存在于代码中；它仅用于展示步进功能。所有拥有这些 typedef 和方法的特定类都被称为步进类；我们也说所有具有这些定义的类都满足步进概念。
现在，我们来看一个具体的方法，欧拉求解器：
template<
class Container ,
class Time = double ,
class Traits = container_traits< Container >
>
class stepper_euler
{
// provide basic typedefs
public:
typedef unsigned short order_type;
typedef Time time_type;
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
// private members
private:
container_type m_dxdt;
// public interface
public:
// return the order of the stepper
order_type order_step() const { return 1; }
// standard constructor, m_dxdt is not adjusted
stepper_euler( void )
{
}
// contructor, which adjusts m_dxdt
stepper_euler( const container_type &x )
{
adjust_size( x );
}
// adjust the size of m_dxdt
void adjust_size( const container_type &x )
{
traits_type::adjust_size( x , m_dxdt );
}
// performs one step with the knowledge of dxdt(t)
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
const container_type &dxdt ,
time_type t ,
time_type dt )
{
detail::it_algebra::increment( traits_type::begin(x) ,
traits_type::end(x) ,
traits_type::begin(dxdt) ,
dt );
}
// performs one step
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
time_type t ,
time_type dt )
{
system( x , m_dxdt , t );
do_step( system , x , m_dxdt , t , dt );
}
};
这个定义很简单，只有几点说明：
大小的调整在 container_traits 类型中完成。标准 traits 类型定义为：
template< class Container >
struct container_traits
{
typedef Container container_type;
typedef typename container_type::value_type value_type;
typedef typename container_type::iterator iterator;
typedef typename container_type::const_iterator const_iterator;
static void resize( const container_type &x , container_type &dxdt )
{ dxdt.resize( x.size() ); }
static bool same_size( const container_type &x1 , const container_type &x2 )
{ return (x1.size() == x2.size()); }
static void adjust_size( const container_type &x1 , container_type &x2 )
{ if( !same_size( x1 , x2 ) ) resize( x1 , x2 ); }
static iterator begin( container_type &x ) { return x.begin(); }
static const_iterator begin( const container_type &x ) { return x.begin(); }
static iterator end( container_type &x ) { return x.end(); }
static const_iterator end( const container_type &x ) { return x.end(); }
};
这将适用于所有满足 序列 和 容器 概念的（STL）容器。但对于其他容器，这些要求并不充分。例如，std::tr1::array 没有 resize 方法。容器 traits 也负责 begin() 和 end() 迭代器。对于 SGI 容器 来说这很简单，但对于像 Blitz++ 或 MTL 矩阵这样的矩阵类型，begin() 和 end() 没有定义，但可以通过不同方式获得。
所有迭代器操作都被放在它们自己的函数中。这样做的好处是代码小巧且清晰，因为在步进函数中不需要定义迭代器。此外，还可以为特定的迭代器或容器专门化和优化这些迭代器函数。
如果系统大小在演化过程中发生变化，则需要 adjust size 功能；例如，考虑一个定义在具有非固定节点数的网络上的 ODE。
这并不是步进函数的全部内容。odeint 还提供了一个高级步进概念，该概念在每一步计算数值误差。通过这种误差计算，可以在迭代的每一步中自适应地调整步长 dt。自适应步长积分背后的概念不会在本文中进行描述。
以下代码演示了步进函数的使用：
#include <iostream>
#include <boost/numeric/odeint.hpp>
#define tab "\t"
using namespace std;
using namespace boost::numeric::odeint;
typedef vector< double > state_type;
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
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];
}
int main( int argc , char **argv )
{
const double dt = 0.01;
state_type x(3);
x[0] = 1.0 ;
x[1] = 0.0 ;
x[2] = 0.0;
stepper_euler< state_type > stepper;
stepper.adjust_size( x );
double t = 0.0;
for( size_t oi=0 ; oi<10000 ; ++oi,t+=dt )
{
stepper.do_step( lorenz , x , t , dt );
cout << x[0] << tab << x[1] << tab << x[2] << endl;
}
}
在某些情况下，希望将 ODE 定义为一个类。在我们的例子中，lorenz 也可以定义为一个类，它具有适当的括号运算符：
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];
}
}
然后通过以下方式调用：
lorenz_class lorenzo;
stepper.do_step( lorenzo , x , t , dt );
通过类可以创建更复杂的动力系统 (ODE)。
积分函数
下一个主要思想是提供一个或多个函数来自动迭代 ODE 的状态，并可以在每一步观察状态。这是通过 integrate_const 函数实现的：
template<
class Stepper ,
class DynamicalSystem ,
class Observer
>
size_t integrate_const(
Stepper &stepper ,
DynamicalSystem &system ,
typename Stepper::container_type &state ,
typename Stepper::time_type start_time ,
typename Stepper::time_type end_time ,
typename Stepper::time_type dt ,
Observer &observer
)
{
stepper.adjust_size( state );
size_t iteration = 0;
while( start_time < end_time )
{
observer( start_time , state , system );
stepper.do_step( system , state , start_time , dt );
start_time += dt;
++iteration;
}
observer( start_time , state , system );
return iteration;
}
该函数以恒定步长 dt 在时间间隔 start_time、end_time 内求解 ODE。还必须提供初始状态。此外，还引入了观察者（observer）的概念，尽管这并非严格意义上的观察者模式。同样，观察者可以是任何类型，它只需要拥有括号运算符 operator()( time , state , system )。例如，可以使用来自 boost 的 lambda 表达式：
integrate_const( stepper , lorenz , x , 0.0 , 10.0 , dt ,
cout << _1 << tab << _2[0] << "\n" );
或
vector< double > lx;
back_insert_iterator< vector<double> > iter( lx );
integrate_const( stepper , lorenz , x , 0.0 , 10.0 , dt , var(*iter++) = _2[0] );
安装
整个库是头文件形式的，意味着不需要进行特殊的安装过程。在类 Unix 系统上，请执行以下操作：
- 解压源代码：
unzip odeint.zip - 进入 odeint 目录：
cd odeint - 设置 boost 环境变量：
export $BOOST_ROOT=path_to_boost - 编译 lorenz 示例：
make
现在，应该会有一个名为 lorenz 的可执行文件。lambda 表达式示例需要 boost 源代码。您可以从 boost.org 下载它们，并且环境变量应该指向 boost 的根目录，在该目录中您将找到子目录 bin.v2、boost、doc、lib 等。
对于 MSVC，odeint 也能工作，但目前我还没有关于如何编译示例的逐步指南。
总结与展望
我向您展示了 odeint 的一些设计方面和用法——一个用于求解常微分方程的 C++ 框架。它非常易于使用且非常灵活。性能也非常出色。我们进行了一些测试运行，将 odeint 与 gsl 进行了比较。作为测试系统，使用了洛伦兹系统。结果表明，odeint 比 gsl 例程快约两倍。
本文提供的源代码是当前开发的一个快照。存在一些用于自适应步长控制的方法，以及一些用于求解哈密顿系统的方程。下表概述了现有方法。S 表示步进器满足上面描述的简单步进器概念，而 E 表示步进器满足此处未介绍但用于自适应步长控制的 Error 步进器概念。
| 方法 | 类名 | Order | 误差（阶） | 步进器概念 |
| 欧拉 | stepper_euler |
1 | 否 | S |
| 龙格-库塔 4 阶 | stepper_rk4 |
4 | 否 | S |
| 龙格-库塔-Cash-Karp | stepper_rk5_ck |
5 | 是 4 | E |
| 龙格-库塔-Fehlberg | stepper_rk78_fehlberg |
7 | 是 8 | S,E |
| 中点法 | stepper_midpoint |
变分。 | 否 | S |
| 辛欧拉法 | hamiltonian_stepper_euler |
1 | 否 | S |
| 辛龙格-库塔法 | hamiltonian_stepper_rk |
6 | 否 | S |
整个库尚未完成，将在不久的将来进行扩展。我们有一个小的待办事项列表，应该在最终发布前完成：
- Blitz++, MTL,
boost::ublas和boost::multi_array的适配器 - 文档
- 求解刚性系统的求解器
- 抽象迭代器代数
如果您想为这个库贡献力量，或者有任何有趣的观点、批评或建议，都热烈欢迎。
历史
- 2009 年 11 月 9 日 - 初始文档
- 2010 年 3 月 16 日 - 更新，引入容器 traits，替换示例
- 2010 年 3 月 23 日 - 更新源代码和文章