Сравнение odeint's runge_kutta4 с ode45 Matlab

Question

Я хотел бы использовать метод runge_kutta4 в odeint C++ library. Я решил проблему в 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. Мой код для использования 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 

---

Update:

Проблема заключается в том, что я забыл включить начальное значение в моем коде 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 

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

Answer 1

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

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

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

Answer 2

Функция 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 не может быть impemented точно как в ode45 Matlab.