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

你可以在最大努力的條件後直接引入你的循環。