Первая специализация - это явные степарды. Explicit означает, что новое состояние од может быть вычислено явно из текущего состояния без решения имплицитных уравнений. Такие степи имеют общее представление о том, что они оценивают систему во времени t таким образом, что результат f(x,t) может быть передан на степку. В odeint у явного степного есть два дополнительных метода

do_step( sys , inout , dxdtin , t , dt>>

do_step( sys , , dxdtin , t , out , dt>>>>>>>>>>>>>

Здесь дополнительным параметром является значение функции f в состоянии x и времени t. Примером может служить Runge-Kutta stepper четвертого порядка:

runge_kutta4< state_type > rk;
rk.do_step( sys1 , inout , t , dt );               // In-place transformation of inout
rk.do_step( sys2 , inout , t , dt );               // call with different system: Ok
rk.do_step( sys1 , in , t , out , dt );            // Out-of-place transformation
rk.do_step( sys1 , inout , dxdtin , t , dt );      // In-place tranformation of inout
rk.do_step( sys1 , in , dxdtin , t , out , dt );   // Out-of-place transformation

На самом деле вам не нужно называть эти два метода. Вы всегда можете использовать более простой do_step( sys , in , t , dt ), но иногда производная состояния необходима внешне для выполнения некоторых внешних вычислений или для проведения статистического анализа.

Особым классом явных степперов являются FSAL (первая, та же, как и последняя оценка функции системы также первая оценка следующего шага. Для таких шагов метод do_step немного отличается:

do_step( sys , inout , dxdtinout , t , dt )

do_step( sys , in , dxdtin , out , dxdtout , t , dt )

Этот метод берет производную на время t, а также сохраняет производную на время t+dt. Называя эти функции впоследствии итерируя вдоль решения, один сохраняет один вызов функции, передав результат для dxdt в следующий вызов функции. Однако при использовании FSAL степперов без поставок производных:

do_step( sys , inout , t , dt )

степпер внутренне удовлетворяет свойства FSAL, что означает, что он помнит последний dxdt и использует его для следующего шага. Примером FSAL stepper является Runge-Kutta-Dopri5 stepper. Фокус FSAL иногда также называют трюком Fehlberg. Примером того, как можно использовать степарды FSAL, является

runge_kutta_dopri5< state_type > rk;
rk.do_step( sys1 , in , t , out , dt );
rk.do_step( sys2 , in , t , out , dt );         // DONT do this, sys1 is assumed
rk.do_step( sys2 , in2 , t , out , dt );
rk.do_step( sys2 , in3 , t , out , dt );        // DONT do this, in2 is assumed
rk.do_step( sys1 , inout , dxdtinout , t , dt );
rk.do_step( sys2 , inout , dxdtinout , t , dt );           // Ok, internal derivative is not used, dxdtinout is updated
rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt );
rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used

Caution

FSAL-шагпы сохраняют производную на время t +dt внутренне, если они называются через do_step( sys , , , t , dt >>>>>>>>>>. Первый вызов do_step инициализирует dxdt, и для всех следующих вызовов предполагается, что используется одна и та же система и одно и то же состояние. Если вы используете FSAL Stepper в рамках интеграционных функций, это делается автоматически. Дополнительные сведения см. в разделе Использование степных шагов или зайдите в таблицу ниже, чтобы увидеть, какой ступенчатый имеет внутреннее состояние.

Как упоминалось выше, для гамильтоновых систем используются симплектические решетки. Симплектические растворители сохраняют объем фазового пространства точно, и если гамильтонская система является энергетической консервативной, они также сохраняют энергию приблизительно. dqdt/dt = f1(p), dpdt/dt = f2(q), dpdt/dt = f2(q), где (q,p) являются состоянием системы. Пространство (q,p) иногда называют фазовым пространством и q и p называют переменными фазового пространства. Симплектические системы в этой специальной форме широко встречаются в природе. Например, полная классическая механика, написанная Ньютоном, Лагранжем и Гамильтоном, может быть сформулирована в этой структуре. Разделимость системы зависит от конкретного выбора координат.

Симплектические системы могут быть решены путем одеинта с помощью симплектического шага и симплектического метода Runge-Kutta-Nystrom четвертого порядка. Эти степеры предполагают, что система автономна, поэтому времени явно не будет. Далее они в принципе выполняют концепцию Stepper по умолчанию, но они ожидают, что система будет парой приемлемых объектов. Первая запись этой пары рассчитывается f1(p), а вторая рассчитывает f2(q). Синтаксис sys.первая(p,dqdt) и sys.секунда(q,dpdt), где первая и вторая часть может быть просто фуна. Примером является гармонический осциллятор:

typedef boost::array< double , 1 > vector_type;
struct harm_osc_f1
{
    void operator()( const vector_type &p , vector_type &dqdt )
    {
        dqdt[0] = p[0];
    }
};
struct harm_osc_f2
{
    void operator()( const vector_type &q , vector_type &dpdt )
    {
        dpdt[0] = -q[0];
    }
};

Состояние такого ODE теперь также состоит из двух частей, часть для q (также называемая координатами) и часть для p (момента). Полный пример гармонического осциллятора сейчас:

pair< vector_type , vector_type > x;
x.first[0] = 1.0; x.second[0] = 0.0;
symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt );

Если вам нравится представлять систему одним классом, вы можете легко связать два общественных метода:

struct harm_osc
{
    void f1( const vector_type &p , vector_type &dqdt ) const
    {
        dqdt[0] = p[0];
    }
    void f2( const vector_type &q , vector_type &dpdt ) const
    {
        dpdt[0] = -q[0];
    }
};

harm_osc h;
rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
        x , t , dt );

Многие гамильтонская система может быть написана как dq/dt=p, dp/dt=f(q), что вычислительно намного проще, чем полная сепарируемая система. Очень часто также возможно трансформировать оригинальные уравнения движения, чтобы привести систему в эту упрощенную форму. Этот вид системы может быть использован в симплектических решателях, просто пройдя f(p) к методу do_step, опять f(p) будет представлен простой C-функцией или фанктором. Здесь приведенный выше пример гармонического осциллятора может быть написан как

pair< vector_type , vector_type > x;
x.first[0] = 1.0; x.second[0] = 0.0;
symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
rkn.do_step( harm_osc_f1() , x , t , dt );

В этом примере функция harm_osc_f1 является точно такой же функцией, как и в приведенных выше примерах.

Обратите внимание, что состояние ODE не должно быть открыто сконструировано через pair< vector_type , vector_type > x. Можно также использовать комбинацию make_pair и ref. Кроме того, существует удобная версия do_step, которая принимает q и p без объединения их в пару:

rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
rkn.do_step( harm_osc_f1() , q , p , t , dt );
rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt );

	Caution
	Этот раздел не обновляется.

