Решение уравнения Лоренца с помощью РК-4 на C ++

Я спросил ОП, могут ли они использовать шаблоны или нет, они ответили "I can use templates, but I wanted it to be simple.". За кулисами предстоит еще немного поработать, чтобы сначала настроить структуры шаблонов; но как только это будет сделано; использование шаблонов действительно немного упрощает. В конце концов, это обычно делает код более общим, а не конкретным для одной задачи.

Вот нереализованный класс, показывающий обобщенный образец интегратора ODE.

#ifndef ODE_GENERALIZED_DEFINITION
#define ODE_GENERALIZED_DEFINITION

// A generalized structure of what an ODE should represent
// This class is not used; but only serves to show the interface
// of what any ODE type integrator - solver needs to contain
template<class Container, class Time = double, class Traits = container_traits<Container> >
class ode_step {
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;

    ode_step() = default;

    order_type order_step() const {};

    void adjust_size( const container_type& x ) {}

    // performs one step
    template<class DynamicSystem>
    void do_step( DynamicSystem& system, container_type& x, const container_type& dxdt, time_type t, time_type dt ) {}

    // performs one step
    template<class DynamicSystem>
    void do_step( DynamicSystem& system, container_type& x, time_type t, time_type dt ) {}
};

#endif // !ODE_GENERALIZED_DEFINITION

Вышеупомянутое ничего не делает, а только показывает, как ODE должно себя вести или действовать. Прежде чем мы сможем смоделировать какое-либо ODE по вышеприведенному шаблону, нам нужно определить несколько вещей.

Сначала нам нужен класс container_traits.

container_traits.h

#ifndef CONTAINER_TRAITS
#define CONTAINER_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();
    }
};

#endif // !CONTAINER_TRAITS

Вышеупомянутый класс довольно прост и довольно прямолинеен.

Затем нам понадобится набор функций типа итератора для итерации этих контейнеров. Функции итератора немного сложнее, но опять же довольно просто относительно их предполагаемых целей. У меня они заключены в пространство имен, поэтому они не будут конфликтовать с любыми другими классами итераторами или функциями, которые могут быть там.

iterator_algebra.h

#ifndef ITERATOR_ALGEBRA
#define ITERATOR_ALGEBRA

#include <cmath>
#include <iostream>

namespace it_algebra { // iterator algebra

// computes y += alpha * x1
template<class InOutIter, class InIter, class T>
void increment( InOutIter first1, InOutIter last1,
                InIter first2, T alpha ) {
    while( first1 != last1 )
        (*first1++) += alpha * (*first2++);
}

// computes y = alpha1 * ( x1 + x2 + alpha2*x3 )
template<class OutIter, class InIter1, class InIter2, class InIter3, class T>
void increment_sum_sum( OutIter first1, OutIter last1,
                        InIter1 first2, InIter2 first3,
                        InIter3 first4, T alpha1, T alpha2 ) {
    while( first1 != last1 )
        (*first1++) += alpha1 *
        ((*first2++) + (*first3++) + alpha2 * (*first4++));
}

// computes y = alpha1*x1 + alpha2*x2
template<class OutIter, class InIter1, class InIter2, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
                       T alpha1, InIter1 x1_begin,
                       T alpha2, InIter2 x2_begin ) {
    while( y_begin != y_end ) {
        (*y_begin++) = alpha1 * (*x1_begin++) +
            alpha2 * (*x2_begin++);
    }
}

