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 - dx/dt = f(x,t), а решение ODEs численно работает, оценивая f(x) и обновляя x с помощью численного метода. Ваш цикл прервет этот алгоритм. В деталях 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); 
} 

Вы можете ввести свой цикл непосредственно после того, как условие для попыток макс.