odeint的runge_kutta4與Matlab的ode45的比較

Question

我想在odeint C++ library中使用runge_kutta4方法。我已經在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 

---

更新：

的問題是，我忘了包括在我的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

在odeint中的runge_kutta4步進器與Matlab的ode45完全不同，它是基於Dormand-Prince方法的自適應方案。要複製Matlab的結果，您應該嘗試使用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可能不會像在Matlab的ode45中一樣受到的限制。