Сравнение odeint’s runge_kutta4 с od45 Матлаба

Я хотел бы использовать runge_kutta4 метод в библиотека odeint C ++. Я решил проблему в Matlab. Мой следующий код в Matlab, чтобы решить x'' = -x - g*x'с начальными значениями x1 = 1, x2 = 0, как следует

main.m

clear all
clc
t = 0:0.1:10;
x0 = [1; 0];
[t, x] = ode45('ODESolver', t, x0);
plot(t, x(:,1));
title('Position');
xlabel('time (sec)');
ylabel('x(t)');

ODESolver.m

function dx = ODESolver(t, x)
dx = zeros(2,1);
g  = 0.15;
dx(1) =  x(2);
dx(2) = -x(1) - g*x(2);
end

Я установил odeint Library. Мой код для использования runge_kutta4 как следует

#include <iostream>

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::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 lorenz( const state_type &x , state_type &dx ,  double t )
{
dx[0] =  x[1];
dx[1] = -x[0] - gam*x[1];
}int main(int argc, char **argv)
{
const double dt = 0.1;
runge_kutta_dopri5<state_type> stepper;
state_type x(2);
x[0] = 1.0;
x[1] = 0.0;double t = 0.0;
cout << x[0] << endl;
for ( size_t i(0); i <= 100; ++i){
stepper.do_step(lorenz, x , t, dt );
t += dt;
cout << x[0] << endl;
}return 0;
}

Результат на следующем рисунке
введите описание изображения здесь

Мой вопрос: почему результат меняется? Что-то не так с моим кодом C ++?

Это первые значения обоих методов

Matlab    C++
-----------------
1.0000    0.9950
0.9950    0.9803
0.9803    0.9560
0.9560    0.9226
0.9226    0.8806
0.8806    0.8304
0.8304    0.7728
0.7728    0.7084
0.7083    0.6379

Обновить:

Проблема в том, что я забыл включить начальное значение в мой код C ++. Спасибо за то, что @horchler это заметил. После включения правильных значений и использования runge_kutta_dopri5 как предложил @horchler, результат

Matlab    C++
-----------------
1.0000    1.0000
0.9950    0.9950
0.9803    0.9803
0.9560    0.9560
0.9226    0.9226
0.8806    0.8806
0.8304    0.8304
0.7728    0.7728
0.7083    0.7084

Я обновил код, чтобы отразить эти изменения.

4

Решение

runge_kutta4 степпер в одеинте не имеет ничего общего с матлабом ode45, которая является адаптивной схемой, основанной на Дорманда-Prince метод. Чтобы повторить результаты Matlab, вы, вероятно, должны попробовать runge_kutta_dopri5 шаговый. Кроме того, убедитесь, что ваш код C ++ использует те же абсолютные и относительные допуски, что и ode45 (по умолчанию 1e-6 а также 1e-3соответственно). Наконец, похоже, что вы не сохраняете / не распечатываете исходное состояние в результатах C ++.

Если вы не уверены, почему ode45 не предпринимает фиксированных шагов, даже если вы указали t = 0:0.1:10;, увидеть мой подробный ответ здесь.

Если вы должны использовать фиксированный шагrunge_kutta4 пошагово, тогда вам нужно будет уменьшить шаг интеграции в вашем коде C ++, чтобы он соответствовал выводу Matlab.

9

Другие решения

Функция Matlab ode45 уже включает в себя контроль ошибок и, я думаю, также интерполяцию (плотный вывод). для сравнения с boost.odeint вы должны использовать ту же функциональность. Boost.odeint обеспечивает integrate функции, которые выполняют управление размером шага и плотный вывод, если используемый пошаговый алгоритм обеспечивает эту функциональность. Следующий фрагмент кода показывает, как это используется с параметрами контроля ошибок по умолчанию из Matlab, предоставленными horchler:

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::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 damped_osc( const state_type &x , state_type &dx , const double t )
{
dx[0] =  x[1];
dx[1] = -x[0] - gam*x[1];
}

void print( const state_type &x, const double t )
{
cout << x[0] << endl;
}

int main(int argc, char **argv)
{
cout.precision(16);  // full precision output
const double dt = 0.1;
typedef runge_kutta_dopri5<state_type> stepper_type;
state_type x(2);
x[0] = 1.0;
x[1] = 0.0;

integrate_const(make_dense_output<stepper_type>( 1E-6, 1E-3 ),
damped_osc, x, 0.0, 10.0, dt , print);

return 0;
}

Обратите внимание, что результаты могут все еще не быть именно так то же самое (как и во всех 16 цифрах), потому что контроль ошибок в Boost.odeint не может быть затруднен именно так как в одлаб Matlab.

3

По вопросам рекламы [email protected]