Boost C++ 库

本教程的下一个示例是外太阳系的模拟，包括太阳、木星、土星、天王星、海王星和冥王星。

每个行星，当然还有太阳，都将由质点表示。每个物体之间的相互作用力是万有引力，可以写成

F_​ij = -γ m_​i m_​j ( q_​i - q_​j ) / | q_​i - q_​j | ³

其中 γ 是引力常数，m_​i 和 m_​j 是质量，q_​i 和 q_​j 是两个物体的位置。运动方程为

dq_​i / dt = p_​i

dp_​i / dt = 1 / m_​i Σ_​ji F_​ij

其中 p_​i 是物体 i 的动量。运动方程也可以从哈密顿量导出

H = Σ_​i p_​i² / ( 2 m_​i ) + Σ_​j V( q_​i , q_​j )

与相互作用势 V(q_​i,q_​j)。哈密顿方程给出运动方程

dq_​i / dt = dH / dp_​i

dp_​i / dt = -dH / dq_​i

在与时间无关的哈密顿系统中，能量和相空间体积是守恒的，必须应用特殊的积分方法以确保这些守恒定律。 odeint 库为可分离的哈密顿系统提供类，这些系统可以写成 H = Σ p_​i² / (2m_​i) + H_​q(q) 的形式，其中 H_​q(q) 仅取决于坐标。虽然这种函数形式可能看起来有点任意，但它几乎涵盖了所有具有惯性且没有耗散的经典力学系统，或者运动方程可以写成 dq_​i / dt = p_​i / m_​i , dp_​i / dt = f( q_​i ) 的形式。

	注意
	一个简短的物理学说明：虽然已知二体问题是可积的，这意味着它可以完全用解析技术解决，但三体问题已经是不可解的了。这是 19 世纪末由 H. Poincare 发现的，它导致了全新的学科 混沌理论。

为了实现这个系统，我们定义一个 3D 点类型，它将表示空间和速度。因此，我们使用来自 Boost.Operators 的运算符

/*the point type */
template< class T , size_t Dim >
class point :
    boost::additive1< point< T , Dim > ,
    boost::additive2< point< T , Dim  > , T ,
    boost::multiplicative2< point< T , Dim > , T
    > > >
    {
    public:

        const static size_t dim = Dim;
        typedef T value_type;
        typedef point< value_type , dim > point_type;

        // ...
        // constructors

        // ...
        // operators

    private:

        T m_val[dim];
    };

    //...
    // more operators

下一步是定义一个容器类型来存储 q 和 p 的值，并定义系统函数。作为容器类型，我们使用 std::array

// we simulate 5 planets and the sun
const size_t n = 6;

typedef point< double , 3 > point_type;
typedef std::array< point_type , n > container_type;
typedef std::array< double , n > mass_type;

container_type 与 ODE 的状态类型不同。 ode 的状态类型只是一个 pair< container_type , container_type >，因为它需要关于坐标和动量的信息。

接下来我们定义系统的方程。由于我们将使用一个考虑系统哈密顿（能量守恒）特性的步进器，我们必须定义与通常情况不同的 rhs，在通常情况下它只是一个函数。步进器将利用可分离的特性，这意味着系统将由两个对象定义，这两个对象分别表示 f(p) = -dH/dq 和 g(q) = dH/dp

const double gravitational_constant = 2.95912208286e-4;

struct solar_system_coor
{
    const mass_type &m_masses;

    solar_system_coor( const mass_type &masses ) : m_masses( masses ) { }

    void operator()( const container_type &p , container_type &dqdt ) const
    {
        for( size_t i=0 ; i<n ; ++i )
            dqdt[i] = p[i] / m_masses[i];
    }
};

struct solar_system_momentum
{
    const mass_type &m_masses;

    solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }

    void operator()( const container_type &q , container_type &dpdt ) const
    {
        const size_t n = q.size();
        for( size_t i=0 ; i<n ; ++i )
        {
            dpdt[i] = 0.0;
            for( size_t j=0 ; j<i ; ++j )
            {
                point_type diff = q[j] - q[i];
                double d = abs( diff );
                diff *= ( gravitational_constant * m_masses[i] * m_masses[j] / d / d / d );
                dpdt[i] += diff;
                dpdt[j] -= diff;

            }
        }
    }
};

通常，三体系统是混沌的，因此我们不能期望系统的任意初始条件都会导致与太阳系动力学相当的解。也就是说，我们必须定义适当的初始条件，这些条件取自 Hairer, Wannier, Lubich 的书 [4] 。

如上所述，我们需要使用一些特殊的积分器来守恒相空间体积。有一系列众所周知的此类积分器，即所谓的 Runge-Kutta-Nystroem 求解器，我们在此以 symplectic_rkn_sb3a_mclachlan 步进器的形式应用它们

typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
const double dt = 100.0;

integrate_const(
        stepper_type() ,
        make_pair( solar_system_coor( masses ) , solar_system_momentum( masses ) ) ,
        make_pair( boost::ref( q ) , boost::ref( p ) ) ,
        0.0 , 200000.0 , dt , streaming_observer( cout ) );

这些积分例程用于生成上述太阳系草图。请注意，此示例中有两个特殊性。首先，辛步进器的状态不是 container_type，而是一对 container_type。因此，我们必须将这样一对传递给 integrate 函数。由于我们想将它们作为引用传递，我们可以简单地将它们打包到 Boost.Ref 中。第二点是观察者，它使用状态类型调用，因此是一对 container_type。引用包装器也被传递，但这根本不是问题

struct streaming_observer
{
    std::ostream& m_out;

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

    template< class State >
    void operator()( const State &x , double t ) const
    {
        container_type &q = x.first;
        m_out << t;
        for( size_t i=0 ; i<q.size() ; ++i ) m_out << "\t" << q[i];
        m_out << "\n";
    }
};

	提示
	您可以使用 C++11 lambda 来创建观察者

完整的示例可以在这里找到： solar_system.cpp