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

9

runge_kutta4odeint 中的步进器与 Matlab 完全不同ode45，它是一种基于Dormand-Prince方法的自适应方案。要复制 Matlab 的结果，您可能应该尝试使用runge_kutta_dopri5步进器。此外，请确保您的 C++ 代码使用与ode45（默认分别为1e-6和1e-3）相同的绝对和相对容差。最后，看起来您可能没有在 C++ 结果中保存/打印初始条件。

如果您对ode45即使您指定了为什么不采取固定步骤感到困惑t = 0:0.1:10;，请在此处查看我的详细答案。

如果您必须使用固定runge_kutta4步进步进器，那么您需要减少 C++ 代码中的集成步长以匹配 Matlab 的输出。

于 2014-11-05T02:41:53.200 回答
3

Matlab ode45 函数已经包括错误控制，我认为也包括插值（密集输出）。要与 boost.odeint 进行比较，您应该在那里使用相同的功能。integrate如果使用的步进算法提供此功能，则Boost.odeint 提供执行步长控制和密集输出的功能。以下代码片段显示了如何将其与 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;
}

请注意，结果可能仍然不完全相同（与所有 16 位数字一样），因为 Boost.odeint 中的错误控制可能与 Matlab 的 ode45 中的错误控制不同。

于 2014-11-05T08:23:22.683 回答