// computes y = x1 + alpha2*x2 + alpha3*x3
template<class OutIter, class InIter1, class InIter2, class InIter3, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
                       T alpha1, InIter1 x1_begin,
                       T alpha2, InIter2 x2_begin,
                       T alpha3, InIter3 x3_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++);
}

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
                       T alpha1, InIter1 x1_begin,
                       T alpha2, InIter2 x2_begin,
                       T alpha3, InIter3 x3_begin,
                       T alpha4, InIter4 x4_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++);
}

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
template<class OutIter, class InIter1, class InIter2,
    class InIter3, class InIter4, class InIter5, class T>
    inline void scale_sum( OutIter y_begin, OutIter y_end,
                           T alpha1, InIter1 x1_begin,
                           T alpha2, InIter2 x2_begin,
                           T alpha3, InIter3 x3_begin,
                           T alpha4, InIter4 x4_begin,
                           T alpha5, InIter5 x5_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++) +
        alpha5 * (*x5_begin++);
}

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6
template<class OutIter, class InIter1, class InIter2,
    class InIter3, class InIter4, class InIter5, class InIter6, class T>
    inline void scale_sum( OutIter y_begin, OutIter y_end,
                           T alpha1, InIter1 x1_begin,
                           T alpha2, InIter2 x2_begin,
                           T alpha3, InIter3 x3_begin,
                           T alpha4, InIter4 x4_begin,
                           T alpha5, InIter5 x5_begin,
                           T alpha6, InIter6 x6_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++) +
        alpha5 * (*x5_begin++) +
        alpha6 * (*x6_begin++);
}    

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7
template<class OutIter, class InIter1, class InIter2, class InIter3,
    class InIter4, class InIter5, class InIter6, class InIter7, class T>
    inline void scale_sum( OutIter y_begin, OutIter y_end,
                           T alpha1, InIter1 x1_begin,
                           T alpha2, InIter2 x2_begin,
                           T alpha3, InIter3 x3_begin,
                           T alpha4, InIter4 x4_begin,
                           T alpha5, InIter5 x5_begin,
                           T alpha6, InIter6 x6_begin,
                           T alpha7, InIter7 x7_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++) +
        alpha5 * (*x5_begin++) +
        alpha6 * (*x6_begin++) +
        alpha7 * (*x7_begin++);
}

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7 + alpha8*x8
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4,
    class InIter5, class InIter6, class InIter7, class InIter8, class T>
    inline void scale_sum( OutIter y_begin, OutIter y_end,
                           T alpha1, InIter1 x1_begin,
                           T alpha2, InIter2 x2_begin,
                           T alpha3, InIter3 x3_begin,
                           T alpha4, InIter4 x4_begin,
                           T alpha5, InIter5 x5_begin,
                           T alpha6, InIter6 x6_begin,
                           T alpha7, InIter7 x7_begin,
                           T alpha8, InIter8 x8_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++) +
        alpha5 * (*x5_begin++) +
        alpha6 * (*x6_begin++) +
        alpha7 * (*x7_begin++) +
        alpha8 * (*x8_begin++);
}

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7 + alpha8*x8 + alpha9*x9
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4,
    class InIter5, class InIter6, class InIter7, class InIter8, class InIter9, class T>
    inline void scale_sum( OutIter y_begin, OutIter y_end,
                           T alpha1, InIter1 x1_begin,
                           T alpha2, InIter2 x2_begin,
                           T alpha3, InIter3 x3_begin,
                           T alpha4, InIter4 x4_begin,
                           T alpha5, InIter5 x5_begin,
                           T alpha6, InIter6 x6_begin,
                           T alpha7, InIter7 x7_begin,
                           T alpha8, InIter8 x8_begin,
                           T alpha9, InIter9 x9_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++) +
        alpha5 * (*x5_begin++) +
        alpha6 * (*x6_begin++) +
        alpha7 * (*x7_begin++) +
        alpha8 * (*x8_begin++) +
        alpha9 * (*x9_begin++);
}

// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7 + alpha8*x8 + alpha9*x9 + alpha10*x10
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4, class InIter5,
    class InIter6, class InIter7, class InIter8, class InIter9, class InIter10, class T>
    inline void scale_sum( OutIter y_begin, OutIter y_end,
                           T alpha1, InIter1 x1_begin,
                           T alpha2, InIter2 x2_begin,
                           T alpha3, InIter3 x3_begin,
                           T alpha4, InIter4 x4_begin,
                           T alpha5, InIter5 x5_begin,
                           T alpha6, InIter6 x6_begin,
                           T alpha7, InIter7 x7_begin,
                           T alpha8, InIter8 x8_begin,
                           T alpha9, InIter9 x9_begin,
                           T alpha10, InIter10 x10_begin ) {
    while( y_begin != y_end )
        (*y_begin++) =
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++) +
        alpha5 * (*x5_begin++) +
        alpha6 * (*x6_begin++) +
        alpha7 * (*x7_begin++) +
        alpha8 * (*x8_begin++) +
        alpha9 * (*x9_begin++) +
        alpha10 * (*x10_begin++);
}

