The next example in this tutorial is a simulation of the outer solar system, consisting of the sun, Jupiter, Saturn, Uranus, Neptune and Pluto.

Each planet and of course the sun will be represented by mass points. The interaction force between each object is the gravitational force which can be written as

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

where γ is the gravitational constant, m_​i and m_​j are the masses and q_​i and q_​j are the locations of the two objects. The equations of motion are then

dq_​i / dt = p_​i

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

where p_​i is the momenta of object i. The equations of motion can also be derived from the Hamiltonian

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

with the interaction potential V(q_​i,q_​j). The Hamiltonian equations give the equations of motion

dq_​i / dt = dH / dp_​i

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

In time independent Hamiltonian system the energy and the phase space volume are conserved and special integration methods have to be applied in order to ensure these conservation laws. The odeint library provides classes for separable Hamiltonian systems, which can be written in the form H = Σ p_​i² / (2m_​i) + H_​q(q), where H_​q(q) only depends on the coordinates. Alltough this functional form might look a bit arbitrary it covers nearly all classical mechanical systems with inertia and without dissipation, or where the equations of motion can be written in the form dq_​i / dt = p_​i / m_​i , dp_​i / dt = f( q_​i ).

To implement this system we define a point type which will represent the space as well as the velocity. Therefore, we use the operators from 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
    

The next step is to define a container type storing the values of q and p and to define system functions. As container type we use std::tr1::array

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

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

The container_type is different from the state type of the ODE. The state type of the ode is simply a pair< container_type , container_type > since it needs the informations about the coordinates and the momenta.

As system function we have to provide f(p) and f(q):

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;

            }
        }
    }
};

In general a three body-system is chaotic, hence we can not expect that arbitray initial conditions of the system will lead to a dynamic which is comparable with the solar system. That is we have to define proper initial conditions, which are taken from the book of Hairer, Wannier, Lubich.

As mentioned above, we need to use some special integrators in order to conserve phase space volume and energy. There is a well known family of such integrators, the so-called Runge-Kutta-Nystroem solvers, which we apply here:

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 ) );

These integration routine was used to produce the above sketch of the solar system. Note, that there are two particularities in this example. First, the state of the symplectic stepper is not container_type but a pair of container_type. Hence, we must pass such a pair to the integrate function. Since, we want to pass them as references we can simply pack them into Boost.Ref. The second point is the observer, which is called with a state type, hence a pair of container_type. The reference wrapper is also passed, but this is not a problem at all:

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";
    }
};

The full example can be found here: solar_system.cpp