我有一个使用“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;
}
}
您是否足够了解该程序?您是否知道如何使它起作用?
谢谢 !
集成代码
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;
}
最佳答案
是的,此循环是无用的,因为它与解决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 );
}
您可以在最大尝试条件之后直接引入循环。
关于c++ - odeint禁止负值,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/26269359/