// generic version for n values
template<class OutIter, class InIter, class InIterIter, class FactorIter, class T>
inline void scale_sum_generic( OutIter y_begin, OutIter y_end,
                               FactorIter alpha_begin, FactorIter alpha_end,
                               T beta, InIter x_begin, InIterIter x_iter_begin ) {
    FactorIter alpha_iter;
    InIterIter x_iter_iter;
    while( y_begin != y_end ) {
        x_iter_iter = x_iter_begin;
        alpha_iter = alpha_begin;
        *y_begin = *x_begin++;
        //std::clog<<(*y_begin);
        while( alpha_iter != alpha_end ) {
            //std::clog<< " + " <<beta<<" * "<<*alpha_iter<<"*"<<(*(*(x_iter_iter)));
            (*y_begin) += beta * (*alpha_iter++) * (*(*x_iter_iter++)++);
        }
        //std::clog<<" = "<<*y_begin<<std::endl;
        y_begin++;
    }
    //std::clog<<std::endl;
}

// computes y += alpha1*x1 + alpha2*x2 + alpha3*x3 + alpha4*x4
template<class OutIter, class InIter1, class InIter2,
    class InIter3, class InIter4, class T>
    inline void scale_sum_inplace( OutIter y_begin, OutIter y_end,
                                   T alpha1, InIter1 x1_begin,
                                   T alpha2, InIter2 x2_begin,
                                   T alpha3, InIter3 x3_begin,
                                   T alpha4, InIter4 x4_begin ) {
    while( y_begin != y_end )
        (*y_begin++) +=
        alpha1 * (*x1_begin++) +
        alpha2 * (*x2_begin++) +
        alpha3 * (*x3_begin++) +
        alpha4 * (*x4_begin++);
}

// calculates tmp = y, y = x1 + alpha*x2, x1 = tmp 
template<class OutIter, class InIter, class T>
inline void scale_sum_swap( OutIter y_begin, OutIter y_end,
                            OutIter x1_begin, T alpha, InIter x2_begin ) {
    T swap;
    while( y_begin != y_end ) {
        swap = (*x1_begin) + alpha * (*x2_begin++);
        *x1_begin++ = *y_begin;
        *y_begin++ = swap;
    }
}

// computes y = x1 + alpha2 * x2 ; x2 += x3
template<class OutIter, class InIter1,
    class InOutIter, class InIter2, class T>
    void assign_sum_increment( OutIter first1, OutIter last1, InIter1 first2,
                               InOutIter first3, InIter2 first4, T alpha ) {
    while( first1 != last1 ) {
        (*first1++) = (*first2++) + alpha * (*first3);
        (*first3++) += (*first4++);
    }
}

template<class OutIter, class InIter1, class InIter2, class T >
void weighted_scale( OutIter y_begin, OutIter y_end, InIter1 x1_begin, InIter2 x2_begin,
                     T eps_abs, T eps_rel, T a_x, T a_dxdt ) {
    using std::abs;

    while( y_begin != y_end ) {
        *y_begin++ = eps_abs +
        eps_rel * (a_x * abs( *x1_begin++ ) +
        a_dxdt * abs( *x2_begin++ ));
    }
}    

template<class InIter1, class InIter2, class T >
T max_ratio( InIter1 x1_begin, InIter1 x1_end, InIter2 x2_begin, T initial_max ) {
    using std::abs;
    while( x1_begin != x1_end ) {
        initial_max = std::max( static_cast<T>(abs( *x1_begin++ ) / abs( *x2_begin++ )), initial_max );
    }
    return initial_max;
}

} // namespace it_algebra

#endif // !ITERATOR_ALGEBRA

