c++ - odeint 禁止负值

Question

我有一个使用“odeint”模拟种群动态的程序。我想设置一个 if 条件来禁止我的颂歌的结果为负。这是我的代码摘要:

class Communities
{
    public :
    typedef boost::array< double , 22 > state_type;

    Communities(...);
    ~Communities();  

    void operator()(state_type &x , state_type &dxdt , double t);
    void operator()(state_type &x , const double t );
    void integration(par::Params parParam);

    private:
    ...
};

void Communities :: operator ()( state_type &x , state_type &dxdt , double t )
{
    for (int i=0; i<nb ; ++i)
    { 
        dxdt[i] = ...
    }
    for (int j=0; j<g ; ++j)
    { 
        dxdt[j] += ...
    }

    for (int k=0; k<nb+g ; ++k)
    { 
        if (x[k] <0)
        {
            x[k] = 0;
        }
    }
 }

...

void Communities :: integration(par::Params parParam)
{
...
    integrate_adaptive(make_controlled( 1E-12 , 1E-12 , runge_kutta_fehlberg78< state_type >()) , boost::ref(*this), B , 0.0 , 10.0 , 0.1 , boost::ref(*this));
}

int main(int argc, const char * argv[])
{
    ...
    Com1.integration(ParamText); 

    return 0;
}

正如所写，下面的循环是无用的:

    for (int k=0; k<nb+g ; ++k)
    { 
        if (x[k] <0)
        {
            x[k] = 0;
        }
    }

你对这个程序有足够的了解吗？您知道我如何让它发挥作用吗？

谢谢!

integrate_adaptive的代码

template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
    Stepper stepper , System system , State &start_state ,
    Time &start_time , Time end_time , Time &dt ,
    Observer observer , controlled_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;

const size_t max_attempts = 1000;
const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
size_t count = 0;
while( less_with_sign( start_time , end_time , dt ) )
{
    obs( start_state , start_time );
    if( less_with_sign( end_time , static_cast<Time>(start_time + dt) , dt ) )
    {
        dt = end_time - start_time;
    }

    size_t trials = 0;
    controlled_step_result res = success;
    do
    {
        res = st.try_step( system , start_state , start_time , dt );
        ++trials;
    }
    while( ( res == fail ) && ( trials < max_attempts ) );
    if( trials == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ) );

    ++count;
}
obs( start_state , start_time );


return count;
}

Answer 1

是的，这个循环是无用的，因为它与求解 ODE 无关。 ODE 是 dx/dt = f(x,t) 并且通过计算 f(x) 和更新 x 以数值方式求解 ODE数值方法。你的循环打破了这个算法。详细地说，odeint 假定 x 是一个输入参数。

您需要的是一个特殊的集成例程。您可以查看 integrate_adaptive 的实现并在那里介绍您的外观。 integrate_adaptive的代码基本是

template< typename Stepper , typename System , typename State , typename Time , typename Obs >
void integrate_adaptive( Stepper stepper , System system , State &x , Time &start_time , Time end_time , Time dt , Obs obs )
{
    const size_t max_attempts = 1000;
    size_t count = 0;
    while( ( start_time + dt ) < end_time )
    {
        obs( start_state , start_time );
        if( ( start_time + dt ) > end_time  )
        {
            dt = end_time - start_time;
        }

        size_t trials = 0;
        controlled_step_result res = success;
        do
        {
            res = st.try_step( system , start_state , start_time , dt );
            ++trials;
        }
        while( ( res == fail ) && ( trials < max_attempts ) );
        if( trials == max_attempts ) throw std::overflow_error( "error" );
    }
    obs( start_state , start_time );
}

您可以在最大尝试次数条件之后直接引入循环。