Следующим примером в этом учебнике является моделирование внешней Солнечной системы, состоящей из Солнца, Юпитера, Сатурна, Урана, Нептуна и Плутона.

Каждая планета и, конечно же, Солнце будут представлены точками масс. Сила взаимодействия между каждым объектом является гравитационной силой, которая может быть записана как

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

гдеи #947;— гравитационная постоянная,m_{& #8203;i}им_{и #8203;j}являются массами иq_{& #8203;i}иq_{& #8203;j}— расположение двух объектов. Тогда уравнения движения

dq_​i / dt = p_​i

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

гдеp_{& #8203;i}является моментом объектаi. Уравнения движения также могут быть получены из Гамильтона.

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

с потенциалом взаимодействияV(q_{и #8203;i},q_{и #8203;j}]. Уравнения Гамильтона дают уравнения движения

dq_​i / dt = dH / dp_​i

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

В независимой гамильтоновой системе сохраняется энергия и объем фазового пространства, и для обеспечения этих законов сохранения необходимо применять специальные методы интеграции. Библиотека odeint предоставляет классы для разделяемых гамильтоновых систем, которые могут быть написаны в форме.H = & #931; p_{& #8203;i}²/ (2m_{& #8203;i}) + H_{& #8203;q}, гдеH_{& #8203;q}(q)зависит только от координат. Хотя эта функциональная форма может выглядеть немного произвольной, она охватывает почти все классические механические системы с инерцией и без рассеивания, или где уравнения движения могут быть написаны в видеdq_{& #8203;i}/ dt = p_{& #8203;i}/ m_{& #8203;i},dp_{& #8203;i}/ dt = f ( q_{& #8203;i}).

	Note
	Короткая физическая нота: Хотя известно, что проблема двух тел интегрируема, это означает, что ее можно решить с помощью чисто аналитических методов, уже проблема трех тел неразрешима. Это было найдено в конце 19-го века Х. Пуанкаре, что привело к совершенно новому вопросу. Теория хаоса.

Для реализации этой системы мы определяем тип точки 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

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 boost::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;

<container_type>отличается от государственного типа ОДЭ. Тип состояния оды просто<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;
            }
        }
    }
};

В целом три системы тела хаотичны, поэтому нельзя ожидать, что произвольные начальные условия системы приведут к решению, сравнимому с динамикой Солнечной системы. То есть мы должны определить надлежащие начальные условия, которые взяты из книги Хайера, Ванье, Любича..

Как уже упоминалось выше, для сохранения объема фазового пространства необходимо использовать специальные интеграторы. Существует хорошо известное семейство таких интеграторов, так называемых растворителей Рунге-Кутта-Нистром, которые мы применяем здесь в терминах степпера<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>. Следовательно, мы должны передать такую пару функции интеграции. Поскольку мы хотим передать их в качестве ссылок, мы можем просто упаковать их в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";
    }
};

	Tip
	Вы можете использовать лямбду C++11 для создания наблюдателей.

The full example can be found here: solar_system.cpp