Теперь, когда у нас есть необходимые структуры и функции, мы можем использовать приведенный выше шаблон для реализации различных типов интеграторов ODE. Я покажу два самых распространенных: euler и rk4. Можно легко реализовать среднюю точку, rk5 или любой другой, если они следуют шаблону.

euler.h

#ifndef EULER_H
#define EULER_H

#include "iterator_algebra.h"
#include "container_traits.h"

template<class Container, class Time = double, class Traits = container_traits<Container> >
class euler_stepper {
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:
    container_type m_dxdt;

public:
    euler_stepper() = default;
    euler_stepper( const container_type& x ) {
        adjust_size( x );
    }

    order_type order_step() const {
        return 1;
    }

    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 DynamicSystem>
    void do_step( DynamicSystem& system, container_type& x, const container_type& dxdt, time_type t, time_type dt ) {
        it_algebra::increment( traits_type::begin( x ),
                               traits_type::end( x ),
                               traits_type::begin( dxdt ),
                               dt );
    }

    // performs one step
    template<class DynamicSystem>
    void do_step( DynamicSystem& system, container_type& x, time_type t, time_type dt ) {
        system( x, m_dxdt, t );
        do_step( system, x, m_dxdt, t, dt );
    }
};

#endif // EULER_H

rk4.h

#ifndef RK4_H
#define RK4_H

#include "iterator_algebra.h"
#include "container_traits.h"

template<class Container, class Time = double, class Traits = container_traits<Container>>
class rk4_stepper {
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;
    // typedef typename traits_type::iterator iterator;
    // typedef typename traits_type::const_iterator const_iterator;

private:
    container_type m_dxdt;
    container_type m_dxt;
    container_type m_dxm;
    container_type m_dxh;
    container_type m_xt;

public:
    rk4_stepper() = default;
    rk4_stepper( const container_type& x ) { adjust_size( x ); }

    order_type order_step() const { return 4; }

    void adjust_size( const container_type& x ) {
        traits_type::adjust_size( x, m_dxdt );
        traits_type::adjust_size( x, m_dxt );
        traits_type::adjust_size( x, m_dxm );
        traits_type::adjust_size( x, m_xt );
        traits_type::adjust_size( x, m_dxh );
    }

    template<class DynamicSystem>
    void do_step( DynamicSystem& system, container_type& x, const container_type& dxdt, time_type t, time_type dt ) {
        using namespace it_algebra;

        const time_type val1 = static_cast<time_type>(1.0);

        time_type dh = static_cast<time_type>(0.5) * dt;
        time_type th = t + dh;

        // dt * dxdt = k1
        // m_xt = x + dh*dxdt
        scale_sum( traits_type::begin( m_xt ), traits_type::end( m_xt ),
                   val1, traits_type::begin( x ), dh, traits_type::begin( dxdt ) );

        // dt * m_dxt = k2
        system( m_xt, m_dxt, th );
        // m_xt = x + dh*m_dxt
        scale_sum( traits_type::begin( m_xt ), traits_type::end( m_xt ),
                   val1, traits_type::begin( x ), dh, traits_type::begin( m_dxt ) );

        // dt * m_dxm = k3
        system( m_xt, m_dxm, th );
        // m_xt = x + dt*m_dxm
        scale_sum( traits_type::begin( m_xt ), traits_type::end( m_xt ),
                   val1, traits_type::begin( x ), dt, traits_type::begin( m_dxm ) );

        // dt * m_dxh = k4
        system( m_xt, m_dxh, t + dt );
        // x += dt / 6 * (m_dxdt + m_dxt + val2*m_dxm)
        time_type dt6 = dt / static_cast<time_type>(6.0);
        time_type dt3 = dt / static_cast<time_type>(3.0);
        scale_sum_inplace( traits_type::begin( x ), traits_type::end( x ),
                           dt6, traits_type::begin( dxdt ),
                           dt3, traits_type::begin( m_dxt ),
                           dt3, traits_type::begin( m_dxm ),
                           dt6, traits_type::begin( m_dxh ) );
    }

