Воображая, вы хотите численно интегрировать гармонический осциллятор с трением. Уравнения движения даныx'' = -x + γ x''. Odeint имеет дело только с ODE первого порядка, которые не имеют более высоких производных, чем x'. Однако любой ODE более высокого порядка может быть преобразован в систему ODE первого порядка путем введения новых переменных.q=xиp=x'так чтоw=(q,p). Для применения численной интеграции сначала нужно спроектировать правую сторону уравненияw' = f(w) = (p,-q+γ p):

/* The type of container used to hold the state vector */
typedef std::vector< double > state_type;
const double gam = 0.15;
/* The rhs of x' = f(x) */
void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
{
    dxdt[0] = x[1];
    dxdt[1] = -x[0] - gam*x[1];
}

Здесь мы выбрали<vector<double>>как тип состояния, но возможны и другие, например<boost::array<double,2>>. odeint разработан таким образом, что вы можете легко использовать свои собственные типы состояния. Далее определяется ОДЭ, являющаяся в данном случае простой функцией вычисленияf(x). Параметрическая подпись этой функции имеет решающее значение: методы интеграции всегда будут называть их в виде<f(x, dxdt, t)>(есть исключения для некоторых специальных процедур). Таким образом, даже если нет явной зависимости времени, необходимо определить<t>как параметр функции.

Теперь мы должны определить начальное состояние, с которого должна начаться интеграция:

state_type x(2);
x[0] = 1.0; // start at x=1.0, p=0.0
x[1] = 0.0;

Для самой интеграции мы используем функцию<integrate>, которая является удобным способом получения быстрых результатов. Он основан на погрешно-управляемом<runge_kutta54_cash_karp>степпере (5-го порядка) и использует адаптивные пошаговые размеры.

size_t steps = integrate( harmonic_oscillator ,
        x , 0.0 , 10.0 , 0.1 );

Интегрирующая функция ожидает в качестве параметров rhs оды, как определено выше, начальное состояние<x>, время начала и окончания интеграции, а также начальный шаг времени = размер. Обратите внимание, что<integrate>использует адаптивный размер шага во время этапов интеграции, поэтому временные точки не будут одинаково разнесены. Интеграция возвращает количество шагов, которые были применены, и обновляет x, который устанавливается на приблизительное решение ODE в конце интеграции.

Также можно представить одную систему как класс. Затем rhs должны быть реализованы в виде функтора - класса с перегруженной функцией вызова оператора:

/* The rhs of x' = f(x) defined as a class */
class harm_osc {
    double m_gam;
public:
    harm_osc( double gam ) : m_gam(gam) { }
    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
    {
        dxdt[0] = x[1];
        dxdt[1] = -x[0] - m_gam*x[1];
    }
};

который можно использовать через

harm_osc ho(0.15);
steps = integrate( ho ,
        x , 0.0 , 10.0 , 0.1 );

Для того, чтобы наблюдать за решением на этапах интеграции, все, что вам нужно сделать, это предоставить разумного наблюдателя. Примером является

struct push_back_state_and_time
{
    std::vector< state_type >& m_states;
    std::vector< double >& m_times;
    push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
    : m_states( states ) , m_times( times ) { }
    void operator()( const state_type &x , double t )
    {
        m_states.push_back( x );
        m_times.push_back( t );
    }
};

который хранит промежуточные ступени в контейнере. Обратите внимание, что структура аргумента ()-оператора: odeint вызывает наблюдателя именно таким образом, обеспечивая текущее состояние и время. Теперь вам нужно только передать этот контейнер функции интеграции:

vector<state_type> x_vec;
vector<double> times;
steps = integrate( harmonic_oscillator ,
        x , 0.0 , 10.0 , 0.1 ,
        push_back_state_and_time( x_vec , times ) );
/* output */
for( size_t i=0; i<=steps; i++ )
{
    cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
}

Вот и все. Вы можете использовать функциональные библиотеки, такие какBoost.LambdaилиBoost.Phoenix, чтобы облегчить создание функций наблюдателя.

Полный файл cpp для этого примера можно найти здесь:harmonic_oscillator.cpp