使用 boost odeint 返回导数

Question

我正在使用 boost odeint 来积分 2n 个耦合方程。我的基本问题很笼统。我如何 return 不仅可以 x_i 使用观察器的每个时间步长的状态 return 还可以 return 导数 d(x_i)/dt？我知道如何 return 给定时间的状态或状态的某些函数，但我不知道如何 return 给定时间和之前的状态函数.

更具体地说，这是我的代码的简化版本。

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/ublas/matrix.hpp>
typedef boost::numeric::ublas::matrix<double> matrixd;
using namespace boost::numeric::odeint;
using namespace std;

struct k_ring
{
    int     m_N;

    k_ring( int N = 5 ) :
    m_N ( N ) {}

    void operator() (const matrixd &x, matrixd &dxdt, const double )
    {
        matrixd loc_coup( m_N, 2, 0 );
    
        for( int i=0 ; i<m_N ; ++i ){
            for( int k=0; k<m_N ; ++k ){
                if(abs(k-i) <= 30 && dist > 0){
                    loc_coup(i,0) += sin( x(i,0) - x(k,0) + 0.45 );
                    loc_coup(i,1) += sin( x(i,1) - x(k,1) + 0.45 );
                }
            }
            dxdt(i,0) = - (0.1/61)*loc_coup(i,0) + 0.2*sin( x(i,1) - x(i,0) );
            dxdt(i,1) = - (0.1/61)*loc_coup(i,1) - 0.2*sin( x(i,1) - x(i,0) );
        }
    }
};

struct observer
{
    const double    m_sim_time;

    observer( const float &sim_time ) :
    m_sim_time( sim_time ) {}

    void operator()( matrixd &x, double t)
    {
        for(int i=0; i < x.size1(); ++i){
            cout << t/m_sim_time << '\t' << x(i,0) << '\t' << '\n';
        }
    }
};


int main()
{
    double dt = 0.01;
    int n = 256;
    double sim_time = 100.0;
    matrixd x( n , 2 , 0. );

    k_ring system(n);

    observer obs(sim_time);

    integrate_const(runge_kutta4< matrixd >(), system, x, 0.0, sim_time, dt, boost::ref( obs ));
    
    return 0;
}

如果我定义 psi(i) = x(i,1)-x(i,0)，我如何计算每个时间步的导数 d psi(i)/dt？

Answer 1

一个简单的例子来说明我的评论

You would have to call the system function from inside the observer to get the derivatives values for the state. That is, pass the system object via the constructor of the observer, provide the return object and just call the function and evaluate the results.

这是通过 m_odefun 和 m_dxdt 字段完成的。该示例使用标准 vanderPol 振荡器，它是一个二阶 DE，给出具有二维状态向量的一阶系统。（观察者不需要模板化，这样做只是为了遵循 odeint 模板实践。）

#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
typedef boost::array<double, 2> vec;
using namespace boost::numeric::odeint;
using namespace std;

struct vanderPol
{
    double m_mu;

    vanderPol( double mu = 3 ) : m_mu ( mu ) {}

    void operator() (const vec &x, vec &dxdt, const double )
    {
        dxdt[0] = x[1];
        dxdt[1] = -m_mu*(x[0]*x[0]-1)*x[1]-x[0];
    }
};

template<class system>
struct observer
{
    vec m_dxdt;
    system m_odefun;
    observer(system odefun) : m_odefun(odefun) { }

    void operator()( vec &x, double t)
    {
        cout << "at t=" << t << "\tx: " << x[0] << '\t'  << x[1];
        m_odefun(x, m_dxdt, t);
        cout << "\tdxdt: " << m_dxdt[0] << '\t'  << m_dxdt[1] << '\n';
        
    }
};


int main()
{
    double dt = 0.01;
    double sim_time = 1.0;
    vec x = {0.0, 1.7};

    vanderPol system(1.3);

    observer< vanderPol > obs(system);

    integrate_const(runge_kutta4< vec >(), system, x, 0.0, sim_time, dt, boost::ref( obs ));
    
    return 0;
}

使用 boost odeint 返回导数

returning derivatives with boost odeint

c++

boost

ode