Прежде всего, необходимо указать тип данных, который представляет состояниеxвашей системы. Математически это обычно n-мерный вектор с реальными числами или комплексными числами в качестве скалярных объектов. Для одеяния наиболее естественным способом является использование<vector<double>>или<vector<complex<double>>>для представления состояния системы. Однако одеинт может иметь дело и с другими типами контейнеров, например.<boost::array<double,N>>До тех пор, пока он удовлетворяет определенным ниже требованиям.

Чтобы численно интегрировать дифференциальное уравнение, нужно также определить rhs уравненияx' = f(x). В odeint вы предоставляете эту функцию с точки зрения объекта, который реализует ()-оператор с определенной структурой параметров. Следовательно, простым способом было бы просто определить функцию, например:

/* 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];
}

Параметры функции должны следовать примеру выше, где<x>является текущим состоянием, здесь двухкомпонентный вектор, содержащий положениеqи импульсpосциллятора,<dxdt>является производнымx'и должен быть заполнен функцией сf(x), а<t>является текущим временем. Обратите внимание, что в этом примереtне требуется вычислятьf, однако одеинт ожидает, что функциональная подпись будет иметь ровно три параметра (есть исключение, обсуждаемое ниже).

Более сложный подход заключается в реализации системы как класса, где функция 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];
    }
};

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

Числовая интеграция работает итеративно, что означает, что вы начинаете с состоянияx(t)и выполняете шаг времени длиныdt, чтобы получить приблизительное состояниеx(t+dt). Существует множество различных способов выполнения такого временного шага, каждый из которых имеет определенный порядокq. Если порядок способаqявляется точным до термина~ dt^q, что означает, что ошибка вx, сделанная таким шагом, является~ dt^{q +1}.

Некоторые из шагов в таблице выше являются особенными: Одни нуждаются в якобиане ОДЭ, другие — в специальных ОДЭ-системах, подобных гамильтоновым. Мы покажем типичные примеры и примеры использования в этом руководстве и какие шаги следует применять.

Базовый степпер просто выполняет один шаг и не дает вам никакой информации об ошибке, которая была сделана (кроме того, что вы знаете, что это порядокq + 1). Такие степеры используются с постоянным размером шага, который должен быть выбран достаточно маленьким, чтобы иметь разумные небольшие ошибки. Тем не менее, вы должны применить некоторую проверку достоверности ваших результатов (например, наблюдение за сохраненными количествами), потому что у вас нет другого контроля над ошибкой. Следующий пример определяет базовый степпер на основе классической схемы Рунге-Кутта 4-го порядка. Декларация степпера требует государственного типа в качестве шаблонного параметра. Интеграция теперь может быть выполнена с помощью функции<integrate_const(Stepper,System,state,start_time,end_time,step_size )>от odeint:

runge_kutta4< state_type > stepper;
integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

Этот вызов интегрирует систему, определенную<harmonic_oscillator>, используя метод RK4 отt=0до10с шагомdt=0,01и начальным условием, данным в<x>. Результатx(t=10)хранится в<x>(на месте). Каждая ступенька определяет метод<do_step>, который также может быть использован непосредственно. Запишите приведенный выше пример как

const double dt = 0.01;
for( double t=0.0 ; t<10.0 ; t+= dt )
    stepper.do_step( harmonic_oscillator , x , t , dt );

	Tip
	Если у вас есть компилятор с поддержкой C++11, вы можете легко использовать лямбда для создания функции системы: { runge_kutta4< state_type > stepper; integrate_const( stepper , []( const state_type &x , state_type &dxdt , double t ) { dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; } , x , 0.0 , 10.0 , 0.01 ); }

Для улучшения численных результатов и дополнительного минимизации вычислительных усилий целесообразно применение контроля размера шага. Контроль размера шага реализуется с помощью алгоритмов шага, которые дополнительно обеспечивают оценку ошибки применяемого шага. Одесса дает такое количество.Ошибки Стэпперови мы покажем их использование на примере<explicit_error_rk54_ck>— метода Рунге-Кутта 5-го порядка с оценкой ошибок 4-го порядка и коэффициентами, введенными Кэшем и Карпом.

typedef runge_kutta_cash_karp54< state_type > error_stepper_type;

Учитывая погрешность ошибки, все еще нужен экземпляр, который проверяет ошибку и соответствующим образом корректирует размер шага. В Одессе это делается.Контролируемые степперы. Для<runge_kutta_cash_karp54>степпера существует<controlled_runge_kutta>степпер, который можно использовать через

typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;
integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

Как указано выше, это интегрирует систему, определенную<harmonic_oscillator>, но теперь используя адаптивный метод размера шага, основанный на схеме Runge-Kutta Cash-Karp 54 отt=0до10с начальным размером шагаdt=0,01(будет отрегулирован) и начальным условием, данным в x. Результатx(t=10)также будет храниться в x (в месте).

В приведенном выше примере степпер ошибки вложен в контролируемую степперу. Это хорошая техника, но один недостаток заключается в том, что всегда нужно определить обоих шагов. Можно также написать инстанциацию контролируемого степпера в вызов функции интеграции, но полное знание основных степперных типов все еще необходимо. Другой момент заключается в том, что допуски ошибок для управления размером шага нелегко включить в контролируемый шаг. Обе проблемы можно решить с помощью<make_controlled>:

integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) ,
                    harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

make_controlled can be used with many of the steppers of odeint. The first parameter is the absolute error tolerance eps_abs and the second is the relative error tolerance eps_rel which is used during the integration. The template parameter determines from which error stepper a controlled stepper should be instantiated. An alternative syntax of make_controlled is

integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) ,
                    harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

Для контроллера Runge-Kutta ошибка, допущенная в течение одного шага, сравнивается сeps_abs + eps_rel * (a_x* |x | + a_dxdt* dt * |dxdt | ). Если ошибка меньше этого значения, принимается текущий шаг, в противном случае он отклоняется и размер шага уменьшается. Обратите внимание, что размер шага также увеличивается, если ошибка становится слишком маленькой по сравнению с rhs вышеупомянутого отношения. Таким образом, полное воплощение<controlled_runge_kutta>со всеми параметрами

double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
controlled_stepper_type controlled_stepper(
    default_error_checker< double , range_algebra , default_operations >( abs_err , rel_err , a_x , a_dxdt ) );
integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );

При использовании<make_controlled>параметрa_xиa_dxdtиспользуются с их стандартными значениями 1.

В таблицах ниже можно найти все степперы, которые работают с<make_controlled>и<make_dense_output>, которые являются аналогом для плотных выходных степперов.

При использовании<make_controlled>или<make_dense_output>следует знать, какой именно тип используется и как работает контроль размера шага.

odeint поддерживает итераторы для решения ODE. То есть вы инстанцируете пару итераторов и вместо того, чтобы использовать интегрированные процедуры с соответствующим наблюдателем, помещаете итераторы в один из алгоритмов из стандартной библиотеки C++ или из Boost. Диапазон. Примером является

std::for_each( make_const_step_time_iterator_begin( stepper , harmonic_oscillator, x , 0.0 , 0.1 , 10.0 ) ,
               make_const_step_time_iterator_end( stepper , harmonic_oscillator, x ) ,
               []( std::pair< const state_type & , const double & > x ) {
                   cout << x.second << " " << x.first[0] << " " << x.first[1] << "\n"; } );