odeint forbid negative values

Question

I've got a program that simulates populations dynamics using "odeint". I would like to set an if condition to forbid the result of my ode from being negative. Here is a summary of my code :

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

As written, the following loop is useless :

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

Do you have enough to understand the program ? Do you have any idea of how I can make it work ?

Thank you !

The code of 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

Yes, this loop is useless, since it has nothing to do with solving ODEs. The ODE is dx/dt = f(x,t) and solving ODEs numerically works by evaluating f(x) and updating x by the numerical method. Your loop break this algorithm. In detail, odeint assumes that x is an input parameter.

What you need is a special integrate routine. You can have a look into the implementation of integrate_adaptive and introduce your look there. The code of integrate_adaptive is basically

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

You could introduce your loop directly after the condition for the max attempts.