одеинт также может быть использован для решения обычных дифференциальных уравнений, определенных на решетках. Ярким примером является система Ферми-Паста-Улам.. Это гамильтоновская система нелинейных сопряженных гармонических осцилляторов. Гамильтониан является

Примечательно, что система Ферми-Паста-Улама была первым численным экспериментом, реализованным на компьютере. Она была изучена в Лос-Аламосе в 1953 году на одном из первых компьютеров (MANIAC I), и это вызвало совершенно новое дерево математической и физической науки.

Как иСолнечная система, FPU снова решается симплектическим решателем, но в этом случае мы можем ускорить вычисления, потому чтоqкомпоненты тривиально уменьшаются доdq_{& #8203;i}/ dt = p_{& #8203;i}. Все, что вам нужно сделать, это позвонить симплектическому решателю с функцией состояния для компонентовp. Вот как выглядит эта функция

typedef vector< double > container_type;
struct fpu
{
    const double m_beta;
    fpu( const double beta = 1.0 ) : m_beta( beta ) { }
    // system function defining the ODE
    void operator()( const container_type &q , container_type &dpdt ) const
    {
        size_t n = q.size();
        double tmp = q[0] - 0.0;
        double tmp2 = tmp + m_beta * tmp * tmp * tmp;
        dpdt[0] = -tmp2;
        for( size_t i=0 ; i<n-1 ; ++i )
        {
            tmp = q[i+1] - q[i];
            tmp2 = tmp + m_beta * tmp * tmp * tmp;
            dpdt[i] += tmp2;
            dpdt[i+1] = -tmp2;
        }
        tmp = - q[n-1];
        tmp2 = tmp + m_beta * tmp * tmp * tmp;
        dpdt[n-1] += tmp2;
    }
    // calculates the energy of the system
    double energy( const container_type &q , const container_type &p ) const
    {
        // ...
    }
    // calculates the local energy of the system
    void local_energy( const container_type &q , const container_type &p , container_type &e ) const
    {
        // ...
    }
};

Вы также можете использовать<boost::array<double,N>>.

Теперь вы должны определить свои начальные значения и выполнить интеграцию:

const size_t n = 64;
container_type q( n , 0.0 ) , p( n , 0.0 );
for( size_t i=0 ; i<n ; ++i )
{
    p[i] = 0.0;
    q[i] = 32.0 * sin( double( i + 1 ) / double( n + 1 ) * M_PI );
}
const double dt = 0.1;
typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
fpu fpu_instance( 8.0 );
integrate_const( stepper_type() , fpu_instance ,
        make_pair( boost::ref( q ) , boost::ref( p ) ) ,
        0.0 , 1000.0 , dt , streaming_observer( cout , fpu_instance , 10 ) );

Наблюдатель использует ссылку на объект системы для расчета локальных энергий:

struct streaming_observer
{
    std::ostream& m_out;
    const fpu &m_fpu;
    size_t m_write_every;
    size_t m_count;
    streaming_observer( std::ostream &out , const fpu &f , size_t write_every = 100 )
    : m_out( out ) , m_fpu( f ) , m_write_every( write_every ) , m_count( 0 ) { }
    template< class State >
    void operator()( const State &x , double t )
    {
        if( ( m_count % m_write_every ) == 0 )
        {
            container_type &q = x.first;
            container_type &p = x.second;
            container_type energy( q.size() );
            m_fpu.local_energy( q , p , energy );
            for( size_t i=0 ; i<q.size() ; ++i )
            {
                m_out << t << "\t" << i << "\t" << q[i] << "\t" << p[i] << "\t" << energy[i] << "\n";
            }
            m_out << "\n";
            clog << t << "\t" << accumulate( energy.begin() , energy.end() , 0.0 ) << "\n";
        }
        ++m_count;
    }
};