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)，通过数值方法评估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); 
} 

你可以在最大努力的条件后直接引入你的循环。