Для некоторых систем стабильность свойств классического Runge-Kutta недостаточна, особенно если система считается жесткой. Жесткая система обладает двумя или более временными масштабами совершенно другого порядка. Растворы для жестких систем обычно имплицитны, что означает, что они решают уравнения, такие как x(t+dt) = x(t) + dt * f(x(t+1)). Эта конкретная схема является неявным методом Euler. Имплицитные методы обычно решают систему уравнений по алгоритму поиска корня, подобно методу Ньютона, и поэтому должны знать якобиана системы J​ij = df​i/dx​j.

Для неявных решателей система снова пара, где первый компонент вычисляет f(x,t) и второй якобиец. sys.первый( x , dxdt , Для неявного решателя state_type является ublas::vector, а якобиец представлен ublas::matrix.

	Important
	Implicit Solutions работает только с ublas::vector как государственный тип. На данный момент никакие другие типы состояний не поддерживаются.

Другой большой класс решающих - многоступенчатый метод. Они сохраняют небольшую часть истории решения и вычисляют следующий шаг с помощью этой истории. Поскольку многоступенчатые методы знают часть своей истории, им не нужно очень часто вычислять функцию системы, обычно она вычисляется только один раз. Это делает многоступенчатые методы предпочтительными, если вызов функции системы стоит дорого. Примерами являются ODE, определенные в сетях, где вычисления взаимодействия обычно стоят дорого (и могут быть порядка O(N^2).

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

adams_bashforth_moulton< 5 , state_type > abm;
abm.initialize( sys , inout , t , dt );
abm.do_step( sys , inout , t , dt );

Инициализация использует пассионер четвертого порядка Runge-Kutta и после вызова иализировать состояние inout изменилось на нынешнее состояние, так что его можно сразу же использовать, передавая его на следующие звонки do_step. Вы также можете использовать собственные степперсы для инициализации внутреннего состояния Adams-Bashforth-Stepper:

abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );

Многие многоступенчатые методы также являются явными шагерами, поэтому параметр do_step метод не отличается от явных шагов.

	Caution
	Многоступенчатые методы имеют некоторые внутренние переменные, которые зависят от явного решения. Следовательно, после любых внешних изменений вашего состояния (например, размера) или системы функция инициализации должна быть снова вызвана, чтобы скорректировать внутреннее состояние степного. Если вы используете интеграционные функции, это будет учтено. Подробнее см. раздел Использование степей.

Многие из приведенных выше степей обладают возможностью использования адаптивного поэтапного контроля. Адаптивная ступенчатая интеграция работает в принципе следующим образом:

Класс контролируемых степеров имеет свою собственную концепцию в odeint - концепцию Controlled Stepper. Они, как правило, построены из лежащих в основе ошибок шагов. Примером является контроллер для явных степных шагов Runge-Kutta. Степперы Runge-Kutta входят в контроллер как аргумент шаблона. Кроме того, можно передать Runge-Kutta stepper конструктору, но этот шаг не нужен; степка по умолчанию построена, если это возможно.

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

Классический способ решить, отклоняется или принимается шаг - это рассчитать

val = || | err​i | / ( ε​abs + ε​rel * ( a​x | x​i | + a​dxdt | | dxdt​i | )||

ε​abs и ε​rel являются абсолютными и относительными допусками ошибок, а || x || является нормой, как правило, ||x||=(Σ​i<3 x&8203;i<3222>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><2>>>>>>>>>><2>>>>>><2>>>>>>>>>>>>>>> Шаг отклоняется, если val больше, чем 1, иначе он принимается. Подробные сведения об использованных нормах и допусках ошибок см. в таблице ниже.

Для контролируемый_runge_kutta шаг за шагом новый размер шага затем рассчитывается через

val > 1: dt​new = dt​current max(0.9 pow(val, -1/ (O​E - 1) ) , 0.2 )

val < 0.5 : dt​new = dt​current min( 0.9 pow( val , -1 / O​S ) , 5 )

else : dt​new = dt​current

Вот. O​S и O​E являются порядком степного и погрешного степного. Эти формулы также содержат некоторые факторы безопасности, избегая того, что размер шага уменьшается или увеличивается до многих. Подробная информация о реализации контролируемых степек в одине содержится в таблице ниже.

Существует четвертый класс степперов, которые являются так называемыми плотными выходными степями. Степейные степи могут предпринимать более масштабные шаги и интерполировать решение между двумя последовательными точками. Эти интерполированные точки обычно имеют тот же порядок, что и порядок степного. Степперы Dense-output часто являются композитными шагерами, которые принимают базовый метод в качестве параметра шаблона. Примером может служить dense_output_runge_kutta stepper, который делает Runge-Kutta stepper с плотно-выходными объектами в качестве аргумента. Не все степные степи Runge-Kutta обеспечивают расчет плотной пропускной способности; на данный момент только степпер Dormand- Prince 5 обеспечивает плотный выход. Пример:

dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
dense.initialize( in , t , dt );
pair< double , double > times = dense.do_step( sys );
(void)times;

У степпера Dense есть своя концепция. Главное отличие от обычных степеров заключается в том, что они управляют состоянием и временем внутренне. Если вы называете do_step, только ODE передается в качестве аргумента. Кроме того, do_step возвращает последний временной интервал: t и t+dt, следовательно, вы можете интерполировать решение между этими двумя точками. Другое отличие заключается в том, что они должны быть инициализированы с помощью интиализировать, в противном случае внутреннее состояние степного устройства по умолчанию строится, что может привести к смешным ошибкам или ошибкам.

Построение плотной выходной степи выглядит немного противно, так как в случае dense_output_runge_kutta на шаге должен быть вложен контролируемый степей и степей ошибки. Для упрощения генерации плотных выходных шаговых генерирующих функций существуют:

typedef boost::numeric::odeint::result_of::make_dense_output<
    runge_kutta_dopri5< state_type > >::type dense_stepper_type;
dense_stepper_type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
(void)dense2;

Это утверждение также длинное; оно демонстрирует, как make_dense_output может быть использован с протоколом result_of. Параметры make_dense_output - абсолютная допуска ошибок, относительная допуска ошибок и ступенчатый. Это явно предполагает, что лежащий в основе степей является контролируемым степи и что этот степей имеет абсолютную и относительную допустимость ошибок. Подробности о функциях генерации см. в Generation functions. Функции поколения были разработаны для простого использования с интегрирующими функциями:

integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );

В этом разделе содержится некоторая общая информация об использовании степных в унитазе.

Степперы копируются по значению

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

Степперы могут иметь внутреннее состояние

	Caution
	Некоторые функции, описанные в этом разделе, еще не реализованы

Некоторые степеры требуют хранить некоторую информацию о состоянии ODE между двумя шагами. Примерами являются многоступенчатые методы, которые хранят часть решения во время эволюции ODE, или FSAL Steppers, которые хранят последнюю производную на момент времени t+dt, которые будут использоваться на следующем этапе. В обоих случаях степеры ожидают, что последовательные звонки do_step происходят из одного и того же решения и одного и того же ODE. В этом случае абсолютно необходимо, чтобы вы назвали do_step с той же системной функцией и тем же состоянием, см. также примеры для FSAL Steppers выше.

Степпер с внутренним государством поддерживает два дополнительных метода: перезагрузка, которые реанимируют государство и инициализируют, которые инициализируют внутреннее состояние. Параметры инициализировать зависят от конкретного ступенчатого. >, , , , > , >>>>> Для случая FSAL степперы, интиализируйте интиализируйте( sys , in , t ) которые просто вычисляют r.h.s. ODE и присваивают ее значение внутреннему производному.

Все эти степарды имеют общее, что изначально заполняют свое внутреннее состояние самостоятельно. Следовательно, вы не обязаны называть инициализацию. Смотрите, как это работает для Adams-Bashforth-Moulton stepper: в примере мы обобщаем четвертый заказ Adams-Bashforth-Moulton stepper, что означает, что он будет хранить 4 внутренних производных решения в разы (t-dt,t<>>><>>>>>><>>>.

adams_bashforth_moulton< 4 , state_type > stepper;
stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the first internal state
                                       // the internal array is now [x(t-dt)]
stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the second internal state
                                       // the internal state array is now [x(t-dt), x(t-2*dt)]
stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the third internal state
                                       // the internal state array is now [x(t-dt), x(t-2*dt), x(t-3*dt)]
stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the fourth internal state
                                       // the internal state array is now [x(t-dt), x(t-2*dt), x(t-3*dt), x(t-4*dt)]
stepper.do_step( sys , x , t , dt );   // make one step with Adam-Bashforth-Moulton, the internal array of states is now rotated

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

Stepper может быть реализуемым

Почти все степи в одине должны хранить некоторые промежуточные результаты типа state_type или deriv_type. Для этого Одеинту необходимо некоторое управление памятью для внутренних современников. Поскольку это управление памятью, как правило, связано с корректировкой размера векторных типов, оно называется смолой в odeint. Таким образом, большинство степперов в одине обеспечивают дополнительный параметр шаблона, который контролирует корректировку размера внутренних переменных - ретранслятора. Подробно odeint предоставляет три политических класса (ресайзеры) always_resizer, initially_resizer, и never_resizer. Кроме того, все степи имеют метод adjust_size, который принимает параметр, представляющий тип состояния и который вручную корректирует размер внутренних переменных, соответствующих размеру данного экземпляра. Прежде чем выполнить фактическую ретрансляцию одина, всегда проверяется, различаются ли размеры состояния и внутренняя переменная и только пересчитывается, если они разные.

	Note
	Вам нужно беспокоиться только о распределении памяти при использовании динамически размерных векторных типов. Если ваш тип штата распределен, как boost::array, никакого распределения памяти не требуется.

По умолчанию параметр первоначально_resizer, что означает, что первый вызов do_step выполняет релизию, отсюда и распределение памяти. Если вы изменили размер своей системы и свое состояние, вы должны назвать adjust_size вручную в этом случае. Второй ретранслятор - always_resizer, который пытается изменить размер внутренних переменных по каждому вызову do_step. Типичные случаи использования для такого рода резизеров - это саморасширяющиеся решетки, как показано в учебнике (Само расширяющиеся решетки) или частичные дифференциальные уравнения с адаптивной сеткой. Здесь не требуется никаких звонков adjust_size, степеры управляют всем собой. Третий класс резизеров - Never_resizer, что означает, что внутренние переменные никогда не корректируются автоматически и всегда должны быть скорректированы вручную.

Существует второй механизм, который влияет на переделку и который контролирует, если тип штата по крайней мере изменен - мета-функция is_resizeable. Эта мета-функция возвращает статическое значение Boolean, если какой-либо тип реализуется. Например, он вернет Правда для std::vector< T > но false для boost::array<T>>>>>>>>>> По умолчанию и для неизвестных типов is_resizeable возвращает false, поэтому, если у вас есть свой тип, вам нужно специализироваться на этой мета-функции. Более подробную информацию о механизме переноса см. в разделе Адаптируйте свои типы состояний.

Частные степи должны быть использованы в какой ситуации

odeint предоставляет довольно большое количество различных степей, таких, что пользователю остается вопрос о том, какой степей соответствует его потребностям. Наши личные рекомендации:

Наконец, можно также написать новые степи, которые полностью совместимы с odeint. Они должны выполнять только один или несколько шаговых Концепций окраски.

Мы будем иллюстрировать, как написать свой собственный степей с примером стохастического метода Euler. Этот метод подходит для решения стохастических дифференциальных уравнений (SDE). A SDE имеет форму

dx/dt = f(x) + g(x) ξ(t)

где ξ - это гусский белый шум с нулевым средним и стандартным отклонением σ(t). f(x) считается детерминированной частью, в то время как g(x) ξ является шумной частью. В случае, если g(x) не зависит от x, SDE, как говорят, имеет аддитивный шум. Невозможно решить SDE с классическими решателями для ODE, поскольку шумная часть SDE должна быть масштабирована по-разному, а затем детерминированной частью по отношению к шагу времени. Но существует много решающих для SDE. Классический и простой метод - это стохастический разработчик Euler. Он работает на итерации

x(t+Δ t) = x(t) + Δ t f(x(t)) + Δ t1/2 g(x) ξ(t)

где ξ (t) является независимой обычной распределенной случайной переменной.

Теперь мы реализуем этот метод. Мы назовем степпер stochastic_euler. Он моделирует концепцию Stepper. Для простоты мы фиксируем тип состояния как array< ДП,N> Определение класса выглядит как

template< size_t N > class stochastic_euler
{
public:
    typedef boost::array< double , N > state_type;
    typedef boost::array< double , N > deriv_type;
    typedef double value_type;
    typedef double time_type;
    typedef unsigned short order_type;
    typedef boost::numeric::odeint::stepper_tag stepper_category;
    static order_type order( void ) { return 1; }
    // ...
};

Типы необходимы для выполнения степной концепции. В качестве внутреннего состояния и типа дерива мы используем простые массивы в стохастическом Евлере, они необходимы для современников. Степпер имеет заказ один, который возвращается из функции order().

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

template< size_t N > class stochastic_euler
{
public:
    // ...
    template< class System >
    void do_step( System system , state_type &x , time_type t , time_type dt ) const
    {
        deriv_type det , stoch ;
        system.first( x , det );
        system.second( x , stoch );
        for( size_t i=0 ; i<x.size() ; ++i )
            x[i] += dt * det[i] + sqrt( dt ) * stoch[i];
    }
};

Вот и все. Это довольно просто, и стохастическая ступенчатая реализация Euler здесь довольно общая. Конечно, это может быть улучшено, например

Теперь давайте посмотрим, как мы используем новый степень. Хороший пример - процесс Орнштейна-Уленбека. Он состоит из простого движения Брауна, совпадающего с процессом расслабления. Его SDE читает

dx/dt = - x + ξ

где ξ является гусийским белым шумом со стандартным отклонением σ. Реализация процесса Ornstein-Uhlenbeck довольно проста. Нам нужны две функции или фанкторы - одна для детерминированной и одна для стохастической части:

const static size_t N = 1;
typedef boost::array< double , N > state_type;
struct ornstein_det
{
    void operator()( const state_type &x , state_type &dxdt ) const
    {
        dxdt[0] = -x[0];
    }
};
struct ornstein_stoch
{
    boost::mt19937 &m_rng;
    boost::normal_distribution<> m_dist;
  ornstein_stoch( boost::mt19937 &rng , double sigma ) : m_rng( rng ) , m_dist( 0.0 , sigma ) { }
    void operator()( const state_type &x , state_type &dxdt )
    {
        dxdt[0] = m_dist( m_rng );
    }
};

В стохастической части мы использовали скручиватель Мерсенна для случайного поколения чисел и гауссийский белый шумогенератор нормальное распределение со стандартным отклонением σ. Теперь мы можем использовать стохастический ступенчатый Euler с интегрирующими функциями:

boost::mt19937 rng;
double dt = 0.1;
state_type x = {{ 1.0 }};
integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( rng , 1.0 ) ),
        x , 0.0 , 10.0 , dt , streaming_observer() );

Обратите внимание, как мы использовали функцию make_pair для генерации функции системы.

odeint предоставляет мета-алгоритм шаблона C++ для создания произвольных схем Runge-Kutta [1]. Некоторые схемы предопределены в odeint, например классический Runge-Kutta четвертого порядка, или Runge-Kutta-Cash-Karp 54 и метод Runge-Kutta-Fehlberg 78. Вы можете использовать этот мета-алгоритм, чтобы построить свои собственные решатели. Это имеет преимущество, что вы можете в полной мере использовать алгебру и операционную систему odeint.

Рассмотрим, например, метод Хена, определенный следующим бакау:

c1 = 0
c2 = 1/3, a21 = 1/3
c3 = 2/3, a31 =  0 , a32 = 2/3
          b1  = 1/4, b2  = 0  , b3 = 3/4

Реализация этого метода очень проста. Сначала нужно определить константы:

template< class Value = double >
struct heun_a1 : boost::array< Value , 1 > {
    heun_a1( void )
    {
        (*this)[0] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
    }
};
template< class Value = double >
struct heun_a2 : boost::array< Value , 2 >
{
    heun_a2( void )
    {
        (*this)[0] = static_cast< Value >( 0 );
        (*this)[1] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
    }
};
template< class Value = double >
struct heun_b : boost::array< Value , 3 >
{
    heun_b( void )
    {
        (*this)[0] = static_cast<Value>( 1 ) / static_cast<Value>( 4 );
        (*this)[1] = static_cast<Value>( 0 );
        (*this)[2] = static_cast<Value>( 3 ) / static_cast<Value>( 4 );
    }
};
template< class Value = double >
struct heun_c : boost::array< Value , 3 >
{
    heun_c( void )
    {
        (*this)[0] = static_cast< Value >( 0 );
        (*this)[1] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
        (*this)[2] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
    }
};

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

const boost::array< double , 1 > heun_a1 = {{ 1.0 / 3.0 }};
const boost::array< double , 2 > heun_a2 = {{ 0.0 , 2.0 / 3.0 }};
const boost::array< double , 3 > heun_b = {{ 1.0 / 4.0 , 0.0 , 3.0 / 4.0 }};
const boost::array< double , 3 > heun_c = {{ 0.0 , 1.0 / 3.0 , 2.0 / 3.0 }};

Но затем вы прибиты, чтобы использовать дубли.

Далее вам нужно определить свой ступенчатый, обратите внимание, что метод Heun имеет 3 этапа и производит приближения порядка 3:

template<
    class State ,
    class Value = double ,
    class Deriv = State ,
    class Time = Value ,
    class Algebra = boost::numeric::odeint::range_algebra ,
    class Operations = boost::numeric::odeint::default_operations ,
    class Resizer = boost::numeric::odeint::initially_resizer
>
class heun : public
boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
                                             Algebra , Operations , Resizer >
{
public:
    typedef boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
                                                         Algebra , Operations , Resizer > stepper_base_type;
    typedef typename stepper_base_type::state_type state_type;
    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
    typedef typename stepper_base_type::value_type value_type;
    typedef typename stepper_base_type::deriv_type deriv_type;
    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
    typedef typename stepper_base_type::time_type time_type;
    typedef typename stepper_base_type::algebra_type algebra_type;
    typedef typename stepper_base_type::operations_type operations_type;
    typedef typename stepper_base_type::resizer_type resizer_type;
    typedef typename stepper_base_type::stepper_type stepper_type;
    heun( const algebra_type &algebra = algebra_type() )
    : stepper_base_type(
            fusion::make_vector(
                heun_a1<Value>() ,
                heun_a2<Value>() ) ,
            heun_b<Value>() , heun_c<Value>() , algebra )
    { }
};

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

typedef boost::array< double , 3 > state_type;
heun< state_type > h;
state_type x = {{ 10.0 , 10.0 , 10.0 }};
integrate_const( h , lorenz() , x , 0.0 , 100.0 , 0.01 ,
                 streaming_observer( std::cout ) );