    template<class DynamicSystem>
    void do_step( DynamicSystem& system, container_type& x, time_type t, time_type dt ) {
        system( x, m_dxdt, t );
        do_step( system, x, m_dxdt, t, dt );
    }
};

#endif // !RK4_H

Теперь, когда у нас есть все, что нам нужно; все, что осталось сделать, - это использовать их.

OP упомянул о решении проблемы Лоренца, поэтому я буду использовать это, чтобы продемонстрировать, как использовать ODE:

main.cpp

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>

#include "euler.h"
#include "rk4.h"

#define tab "\t"

typedef std::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() {
    const double dt = 0.01;
    state_type x1( 3 );
    x1[0] = 1.0;
    x1[1] = 0.0;
    x1[2] = 0.0;

    state_type x2( 3 );
    x2[0] = 1.0;
    x2[1] = 0.0;
    x2[2] = 0.0;

    euler_stepper<state_type> stepper_euler;
    stepper_euler.adjust_size( x1 );

    rk4_stepper<state_type> stepper_rk4;
    stepper_rk4.adjust_size( x2 );

    std::fstream file( "compare.txt", std::ios::out );

    file << tab << "Euler Stepper to Solve Lorenz" << tab << tab << tab << tab << "RK4 Stepper to Solve Lorenz\n"
         << tab << "========================" << tab << tab << tab << "=======================\n";

    double t1 = 0.0;
    double t2 = 0.0;
    for( size_t oi = 0; oi < 10000; ++oi, t1 += dt, t2 += dt ) {
        stepper_euler.do_step( lorenz, x1, t1, dt );
        stepper_rk4.do_step( lorenz, x2, t2, dt );
        file << " " << std::setw( 15 ) << std::setprecision( 10 ) << x1[0]
             << tab << std::setw( 15 ) << std::setprecision( 10 ) << x1[1]
             << tab << std::setw( 15 ) << std::setprecision( 10 ) << x1[2]
             << tab << "||"
             << " " << std::setw( 15 ) << std::setprecision( 10 ) << x2[0]
             << tab << std::setw( 15 ) << std::setprecision( 10 ) << x2[1]
             << tab << std::setw( 15 ) << std::setprecision( 10 ) << x2[2]
             << '\n';
    }

    file.close();

    std::cout << "\nPress any key and enter to quit.\n";
    std::cin.get();
    return 0;
}

Теперь вы можете видеть, насколько просто определить функцию Лоренца и просто передать ее тому типу интегратора ODE, который вы хотите использовать.

В качестве дополнительного бонуса есть небольшая приятная функция, которую можно использовать с этими структурами для интеграции констант.

integrate_const.h

#ifndef INTEGRATE_CONST_H
#define INTEGRATE_CONST_H

// This function will iterate the state of the ODE, 
// with the possibility to observe that state in each step

template<class Stepper, class DynamicSystem, class Observer>
size_t integrate_const( Stepper& stepper, DynamicSystem& 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;
}

#endif // !INTEGRATE_CONST_H

Резюме. Вместо использования базовых массивов или указателей теперь можно использовать любой контейнер, который имеет функции начала, конца и изменения размера, что делает код универсальным, простым в реализации и использовании. Приведенный выше код является рабочим кодом, поскольку я не видел ошибок компиляции или компоновки. Я только умеренно протестировал на наличие ошибок времени выполнения. Если у кого-нибудь найдутся баги; Пожалуйста, не голосуйте против из-за ошибки, но не стесняйтесь оставлять комментарий с описанием обнаруженной вами ошибки, и я с радостью исправлю соответствующие ошибки, чтобы сделать эту базу кода более точной.

Примечание. Эта работа была вдохновлена ​​проектом, который я прочитал из code project. Я внес несколько изменений в исходную версию, которую нашел, и удалил все зависимости из библиотеки boost. Вы можете найти его здесь. Также я считаю, что библиотека odeint теперь официально является частью новых версий библиотеки boost.

Вот еще несколько ссылок о работе с ODE, которые могут быть полезны: одни о программировании, другие о математической нотации, стоящей за ним.

• Мечта в коде
• YouTube
• Odeint на Github
• Вопросы и ответы по стеку
• YouTube