4

odeint C++ ライブラリrunge_kutta4のメソッドを使用したいと考えています。Matlabで問題を解決しました。を初期値,で解決するための Matlab の私の次のコードは次のとおりです。x'' = -x - g*x'x1 = 1x2 = 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

これらの変更を反映するようにコードを更新しました。

4

2 に答える 2

3

Matlab ode45 関数には既にエラー制御が含まれており、補間 (密な出力) も含まれていると思います。boost.odeint と比較するには、そこで同じ機能を使用する必要があります。Boost.odeint はintegrate、使用されるステッパー アルゴリズムがこの機能を提供する場合、ステップ サイズ制御と高密度出力を実行する関数を提供します。次のコード部分は、horchler によって提供された Matlab の既定のエラー制御パラメーターでこれを使用する方法を示しています。

#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;
}

Boost.odeint のエラー制御は、Matlab の ode45 のように正確に実装されていない可能性があるため、結果が(16 桁すべてのように)まったく同じではない可能性があることに注意してください。

于 2014-11-05T08:23:22.683 に答える