Changeset 81754


Ignore:
Timestamp:
Dec 7, 2012, 12:26:55 PM (6 years ago)
Author:
karsten
Message:

removing iterators from odeint, updating the documentation for odeint

Location:
trunk
Files:
1 added
9 deleted
37 edited

Legend:

Unmodified
Added
Removed
  • trunk/boost/numeric/odeint/integrate/integrate.hpp

    r81491 r81754  
    3636 * determine type of dxdt for units
    3737 *
    38  * the two overloads are needed in order to solve the forwarding problem
    3938 */
    4039template< class System , class State , class Time , class Observer >
     
    4342    return integrate_adaptive( controlled_runge_kutta< runge_kutta_dopri5< State > >() , system , start_state , start_time , end_time , dt , observer );
    4443}
    45 
    46 template< class System , class State , class Time , class Observer >
    47 size_t integrate( System system , const State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
    48 {
    49     return integrate_adaptive( controlled_runge_kutta< runge_kutta_dopri5< State > >() , system , start_state , start_time , end_time , dt , observer );
    50 }
    51 
    52 
    53 
    54 
    5544
    5645
     
    6554}
    6655
    67 template< class System , class State , class Time >
    68 size_t integrate( System system , const State &start_state , Time start_time , Time end_time , Time dt )
    69 {
    70     return integrate( system , start_state , start_time , end_time , dt , null_observer() );
    71 }
     56
     57/**
     58 * \fn integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
     59 * \brief Integrates the ODE.
     60 *
     61 * Integrates the ODE given by system from start_time to end_time starting
     62 * with start_state as initial condtition and dt as initial time step.
     63 * This function uses a dense output dopri5 stepper and performs an adaptive
     64 * integration with step size control, thus dt changes during the integration.
     65 * This method uses standard error bounds of 1E-6.
     66 * After each step, the observer is called.
     67 *
     68 * \param system The system function to solve, hence the r.h.s. of the
     69 * ordinary differential equation.
     70 * \param start_state The initial state.
     71 * \param start_time Start time of the integration.
     72 * \param end_time End time of the integration.
     73 * \param dt Initial step size, will be adjusted durinf the integration.
     74 * \param observer Observer that will be called after each time step.
     75 * \return The number of steps performed.
     76 */
    7277
    7378
    74 
     79/**
     80 * \fn integrate( System system , State &start_state , Time start_time , Time end_time , Time dt )
     81 * \brief Integrates the ODE without observer calls.
     82 *
     83 * Integrates the ODE given by system from start_time to end_time starting
     84 * with start_state as initial condtition and dt as initial time step.
     85 * This function uses a dense output dopri5 stepper and performs an adaptive
     86 * integration with step size control, thus dt changes during the integration.
     87 * This method uses standard error bounds of 1E-6.
     88 * No observer is called.
     89 *
     90 * \param system The system function to solve, hence the r.h.s. of the
     91 * ordinary differential equation.
     92 * \param start_state The initial state.
     93 * \param start_time Start time of the integration.
     94 * \param end_time End time of the integration.
     95 * \param dt Initial step size, will be adjusted durinf the integration.
     96 * \return The number of steps performed.
     97 */
    7598
    7699} // namespace odeint
  • trunk/boost/numeric/odeint/integrate/integrate_adaptive.hpp

    r81491 r81754  
    5252}
    5353
     54/**
     55 * \brief Second version to solve the forwarding problem,
     56 * can be called with Boost.Range as start_state.
     57 */
    5458template< class Stepper , class System , class State , class Time , class Observer >
    5559size_t integrate_adaptive(
     
    6771
    6872
    69 /*
    70  * the two overloads are needed in order to solve the forwarding problem
     73/**
     74 * \brief integrate_adaptive without an observer.
    7175 */
    7276template< class Stepper , class System , class State , class Time >
     
    7882}
    7983
     84/**
     85 * \brief Second version to solve the forwarding problem,
     86 * can be called with Boost.Range as start_state.
     87 */
    8088template< class Stepper , class System , class State , class Time >
    8189size_t integrate_adaptive(
     
    8795
    8896
     97/************* DOXYGEN ************/
     98
     99    /**
     100     * \fn integrate_adaptive( Stepper stepper , System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
     101     * \brief Integrates the ODE with adaptive step size.
     102     *
     103     * This function integrates the ODE given by system with the given stepper.
     104     * The observer is called after each step. If the stepper has no error
     105     * control, the step size remains constant and the observer is calleda at
     106     * equidistant time points t0+n*dt. If the stepper is a ControlledStepper,
     107     * the step size is adjusted and the observer is called in non-equidistant
     108     * intervals.
     109     *
     110     * \param stepper The stepper to be used for numerical integration.
     111     * \param system Function/Functor defining the rhs of the ODE.
     112     * \param start_state The initial condition x0.
     113     * \param start_time The initial time t0.
     114     * \param end_time The final integration time tend.
     115     * \param dt The time step between observer calls, _not_ neccessarily the
     116     * time step of the integration.
     117     * \param observer Function/Functor called at equidistant time intervals.
     118     * \return The number of steps performed.
     119     */
     120
    89121} // namespace odeint
    90122} // namespace numeric
  • trunk/boost/numeric/odeint/integrate/integrate_const.hpp

    r81491 r81754  
    6161}
    6262
    63 
     63/**
     64 * \brief Second version to solve the forwarding problem,
     65 * can be called with Boost.Range as start_state.
     66 */
    6467template< class Stepper , class System , class State , class Time , class Observer >
    6568size_t integrate_const(
     
    8992
    9093
    91 /*
    92  * Without observers
     94/**
     95 * \brief integrate_const without observer calls
    9396 */
    9497template< class Stepper , class System , class State , class Time >
     
    101104}
    102105
     106/**
     107 * \brief Second version to solve the forwarding problem,
     108 * can be called with Boost.Range as start_state.
     109 */
    103110template< class Stepper , class System , class State , class Time >
    104111size_t integrate_const(
     
    115122
    116123
     124/********* DOXYGEN *********/
     125    /**
     126     * \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
     127     * \brief Integrates the ODE with constant step size.
     128     *
     129     * Integrates the ODE defined by system using the given stepper.
     130     * This method ensures that the observer is called at constant intervals dt.
     131     * If the Stepper is a normal stepper without step size control, dt is also
     132     * used for the numerical scheme. If a ControlledStepper is provided, the
     133     * algorithm might reduce the step size to meet the error bounds, but it is
     134     * ensured that the observer is always called at equidistant time points
     135     * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
     136     * and the dense output is used to call the observer at equidistant time
     137     * points.
     138     *
     139     * \param stepper The stepper to be used for numerical integration.
     140     * \param system Function/Functor defining the rhs of the ODE.
     141     * \param start_state The initial condition x0.
     142     * \param start_time The initial time t0.
     143     * \param end_time The final integration time tend.
     144     * \param dt The time step between observer calls, _not_ neccessarily the
     145     * time step of the integration.
     146     * \param observer Function/Functor called at equidistant time intervals.
     147     * \return The number of steps performed.
     148     */
    117149
    118150} // namespace odeint
  • trunk/boost/numeric/odeint/integrate/integrate_n_steps.hpp

    r81491 r81754  
    4848}
    4949
     50/**
     51 * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
     52 */
    5053template< class Stepper , class System , class State , class Time , class Observer >
    5154Time integrate_n_steps(
     
    6164
    6265
    63 
    64 
    65 
    66 
    67 
    68 
    69 /*
    70  * the two overloads are needed in order to solve the forwarding problem
     66/**
     67 * \brief The same function as above, but without observer calls.
    7168 */
    7269template< class Stepper , class System , class State , class Time >
     
    7875}
    7976
     77/**
     78 * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
     79 */
    8080template< class Stepper , class System , class State , class Time >
    8181Time integrate_n_steps(
     
    8888
    8989
     90/************* DOXYGEN *************/
     91    /**
     92     * \fn Time integrate_n_steps( Stepper stepper , System system , State &start_state , Time start_time , Time dt , size_t num_of_steps , Observer observer )
     93     * \brief Integrates the ODE with constant step size.
     94     *
     95     * This function is similar to integrate_const. The observer is called at
     96     * equidistant time intervals t0 + n*dt.
     97     * If the Stepper is a normal stepper without step size control, dt is also
     98     * used for the numerical scheme. If a ControlledStepper is provided, the
     99     * algorithm might reduce the step size to meet the error bounds, but it is
     100     * ensured that the observer is always called at equidistant time points
     101     * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
     102     * and the dense output is used to call the observer at equidistant time
     103     * points. The final integration time is always t0 + num_of_steps*dt.
     104     *
     105     * \param stepper The stepper to be used for numerical integration.
     106     * \param system Function/Functor defining the rhs of the ODE.
     107     * \param start_state The initial condition x0.
     108     * \param start_time The initial time t0.
     109     * \param dt The time step between observer calls, _not_ neccessarily the
     110     * time step of the integration.
     111     * \param num_of_steps Number of steps to be performed
     112     * \param observer Function/Functor called at equidistant time intervals.
     113     * \return The number of steps performed.
     114     */
     115
     116
    90117
    91118} // namespace odeint
  • trunk/boost/numeric/odeint/integrate/integrate_times.hpp

    r81491 r81754  
    4747}
    4848
     49/**
     50 * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
     51 */
    4952template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
    5053size_t integrate_times(
     
    5962}
    6063
    61 /*
    62  * Range Version:
    63  * the two overloads are needed in order to solve the forwarding problem
     64/**
     65 * \brief The same function as above, but without observer calls.
    6466 */
    6567template< class Stepper , class System , class State , class TimeRange , class Time , class Observer >
     
    7476}
    7577
     78/**
     79 * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
     80 */
    7681template< class Stepper , class System , class State , class TimeRange , class Time , class Observer >
    7782size_t integrate_times(
     
    8590}
    8691
     92
     93
     94
     95/********* DOXYGEN ***********/
     96
     97    /**
     98     * \fn size_t integrate_times( Stepper stepper , System system , State &start_state , TimeIterator times_start , TimeIterator times_end , Time dt , Observer observer )
     99     * \brief Integrates the ODE with observer calls at given time points.
     100     *
     101     * Integrates the ODE given by system using the given stepper. This function
     102     * does observer calls at the subsequent time points given by the range
     103     * times_start, times_end. If the stepper has not step size control, the
     104     * step size might be reduced occasionally to ensure observer calls exactly
     105     * at the time points from the given sequence. If the stepper is a
     106     * ControlledStepper, the step size is adjusted to meet the error bounds,
     107     * but also might be reduced occasionally to ensure correct observer calls.
     108     * If a DenseOutputStepper is provided, the dense output functionality is
     109     * used to call the observer at the given times. The end time of the
     110     * integration is always *(end_time-1).
     111     *
     112     * \param stepper The stepper to be used for numerical integration.
     113     * \param system Function/Functor defining the rhs of the ODE.
     114     * \param start_state The initial condition x0.
     115     * \param times_start Iterator to the start time
     116     * \param times_end Iterator to the end time
     117     * \param dt The time step between observer calls, _not_ neccessarily the
     118     * time step of the integration.
     119     * \param observer Function/Functor called at equidistant time intervals.
     120     * \return The number of steps performed.
     121     */
     122
     123
     124
    87125} // namespace odeint
    88126} // namespace numeric
  • trunk/boost/numeric/odeint/stepper/adams_bashforth.hpp

    r81491 r81754  
    4848
    4949
    50 /*
    51  * Static explicit Adams-Bashforth multistep-solver without step size control and without dense output.
    52  */
    5350template<
    5451size_t Steps ,
     
    6461class adams_bashforth : public algebra_stepper_base< Algebra , Operations >
    6562{
     63
     64#ifndef DOXYGEN_SKIP
    6665    BOOST_STATIC_ASSERT(( Steps > 0 ));
     66    BOOST_STATIC_ASSERT(( Steps < 9 ));
     67#endif
    6768
    6869public :
     
    7778    typedef stepper_tag stepper_category;
    7879
     80    typedef InitializingStepper initializing_stepper_type;
     81
    7982    typedef typename algebra_stepper_base< Algebra , Operations >::algebra_type algebra_type;
    8083    typedef typename algebra_stepper_base< Algebra , Operations >::operations_type operations_type;
     84#ifndef DOXYGEN_SKIP
    8185    typedef adams_bashforth< Steps , State , Value , Deriv , Time , Algebra , Operations , Resizer , InitializingStepper > stepper_type;
    82     typedef InitializingStepper initializing_stepper_type;
    83 
     86#endif
    8487    static const size_t steps = Steps;
    8588
     
    9295
    9396
     97   
    9498    order_type order( void ) const { return order_value; }
    95 
    9699
    97100    adams_bashforth( const algebra_type &algebra = algebra_type() )
     
    127130    }
    128131
     132    /**
     133     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     134     */
    129135    template< class System , class StateInOut >
    130136    void do_step( System system , const StateInOut &x , time_type t , time_type dt )
     
    140146     * solves the forwarding problem
    141147     */
     148
    142149    template< class System , class StateIn , class StateOut >
    143150    void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     
    146153    }
    147154
     155    /**
     156     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateOut.
     157     */
    148158    template< class System , class StateIn , class StateOut >
    149159    void do_step( System system , const StateIn &in , time_type t , const StateOut &out , time_type dt )
     
    151161        do_step_impl( system , in , t , out , dt );
    152162    }
    153 
    154 
    155 
    156 
    157 
    158 
    159 
    160 
    161 
    162 
    163 
    164 
    165163
    166164
     
    275273
    276274
    277 
     275/***** DOXYGEN *****/
     276
     277/**
     278 * \class adams_bashforth
     279 * \brief The Adams-Bashforth multistep algorithm.
     280 *
     281 * The Adams-Bashforth method is a multi-step algorithm with configurable step
     282 * number. The step number is specified as template parameter Steps and it
     283 * then uses the result from the previous Steps steps. See also
     284 * <a href="http://en.wikipedia.org/wiki/Linear_multistep_method">en.wikipedia.org/wiki/Linear_multistep_method</a>.
     285 * Currently, a maximum of Steps=8 is supported.
     286 * The method is explicit and fullfils the Stepper concept. Step size control
     287 * or continous output are not provided.
     288 *
     289 * This class derives from algebra_base and inherits its interface via
     290 * CRTP (current recurring template pattern). For more details see
     291 * algebra_stepper_base.
     292 *
     293 * \tparam Steps The number of steps (maximal 8).
     294 * \tparam State The state type.
     295 * \tparam Value The value type.
     296 * \tparam Deriv The type representing the time derivative of the state.
     297 * \tparam Time The time representing the independent variable - the time.
     298 * \tparam Algebra The algebra type.
     299 * \tparam Operations The operations type.
     300 * \tparam Resizer The resizer policy type.
     301 * \tparam InitializingStepper The stepper for the first two steps.
     302 */
     303
     304    /**
     305     * \fn adams_bashforth::adams_bashforth( const algebra_type &algebra )
     306     * \brief Constructs the adams_bashforth class. This constructor can be used as a default
     307     * constructor if the algebra has a default constructor.
     308     * \param algebra A copy of algebra is made and stored.
     309     */
     310
     311    /**
     312     * \fn order_type adams_bashforth::order( void ) const
     313     * \brief Returns the order of the algorithm, which is equal to the number of steps.
     314     * \return order of the method.
     315     */
     316
     317    /**
     318     * \fn void adams_bashforth::do_step( System system , StateInOut &x , time_type t , time_type dt )
     319     * \brief This method performs one step. It transforms the result in-place.
     320     *
     321     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     322     *               Simple System concept.
     323     * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
     324     * \param t The value of the time, at which the step should be performed.
     325     * \param dt The step size.
     326     */
     327
     328    /**
     329     * \fn void adams_bashforth::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     330     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
     331     *
     332     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     333     *               Simple System concept.
     334     * \param in The state of the ODE which should be solved. in is not modified in this method
     335     * \param t The value of the time, at which the step should be performed.
     336     * \param out The result of the step is written in out.
     337     * \param dt The step size.
     338     */
     339
     340    /**
     341     * \fn void adams_bashforth::adjust_size( const StateType &x )
     342     * \brief Adjust the size of all temporaries in the stepper manually.
     343     * \param x A state from which the size of the temporaries to be resized is deduced.
     344     */
     345
     346
     347    /**
     348     * \fn const step_storage_type& adams_bashforth::step_storage( void ) const
     349     * \brief Returns the storage of intermediate results.
     350     * \return The storage of intermediate results.
     351     */
     352
     353    /**
     354     * \fn step_storage_type& adams_bashforth::step_storage( void )
     355     * \brief Returns the storage of intermediate results.
     356     * \return The storage of intermediate results.
     357     */
     358
     359    /**
     360     * \fn void adams_bashforth::initialize( ExplicitStepper explicit_stepper , System system , StateIn &x , time_type &t , time_type dt )
     361     * \brief Initialized the stepper. Does Steps-1 steps with the explicit_stepper to fill the buffer.
     362     * \param explicit_stepper the stepper used to fill the buffer of previous step results
     363     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     364     *               Simple System concept.
     365     * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
     366     * \param t The value of the time, at which the step should be performed.
     367     * \param dt The step size.
     368     */
     369
     370    /**
     371     * \fn void adams_bashforth::initialize( System system , StateIn &x , time_type &t , time_type dt )
     372     * \brief Initialized the stepper. Does Steps-1 steps with an internal instance of InitializingStepper to fill the buffer.
     373     * \note The state x and time t are updated to the values after Steps-1 initial steps.
     374     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     375     *               Simple System concept.
     376     * \param x The initial state of the ODE which should be solved, updated in this method.
     377     * \param t The initial value of the time, updated in this method.
     378     * \param dt The step size.
     379     */
     380
     381    /**
     382     * \fn void adams_bashforth::reset( void )
     383     * \brief Resets the internal buffer of the stepper.
     384     */
     385
     386    /**
     387     * \fn bool adams_bashforth::is_initialized( void ) const
     388     * \brief Returns true if the stepper has been initialized.
     389     * \return bool true if stepper is initialized, false otherwise
     390     */
     391
     392    /**
     393     * \fn const initializing_stepper_type& adams_bashforth::initializing_stepper( void ) const
     394     * \brief Returns the internal initializing stepper instance.
     395     * \return initializing_stepper
     396     */
     397
     398    /**
     399     * \fn const initializing_stepper_type& adams_bashforth::initializing_stepper( void ) const
     400     * \brief Returns the internal initializing stepper instance.
     401     * \return initializing_stepper
     402     */
     403
     404    /**
     405     * \fn initializing_stepper_type& adams_bashforth::initializing_stepper( void )
     406     * \brief Returns the internal initializing stepper instance.
     407     * \return initializing_stepper
     408     */
    278409
    279410} // odeint
  • trunk/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp

    r81491 r81754  
    3838
    3939
    40 /*
    41  * Static Adams-Bashforth-Moulton multistep-solver without step size control and without dense output.
    42  */
    4340template<
    4441size_t Steps ,
     
    5451{
    5552
     53#ifndef DOXYGEN_SKIP
     54    BOOST_STATIC_ASSERT(( Steps > 0 ));
     55    BOOST_STATIC_ASSERT(( Steps < 9 ));
     56#endif
     57
    5658public :
    5759
     
    6870
    6971    static const size_t steps = Steps;
    70 
     72#ifndef DOXYGEN_SKIP
    7173    typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_bashforth_type;
    7274    typedef adams_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_moulton_type;
    73 
     75#endif //DOXYGEN_SKIP
    7476    typedef unsigned short order_type;
    7577    static const order_type order_value = steps + 1;
    7678
     79    /** \brief Constructs the adams_bashforth class. */
    7780    adams_bashforth_moulton( void )
    7881    : m_adams_bashforth() , m_adams_moulton( m_adams_bashforth.algebra() )
     
    9295    }
    9396
     97    /**
     98     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     99     */
    94100    template< class System , class StateInOut >
    95101    void do_step( System system , const StateInOut &x , time_type t , time_type dt )
     
    106112    }
    107113
     114    /**
     115     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateOut.
     116     */
    108117    template< class System , class StateIn , class StateOut >
    109118    void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     
    113122    }
    114123
     124
    115125    template< class StateType >
    116126    void adjust_size( const StateType &x )
     
    126136        m_adams_bashforth.initialize( explicit_stepper , system , x , t , dt );
    127137    }
     138
    128139
    129140    template< class System , class StateIn >
     
    142153
    143154
     155/********* DOXYGEN ********/
     156
     157/**
     158 * \class adams_bashforth_moulton
     159 * \brief The Adams-Bashforth-Moulton multistep algorithm.
     160 *
     161 * The Adams-Bashforth method is a multi-step predictor-corrector algorithm
     162 * with configurable step number. The step number is specified as template
     163 * parameter Steps and it then uses the result from the previous Steps steps.
     164 * See also
     165 * <a href="http://en.wikipedia.org/wiki/Linear_multistep_method">en.wikipedia.org/wiki/Linear_multistep_method</a>.
     166 * Currently, a maximum of Steps=8 is supported.
     167 * The method is explicit and fullfils the Stepper concept. Step size control
     168 * or continous output are not provided.
     169 *
     170 * This class derives from algebra_base and inherits its interface via
     171 * CRTP (current recurring template pattern). For more details see
     172 * algebra_stepper_base.
     173 *
     174 * \tparam Steps The number of steps (maximal 8).
     175 * \tparam State The state type.
     176 * \tparam Value The value type.
     177 * \tparam Deriv The type representing the time derivative of the state.
     178 * \tparam Time The time representing the independent variable - the time.
     179 * \tparam Algebra The algebra type.
     180 * \tparam Operations The operations type.
     181 * \tparam Resizer The resizer policy type.
     182 * \tparam InitializingStepper The stepper for the first two steps.
     183 */
     184
     185    /**
     186     * \fn adams_bashforth_moulton::adams_bashforth_moulton( const algebra_type &algebra )
     187     * \brief Constructs the adams_bashforth class. This constructor can be used as a default
     188     * constructor if the algebra has a default constructor.
     189     * \param algebra A copy of algebra is made and stored.
     190     */
     191
     192    /**
     193     * \fn adams_bashforth_moulton::order( void ) const
     194     * \brief Returns the order of the algorithm, which is equal to the number of steps+1.
     195     * \return order of the method.
     196     */
     197
     198    /**
     199     * \fn adams_bashforth_moulton::do_step( System system , StateInOut &x , time_type t , time_type dt )
     200     * \brief This method performs one step. It transforms the result in-place.
     201     *
     202     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     203     *               Simple System concept.
     204     * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
     205     * \param t The value of the time, at which the step should be performed.
     206     * \param dt The step size.
     207     */
     208
     209
     210    /**
     211     * \fn adams_bashforth_moulton::do_step( System system , const StateIn &in , time_type t , const StateOut &out , time_type dt )
     212     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
     213     *
     214     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     215     *               Simple System concept.
     216     * \param in The state of the ODE which should be solved. in is not modified in this method
     217     * \param t The value of the time, at which the step should be performed.
     218     * \param out The result of the step is written in out.
     219     * \param dt The step size.
     220     */
     221
     222    /**
     223     * \fn adams_bashforth_moulton::adjust_size( const StateType &x )
     224     * \brief Adjust the size of all temporaries in the stepper manually.
     225     * \param x A state from which the size of the temporaries to be resized is deduced.
     226     */
     227
     228    /**
     229     * \fn adams_bashforth_moulton::initialize( ExplicitStepper explicit_stepper , System system , StateIn &x , time_type &t , time_type dt )
     230     * \brief Initialized the stepper. Does Steps-1 steps with the explicit_stepper to fill the buffer.
     231     * \note The state x and time t are updated to the values after Steps-1 initial steps.
     232     * \param explicit_stepper the stepper used to fill the buffer of previous step results
     233     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     234     *               Simple System concept.
     235     * \param x The initial state of the ODE which should be solved, updated after in this method.
     236     * \param t The initial time, updated in this method.
     237     * \param dt The step size.
     238     */
     239
     240    /**
     241     * \fn adams_bashforth_moulton::initialize( System system , StateIn &x , time_type &t , time_type dt )
     242     * \brief Initialized the stepper. Does Steps-1 steps using the standard initializing stepper
     243     * of the underlying adams_bashforth stepper.
     244     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     245     *               Simple System concept.
     246     * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
     247     * \param t The value of the time, at which the step should be performed.
     248     * \param dt The step size.
     249     */
    144250
    145251
  • trunk/boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp

    r81491 r81754  
    2424namespace odeint {
    2525
     26template< class Algebra , class Operations >
     27class algebra_stepper_base
     28{
     29public:
     30
     31    typedef Algebra algebra_type;
     32    typedef Operations operations_type;
     33
     34    algebra_stepper_base( const algebra_type &algebra = algebra_type() )
     35    : m_algebra( algebra ) { }
     36
     37    algebra_type& algebra()
     38    {
     39        return m_algebra;
     40    }
     41
     42    const algebra_type& algebra() const
     43    {
     44        return m_algebra;
     45    }
     46
     47protected:
     48
     49    algebra_type m_algebra;
     50};
     51
     52
     53/******* DOXYGEN *******/
     54
    2655/**
    2756 * \class algebra_stepper_base
     
    3766 * to work with the stepper.
    3867 */
    39 template< class Algebra , class Operations >
    40 class algebra_stepper_base
    41 {
    42 public:
    43 
    44     typedef Algebra algebra_type;
    45     typedef Operations operations_type;
    4668
    4769    /**
     70     * \fn algebra_stepper_base::algebra_stepper_base( const algebra_type &algebra = algebra_type() )
    4871     * \brief Constructs a algebra_stepper_base and creates the algebra. This constructor can be used as a default
    4972     * constructor if the algebra has a default constructor.
    5073     * \param algebra The algebra_stepper_base stores and uses a copy of algebra.
    5174     */
    52     algebra_stepper_base( const algebra_type &algebra = algebra_type() )
    53     : m_algebra( algebra ) { }
    5475
    5576    /**
     77     * \fn algebra_type& algebra_stepper_base::algebra()
    5678     * \return A reference to the algebra which is held by this class.
    5779     */
    58     algebra_type& algebra()
    59     {
    60         return m_algebra;
    61     }
    6280
    6381    /**
     82     * \fn const algebra_type& algebra_stepper_base::algebra() const
    6483     * \return A const reference to the algebra which is held by this class.
    6584     */
    66     const algebra_type& algebra() const
    67     {
    68         return m_algebra;
    69     }
    70 
    71 protected:
    72 
    73     algebra_type m_algebra;
    74 };
    7585
    7686} // odeint
  • trunk/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp

    r81491 r81754  
    5252    * do_step( sys , in , dxdt , t , out , dt , xerr )
    5353 */
     54template<
     55class Stepper ,
     56unsigned short Order ,
     57unsigned short StepperOrder ,
     58unsigned short ErrorOrder ,
     59class State ,
     60class Value ,
     61class Deriv ,
     62class Time ,
     63class Algebra ,
     64class Operations ,
     65class Resizer
     66>
     67class explicit_error_stepper_base : public algebra_stepper_base< Algebra , Operations >
     68{
     69public:
     70
     71    typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
     72    typedef typename algebra_stepper_base_type::algebra_type algebra_type;
     73
     74
     75    typedef State state_type;
     76    typedef Value value_type;
     77    typedef Deriv deriv_type;
     78    typedef Time time_type;
     79    typedef Resizer resizer_type;
     80    typedef Stepper stepper_type;
     81    typedef explicit_error_stepper_tag stepper_category;
     82    #ifndef DOXYGEN_SKIP
     83    typedef state_wrapper< state_type > wrapped_state_type;
     84    typedef state_wrapper< deriv_type > wrapped_deriv_type;
     85    typedef explicit_error_stepper_base< Stepper , Order , StepperOrder , ErrorOrder ,
     86            State , Value , Deriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
     87    #endif
     88
     89    typedef unsigned short order_type;
     90    static const order_type order_value = Order;
     91    static const order_type stepper_order_value = StepperOrder;
     92    static const order_type error_order_value = ErrorOrder;
     93
     94
     95    explicit_error_stepper_base( const algebra_type &algebra = algebra_type() )
     96    : algebra_stepper_base_type( algebra )
     97    { }
     98
     99    order_type order( void ) const
     100    {
     101        return order_value;
     102    }
     103
     104    order_type stepper_order( void ) const
     105    {
     106        return stepper_order_value;
     107    }
     108
     109    order_type error_order( void ) const
     110    {
     111        return error_order_value;
     112    }
     113
     114
     115
     116    /*
     117     * Version 1 : do_step( sys , x , t , dt )
     118     *
     119     * the two overloads are needed in order to solve the forwarding problem
     120     */
     121    template< class System , class StateInOut >
     122    void do_step( System system , StateInOut &x , time_type t , time_type dt )
     123    {
     124        do_step_v1( system , x , t , dt );
     125    }
     126
     127    /**
     128     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     129     */
     130    template< class System , class StateInOut >
     131    void do_step( System system , const StateInOut &x , time_type t , time_type dt )
     132    {
     133        do_step_v1( system , x , t , dt );
     134    }
     135
     136
     137
     138    /*
     139     * Version 2 : do_step( sys , x , dxdt , t , dt )
     140     *
     141     * this version does not solve the forwarding problem, boost.range can not be used
     142     *
     143     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     144     */
     145    template< class System , class StateInOut , class DerivIn >
     146    typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
     147    do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt )
     148    {
     149        this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
     150    }
     151
     152
     153    /*
     154     * Version 3 : do_step( sys , in , t , out , dt )
     155     *
     156     * this version does not solve the forwarding problem, boost.range can not be used
     157     *
     158     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     159     */
     160    template< class System , class StateIn , class StateOut >
     161    typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
     162    do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     163    {
     164        typename odeint::unwrap_reference< System >::type &sys = system;
     165        m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
     166        sys( in , m_dxdt.m_v ,t );
     167        this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , dt );
     168    }
     169
     170    /*
     171     * Version 4 :do_step( sys , in , dxdt , t , out , dt )
     172     *
     173     * this version does not solve the forwarding problem, boost.range can not be used
     174     *
     175     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     176     */
     177    template< class System , class StateIn , class DerivIn , class StateOut >
     178    typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
     179    do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
     180    {
     181        this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
     182    }
     183
     184
     185
     186
     187
     188    /*
     189     * Version  5 :do_step( sys , x , t , dt , xerr )
     190     *
     191     * the two overloads are needed in order to solve the forwarding problem
     192     */
     193    template< class System , class StateInOut , class Err >
     194    void do_step( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
     195    {
     196        do_step_v5( system , x , t , dt , xerr );
     197    }
     198
     199    /**
     200     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     201     */
     202    template< class System , class StateInOut , class Err >
     203    void do_step( System system , const StateInOut &x , time_type t , time_type dt , Err &xerr )
     204    {
     205        do_step_v5( system , x , t , dt , xerr );
     206    }
     207
     208
     209    /*
     210     * Version 6 :do_step( sys , x , dxdt , t , dt , xerr )
     211     *
     212     * this version does not solve the forwarding problem, boost.range can not be used
     213     *
     214     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     215     */
     216    template< class System , class StateInOut , class DerivIn , class Err >
     217    typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
     218    do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt , Err &xerr )
     219    {
     220        this->stepper().do_step_impl( system , x , dxdt , t , x , dt , xerr );
     221    }
     222
     223
     224    /*
     225     * Version 7 : do_step( sys , in , t , out , dt , xerr )
     226     *
     227     * this version does not solve the forwarding problem, boost.range can not be used
     228     */
     229    template< class System , class StateIn , class StateOut , class Err >
     230    void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , Err &xerr )
     231    {
     232        typename odeint::unwrap_reference< System >::type &sys = system;
     233        m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
     234        sys( in , m_dxdt.m_v ,t );
     235        this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , dt , xerr );
     236    }
     237
     238
     239    /*
     240     * Version 8 : do_step( sys , in , dxdt , t , out , dt , xerr )
     241     *
     242     * this version does not solve the forwarding problem, boost.range can not be used
     243     */
     244    template< class System , class StateIn , class DerivIn , class StateOut , class Err >
     245    void do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr )
     246    {
     247        this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr );
     248    }
     249
     250    template< class StateIn >
     251    void adjust_size( const StateIn &x )
     252    {
     253        resize_impl( x );
     254    }
     255
     256
     257
     258private:
     259
     260    template< class System , class StateInOut >
     261    void do_step_v1( System system , StateInOut &x , time_type t , time_type dt )
     262    {
     263        typename odeint::unwrap_reference< System >::type &sys = system;
     264        m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl<StateInOut> , detail::ref( *this ) , detail::_1 ) );
     265        sys( x , m_dxdt.m_v ,t );
     266        this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , dt );
     267    }
     268
     269    template< class System , class StateInOut , class Err >
     270    void do_step_v5( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
     271    {
     272        typename odeint::unwrap_reference< System >::type &sys = system;
     273        m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl<StateInOut> , detail::ref( *this ) , detail::_1 ) );
     274        sys( x , m_dxdt.m_v ,t );
     275        this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , dt , xerr );
     276    }
     277
     278    template< class StateIn >
     279    bool resize_impl( const StateIn &x )
     280    {
     281        return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
     282    }
     283
     284    stepper_type& stepper( void )
     285    {
     286        return *static_cast< stepper_type* >( this );
     287    }
     288
     289    const stepper_type& stepper( void ) const
     290    {
     291        return *static_cast< const stepper_type* >( this );
     292    }
     293
     294
     295    resizer_type m_resizer;
     296
     297protected:
     298
     299    wrapped_deriv_type m_dxdt;
     300};
     301
     302
     303
     304
     305/******** DOXYGEN *******/
     306
    54307/**
    55308 * \class explicit_error_stepper_base
     
    118371 * \tparam Resizer The resizer policy class.
    119372 */
    120 template<
    121 class Stepper ,
    122 unsigned short Order ,
    123 unsigned short StepperOrder ,
    124 unsigned short ErrorOrder ,
    125 class State ,
    126 class Value ,
    127 class Deriv ,
    128 class Time ,
    129 class Algebra ,
    130 class Operations ,
    131 class Resizer
    132 >
    133 class explicit_error_stepper_base : public algebra_stepper_base< Algebra , Operations >
    134 {
    135 public:
    136 
    137     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
    138     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
    139 
    140 
    141     typedef State state_type;
    142     typedef Value value_type;
    143     typedef Deriv deriv_type;
    144     typedef Time time_type;
    145     typedef Resizer resizer_type;
    146     typedef Stepper stepper_type;
    147     typedef explicit_error_stepper_tag stepper_category;
    148     #ifndef DOXYGEN_SKIP
    149     typedef state_wrapper< state_type > wrapped_state_type;
    150     typedef state_wrapper< deriv_type > wrapped_deriv_type;
    151     typedef explicit_error_stepper_base< Stepper , Order , StepperOrder , ErrorOrder ,
    152             State , Value , Deriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
    153     #endif
    154 
    155     typedef unsigned short order_type;
    156     static const order_type order_value = Order;
    157     static const order_type stepper_order_value = StepperOrder;
    158     static const order_type error_order_value = ErrorOrder;
    159 
    160 
    161     /**
    162      * \brief Constructs a explicit_stepper_base class. This constructor can be used as a default
     373
     374
     375    /**
     376     * \fn explicit_error_stepper_base::explicit_error_stepper_base( const algebra_type &algebra = algebra_type() )
     377     *
     378     * \brief Constructs a explicit_error_stepper_base class. This constructor can be used as a default
    163379     * constructor if the algebra has a default constructor.
    164380     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    165381     */
    166     explicit_error_stepper_base( const algebra_type &algebra = algebra_type() )
    167     : algebra_stepper_base_type( algebra )
    168     { }
    169 
    170     /**
     382
     383    /**
     384     * \fn explicit_error_stepper_base::order( void ) const
    171385     * \return Returns the order of the stepper if it used without error estimation.
    172386     */
    173     order_type order( void ) const
    174     {
    175         return order_value;
    176     }
    177 
    178     /**
     387
     388    /**
     389     * \fn explicit_error_stepper_base::stepper_order( void ) const
    179390     * \return Returns the order of a step if the stepper is used without error estimation.
    180391     */
    181     order_type stepper_order( void ) const
    182     {
    183         return stepper_order_value;
    184     }
    185 
    186     /**
     392
     393    /**
     394     * \fn explicit_error_stepper_base::error_order( void ) const
    187395     * \return Returns the order of an error step if the stepper is used without error estimation.
    188396     */
    189     order_type error_order( void ) const
    190     {
    191         return error_order_value;
    192     }
    193 
    194 
    195 
    196     /*
    197      * Version 1 : do_step( sys , x , t , dt )
    198      *
    199      * the two overloads are needed in order to solve the forwarding problem
    200      */
    201     /**
     397
     398    /**
     399     * \fn explicit_error_stepper_base::do_step( System system , StateInOut &x , time_type t , time_type dt )
    202400     * \brief This method performs one step. It transforms the result in-place.
    203401     *
     
    208406     * \param dt The step size.
    209407     */
    210     template< class System , class StateInOut >
    211     void do_step( System system , StateInOut &x , time_type t , time_type dt )
    212     {
    213         do_step_v1( system , x , t , dt );
    214     }
    215 
    216     /**
    217      * \brief This method performs one step with the stepper passed by Stepper.
    218      * It transforms the result in-place. This method is needed in order to solve the forwarding problem.
    219      * The difference to the other version is that it can be used like
    220      * `stepper.do_step( sys , make_range( iter1 , iter2 ) , t , dt )`
    221      *
    222      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
    223      *               Simple System concept.
    224      * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
    225      * \param t The value of the time, at which the step should be performed.
    226      * \param dt The step size.
    227      */
    228     template< class System , class StateInOut >
    229     void do_step( System system , const StateInOut &x , time_type t , time_type dt )
    230     {
    231         do_step_v1( system , x , t , dt );
    232     }
    233 
    234 
    235 
    236     /*
    237      * Version 2 : do_step( sys , x , dxdt , t , dt )
    238      *
    239      * this version does not solve the forwarding problem, boost.range can not be used
    240      *
    241      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    242      */
    243     /**
     408
     409    /**
     410     * \fn explicit_error_stepper_base::do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt )
    244411     * \brief The method performs one step with the stepper passed by Stepper. Additionally to the other method
    245      * the derivative of x is also passed to this method. It is equivalent to
     412     * the derivative of x is also passed to this method. It is supposed to be used in the following way:
    246413     *
    247414     * \code
     
    262429     * \param dt The step size.
    263430     */
    264     template< class System , class StateInOut , class DerivIn >
    265     typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
    266     do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt )
    267     {
    268         this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
    269     }
    270 
    271 
    272     /*
    273      * Version 3 : do_step( sys , in , t , out , dt )
    274      *
    275      * this version does not solve the forwarding problem, boost.range can not be used
    276      *
    277      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    278      */
    279     /**
     431
     432    /**
     433     * \fn explicit_error_stepper_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    280434     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    281435     * This method is disabled if StateIn and Time are the same type. In this case the method can not be distinguished from
     
    290444     * \param dt The step size.
    291445     */
    292     template< class System , class StateIn , class StateOut >
    293     typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
    294     do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    295     {
    296         typename odeint::unwrap_reference< System >::type &sys = system;
    297         m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
    298         sys( in , m_dxdt.m_v ,t );
    299         this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , dt );
    300     }
    301 
    302     /*
    303      * Version 4 :do_step( sys , in , dxdt , t , out , dt )
    304      *
    305      * this version does not solve the forwarding problem, boost.range can not be used
    306      *
    307      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    308      */
    309     /**
     446
     447
     448    /**
     449     * \fn explicit_error_stepper_base::do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    310450     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    311      * Furthermore, the derivative of x at t is passed to the stepper. It is equivalent to:
     451     * Furthermore, the derivative of x at t is passed to the stepper. It is supposed to be used in the following way:
    312452     *
    313453     * \code
     
    328468     * \param dt The step size.
    329469     */
    330     template< class System , class StateIn , class DerivIn , class StateOut >
    331     typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
    332     do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    333     {
    334         this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
    335     }
    336 
    337 
    338 
    339 
    340 
    341     /*
    342      * Version  5 :do_step( sys , x , t , dt , xerr )
    343      *
    344      * the two overloads are needed in order to solve the forwarding problem
    345      */
    346     /**
     470
     471    /**
     472     * \fn explicit_error_stepper_base::do_step( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
    347473     * \brief The method performs one step with the stepper passed by Stepper and estimates the error. The state of the ODE
    348474     * is updated in-place.
     
    355481     * \param xerr The estimation of the error is stored in xerr.
    356482     */
    357     template< class System , class StateInOut , class Err >
    358     void do_step( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
    359     {
    360         do_step_v5( system , x , t , dt , xerr );
    361     }
    362 
    363     /**
    364      * \brief The method performs one step with the stepper passed by Stepper and estimates the error. The state of the ODE
    365      * is updated in-place. This method is needed in order to solve the forwarding problem.
    366      * The difference to the other version is that it can be used like
    367      * `stepper.do_step( sys , make_range( iter1 , iter2 ) , t , dt )`
    368      *
    369      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
    370      *               Simple System concept.
    371      * \param x The state of the ODE which should be solved. x is updated by this method.
    372      * \param t The value of the time, at which the step should be performed.
    373      * \param dt The step size.
    374      * \param xerr The estimation of the error is stored in xerr.
    375      */
    376     template< class System , class StateInOut , class Err >
    377     void do_step( System system , const StateInOut &x , time_type t , time_type dt , Err &xerr )
    378     {
    379         do_step_v5( system , x , t , dt , xerr );
    380     }
    381 
    382 
    383     /*
    384      * Version 6 :do_step( sys , x , dxdt , t , dt , xerr )
    385      *
    386      * this version does not solve the forwarding problem, boost.range can not be used
    387      *
    388      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    389      */
    390     /**
     483
     484    /**
     485     * \fn explicit_error_stepper_base::do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt , Err &xerr )
    391486     * \brief The method performs one step with the stepper passed by Stepper. Additionally to the other method
    392      * the derivative of x is also passed to this method. It is equivalent to
     487     * the derivative of x is also passed to this method. It is supposed to be used in the following way:
    393488     *
    394489     * \code
     
    410505     * \param xerr The error estimate is stored in xerr.
    411506     */
    412     template< class System , class StateInOut , class DerivIn , class Err >
    413     typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
    414     do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt , Err &xerr )
    415     {
    416         this->stepper().do_step_impl( system , x , dxdt , t , x , dt , xerr );
    417     }
    418 
    419 
    420     /*
    421      * Version 7 : do_step( sys , in , t , out , dt , xerr )
    422      *
    423      * this version does not solve the forwarding problem, boost.range can not be used
    424      */
    425     /**
     507
     508    /**
     509     * \fn explicit_error_stepper_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , Err &xerr )
    426510     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    427511     * Furthermore, the error is estimated.
     
    437521     * \param xerr The error estimate.
    438522     */
    439     template< class System , class StateIn , class StateOut , class Err >
    440     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , Err &xerr )
    441     {
    442         typename odeint::unwrap_reference< System >::type &sys = system;
    443         m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
    444         sys( in , m_dxdt.m_v ,t );
    445         this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , dt , xerr );
    446     }
    447 
    448 
    449     /*
    450      * Version 8 : do_step( sys , in , dxdt , t , out , dt , xerr )
    451      *
    452      * this version does not solve the forwarding problem, boost.range can not be used
    453      */
    454     /**
     523
     524
     525    /**
     526     * \fn explicit_error_stepper_base::do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr )
    455527     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    456      * Furthermore, the derivative of x at t is passed to the stepper and the error is estimated. It is equivalent to:
     528     * Furthermore, the derivative of x at t is passed to the stepper and the error is estimated. It is supposed to be used in the following way:
    457529     *
    458530     * \code
     
    474546     * \param xerr The error estimate.
    475547     */
    476     template< class System , class StateIn , class DerivIn , class StateOut , class Err >
    477     void do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr )
    478     {
    479         this->stepper().do_step_impl( system , in , dxdt , t , out , dt , xerr );
    480     }
    481 
    482     /**
     548
     549
     550    /**
     551     * \fn explicit_error_stepper_base::adjust_size( const StateIn &x )
    483552     * \brief Adjust the size of all temporaries in the stepper manually.
    484553     * \param x A state from which the size of the temporaries to be resized is deduced.
    485554     */
    486     template< class StateIn >
    487     void adjust_size( const StateIn &x )
    488     {
    489         resize_impl( x );
    490     }
    491 
    492 
    493 
    494 private:
    495 
    496     template< class System , class StateInOut >
    497     void do_step_v1( System system , StateInOut &x , time_type t , time_type dt )
    498     {
    499         typename odeint::unwrap_reference< System >::type &sys = system;
    500         m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl<StateInOut> , detail::ref( *this ) , detail::_1 ) );
    501         sys( x , m_dxdt.m_v ,t );
    502         this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , dt );
    503     }
    504 
    505     template< class System , class StateInOut , class Err >
    506     void do_step_v5( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
    507     {
    508         typename odeint::unwrap_reference< System >::type &sys = system;
    509         m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl<StateInOut> , detail::ref( *this ) , detail::_1 ) );
    510         sys( x , m_dxdt.m_v ,t );
    511         this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , dt , xerr );
    512     }
    513 
    514     template< class StateIn >
    515     bool resize_impl( const StateIn &x )
    516     {
    517         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
    518     }
    519 
    520     stepper_type& stepper( void )
    521     {
    522         return *static_cast< stepper_type* >( this );
    523     }
    524 
    525     const stepper_type& stepper( void ) const
    526     {
    527         return *static_cast< const stepper_type* >( this );
    528     }
    529 
    530 
    531     resizer_type m_resizer;
    532 
    533 protected:
    534 
    535     wrapped_deriv_type m_dxdt;
    536 };
    537 
    538555
    539556} // odeint
  • trunk/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp

    r81491 r81754  
    5252    * do_step( sys , in , dxdt_in , t , out , dxdt_out , dt , xerr )
    5353 */
     54template<
     55class Stepper ,
     56unsigned short Order ,
     57unsigned short StepperOrder ,
     58unsigned short ErrorOrder ,
     59class State ,
     60class Value ,
     61class Deriv ,
     62class Time ,
     63class Algebra ,
     64class Operations ,
     65class Resizer
     66>
     67class explicit_error_stepper_fsal_base : public algebra_stepper_base< Algebra , Operations >
     68{
     69public:
     70
     71    typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
     72    typedef typename algebra_stepper_base_type::algebra_type algebra_type;
     73
     74    typedef State state_type;
     75    typedef Value value_type;
     76    typedef Deriv deriv_type;
     77    typedef Time time_type;
     78    typedef Resizer resizer_type;
     79    typedef Stepper stepper_type;
     80    typedef explicit_error_stepper_fsal_tag stepper_category;
     81
     82    #ifndef DOXYGEN_SKIP
     83    typedef state_wrapper< state_type > wrapped_state_type;
     84    typedef state_wrapper< deriv_type > wrapped_deriv_type;
     85    typedef explicit_error_stepper_fsal_base< Stepper , Order , StepperOrder , ErrorOrder ,
     86            State , Value , Deriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
     87    #endif
     88
     89
     90    typedef unsigned short order_type;
     91    static const order_type order_value = Order;
     92    static const order_type stepper_order_value = StepperOrder;
     93    static const order_type error_order_value = ErrorOrder;
     94
     95    explicit_error_stepper_fsal_base( const algebra_type &algebra = algebra_type() )
     96    : algebra_stepper_base_type( algebra ) , m_first_call( true )
     97    { }
     98
     99    order_type order( void ) const
     100    {
     101        return order_value;
     102    }
     103
     104    order_type stepper_order( void ) const
     105    {
     106        return stepper_order_value;
     107    }
     108
     109    order_type error_order( void ) const
     110    {
     111        return error_order_value;
     112    }
     113
     114
     115    /*
     116     * version 1 : do_step( sys , x , t , dt )
     117     *
     118     * the two overloads are needed in order to solve the forwarding problem
     119     */
     120    template< class System , class StateInOut >
     121    void do_step( System system , StateInOut &x , time_type t , time_type dt )
     122    {
     123        do_step_v1( system , x , t , dt );
     124    }
     125
     126    /**
     127     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     128     */
     129    template< class System , class StateInOut >
     130    void do_step( System system , const StateInOut &x , time_type t , time_type dt )
     131    {
     132        do_step_v1( system , x , t , dt );
     133    }
     134
     135
     136    /*
     137     * version 2 : do_step( sys , x , dxdt , t , dt )
     138     *
     139     * this version does not solve the forwarding problem, boost.range can not be used
     140     *
     141     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     142     */
     143    template< class System , class StateInOut , class DerivInOut >
     144    typename boost::disable_if< boost::is_same< StateInOut , time_type > , void >::type
     145    do_step( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt )
     146    {
     147        m_first_call = true;
     148        this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt );
     149    }
     150
     151
     152    /*
     153     * version 3 : do_step( sys , in , t , out , dt )
     154     *
     155     * this version does not solve the forwarding problem, boost.range can not be used
     156     *
     157     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     158     */
     159    template< class System , class StateIn , class StateOut >
     160    typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
     161    do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     162    {
     163        if( m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
     164        {
     165            initialize( system , in , t );
     166        }
     167        this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , m_dxdt.m_v , dt );
     168    }
     169
     170
     171    /*
     172     * version 4 : do_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
     173     *
     174     * this version does not solve the forwarding problem, boost.range can not be used
     175     */
     176    template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
     177    void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
     178            StateOut &out , DerivOut &dxdt_out , time_type dt )
     179    {
     180        m_first_call = true;
     181        this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
     182    }
     183
     184
     185
     186
     187
     188    /*
     189     * version 5 : do_step( sys , x , t , dt , xerr )
     190     *
     191     * the two overloads are needed in order to solve the forwarding problem
     192     */
     193    template< class System , class StateInOut , class Err >
     194    void do_step( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
     195    {
     196        do_step_v5( system , x , t , dt , xerr );
     197    }
     198
     199    /**
     200     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     201     */
     202    template< class System , class StateInOut , class Err >
     203    void do_step( System system , const StateInOut &x , time_type t , time_type dt , Err &xerr )
     204    {
     205        do_step_v5( system , x , t , dt , xerr );
     206    }
     207
     208
     209    /*
     210     * version 6 : do_step( sys , x , dxdt , t , dt , xerr )
     211     *
     212     * this version does not solve the forwarding problem, boost.range can not be used
     213     *
     214     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     215     */
     216    template< class System , class StateInOut , class DerivInOut , class Err >
     217    typename boost::disable_if< boost::is_same< StateInOut , time_type > , void >::type
     218    do_step( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt , Err &xerr )
     219    {
     220        m_first_call = true;
     221        this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt , xerr );
     222    }
     223
     224
     225
     226
     227    /*
     228     * version 7 : do_step( sys , in , t , out , dt , xerr )
     229     *
     230     * this version does not solve the forwarding problem, boost.range can not be used
     231     */
     232    template< class System , class StateIn , class StateOut , class Err >
     233    void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , Err &xerr )
     234    {
     235        if( m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
     236        {
     237            initialize( system , in , t );
     238        }
     239        this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , m_dxdt.m_v , dt , xerr );
     240    }
     241
     242
     243    /*
     244     * version 8 : do_step( sys , in , dxdt_in , t , out , dxdt_out , dt , xerr )
     245     *
     246     * this version does not solve the forwarding problem, boost.range can not be used
     247     */
     248    template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut , class Err >
     249    void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
     250            StateOut &out , DerivOut &dxdt_out , time_type dt , Err &xerr )
     251    {
     252        m_first_call = true;
     253        this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt , xerr );
     254    }
     255
     256    template< class StateIn >
     257    void adjust_size( const StateIn &x )
     258    {
     259        resize_impl( x );
     260    }
     261
     262    void reset( void )
     263    {
     264        m_first_call = true;
     265    }
     266
     267    template< class DerivIn >
     268    void initialize( const DerivIn &deriv )
     269    {
     270        boost::numeric::odeint::copy( deriv , m_dxdt.m_v );
     271        m_first_call = false;
     272    }
     273
     274    template< class System , class StateIn >
     275    void initialize( System system , const StateIn &x , time_type t )
     276    {
     277        typename odeint::unwrap_reference< System >::type &sys = system;
     278        sys( x , m_dxdt.m_v , t );
     279        m_first_call = false;
     280    }
     281
     282    bool is_initialized( void ) const
     283    {
     284        return ! m_first_call;
     285    }
     286
     287
     288
     289private:
     290
     291    template< class System , class StateInOut >
     292    void do_step_v1( System system , StateInOut &x , time_type t , time_type dt )
     293    {
     294        if( m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
     295        {
     296            initialize( system , x , t );
     297        }
     298        this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , m_dxdt.m_v , dt );
     299    }
     300
     301    template< class System , class StateInOut , class Err >
     302    void do_step_v5( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
     303    {
     304        if( m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
     305        {
     306            initialize( system , x , t );
     307        }
     308        this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , m_dxdt.m_v , dt , xerr );
     309    }
     310
     311    template< class StateIn >
     312    bool resize_impl( const StateIn &x )
     313    {
     314        return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
     315    }
     316
     317
     318    stepper_type& stepper( void )
     319    {
     320        return *static_cast< stepper_type* >( this );
     321    }
     322
     323    const stepper_type& stepper( void ) const
     324    {
     325        return *static_cast< const stepper_type* >( this );
     326    }
     327
     328
     329    resizer_type m_resizer;
     330    bool m_first_call;
     331
     332protected:
     333
     334
     335    wrapped_deriv_type m_dxdt;
     336};
     337
     338
     339/******* DOXYGEN *******/
     340
    54341/**
    55342 * \class explicit_error_stepper_fsal_base
     
    135422 * \tparam Resizer The resizer policy class.
    136423 */
    137 template<
    138 class Stepper ,
    139 unsigned short Order ,
    140 unsigned short StepperOrder ,
    141 unsigned short ErrorOrder ,
    142 class State ,
    143 class Value ,
    144 class Deriv ,
    145 class Time ,
    146 class Algebra ,
    147 class Operations ,
    148 class Resizer
    149 >
    150 class explicit_error_stepper_fsal_base : public algebra_stepper_base< Algebra , Operations >
    151 {
    152 public:
    153 
    154     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
    155     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
    156 
    157     typedef State state_type;
    158     typedef Value value_type;
    159     typedef Deriv deriv_type;
    160     typedef Time time_type;
    161     typedef Resizer resizer_type;
    162     typedef Stepper stepper_type;
    163     typedef explicit_error_stepper_fsal_tag stepper_category;
    164 
    165     #ifndef DOXYGEN_SKIP
    166     typedef state_wrapper< state_type > wrapped_state_type;
    167     typedef state_wrapper< deriv_type > wrapped_deriv_type;
    168     typedef explicit_error_stepper_fsal_base< Stepper , Order , StepperOrder , ErrorOrder ,
    169             State , Value , Deriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
    170     #endif
    171 
    172 
    173     typedef unsigned short order_type;
    174     static const order_type order_value = Order;
    175     static const order_type stepper_order_value = StepperOrder;
    176     static const order_type error_order_value = ErrorOrder;
    177 
    178     /**
     424
     425
     426
     427    /**
     428     * \fn explicit_error_stepper_fsal_base::explicit_error_stepper_fsal_base( const algebra_type &algebra )
    179429     * \brief Constructs a explicit_stepper_fsal_base class. This constructor can be used as a default
    180430     * constructor if the algebra has a default constructor.
    181431     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    182432     */
    183     explicit_error_stepper_fsal_base( const algebra_type &algebra = algebra_type() )
    184     : algebra_stepper_base_type( algebra ) , m_first_call( true )
    185     { }
    186 
    187     /**
     433
     434
     435    /**
     436     * \fn explicit_error_stepper_fsal_base::order( void ) const
    188437     * \return Returns the order of the stepper if it used without error estimation.
    189438     */
    190     order_type order( void ) const
    191     {
    192         return order_value;
    193     }
    194 
    195     /**
     439
     440    /**
     441     * \fn explicit_error_stepper_fsal_base::stepper_order( void ) const
    196442     * \return Returns the order of a step if the stepper is used without error estimation.
    197443     */
    198     order_type stepper_order( void ) const
    199     {
    200         return stepper_order_value;
    201     }
    202 
    203     /**
     444
     445
     446    /**
     447     * \fn explicit_error_stepper_fsal_base::error_order( void ) const
    204448     * \return Returns the order of an error step if the stepper is used without error estimation.
    205449     */
    206     order_type error_order( void ) const
    207     {
    208         return error_order_value;
    209     }
    210 
    211 
    212     /*
    213      * version 1 : do_step( sys , x , t , dt )
    214      *
    215      * the two overloads are needed in order to solve the forwarding problem
    216      */
    217     /**
     450
     451    /**
     452     * \fn explicit_error_stepper_fsal_base::do_step( System system , StateInOut &x , time_type t , time_type dt )
    218453     * \brief This method performs one step. It transforms the result in-place.
    219454     *
     
    226461     * \param dt The step size.
    227462     */
    228     template< class System , class StateInOut >
    229     void do_step( System system , StateInOut &x , time_type t , time_type dt )
    230     {
    231         do_step_v1( system , x , t , dt );
    232     }
    233 
    234     /**
    235      * \brief This method performs one step with the stepper passed by Stepper.
    236      * It transforms the result in-place. This method is needed in order to solve the forwarding problem.
    237      * The difference to the other version is that it can be used like
    238      * `stepper.do_step( sys , make_range( iter1 , iter2 ) , t , dt )`
    239      *
    240      * \note This method uses the internal state of the stepper.
    241      *
    242      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
    243      *               Simple System concept.
    244      * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
    245      * \param t The value of the time, at which the step should be performed.
    246      * \param dt The step size.
    247      */
    248     template< class System , class StateInOut >
    249     void do_step( System system , const StateInOut &x , time_type t , time_type dt )
    250     {
    251         do_step_v1( system , x , t , dt );
    252     }
    253 
    254 
    255     /*
    256      * version 2 : do_step( sys , x , dxdt , t , dt )
    257      *
    258      * this version does not solve the forwarding problem, boost.range can not be used
    259      *
    260      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    261      */
    262     /**
     463
     464
     465    /**
     466     * \fn explicit_error_stepper_fsal_base::do_step( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt )
    263467     * \brief The method performs one step with the stepper passed by Stepper. Additionally to the other methods
    264468     * the derivative of x is also passed to this method. Therefore, dxdt must be evaluated initially:
     
    288492     * \param dt The step size.
    289493     */
    290     template< class System , class StateInOut , class DerivInOut >
    291     typename boost::disable_if< boost::is_same< StateInOut , time_type > , void >::type
    292     do_step( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt )
    293     {
    294         m_first_call = true;
    295         this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt );
    296     }
    297 
    298 
    299     /*
    300      * version 3 : do_step( sys , in , t , out , dt )
    301      *
    302      * this version does not solve the forwarding problem, boost.range can not be used
    303      *
    304      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    305      */
    306     /**
     494
     495    /**
     496     * \fn explicit_error_stepper_fsal_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    307497     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    308498     * This method is disabled if StateIn and Time are the same type. In this case the method can not be distinguished from
     
    320510     * \param dt The step size.
    321511     */
    322     template< class System , class StateIn , class StateOut >
    323     typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
    324     do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    325     {
    326         if( m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
    327         {
    328             initialize( system , in , t );
    329         }
    330         this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , m_dxdt.m_v , dt );
    331     }
    332 
    333 
    334     /*
    335      * version 4 : do_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
    336      *
    337      * this version does not solve the forwarding problem, boost.range can not be used
    338      */
    339     /**
     512
     513    /**
     514     * \fn explicit_error_stepper_fsal_base::do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t , StateOut &out , DerivOut &dxdt_out , time_type dt )
    340515     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    341516     * Furthermore, the derivative of x at t is passed to the stepper and updated by the stepper to its new value at
     
    355530     * \param dt The step size.
    356531     */
    357     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
    358     void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
    359             StateOut &out , DerivOut &dxdt_out , time_type dt )
    360     {
    361         m_first_call = true;
    362         this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
    363     }
    364 
    365 
    366 
    367 
    368 
    369     /*
    370      * version 5 : do_step( sys , x , t , dt , xerr )
    371      *
    372      * the two overloads are needed in order to solve the forwarding problem
    373      */
    374     /**
     532
     533    /**
     534     * \fn explicit_error_stepper_fsal_base::do_step( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
    375535     * \brief The method performs one step with the stepper passed by Stepper and estimates the error. The state of the ODE
    376536     * is updated in-place.
     
    386546     * \param xerr The estimation of the error is stored in xerr.
    387547     */
    388     template< class System , class StateInOut , class Err >
    389     void do_step( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
    390     {
    391         do_step_v5( system , x , t , dt , xerr );
    392     }
    393 
    394     /**
    395      * \brief The method performs one step with the stepper passed by Stepper and estimates the error. The state of the ODE
    396      * is updated in-place. This method is needed in order to solve the forwarding problem.
    397      * The difference to the other version is that it can be used like
    398      * `stepper.do_step( sys , make_range( iter1 , iter2 ) , t , dt )`
    399      *
    400      * \note This method uses the internal state of the stepper.
    401      *
    402      * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
    403      *               Simple System concept.
    404      * \param x The state of the ODE which should be solved. x is updated by this method.
    405      * \param t The value of the time, at which the step should be performed.
    406      * \param dt The step size.
    407      * \param xerr The estimation of the error is stored in xerr.
    408      */
    409     template< class System , class StateInOut , class Err >
    410     void do_step( System system , const StateInOut &x , time_type t , time_type dt , Err &xerr )
    411     {
    412         do_step_v5( system , x , t , dt , xerr );
    413     }
    414 
    415 
    416     /*
    417      * version 6 : do_step( sys , x , dxdt , t , dt , xerr )
    418      *
    419      * this version does not solve the forwarding problem, boost.range can not be used
    420      *
    421      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    422          *
    423      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    424      */
    425         /**
     548
     549    /**
     550     * \fn explicit_error_stepper_fsal_base::do_step( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt , Err &xerr )
    426551     * \brief The method performs one step with the stepper passed by Stepper. Additionally to the other method
    427552     * the derivative of x is also passed to this method and updated by this method.
     
    445570     * \param xerr The error estimate is stored in xerr.
    446571     */
    447     template< class System , class StateInOut , class DerivInOut , class Err >
    448     typename boost::disable_if< boost::is_same< StateInOut , time_type > , void >::type
    449     do_step( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt , Err &xerr )
    450     {
    451         m_first_call = true;
    452         this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt , xerr );
    453     }
    454 
    455 
    456 
    457 
    458     /*
    459      * version 7 : do_step( sys , in , t , out , dt , xerr )
    460      *
    461      * this version does not solve the forwarding problem, boost.range can not be used
    462      */
    463     /**
     572
     573
     574    /**
     575     * \fn explicit_error_stepper_fsal_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , Err &xerr )
    464576     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    465577     * Furthermore, the error is estimated.
     
    477589     * \param xerr The error estimate.
    478590     */
    479     template< class System , class StateIn , class StateOut , class Err >
    480     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , Err &xerr )
    481     {
    482         if( m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
    483         {
    484             initialize( system , in , t );
    485         }
    486         this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , m_dxdt.m_v , dt , xerr );
    487     }
    488 
    489 
    490     /*
    491      * version 8 : do_step( sys , in , dxdt_in , t , out , dxdt_out , dt , xerr )
    492      *
    493      * this version does not solve the forwarding problem, boost.range can not be used
    494      */
    495     /**
     591
     592    /**
     593     * \fn explicit_error_stepper_fsal_base::do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t , StateOut &out , DerivOut &dxdt_out , time_type dt , Err &xerr )
    496594     * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    497595     * Furthermore, the derivative of x at t is passed to the stepper and the error is estimated.
     
    511609     * \param xerr The error estimate.
    512610     */
    513     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut , class Err >
    514     void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
    515             StateOut &out , DerivOut &dxdt_out , time_type dt , Err &xerr )
    516     {
    517         m_first_call = true;
    518         this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt , xerr );
    519     }
    520 
    521     /**
     611
     612    /**
     613     * \fn explicit_error_stepper_fsal_base::adjust_size( const StateIn &x )
    522614     * \brief Adjust the size of all temporaries in the stepper manually.
    523615     * \param x A state from which the size of the temporaries to be resized is deduced.
    524616     */
    525     template< class StateIn >
    526     void adjust_size( const StateIn &x )
    527     {
    528         resize_impl( x );
    529     }
    530 
    531 
    532     /**
     617
     618    /**
     619     * \fn explicit_error_stepper_fsal_base::reset( void )
    533620     * \brief Resets the internal state of this stepper. After calling this method it is safe to use all
    534      *        `do_step` method without explicitly initializing the stepper.
    535      */
    536     void reset( void )
    537     {
    538         m_first_call = true;
    539     }
    540 
    541     /**
     621     * `do_step` method without explicitly initializing the stepper.
     622     */
     623
     624    /**
     625     * \fn explicit_error_stepper_fsal_base::initialize( const DerivIn &deriv )
    542626     * \brief Initializes the interal state of the stepper.
    543627     * \param deriv The derivative of x. The next call of `do_step` expects that the derivative of `x` passed to `do_step`
    544628     *              has the value of `deriv`.
    545629     */
    546     template< class DerivIn >
    547     void initialize( const DerivIn &deriv )
    548     {
    549         boost::numeric::odeint::copy( deriv , m_dxdt.m_v );
    550         m_first_call = false;
    551     }
    552 
    553     /**
     630
     631    /**
     632     * \fn explicit_error_stepper_fsal_base::initialize( System system , const StateIn &x , time_type t )
    554633     * \brief Initializes the interal state of the stepper.
    555634     *
     
    565644     * \param t The current time of the ODE.
    566645     */
    567     template< class System , class StateIn >
    568     void initialize( System system , const StateIn &x , time_type t )
    569     {
    570         typename odeint::unwrap_reference< System >::type &sys = system;
    571         sys( x , m_dxdt.m_v , t );
    572         m_first_call = false;
    573     }
    574 
    575     /**
     646
     647    /**
     648     * \fn explicit_error_stepper_fsal_base::is_initialized( void ) const
    576649     * \brief Returns if the stepper is already initialized. If the stepper is not initialized, the first
    577650     * call of `do_step` will initialize the state of the stepper. If the stepper is already initialized
    578651     * the system function can not be safely exchanged between consecutive `do_step` calls.
    579652     */
    580     bool is_initialized( void ) const
    581     {
    582         return ! m_first_call;
    583     }
    584 
    585 
    586 
    587 private:
    588 
    589     template< class System , class StateInOut >
    590     void do_step_v1( System system , StateInOut &x , time_type t , time_type dt )
    591     {
    592         if( m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
    593         {
    594             initialize( system , x , t );
    595         }
    596         this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , m_dxdt.m_v , dt );
    597     }
    598 
    599     template< class System , class StateInOut , class Err >
    600     void do_step_v5( System system , StateInOut &x , time_type t , time_type dt , Err &xerr )
    601     {
    602         if( m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ) || m_first_call )
    603         {
    604             initialize( system , x , t );
    605         }
    606         this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , m_dxdt.m_v , dt , xerr );
    607     }
    608 
    609     template< class StateIn >
    610     bool resize_impl( const StateIn &x )
    611     {
    612         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
    613     }
    614 
    615 
    616     stepper_type& stepper( void )
    617     {
    618         return *static_cast< stepper_type* >( this );
    619     }
    620 
    621     const stepper_type& stepper( void ) const
    622     {
    623         return *static_cast< const stepper_type* >( this );
    624     }
    625 
    626 
    627     resizer_type m_resizer;
    628     bool m_first_call;
    629 
    630 protected:
    631 
    632 
    633     wrapped_deriv_type m_dxdt;
    634 };
    635 
    636653
    637654} // odeint
  • trunk/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp

    r81491 r81754  
    4848    * do_step( sys , in , dxdt_in , t , out , dt )
    4949 */
     50
     51template<
     52class Stepper ,
     53unsigned short Order ,
     54class State ,
     55class Value ,
     56class Deriv ,
     57class Time ,
     58class Algebra ,
     59class Operations ,
     60class Resizer
     61>
     62class explicit_stepper_base : public algebra_stepper_base< Algebra , Operations >
     63{
     64public:
     65
     66    #ifndef DOXYGEN_SKIP
     67    typedef explicit_stepper_base< Stepper , Order , State , Value , Deriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
     68    #endif // DOXYGEN_SKIP
     69
     70
     71    typedef State state_type;
     72    typedef Value value_type;
     73    typedef Deriv deriv_type;
     74    typedef Time time_type;
     75    typedef Resizer resizer_type;
     76    typedef Stepper stepper_type;
     77    typedef stepper_tag stepper_category;
     78    typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
     79    typedef typename algebra_stepper_base_type::algebra_type algebra_type;
     80    typedef typename algebra_stepper_base_type::operations_type operations_type;
     81    typedef unsigned short order_type;
     82
     83    #ifndef DOXYGEN_SKIP
     84    typedef state_wrapper< state_type > wrapped_state_type;
     85    typedef state_wrapper< deriv_type > wrapped_deriv_type;
     86    #endif // DOXYGEN_SKIP
     87
     88
     89    static const order_type order_value = Order;
     90
     91
     92    explicit_stepper_base( const algebra_type &algebra = algebra_type() )
     93    : algebra_stepper_base_type( algebra )
     94    { }
     95
     96    /**
     97     * \return Returns the order of the stepper.
     98     */
     99    order_type order( void ) const
     100    {
     101        return order_value;
     102    }
     103
     104
     105    /*
     106     * Version 1 : do_step( sys , x , t , dt )
     107     *
     108     * the two overloads are needed in order to solve the forwarding problem
     109     */
     110    template< class System , class StateInOut >
     111    void do_step( System system , StateInOut &x , time_type t , time_type dt )
     112    {
     113        do_step_v1( system , x , t , dt );
     114    }
     115
     116    /**
     117     * \brief Second version to solve the forwarding problem, can be called with Boost.Range as StateInOut.
     118     */
     119    template< class System , class StateInOut >
     120    void do_step( System system , const StateInOut &x , time_type t , time_type dt )
     121    {
     122        do_step_v1( system , x , t , dt );
     123    }
     124
     125    /*
     126     * Version 2 : do_step( sys , x , dxdt , t , dt )
     127     *
     128      * this version does not solve the forwarding problem, boost.range can not be used
     129     *
     130     * the disable is needed to avoid ambiguous overloads if state_type = time_type
     131     */
     132    template< class System , class StateInOut , class DerivIn >
     133    typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
     134    do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt )
     135    {
     136        this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
     137    }
     138
     139
     140    /*
     141     * Version 3 : do_step( sys , in , t , out , dt )
     142     *
     143     * this version does not solve the forwarding problem, boost.range can not be used
     144     */
     145    template< class System , class StateIn , class StateOut >
     146    void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     147    {
     148        typename odeint::unwrap_reference< System >::type &sys = system;
     149        m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
     150        sys( in , m_dxdt.m_v ,t );
     151        this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , dt );
     152    }
     153
     154
     155    /*
     156     * Version 4 : do_step( sys , in , dxdt , t , out , dt )
     157     *
     158     * this version does not solve the forwarding problem, boost.range can not be used
     159     */
     160    template< class System , class StateIn , class DerivIn , class StateOut >
     161    void do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
     162    {
     163        this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
     164    }
     165
     166    template< class StateIn >
     167    void adjust_size( const StateIn &x )
     168    {
     169        resize_impl( x );
     170    }
     171
     172private:
     173
     174    stepper_type& stepper( void )
     175    {
     176        return *static_cast< stepper_type* >( this );
     177    }
     178
     179    const stepper_type& stepper( void ) const
     180    {
     181        return *static_cast< const stepper_type* >( this );
     182    }
     183
     184
     185    template< class StateIn >
     186    bool resize_impl( const StateIn &x )
     187    {
     188        return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
     189    }
     190
     191
     192    template< class System , class StateInOut >
     193    void do_step_v1( System system , StateInOut &x , time_type t , time_type dt )
     194    {
     195        typename odeint::unwrap_reference< System >::type &sys = system;
     196        m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
     197        sys( x , m_dxdt.m_v ,t );
     198        this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , dt );
     199    }
     200
     201
     202    resizer_type m_resizer;
     203
     204protected:
     205
     206    wrapped_deriv_type m_dxdt;
     207};
     208
     209
     210/******* DOXYGEN *********/
    50211
    51212/**
     
    125286 * \tparam Resizer The resizer policy class.
    126287 */
    127 template<
    128 class Stepper ,
    129 unsigned short Order ,
    130 class State ,
    131 class Value ,
    132 class Deriv ,
    133 class Time ,
    134 class Algebra ,
    135 class Operations ,
    136 class Resizer
    137 >
    138 class explicit_stepper_base : public algebra_stepper_base< Algebra , Operations >
    139 {
    140 public:
    141 
    142     #ifndef DOXYGEN_SKIP
    143     typedef explicit_stepper_base< Stepper , Order , State , Value , Deriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
    144     #endif // DOXYGEN_SKIP
    145 
    146 
    147     typedef State state_type;
    148     typedef Value value_type;
    149     typedef Deriv deriv_type;
    150     typedef Time time_type;
    151     typedef Resizer resizer_type;
    152     typedef Stepper stepper_type;
    153     typedef stepper_tag stepper_category;
    154     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
    155     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
    156     typedef typename algebra_stepper_base_type::operations_type operations_type;
    157     typedef unsigned short order_type;
    158 
    159     #ifndef DOXYGEN_SKIP
    160     typedef state_wrapper< state_type > wrapped_state_type;
    161     typedef state_wrapper< deriv_type > wrapped_deriv_type;
    162     #endif // DOXYGEN_SKIP
    163 
    164 
    165     static const order_type order_value = Order;
    166 
    167 
    168     /**
     288
     289
     290    /**
     291     * \fn explicit_stepper_base::explicit_stepper_base( const algebra_type &algebra )
    169292     * \brief Constructs a explicit_stepper_base class. This constructor can be used as a default
    170293     * constructor if the algebra has a default constructor.
    171294     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    172295     */
    173     explicit_stepper_base( const algebra_type &algebra = algebra_type() )
    174     : algebra_stepper_base_type( algebra )
    175     { }
    176 
    177     /**
     296
     297    /**
     298     * \fn explicit_stepper_base::order_type order( void ) const
    178299     * \return Returns the order of the stepper.
    179300     */
    180     order_type order( void ) const
    181     {
    182         return order_value;
    183     }
    184 
    185 
    186     /*
    187      * Version 1 : do_step( sys , x , t , dt )
    188      *
    189      * the two overloads are needed in order to solve the forwarding problem
    190      */
    191 
    192     /**
     301
     302    /**
     303     * \fn explicit_stepper_base::do_step( System system , StateInOut &x , time_type t , time_type dt )
    193304     * \brief This method performs one step. It transforms the result in-place.
    194305     *
     
    199310     * \param dt The step size.
    200311     */
    201     template< class System , class StateInOut >
    202     void do_step( System system , StateInOut &x , time_type t , time_type dt )
    203     {
    204         do_step_v1( system , x , t , dt );
    205     }
    206 
    207     /**
    208      * \brief This method performs one step with the stepper passed by Stepper.
    209      * It transforms the result in-place. This method is needed in order to solve the forwarding problem.
    210      * The difference to the other version is that it can be used like
    211      * `stepper.do_step( sys , make_range( iter1 , iter2 ) , t , dt )`
    212      *
    213      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
    214      *               Simple System concept.
    215      * \param x The state of the ODE which should be solved. After calling do_step the result is updated in x.
    216      * \param t The value of the time, at which the step should be performed.
    217      * \param dt The step size.
    218      */
    219     template< class System , class StateInOut >
    220     void do_step( System system , const StateInOut &x , time_type t , time_type dt )
    221     {
    222         do_step_v1( system , x , t , dt );
    223     }
    224 
    225 
    226     /*
    227      * Version 2 : do_step( sys , x , dxdt , t , dt )
    228      *
    229       * this version does not solve the forwarding problem, boost.range can not be used
    230      *
    231      * the disable is needed to avoid ambiguous overloads if state_type = time_type
    232      */
    233     /**
    234      * \brief The method performs one step with the stepper passed by Stepper. Additionally to the other method
    235      * the derivative of x is also passed to this method. It is equivalent to
     312
     313
     314    /**
     315     * \fn explicit_stepper_base::do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt )
     316
     317     * \brief The method performs one step. Additionally to the other method
     318     * the derivative of x is also passed to this method. It is supposed to be used in the following way:
    236319     *
    237320     * \code
     
    252335     * \param dt The step size.
    253336     */
    254     template< class System , class StateInOut , class DerivIn >
    255     typename boost::disable_if< boost::is_same< DerivIn , time_type > , void >::type
    256     do_step( System system , StateInOut &x , const DerivIn &dxdt , time_type t , time_type dt )
    257     {
    258         this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
    259     }
    260 
    261 
    262     /*
    263      * Version 3 : do_step( sys , in , t , out , dt )
    264      *
    265      * this version does not solve the forwarding problem, boost.range can not be used
    266      */
    267     /**
    268      * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
     337
     338    /**
     339     * \fn void explicit_stepper_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     340     * \brief The method performs one step. The state of the ODE is updated out-of-place.
    269341     * \note This method does not solve the forwarding problem.
    270342     *
     
    276348     * \param dt The step size.
    277349     */
    278     template< class System , class StateIn , class StateOut >
    279     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    280     {
    281         typename odeint::unwrap_reference< System >::type &sys = system;
    282         m_resizer.adjust_size( in , detail::bind( &internal_stepper_base_type::template resize_impl<StateIn> , detail::ref( *this ) , detail::_1 ) );
    283         sys( in , m_dxdt.m_v ,t );
    284         this->stepper().do_step_impl( system , in , m_dxdt.m_v , t , out , dt );
    285     }
    286 
    287 
    288     /*
    289      * Version 4 : do_step( sys , in , dxdt , t , out , dt )
    290      *
    291      * this version does not solve the forwarding problem, boost.range can not be used
    292      */
    293     /**
    294      * \brief The method performs one step with the stepper passed by Stepper. The state of the ODE is updated out-of-place.
    295      * Furthermore, the derivative of x at t is passed to the stepper. It is equivalent to:
     350
     351    /**
     352     * \fn void explicit_stepper_base::do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
     353     * \brief The method performs one step. The state of the ODE is updated out-of-place.
     354     * Furthermore, the derivative of x at t is passed to the stepper.
     355     * It is supposed to be used in the following way:
    296356     *
    297357     * \code
     
    310370     * \param dt The step size.
    311371     */
    312     template< class System , class StateIn , class DerivIn , class StateOut >
    313     void do_step( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    314     {
    315         this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
    316     }
    317 
    318 
    319     /**
     372
     373    /**
     374     * \fn void explicit_stepper_base::adjust_size( const StateIn &x )
    320375     * \brief Adjust the size of all temporaries in the stepper manually.
    321376     * \param x A state from which the size of the temporaries to be resized is deduced.
    322377     */
    323     template< class StateIn >
    324     void adjust_size( const StateIn &x )
    325     {
    326         resize_impl( x );
    327     }
    328 
    329 private:
    330 
    331     stepper_type& stepper( void )
    332     {
    333         return *static_cast< stepper_type* >( this );
    334     }
    335 
    336     const stepper_type& stepper( void ) const
    337     {
    338         return *static_cast< const stepper_type* >( this );
    339     }
    340 
    341 
    342     template< class StateIn >
    343     bool resize_impl( const StateIn &x )
    344     {
    345         return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
    346     }
    347 
    348 
    349     template< class System , class StateInOut >
    350     void do_step_v1( System system , StateInOut &x , time_type t , time_type dt )
    351     {
    352         typename odeint::unwrap_reference< System >::type &sys = system;
    353         m_resizer.adjust_size( x , detail::bind( &internal_stepper_base_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
    354         sys( x , m_dxdt.m_v ,t );
    355         this->stepper().do_step_impl( system , x , m_dxdt.m_v , t , x , dt );
    356     }
    357 
    358 
    359     resizer_type m_resizer;
    360 
    361 protected:
    362 
    363     wrapped_deriv_type m_dxdt;
    364 };
    365 
    366378
    367379} // odeint
  • trunk/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp

    r81491 r81754  
    4141namespace odeint {
    4242
     43
     44template<
     45size_t NumOfStages ,
     46unsigned short Order ,
     47class Coor ,
     48class Momentum ,
     49class Value ,
     50class CoorDeriv ,
     51class MomentumDeriv ,
     52class Time ,
     53class Algebra ,
     54class Operations ,
     55class Resizer
     56>
     57class symplectic_nystroem_stepper_base : public algebra_stepper_base< Algebra , Operations >
     58{
     59
     60public:
     61
     62    typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
     63    typedef typename algebra_stepper_base_type::algebra_type algebra_type;
     64    typedef typename algebra_stepper_base_type::operations_type operations_type;
     65
     66    const static size_t num_of_stages = NumOfStages;
     67    typedef Coor coor_type;
     68    typedef Momentum momentum_type;
     69    typedef std::pair< coor_type , momentum_type > state_type;
     70    typedef CoorDeriv coor_deriv_type;
     71    typedef state_wrapper< coor_deriv_type> wrapped_coor_deriv_type;
     72    typedef MomentumDeriv momentum_deriv_type;
     73    typedef state_wrapper< momentum_deriv_type > wrapped_momentum_deriv_type;
     74    typedef std::pair< coor_deriv_type , momentum_deriv_type > deriv_type;
     75    typedef Value value_type;
     76    typedef Time time_type;
     77    typedef Resizer resizer_type;
     78    typedef stepper_tag stepper_category;
     79   
     80    #ifndef DOXYGEN_SKIP
     81    typedef symplectic_nystroem_stepper_base< NumOfStages , Order , Coor , Momentum , Value ,
     82            CoorDeriv , MomentumDeriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
     83    #endif
     84    typedef unsigned short order_type;
     85
     86    static const order_type order_value = Order;
     87
     88    typedef boost::array< value_type , num_of_stages > coef_type;
     89
     90    symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra = algebra_type() )
     91        : algebra_stepper_base_type( algebra ) , m_coef_a( coef_a ) , m_coef_b( coef_b ) ,
     92          m_dqdt_resizer() , m_dpdt_resizer() , m_dqdt() , m_dpdt()
     93    { }
     94
     95
     96    order_type order( void ) const
     97    {
     98        return order_value;
     99    }
     100
     101    /*
     102     * Version 1 : do_step( system , x , t , dt )
     103     *
     104     * This version does not solve the forwarding problem, boost.range can not be used.
     105     */
     106    template< class System , class StateInOut >
     107    void do_step( System system , const StateInOut &state , time_type t , time_type dt )
     108    {
     109        typedef typename odeint::unwrap_reference< System >::type system_type;
     110        do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
     111    }
     112
     113    /**
     114     * \brief Same function as above. It differs only in a different const specifier in order
     115     * to solve the forwarding problem, can be used with Boost.Range.
     116     */
     117    template< class System , class StateInOut >
     118    void do_step( System system , StateInOut &state , time_type t , time_type dt )
     119    {
     120        typedef typename odeint::unwrap_reference< System >::type system_type;
     121        do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
     122    }
     123
     124
     125
     126
     127    /*
     128     * Version 2 : do_step( system , q , p , t , dt );
     129     *
     130     * For Convenience
     131     *
     132     * The two overloads are needed in order to solve the forwarding problem.
     133     */
     134    template< class System , class CoorInOut , class MomentumInOut >
     135    void do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
     136    {
     137        do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
     138    }
     139
     140    /**
     141     * \brief Same function as do_step( system , q , p , t , dt ). It differs only in a different const specifier in order
     142     * to solve the forwarding problem, can be called with Boost.Range.
     143     */
     144    template< class System , class CoorInOut , class MomentumInOut >
     145    void do_step( System system , const CoorInOut &q , const MomentumInOut &p , time_type t , time_type dt )
     146    {
     147        do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
     148    }
     149
     150
     151
     152
     153
     154    /*
     155     * Version 3 : do_step( system , in , t , out , dt )
     156     *
     157     * The forwarding problem is not solved in this version
     158     */
     159    template< class System , class StateIn , class StateOut >
     160    void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
     161    {
     162        typedef typename odeint::unwrap_reference< System >::type system_type;
     163        do_step_impl( system , in , t , out , dt , typename is_pair< system_type >::type() );
     164    }
     165
     166
     167    template< class StateType >
     168    void adjust_size( const StateType &x )
     169    {
     170        resize_dqdt( x );
     171        resize_dpdt( x );
     172    }
     173
     174    /** \brief Returns the coefficients a. */
     175    const coef_type& coef_a( void ) const { return m_coef_a; }
     176
     177    /** \brief Returns the coefficients b. */
     178    const coef_type& coef_b( void ) const { return m_coef_b; }
     179
     180private:
     181
     182    // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
     183    template< class System , class StateIn , class StateOut >
     184    void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ )
     185    {
     186        typedef typename odeint::unwrap_reference< System >::type system_type;
     187        typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
     188        typedef typename odeint::unwrap_reference< typename system_type::second_type >::type momentum_deriv_func_type;
     189        system_type &sys = system;
     190        coor_deriv_func_type &coor_func = sys.first;
     191        momentum_deriv_func_type &momentum_func = sys.second;
     192
     193        typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
     194        typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
     195        typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
     196        const state_in_type &state_in = in;
     197        const coor_in_type &coor_in = state_in.first;
     198        const momentum_in_type &momentum_in = state_in.second;
     199
     200        typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
     201        typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
     202        typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
     203        state_out_type &state_out = out;
     204        coor_out_type &coor_out = state_out.first;
     205        momentum_out_type &momentum_out = state_out.second;
     206
     207        m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
     208        m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
     209
     210        // ToDo: check sizes?
     211
     212        for( size_t l=0 ; l<num_of_stages ; ++l )
     213        {
     214            if( l == 0 )
     215            {
     216                coor_func( momentum_in , m_dqdt.m_v );
     217                this->m_algebra.for_each3( coor_out , coor_in , m_dqdt.m_v ,
     218                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
     219                momentum_func( coor_out , m_dpdt.m_v );
     220                this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
     221                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
     222            }
     223            else
     224            {
     225                coor_func( momentum_out , m_dqdt.m_v );
     226                this->m_algebra.for_each3( coor_out , coor_out , m_dqdt.m_v ,
     227                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
     228                momentum_func( coor_out , m_dpdt.m_v );
     229                this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
     230                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
     231            }
     232        }
     233    }
     234
     235
     236    // stepper for systems with only function dp /dt = -f(q), dq/dt = p, time not required but still expected for compatibility reasons
     237    template< class System , class StateIn , class StateOut >
     238    void do_step_impl( System system , const StateIn &in , time_type  /* t */ , StateOut &out , time_type dt , boost::mpl::false_ )
     239    {
     240        typedef typename odeint::unwrap_reference< System >::type momentum_deriv_func_type;
     241        momentum_deriv_func_type &momentum_func = system;
     242
     243        typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
     244        typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
     245        typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
     246        const state_in_type &state_in = in;
     247        const coor_in_type &coor_in = state_in.first;
     248        const momentum_in_type &momentum_in = state_in.second;
     249
     250        typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
     251        typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
     252        typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
     253        state_out_type &state_out = out;
     254        coor_out_type &coor_out = state_out.first;
     255        momentum_out_type &momentum_out = state_out.second;
     256
     257
     258        m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
     259        m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
     260
     261
     262        // ToDo: check sizes?
     263
     264        for( size_t l=0 ; l<num_of_stages ; ++l )
     265        {
     266            if( l == 0 )
     267            {
     268                this->m_algebra.for_each3( coor_out  , coor_in , momentum_in ,
     269                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
     270                momentum_func( coor_out , m_dpdt.m_v );
     271                this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
     272                                           typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
     273            }
     274            else
     275            {
     276                this->m_algebra.for_each3( coor_out , coor_out , momentum_out ,
     277                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
     278                momentum_func( coor_out , m_dpdt.m_v );
     279                this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
     280                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
     281            }
     282        }
     283    }
     284
     285    template< class StateIn >
     286    bool resize_dqdt( const StateIn &x )
     287    {
     288        return adjust_size_by_resizeability( m_dqdt , x , typename is_resizeable<coor_deriv_type>::type() );
     289    }
     290
     291    template< class StateIn >
     292    bool resize_dpdt( const StateIn &x )
     293    {
     294        return adjust_size_by_resizeability( m_dpdt , x , typename is_resizeable<momentum_deriv_type>::type() );
     295    }
     296
     297
     298    const coef_type m_coef_a;
     299    const coef_type m_coef_b;
     300
     301    resizer_type m_dqdt_resizer;
     302    resizer_type m_dpdt_resizer;
     303    wrapped_coor_deriv_type m_dqdt;
     304    wrapped_momentum_deriv_type m_dpdt;
     305
     306};
     307
     308/********* DOXYGEN *********/
    43309
    44310/**
     
    79345 * \tparam Resizer The resizer policy.
    80346 */
    81 template<
    82 size_t NumOfStages ,
    83 unsigned short Order ,
    84 class Coor ,
    85 class Momentum ,
    86 class Value ,
    87 class CoorDeriv ,
    88 class MomentumDeriv ,
    89 class Time ,
    90 class Algebra ,
    91 class Operations ,
    92 class Resizer
    93 >
    94 class symplectic_nystroem_stepper_base : public algebra_stepper_base< Algebra , Operations >
    95 {
    96 
    97 public:
    98 
    99     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
    100     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
    101     typedef typename algebra_stepper_base_type::operations_type operations_type;
    102 
    103     const static size_t num_of_stages = NumOfStages;
    104     typedef Coor coor_type;
    105     typedef Momentum momentum_type;
    106     typedef std::pair< coor_type , momentum_type > state_type;
    107     typedef CoorDeriv coor_deriv_type;
    108     typedef state_wrapper< coor_deriv_type> wrapped_coor_deriv_type;
    109     typedef MomentumDeriv momentum_deriv_type;
    110     typedef state_wrapper< momentum_deriv_type > wrapped_momentum_deriv_type;
    111     typedef std::pair< coor_deriv_type , momentum_deriv_type > deriv_type;
    112     typedef Value value_type;
    113     typedef Time time_type;
    114     typedef Resizer resizer_type;
    115     typedef stepper_tag stepper_category;
    116    
    117     #ifndef DOXYGEN_SKIP
    118     typedef symplectic_nystroem_stepper_base< NumOfStages , Order , Coor , Momentum , Value ,
    119             CoorDeriv , MomentumDeriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
    120     #endif
    121     typedef unsigned short order_type;
    122 
    123     static const order_type order_value = Order;
    124 
    125     typedef boost::array< value_type , num_of_stages > coef_type;
    126 
    127      /**
     347
     348    /**
     349     * \fn symplectic_nystroem_stepper_base::symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra )
    128350     * \brief Constructs a symplectic_nystroem_stepper_base class. The parameters of the specific Nystroem method and the
    129351     * algebra have to be passed.
     
    132354     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    133355     */
    134     symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra = algebra_type() )
    135         : algebra_stepper_base_type( algebra ) , m_coef_a( coef_a ) , m_coef_b( coef_b ) ,
    136           m_dqdt_resizer() , m_dpdt_resizer() , m_dqdt() , m_dpdt()
    137     { }
    138 
    139 
    140     /**
     356
     357    /**
     358     * \fn symplectic_nystroem_stepper_base::order( void ) const
    141359     * \return Returns the order of the stepper.
    142360     */
    143     order_type order( void ) const
    144     {
    145         return order_value;
    146     }
    147 
    148     /*
    149      * Version 1 : do_step( system , x , t , dt )
    150      *
    151      * This version does not solve the forwarding problem, boost.range can not be used.
    152      */
    153     /**
     361
     362    /**
     363     * \fn symplectic_nystroem_stepper_base::do_step( System system , const StateInOut &state , time_type t , time_type dt )
    154364     * \brief This method performs one step. The system can be either a pair of two function object
    155365     * decribing the momentum part and the coordinate part or one function object describing only
     
    168378     * \param dt The time step.
    169379     */
    170     template< class System , class StateInOut >
    171     void do_step( System system , const StateInOut &state , time_type t , time_type dt )
    172     {
    173         typedef typename odeint::unwrap_reference< System >::type system_type;
    174         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
    175     }
    176 
    177     /**
    178      * \brief Same function as do_step( system , x , t , dt ). It differs only in a different const specifier in order
    179      * to solve the forwarding problem.
    180      */
    181     template< class System , class StateInOut >
    182     void do_step( System system , StateInOut &state , time_type t , time_type dt )
    183     {
    184         typedef typename odeint::unwrap_reference< System >::type system_type;
    185         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
    186     }
    187 
    188 
    189 
    190 
    191     /*
    192      * Version 2 : do_step( system , q , p , t , dt );
    193      *
    194      * For Convenience
    195      *
    196      * The two overloads are needed in order to solve the forwarding problem.
    197      */
    198     /**
     380
     381    /**
     382     * \fn symplectic_nystroem_stepper_base::do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
    199383     * \brief This method performs one step. The system can be either a pair of two function object
    200384     * decribing the momentum part and the coordinate part or one function object describing only
     
    215399     * \param dt The time step.
    216400     */
    217     template< class System , class CoorInOut , class MomentumInOut >
    218     void do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
    219     {
    220         do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
    221     }
    222 
    223     /**
    224      * \brief Same function as do_step( system , q , p , t , dt ). It differs only in a different const specifier in order
    225      * to solve the forwarding problem.
    226      */
    227     template< class System , class CoorInOut , class MomentumInOut >
    228     void do_step( System system , const CoorInOut &q , const MomentumInOut &p , time_type t , time_type dt )
    229     {
    230         do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
    231     }
    232 
    233 
    234 
    235 
    236 
    237     /*
    238      * Version 3 : do_step( system , in , t , out , dt )
    239      *
    240      * The forwarding problem is not solved in this version
    241      */
    242     /**
     401
     402    /**
     403     * \fn symplectic_nystroem_stepper_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    243404     * \brief This method performs one step. The system can be either a pair of two function object
    244405     * decribing the momentum part and the coordinate part or one function object describing only
     
    258419     * \param dt The time step.
    259420     */
    260     template< class System , class StateIn , class StateOut >
    261     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
    262     {
    263         typedef typename odeint::unwrap_reference< System >::type system_type;
    264         do_step_impl( system , in , t , out , dt , typename is_pair< system_type >::type() );
    265     }
    266 
    267 
    268     /**
     421
     422    /**
     423     * \fn symplectic_nystroem_stepper_base::adjust_size( const StateType &x )
    269424     * \brief Adjust the size of all temporaries in the stepper manually.
    270425     * \param x A state from which the size of the temporaries to be resized is deduced.
    271426     */
    272     template< class StateType >
    273     void adjust_size( const StateType &x )
    274     {
    275         resize_dqdt( x );
    276         resize_dpdt( x );
    277     }
    278 
    279     /** \brief Returns the coefficients a. */
    280     const coef_type& coef_a( void ) const { return m_coef_a; }
    281 
    282     /** \brief Returns the coefficients b. */
    283     const coef_type& coef_b( void ) const { return m_coef_b; }
    284 
    285 private:
    286 
    287     // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
    288     template< class System , class StateIn , class StateOut >
    289     void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ )
    290     {
    291         typedef typename odeint::unwrap_reference< System >::type system_type;
    292         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
    293         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type momentum_deriv_func_type;
    294         system_type &sys = system;
    295         coor_deriv_func_type &coor_func = sys.first;
    296         momentum_deriv_func_type &momentum_func = sys.second;
    297 
    298         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
    299         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
    300         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
    301         const state_in_type &state_in = in;
    302         const coor_in_type &coor_in = state_in.first;
    303         const momentum_in_type &momentum_in = state_in.second;
    304 
    305         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
    306         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
    307         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
    308         state_out_type &state_out = out;
    309         coor_out_type &coor_out = state_out.first;
    310         momentum_out_type &momentum_out = state_out.second;
    311 
    312         m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
    313         m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
    314 
    315         // ToDo: check sizes?
    316 
    317         for( size_t l=0 ; l<num_of_stages ; ++l )
    318         {
    319             if( l == 0 )
    320             {
    321                 coor_func( momentum_in , m_dqdt.m_v );
    322                 this->m_algebra.for_each3( coor_out , coor_in , m_dqdt.m_v ,
    323                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
    324                 momentum_func( coor_out , m_dpdt.m_v );
    325                 this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
    326                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
    327             }
    328             else
    329             {
    330                 coor_func( momentum_out , m_dqdt.m_v );
    331                 this->m_algebra.for_each3( coor_out , coor_out , m_dqdt.m_v ,
    332                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
    333                 momentum_func( coor_out , m_dpdt.m_v );
    334                 this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
    335                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
    336             }
    337         }
    338     }
    339 
    340 
    341     // stepper for systems with only function dp /dt = -f(q), dq/dt = p, time not required but still expected for compatibility reasons
    342     template< class System , class StateIn , class StateOut >
    343     void do_step_impl( System system , const StateIn &in , time_type  /* t */ , StateOut &out , time_type dt , boost::mpl::false_ )
    344     {
    345         typedef typename odeint::unwrap_reference< System >::type momentum_deriv_func_type;
    346         momentum_deriv_func_type &momentum_func = system;
    347 
    348         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
    349         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
    350         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
    351         const state_in_type &state_in = in;
    352         const coor_in_type &coor_in = state_in.first;
    353         const momentum_in_type &momentum_in = state_in.second;
    354 
    355         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
    356         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
    357         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
    358         state_out_type &state_out = out;
    359         coor_out_type &coor_out = state_out.first;
    360         momentum_out_type &momentum_out = state_out.second;
    361 
    362 
    363         m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
    364         m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
    365 
    366 
    367         // ToDo: check sizes?
    368 
    369         for( size_t l=0 ; l<num_of_stages ; ++l )
    370         {
    371             if( l == 0 )
    372             {
    373                 this->m_algebra.for_each3( coor_out  , coor_in , momentum_in ,
    374                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
    375                 momentum_func( coor_out , m_dpdt.m_v );
    376                 this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
    377                                            typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
    378             }
    379             else
    380             {
    381                 this->m_algebra.for_each3( coor_out , coor_out , momentum_out ,
    382                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
    383                 momentum_func( coor_out , m_dpdt.m_v );
    384                 this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
    385                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
    386             }
    387         }
    388     }
    389 
    390     template< class StateIn >
    391     bool resize_dqdt( const StateIn &x )
    392     {
    393         return adjust_size_by_resizeability( m_dqdt , x , typename is_resizeable<coor_deriv_type>::type() );
    394     }
    395 
    396     template< class StateIn >
    397     bool resize_dpdt( const StateIn &x )
    398     {
    399         return adjust_size_by_resizeability( m_dpdt , x , typename is_resizeable<momentum_deriv_type>::type() );
    400     }
    401 
    402 
    403     const coef_type m_coef_a;
    404     const coef_type m_coef_b;
    405 
    406     resizer_type m_dqdt_resizer;
    407     resizer_type m_dpdt_resizer;
    408     wrapped_coor_deriv_type m_dqdt;
    409     wrapped_momentum_deriv_type m_dpdt;
    410 
    411 };
    412427
    413428} // namespace odeint
  • trunk/boost/numeric/odeint/stepper/bulirsch_stoer.hpp

    r81491 r81754  
    44
    55  [begin_description]
    6   Implementaiton of the Burlish-Stoer method. As described in
     6  Implementation of the Burlish-Stoer method. As described in
    77  Ernst Hairer, Syvert Paul Norsett, Gerhard Wanner
    88  Solving Ordinary Differential Equations I. Nonstiff Problems.
     
    4747namespace numeric {
    4848namespace odeint {
    49 
    50 /** ToDo try_step stepsize changed return values doesn't make too much sense here as we have order control as well */
    5149
    5250template<
     
    7068    typedef Operations operations_type;
    7169    typedef Resizer resizer_type;
     70#ifndef DOXYGEN_SKIP
    7271    typedef state_wrapper< state_type > wrapped_state_type;
    7372    typedef state_wrapper< deriv_type > wrapped_deriv_type;
     
    8483    typedef std::vector< size_t > int_vector;
    8584    typedef std::vector< wrapped_state_type > state_table_type;
    86 
     85#endif //DOXYGEN_SKIP
    8786    const static size_t m_k_max = 8;
    88 
    8987
    9088    bulirsch_stoer(
     
    9391        : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
    9492          m_last_step_rejected( false ) , m_first( true ) ,
    95           /* , m_t_last() ,
    96              m_current_k_opt() ,
    97              m_algebra() ,
    98              m_dxdt_resizer() , m_xnew_resizer() , m_resizer() ,
    99              m_xnew() , m_err() , m_dxdt() ,*/
    10093          m_interval_sequence( m_k_max+1 ) ,
    10194          m_coeff( m_k_max+1 ) ,
     
    10699        BOOST_USING_STD_MIN();
    107100        BOOST_USING_STD_MAX();
    108         //m_dt_last = 1.0E30;
     101        /* initialize sequence of stage numbers and work */
    109102        for( unsigned short i = 0; i < m_k_max+1; i++ )
    110103        {
     
    119112                const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
    120113                m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
    121                 //std::cout << i << "," << k << " " << m_coeff[i][k] << '\t' ;
    122114            }
    123             //std ::cout << std::endl;
     115
    124116            // crude estimate of optimal order
    125             const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , 1.0E-12 ) ) * 0.6 + 0.5 );
    126             m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 1 , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>( m_k_max-1 ) , static_cast<int>( logfact ) ));
    127             //m_current_k_opt = m_k_max - 1;
    128             //std::cout << m_cost[i] << std::endl;
     117
     118            m_current_k_opt = 4;
     119            /* no calculation because log10 might not exist for value_type!
     120            const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
     121            m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
     122            */
    129123        }
    130124
     
    143137    }
    144138
     139    /**
     140     * \brief Second version to solve the forwarding problem, can be used with Boost.Range as StateInOut.
     141     */
    145142    template< class System , class StateInOut >
    146143    controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
     
    172169     */
    173170    template< class System , class StateIn , class StateOut >
    174     controlled_step_result try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
     171    typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
     172    try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
    175173    {
    176174        typename odeint::unwrap_reference< System >::type &sys = system;
     
    181179
    182180
     181    /*
     182     * Full version : try_step( sys , in , dxdt_in , t , out , dt )
     183     *
     184     * contains the actual implementation
     185     */
    183186    template< class System , class StateIn , class DerivIn , class StateOut >
    184187    controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
     
    205208        inv_time_vector work( m_k_max+1 );
    206209
    207         //std::cout << "t=" << t <<", dt=" << dt << "(" << m_dt_last << ")" << ", k_opt=" << m_current_k_opt << std::endl;
    208 
    209210        time_type new_h = dt;
    210211
     212        /* m_current_k_opt is the estimated current optimal stage number */
    211213        for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
    212214        {
    213             //std::cout << "  k=" << k; //<<": " << ", first: " << m_first << std::endl;
     215            /* the stage counts are stored in m_interval_sequence */
    214216            m_midpoint.set_steps( m_interval_sequence[k] );
    215217            if( k == 0 )
    216218            {
    217219                m_midpoint.do_step( sys , in , dxdt , t , out , dt );
     220                /* the first step, nothing more to do */
    218221            }
    219222            else
     
    227230                h_opt[k] = calc_h_opt( dt , error , k );
    228231                work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
    229                 //std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl;
    230                 //std::cout << '\t' << "error: " << error << std::endl;
    231232
    232233                if( (k == m_current_k_opt-1) || m_first )
     
    271272                            new_h = h_opt[k];
    272273                            new_h *= m_cost[m_current_k_opt]/m_cost[k];
    273                             //std::cout << new_h << std::endl;
    274274                        } else
    275275                            new_h = h_opt[m_current_k_opt];
     
    285285                if( k == m_current_k_opt+1 )
    286286                { // convergence at k_opt+1 ?
    287                     //std::cout << "convergence at k_opt+1 ?" << std::endl;
    288287                    if( error < 1.0 )
    289288                    {   //convergence
     
    307306        {
    308307            t += dt;
    309         }// else
    310         //   std::cout << "REJECT!" << std::endl;
     308        }
    311309
    312310        if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) )
     
    325323    }
    326324
     325    /** \brief Resets the internal state of the stepper */
    327326    void reset()
    328327    {
    329         //std::cout << "reset" << std::endl;
    330328        m_first = true;
    331329        m_last_step_rejected = false;
     
    334332
    335333    /* Resizer methods */
    336 
    337334
    338335    template< class StateIn >
     
    383380    template< class StateInOut >
    384381    void extrapolate( size_t k , state_table_type &table , const value_matrix &coeff , StateInOut &xest )
    385     //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
    386     {
    387         //std::cout << "extrapolate k=" << k << ":" << std::endl;
     382    /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
     383       uses the obtained intermediate results to extrapolate to dt->0
     384    */
     385    {
    388386        static const value_type val1 = static_cast< value_type >( 1.0 );
    389387        for( int j=k-1 ; j>0 ; --j )
    390388        {
    391             //std::cout << '\t' << m_coeff[k][j];
    392389            m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
    393390                                 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) );
    394391        }
    395         //std::cout << std::endl << m_coeff[k][0] << std::endl;
    396392        m_algebra.for_each3( xest , table[0].m_v , xest ,
    397393                             typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) );
     
    399395
    400396    time_type calc_h_opt( time_type h , value_type error , size_t k ) const
     397    /* calculates the optimal step size for a given error and stage number */
    401398    {
    402399        BOOST_USING_STD_MIN();
    403400        BOOST_USING_STD_MAX();
    404401        using std::pow;
    405         value_type expo=1.0/(2*k+1);
     402        value_type expo( 1.0/(2*k+1) );
    406403        value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
    407404        value_type fac;
     
    413410            fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( facmin/STEPFAC4 , min BOOST_PREVENT_MACRO_SUBSTITUTION( 1.0/facmin , fac ) );
    414411        }
    415         //return std::abs(h*fac);
    416412        return h*fac;
    417413    }
    418414
    419415    controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt )
    420     {
    421         //std::cout << "finding k_opt..." << std::endl;
     416    /* calculates the optimal stage number */
     417    {
    422418        if( k == 1 )
    423419        {
    424420            m_current_k_opt = 2;
    425             //dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ;
    426421            return success;
    427422        }
     
    501496
    502497
     498/******** DOXYGEN ********/
     499/**
     500 * \class bulirsch_stoer
     501 * \brief The Bulirsch-Stoer algorithm.
     502 *
     503 * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
     504 * and order of the method. The algorithm uses the modified midpoint and
     505 * a polynomial extrapolation compute the solution.
     506 *
     507 * \tparam State The state type.
     508 * \tparam Value The value type.
     509 * \tparam Deriv The type representing the time derivative of the state.
     510 * \tparam Time The time representing the independent variable - the time.
     511 * \tparam Algebra The algebra type.
     512 * \tparam Operations The operations type.
     513 * \tparam Resizer The resizer policy type.
     514 */
     515
     516    /**
     517     * \fn bulirsch_stoer::bulirsch_stoer( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt )
     518     * \brief Constructs the bulirsch_stoer class, including initialization of
     519     * the error bounds.
     520     *
     521     * \param eps_abs Absolute tolerance level.
     522     * \param eps_rel Relative tolerance level.
     523     * \param factor_x Factor for the weight of the state.
     524     * \param factor_dxdt Factor for the weight of the derivative.
     525     */
     526
     527    /**
     528     * \fn bulirsch_stoer::try_step( System system , StateInOut &x , time_type &t , time_type &dt )
     529     * \brief Tries to perform one step.
     530     *
     531     * This method tries to do one step with step size dt. If the error estimate
     532     * is to large, the step is rejected and the method returns fail and the
     533     * step size dt is reduced. If the error estimate is acceptably small, the
     534     * step is performed, success is returned and dt might be increased to make
     535     * the steps as large as possible. This method also updates t if a step is
     536     * performed. Also, the internal order of the stepper is adjusted if requried.
     537     *
     538     * \param system The system function to solve, hence the r.h.s. of the ODE.
     539     * It must fullfil the Simple System concept.
     540     * \param x The state of the ODE which should be solved. Overwritten if
     541     * the step is successfull.
     542     * \param t The value of the time. Updated if the step is successfull.
     543     * \param dt The step size. Updated.
     544     * \return success if the step was accepted, fail otherwise.
     545     */
     546
     547    /**
     548     * \fn bulirsch_stoer::try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
     549     * \brief Tries to perform one step.
     550     *
     551     * This method tries to do one step with step size dt. If the error estimate
     552     * is to large, the step is rejected and the method returns fail and the
     553     * step size dt is reduced. If the error estimate is acceptably small, the
     554     * step is performed, success is returned and dt might be increased to make
     555     * the steps as large as possible. This method also updates t if a step is
     556     * performed. Also, the internal order of the stepper is adjusted if requried.
     557     *
     558     * \param system The system function to solve, hence the r.h.s. of the ODE.
     559     * It must fullfil the Simple System concept.
     560     * \param x The state of the ODE which should be solved. Overwritten if
     561     * the step is successfull.
     562     * \param dxdt The derivative of state.
     563     * \param t The value of the time. Updated if the step is successfull.
     564     * \param dt The step size. Updated.
     565     * \return success if the step was accepted, fail otherwise.
     566     */
     567
     568    /**
     569     * \fn bulirsch_stoer::try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
     570     * \brief Tries to perform one step.
     571     *
     572     * \note This method is disabled if state_type=time_type to avoid ambiguity.
     573     *
     574     * This method tries to do one step with step size dt. If the error estimate
     575     * is to large, the step is rejected and the method returns fail and the
     576     * step size dt is reduced. If the error estimate is acceptably small, the
     577     * step is performed, success is returned and dt might be increased to make
     578     * the steps as large as possible. This method also updates t if a step is
     579     * performed. Also, the internal order of the stepper is adjusted if requried.
     580     *
     581     * \param system The system function to solve, hence the r.h.s. of the ODE.
     582     * It must fullfil the Simple System concept.
     583     * \param in The state of the ODE which should be solved.
     584     * \param t The value of the time. Updated if the step is successfull.
     585     * \param out Used to store the result of the step.
     586     * \param dt The step size. Updated.
     587     * \return success if the step was accepted, fail otherwise.
     588     */
     589
     590
     591    /**
     592     * \fn bulirsch_stoer::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
     593     * \brief Tries to perform one step.
     594     *
     595     * This method tries to do one step with step size dt. If the error estimate
     596     * is to large, the step is rejected and the method returns fail and the
     597     * step size dt is reduced. If the error estimate is acceptably small, the
     598     * step is performed, success is returned and dt might be increased to make
     599     * the steps as large as possible. This method also updates t if a step is
     600     * performed. Also, the internal order of the stepper is adjusted if requried.
     601     *
     602     * \param system The system function to solve, hence the r.h.s. of the ODE.
     603     * It must fullfil the Simple System concept.
     604     * \param in The state of the ODE which should be solved.
     605     * \param dxdt The derivative of state.
     606     * \param t The value of the time. Updated if the step is successfull.
     607     * \param out Used to store the result of the step.
     608     * \param dt The step size. Updated.
     609     * \return success if the step was accepted, fail otherwise.
     610     */
     611
     612
     613    /**
     614     * \fn bulirsch_stoer::adjust_size( const StateIn &x )
     615     * \brief Adjust the size of all temporaries in the stepper manually.
     616     * \param x A state from which the size of the temporaries to be resized is deduced.
     617     */
     618
    503619}
    504620}
  • trunk/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp

    r81491 r81754  
    4747namespace numeric {
    4848namespace odeint {
    49 
    50 /*
    51 
    52 template< class T , size_t N >
    53 std::ostream& operator<<( std::ostream& output , const boost::array< T , N >& a ) {
    54     output << "[ " << a[0] ;
    55     for( size_t n = 1 ; n<N ; ++n )
    56         output << " , " << a[n];
    57     output << " ]";
    58     return output;  // for multiple << operators.
    59 }
    60 
    61 */
    62 
    63 
    64 
    6549
    6650template<
     
    8569    typedef Operations operations_type;
    8670    typedef Resizer resizer_type;
     71    typedef dense_output_stepper_tag stepper_category;
     72#ifndef DOXYGEN_SKIP
    8773    typedef state_wrapper< state_type > wrapped_state_type;
    8874    typedef state_wrapper< deriv_type > wrapped_deriv_type;
    89     typedef dense_output_stepper_tag stepper_category;
    9075
    9176    typedef bulirsch_stoer_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
     
    10186    typedef std::vector< wrapped_deriv_type > deriv_vector_type;
    10287    typedef std::vector< deriv_vector_type > deriv_table_type;
     88#endif //DOXYGEN_SKIP
    10389
    10490    const static size_t m_k_max = 8;
    105 
    10691
    10792
     
    130115        for( unsigned short i = 0; i < m_k_max+1; i++ )
    131116        {
     117            /* only this specific sequence allows for dense ouput */
    132118            m_interval_sequence[i] = 2 + 4*i;  // 2 6 10 14 ...
    133119            m_derivs[i].resize( m_interval_sequence[i] );
     
    141127                const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
    142128                m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
    143                 //std::cout << i << "," << k << " " << m_coeff[i][k] << '\t' ;
    144129            }
    145             //std ::cout << std::endl;
    146130            // crude estimate of optimal order
     131
     132            m_current_k_opt = 4;
     133            /* no calculation because log10 might not exist for value_type!
    147134            const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >( 1.0E-12 ) ) ) * 0.6 + 0.5 );
    148135            m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 1 , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>( m_k_max-1 ) , static_cast<int>( logfact ) ));
    149             //m_current_k_opt = m_k_max - 1;
    150             //std::cout << m_cost[i] << std::endl;
     136            */
    151137        }
    152138        int num = 1;
     
    154140        {
    155141            m_diffs[i].resize( num );
    156             //std::cout << "m_diffs[" << i << "] size: " << num << std::endl;
    157142            num += (i+1)%2;
    158143        }
    159144    }
    160 
    161 
    162 
    163 /*
    164   template< class System , class StateInOut >
    165   controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
    166   {
    167   m_xnew_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , detail::ref( *this ) , detail::_1 ) );
    168   controlled_step_result res = try_step( system , x , t , m_xnew.m_v , dt );
    169   if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged ) )
    170   {
    171   boost::numeric::odeint::copy( m_xnew.m_v , x );
    172   }
    173   return res;
    174   }
    175 */
    176145
    177146    template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
     
    185154
    186155        typename odeint::unwrap_reference< System >::type &sys = system;
    187 //        if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
    188 //            reset(); // system resized -> reset
    189 //        if( dt != m_dt_last )
    190 //            reset(); // step size changed from outside -> reset
    191156
    192157        bool reject( true );
     
    202167        for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
    203168        {
    204             //std::cout << "k=" << k <<" (steps=" << m_interval_sequence[k] << "): " << std::endl;
    205169            m_midpoint.set_steps( m_interval_sequence[k] );
    206170            if( k == 0 )
     
    218182                h_opt[k] = calc_h_opt( dt , error , k );
    219183                work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
    220                 //std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl;
    221                 //std::cout << '\t' << "error: " << error << std::endl;
    222184
    223185                m_k_final = k;
     
    306268                reject = true;
    307269                new_h = dt * pow BOOST_PREVENT_MACRO_SUBSTITUTION( error , static_cast<value_type>(-1)/(2*m_k_final+2) );
    308                 //std::cout << "####### rejected #######! (t=" << t << ") interpolation error: " << error << " , new dt: " << new_h << std::endl;
    309270            } else {
    310                 //std::cout << "####### accepted ####### - new k: " << m_current_k_opt << ", new stepsize: " << new_h << std::endl;
    311                 //increase time
    312271                t += dt;
    313272            }
    314         }// else  std::cout << "####### rejected #######!" << std::endl;
    315 
     273        }
    316274        //set next stepsize
    317275        if( !m_last_step_rejected || (new_h < dt) )
     
    368326    void calc_state( time_type t , StateOut &x ) const
    369327    {
    370         //std::cout << "===========" << std::endl << "doing interpolation for t=" << t << std::endl;
    371328        do_interpolation( t , x );
    372         //std::cout << "===========" << std::endl;
    373     }
    374 
     329    }
    375330
    376331    const state_type& current_state( void ) const
     
    399354    }
    400355
    401 
     356    /** \brief Resets the internal state of the stepper. */
    402357    void reset()
    403358    {
     
    405360        m_last_step_rejected = false;
    406361    }
    407 
    408362
    409363    template< class StateIn >
     
    439393    {
    440394        // result is written into table[0]
    441         //std::cout << "extrapolate k=" << k << ":" << std::endl;
    442395        static const value_type val1( 1.0 );
    443396        for( int j=k ; j>1 ; --j )
    444397        {
    445             //std::cout << '\t' << coeff[k + order_start_index][j + order_start_index - 1];
    446398            m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
    447399                                 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] ,
    448400                                                                                                           -coeff[k + order_start_index][j + order_start_index - 1] ) );
    449401        }
    450         //std::cout << std::endl << coeff[k + order_start_index][order_start_index] << std::endl;
    451402        m_algebra.for_each3( table[0].m_v , table[1].m_v , table[0].m_v ,
    452403                             typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][order_start_index] ,
     
    470421            fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( facmin/STEPFAC4 , min BOOST_PREVENT_MACRO_SUBSTITUTION( 1.0/facmin , fac ) );
    471422        }
    472         //return std::abs(h*fac); //std::min( 0.1 , std::abs(h*fac) );
    473423        return h*fac;
    474424    }
     
    500450    template< class StateIn1 , class DerivIn1 , class StateIn2 , class DerivIn2 >
    501451    value_type prepare_dense_output( int k , const StateIn1 &x_start , const DerivIn1 &dxdt_start ,
    502                                      const StateIn2 & /* x_end */ , const DerivIn2 & /*dxdt_end */ , time_type dt )  // k is the order to which the result was approximated
    503     {
    504 
    505         // compute the coefficients of the interpolation polynomial
    506         // we parametrize the interval t .. t+dt by theta = -1 .. 1
    507         // we use 2k+3 values at the interval center theta=0 to obtain the interpolation coefficients
    508         // the values are x(t+dt/2) and the derivatives dx/dt , ... d^(2k+2) x / dt^(2k+2) at the midpoints
    509         // the derivatives are approximated via finite differences
    510         // all values are obtained from interpolation of the results from the increasing orders of the midpoint calls
     452                                     const StateIn2 & /* x_end */ , const DerivIn2 & /*dxdt_end */ , time_type dt ) 
     453    /* k is the order to which the result was approximated */
     454    {
     455
     456        /* compute the coefficients of the interpolation polynomial
     457         * we parametrize the interval t .. t+dt by theta = -1 .. 1
     458         * we use 2k+3 values at the interval center theta=0 to obtain the interpolation coefficients
     459         * the values are x(t+dt/2) and the derivatives dx/dt , ... d^(2k+2) x / dt^(2k+2) at the midpoints
     460         * the derivatives are approximated via finite differences
     461         * all values are obtained from interpolation of the results from the increasing orders of the midpoint calls
     462         */
    511463
    512464        // calculate finite difference approximations to derivatives at the midpoint
     
    521473                f *= d;
    522474            }
    523             //std::cout << "x_mp[" << j << "] = " << m_mp_states[j].m_v << std::endl;
    524475
    525476            if( j > 0 )
     
    527478        }
    528479
    529         //std::cout << "a_0 = " << m_mp_states[0].m_v << std::endl;
    530 
    531480        time_type d = dt/2;
    532481
     
    538487
    539488            // extrapolation results are now stored in m_diffs[kappa][0]
    540             //std::cout << "extrapolation result: " << m_diffs[kappa][0].m_v << std::endl;
    541489
    542490            // divide kappa-th derivative by kappa because we need these terms for dense output interpolation
     
    544492
    545493            d *= dt/(2*(kappa+2));
    546 
    547             //std::cout << "a_" << kappa+1 << " = " << m_diffs[kappa][0].m_v << std::endl;
    548494        }
    549495
     
    561507
    562508        return error;
    563 
    564         // calculate coefficient a_{2k+1} = (2k+5)/4 * x_end - (2k+5)/4 * x_start - 1/4 * dxdt_end - 1/4 * dxdt_start + sum_i=0...k-1 (i-k-2)*a_{2i+1}
    565 
    566 
    567         //std::cout << std::endl << x_start << std::endl << x_end << std::endl;
    568         //std::cout << std::endl << dxdt_start << std::endl << dxdt_end << std::endl << std::endl;
    569 
    570         // we don't use additional terms in the polynomial, the following calculations are thus not required
    571 
    572 /*
    573   m_algebra.for_each5( m_a1.m_v , x_end , x_start , dxdt_end , dxdt_start ,
    574   typename operations_type::template scale_sum4< value_type >( static_cast<value_type>(2*k+5)/static_cast<value_type>(4),
    575   static_cast<value_type>(-2*k-5)/static_cast<value_type>(4),
    576   static_cast<value_type>(-dt)/static_cast<value_type>(8) ,
    577   static_cast<value_type>(-dt)/static_cast<value_type>(8) ) );
    578   for( int i = 0 ; i<k ; ++i )
    579   m_algebra.for_each3( m_a1.m_v , m_a1.m_v , m_diffs[2*i][0].m_v ,
    580   typename operations_type::template scale_sum2< value_type >( 1 , i-k-2 ) );
    581 
    582   //std::cout << "a_" << 2*k+1 << " = " << m_a1.m_v << std::endl;
    583 
    584   // calculate coefficient a_{2k+2} = (k+2)/2 * x_end + (k+2)/2 * x_start - 1/4 * dxdt_end + 1/4 * dxdt_start + (k+2)/2 * x_mp + sum_i=1...k (i-k-2)*a_{2i}
    585   m_algebra.for_each6( m_a2.m_v , x_end , x_start , dxdt_end , dxdt_start , m_mp_states[0].m_v ,
    586   typename operations_type::template scale_sum5< value_type >( static_cast<value_type>(k+2)/static_cast<value_type>(2),
    587   static_cast<value_type>(k+2)/static_cast<value_type>(2),
    588   static_cast<value_type>(-dt)/static_cast<value_type>(8) ,
    589   static_cast<value_type>(dt)/static_cast<value_type>(8) ,
    590   static_cast<value_type>(-k-2) ) );
    591   for( int i = 1 ; i<=k ; ++i )
    592   m_algebra.for_each3( m_a2.m_v , m_a2.m_v , m_diffs[2*i-1][0].m_v ,
    593   typename operations_type::template scale_sum2< value_type >( 1 , i-k-2 ) );
    594 
    595   //std::cout << "a_" << 2*k+2 << " = " << m_a2.m_v << std::endl;
    596 
    597   // calculate coefficient a_{2k+3} = -(2k+3)/4 * x_end + (2k+3)/4 * x_start + 1/4 * dxdt_end + 1/4 * dxdt_start + sum_i=0...k-1 (k+1-i)*a_{2i+1}
    598   m_algebra.for_each5( m_a3.m_v , x_end , x_start , dxdt_end , dxdt_start ,
    599   typename operations_type::template scale_sum4< value_type >( static_cast<value_type>(-2*k-3)/static_cast<value_type>(4),
    600   static_cast<value_type>(2*k+3)/static_cast<value_type>(4),
    601   static_cast<value_type>(dt)/static_cast<value_type>(8) ,
    602   static_cast<value_type>(dt)/static_cast<value_type>(8) ) );
    603   for( int i = 0 ; i<k ; ++i )
    604   m_algebra.for_each3( m_a3.m_v , m_a3.m_v , m_diffs[2*i][0].m_v ,
    605   typename operations_type::template scale_sum2< value_type >( 1 , k+1-i ) );
    606 
    607   //std::cout << "a_" << 2*k+3 << " = " << m_a3.m_v << std::endl;
    608 
    609   // calculate coefficient a_{2k+4} = -(k+1)/2 * x_end - (k+1)/2 * x_start + 1/4 * dxdt_end - 1/4 * dxdt_start - (k+1)/2 * x_mp + sum_i=0...k-1 (k+1-i)*a_{2i}
    610   m_algebra.for_each6( m_a4.m_v , x_end , x_start , dxdt_end , dxdt_start , m_mp_states[0].m_v ,
    611   typename operations_type::template scale_sum5< value_type >( static_cast<value_type>(-k-1)/static_cast<value_type>(2),
    612   static_cast<value_type>(-k-1)/static_cast<value_type>(2),
    613   static_cast<value_type>(dt)/static_cast<value_type>(8) ,
    614   static_cast<value_type>(-dt)/static_cast<value_type>(8),
    615   static_cast<value_type>(k+1) ) );
    616   for( int i = 1 ; i<=k ; ++i )
    617   m_algebra.for_each3( m_a4.m_v , m_a4.m_v , m_diffs[2*i-1][0].m_v ,
    618   typename operations_type::template scale_sum2< value_type >( 1 , k+1-i ) );
    619 */
    620         //std::cout << "a_" << 2*k+4 << " = " << m_a4.m_v << std::endl;
    621509    }
    622510
     
    629517            m_algebra.for_each2( m_diffs[0][j].m_v , m_derivs[j][m].m_v ,
    630518                                 typename operations_type::template scale_sum1< value_type >( fac ) );
    631             //std::cout << "j=" << j << ", kappa=" << kappa << ", m=" << m;
    632             //std::cout << ": m_diffs[" << kappa << "][" << j << "] = " << fac << " * f[" << m << "]  ";
    633             //std::cout << "(size(f)=" << m_derivs[j].size() << ") = " << m_diffs[0][j].m_v << std::endl;
    634 
    635519        }
    636520        else
    637521        {
    638             //std::cout << m_derivs[j][1].m_v << " , " << dxdt << std::endl;
    639 
    640522            // calculate the index of m_diffs for this kappa-j-combination
    641523            const int j_diffs = j - kappa/2;
    642 
    643             //std::cout << "j=" << j << ", kappa=" << kappa << ", m=" << m << ": m_diffs[" << kappa << "][" << j_diffs << "] = " << fac << " ( 1*f[" << m+kappa << "]";
    644524
    645525            m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v ,
     
    655535                                         typename operations_type::template scale_sum2< value_type , value_type >( 1.0 ,
    656536                                                                                                                   sign * fac * boost::math::binomial_coefficient< value_type >( kappa , c ) ) );
    657                     //std::cout << ( (sign > 0.0) ? " + " : " - " ) <<
    658                     //        boost::math::binomial_coefficient< double >( kappa , c ) << "*f[" << i << "]";
    659537                }
    660538                else
     
    662540                    m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , dxdt ,
    663541                                         typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , sign * fac ) );
    664                     //std::cout << ( (sign > 0.0) ? " + " : " - " ) << "dxdt";
    665542                }
    666543                sign *= -1;
    667544                ++c;
    668545            }
    669             //std::cout << " ) = " << m_diffs[kappa][j_diffs].m_v << std::endl;
    670546        }
    671547    }
     
    677553        // m_k_final is the number of order-iterations done for the last step - it governs the order of the interpolation polynomial
    678554        const value_type theta = 2 * get_unit_value( (t - m_t_last) / (m_t - m_t_last) ) - 1;
    679         //std::cout << "theta=" << theta << std::endl;
    680         //start with x = a0 + a_{2k+1} theta^{2k+1} + a_{2k+2} theta^{2k+2} + a_{2k+3} theta^{2k+3} + a_{2k+4} theta^{2k+4}
    681         //std::cout << "x = a_0 + ";
    682 
    683 /*        m_algebra.for_each6( out , m_mp_states[0].m_v , m_a1.m_v , m_a2.m_v , m_a3.m_v , m_a4.m_v ,
    684           typename operations_type::template scale_sum5< time_type >(
    685           static_cast<time_type>( 1 ) ,
    686           std::pow( theta , 2*m_k_final+1 ) ,
    687           std::pow( theta , 2*m_k_final+2 ) ,
    688           std::pow( theta , 2*m_k_final+3 ) ,
    689           std::pow( theta , 2*m_k_final+4 ) ) );
    690 */
    691 
    692555        // we use only values at interval center, that is theta=0, for interpolation
    693556        // our interpolation polynomial is thus of order 2k+2, hence we have 2k+3 terms
     
    698561        for( size_t i=0 ; i<=2*m_k_final+1 ; ++i )
    699562        {
    700             //std::cout << "a_" << i+1 << " theta^" << i+1 << " = " << m_diffs[i][0].m_v[0] * std::pow( theta , i+1 ) << std::endl;
    701563            m_algebra.for_each3( out , out , m_diffs[i][0].m_v ,
    702564                                 typename operations_type::template scale_sum2< value_type >( static_cast<value_type>(1) , theta_pow ) );
     
    824686};
    825687
     688
     689
     690/********** DOXYGEN **********/
     691
     692/**
     693 * \class bulirsch_stoer_dense_out
     694 * \brief The Bulirsch-Stoer algorithm.
     695 *
     696 * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
     697 * and order of the method. The algorithm uses the modified midpoint and
     698 * a polynomial extrapolation compute the solution. This class also provides
     699 * dense outout facility.
     700 *
     701 * \tparam State The state type.
     702 * \tparam Value The value type.
     703 * \tparam Deriv The type representing the time derivative of the state.
     704 * \tparam Time The time representing the independent variable - the time.
     705 * \tparam Algebra The algebra type.
     706 * \tparam Operations The operations type.
     707 * \tparam Resizer The resizer policy type.
     708 */
     709
     710    /**
     711     * \fn bulirsch_stoer_dense_out::bulirsch_stoer_dense_out( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt , bool control_interpolation )
     712     * \brief Constructs the bulirsch_stoer class, including initialization of
     713     * the error bounds.
     714     *
     715     * \param eps_abs Absolute tolerance level.
     716     * \param eps_rel Relative tolerance level.
     717     * \param factor_x Factor for the weight of the state.
     718     * \param factor_dxdt Factor for the weight of the derivative.
     719     * \param control_interpolation Set true to additionally control the error of
     720     * the interpolation.
     721     */
     722
     723    /**
     724     * \fn bulirsch_stoer_dense_out::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
     725     * \brief Tries to perform one step.
     726     *
     727     * This method tries to do one step with step size dt. If the error estimate
     728     * is to large, the step is rejected and the method returns fail and the
     729     * step size dt is reduced. If the error estimate is acceptably small, the
     730     * step is performed, success is returned and dt might be increased to make
     731     * the steps as large as possible. This method also updates t if a step is
     732     * performed. Also, the internal order of the stepper is adjusted if requried.
     733     *
     734     * \param system The system function to solve, hence the r.h.s. of the ODE.
     735     * It must fullfil the Simple System concept.
     736     * \param in The state of the ODE which should be solved.
     737     * \param dxdt The derivative of state.
     738     * \param t The value of the time. Updated if the step is successfull.
     739     * \param out Used to store the result of the step.
     740     * \param dt The step size. Updated.
     741     * \return success if the step was accepted, fail otherwise.
     742     */
     743
     744    /**
     745     * \fn bulirsch_stoer_dense_out::initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
     746     * \brief Initializes the dense output stepper.
     747     *
     748     * \param x0 The initial state.
     749     * \param t0 The initial time.
     750     * \param dt0 The initial time step.
     751     */
     752
     753    /**
     754     * \fn bulirsch_stoer_dense_out::do_step( System system )
     755     * \brief Does one time step. This is the main method that should be used to
     756     * integrate an ODE with this stepper.
     757     * \note initialize has to be called before using this method to set the
     758     * inital conditions x,t and the stepsize.
     759     * \param system The system function to solve, hence the r.h.s. of the
     760     * ordinary differential equation. It must fullfil the Simple System concept.
     761     * \return Pair with start and end time of the integration step.
     762     */
     763
     764    /**
     765     * \fn bulirsch_stoer_dense_out::calc_state( time_type t , StateOut &x ) const
     766     * \brief Calculates the solution at an intermediate point within the las step
     767     * \param t The time at which the solution should be calculated, has to be
     768     * in the current time interval.
     769     * \param x The output variable where the result is written into.
     770     */
     771
     772    /**
     773     * \fn bulirsch_stoer_dense_out::current_state( void ) const
     774     * \brief Returns the current state of the solution.
     775     * \return The current state of the solution x(t).
     776     */
     777
     778    /**
     779     * \fn bulirsch_stoer_dense_out::current_time( void ) const
     780     * \brief Returns the current time of the solution.
     781     * \return The current time of the solution t.
     782     */
     783
     784    /**
     785     * \fn bulirsch_stoer_dense_out::previous_state( void ) const
     786     * \brief Returns the last state of the solution.
     787     * \return The last state of the solution x(t-dt).
     788     */
     789
     790    /**
     791     * \fn bulirsch_stoer_dense_out::previous_time( void ) const
     792     * \brief Returns the last time of the solution.
     793     * \return The last time of the solution t-dt.
     794     */
     795
     796    /**
     797     * \fn bulirsch_stoer_dense_out::current_time_step( void ) const
     798     * \brief Returns the current step size.
     799     * \return The current step size.
     800     */
     801
     802    /**
     803     * \fn bulirsch_stoer_dense_out::adjust_size( const StateIn &x )
     804     * \brief Adjust the size of all temporaries in the stepper manually.
     805     * \param x A state from which the size of the temporaries to be resized is deduced.
     806     */
     807
    826808}
    827809}
  • trunk/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp

    r81491 r81754  
    4545
    4646
    47 /*
    48  * Error checker for controlled_error_stepper
    49  */
    5047template
    5148<
     
    6158    typedef Algebra algebra_type;
    6259    typedef Operations operations_type;
    63 
    6460
    6561    default_error_checker(
     
    121117
    122118
    123 
    124119/*
    125120 * explicit stepper version
     
    130125    * try_step( sys , in , t , out , dt )
    131126    * try_step( sys , in , dxdt , t , out , dt )
     127 */
     128/**
     129 * \brief Implements step size control for Runge-Kutta steppers with error
     130 * estimation.
     131 *
     132 * This class implements the step size control for standard Runge-Kutta
     133 * steppers with error estimation.
     134 *
     135 * \tparam ErrorStepper The stepper type with error estimation, has to fullfill the ErrorStepper concept.
     136 * \tparam ErrorChecker The error checker
     137 * \tparam Resizer The resizer policy type.
    132138 */
    133139template<
     
    151157    typedef ErrorChecker error_checker_type;
    152158    typedef explicit_controlled_stepper_tag stepper_category;
     159
     160#ifndef DOXYGEN_SKIP
    153161    typedef typename stepper_type::wrapped_state_type wrapped_state_type;
    154162    typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
    155163
    156164    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
    157 
    158 
    159 
     165#endif //DOXYGEN_SKIP
     166
     167
     168    /**
     169     * \brief Constructs the controlled Runge-Kutta stepper.
     170     * \param error_checker An instance of the error checker.
     171     * \param stepper An instance of the underlying stepper.
     172     */
    160173    controlled_runge_kutta(
    161174            const error_checker_type &error_checker = error_checker_type( ) ,
     
    172185     * The overloads are needed to solve the forwarding problem
    173186     */
     187    /**
     188     * \brief Tries to perform one step.
     189     *
     190     * This method tries to do one step with step size dt. If the error estimate
     191     * is to large, the step is rejected and the method returns fail and the
     192     * step size dt is reduced. If the error estimate is acceptably small, the
     193     * step is performed, success is returned and dt might be increased to make
     194     * the steps as large as possible. This method also updates t if a step is
     195     * performed.
     196     *
     197     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     198     *               Simple System concept.
     199     * \param x The state of the ODE which should be solved. Overwritten if
     200     * the step is successfull.
     201     * \param t The value of the time. Updated if the step is successfull.
     202     * \param dt The step size. Updated.
     203     * \return success if the step was accepted, fail otherwise.
     204     */
    174205    template< class System , class StateInOut >
    175206    controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
     
    178209    }
    179210
     211    /**
     212     * \brief Tries to perform one step. Solves the forwarding problem and
     213     * allows for using boost range as state_type.
     214     *
     215     * This method tries to do one step with step size dt. If the error estimate
     216     * is to large, the step is rejected and the method returns fail and the
     217     * step size dt is reduced. If the error estimate is acceptably small, the
     218     * step is performed, success is returned and dt might be increased to make
     219     * the steps as large as possible. This method also updates t if a step is
     220     * performed.
     221     *
     222     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     223     *               Simple System concept.
     224     * \param x The state of the ODE which should be solved. Overwritten if
     225     * the step is successfull. Can be a boost range.
     226     * \param t The value of the time. Updated if the step is successfull.
     227     * \param dt The step size. Updated.
     228     * \return success if the step was accepted, fail otherwise.
     229     */
    180230    template< class System , class StateInOut >
    181231    controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
     
    190240     *
    191241     * this version does not solve the forwarding problem, boost.range can not be used
     242     */
     243    /**
     244     * \brief Tries to perform one step.
     245     *
     246     * This method tries to do one step with step size dt. If the error estimate
     247     * is to large, the step is rejected and the method returns fail and the
     248     * step size dt is reduced. If the error estimate is acceptably small, the
     249     * step is performed, success is returned and dt might be increased to make
     250     * the steps as large as possible. This method also updates t if a step is
     251     * performed.
     252     *
     253     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     254     *               Simple System concept.
     255     * \param x The state of the ODE which should be solved. Overwritten if
     256     * the step is successfull.
     257     * \param dxdt The derivative of state.
     258     * \param t The value of the time. Updated if the step is successfull.
     259     * \param dt The step size. Updated.
     260     * \return success if the step was accepted, fail otherwise.
    192261     */
    193262    template< class System , class StateInOut , class DerivIn >
     
    210279     * the disable is needed to avoid ambiguous overloads if state_type = time_type
    211280     */
     281    /**
     282     * \brief Tries to perform one step.
     283     *
     284     * \note This method is disabled if state_type=time_type to avoid ambiguity.
     285     *
     286     * This method tries to do one step with step size dt. If the error estimate
     287     * is to large, the step is rejected and the method returns fail and the
     288     * step size dt is reduced. If the error estimate is acceptably small, the
     289     * step is performed, success is returned and dt might be increased to make
     290     * the steps as large as possible. This method also updates t if a step is
     291     * performed.
     292     *
     293     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     294     *               Simple System concept.
     295     * \param in The state of the ODE which should be solved.
     296     * \param t The value of the time. Updated if the step is successfull.
     297     * \param out Used to store the result of the step.
     298     * \param dt The step size. Updated.
     299     * \return success if the step was accepted, fail otherwise.
     300     */
    212301    template< class System , class StateIn , class StateOut >
    213302    typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
     
    225314     *
    226315     * this version does not solve the forwarding problem, boost.range can not be used
     316     */
     317    /**
     318     * \brief Tries to perform one step.
     319     *
     320     * This method tries to do one step with step size dt. If the error estimate
     321     * is to large, the step is rejected and the method returns fail and the
     322     * step size dt is reduced. If the error estimate is acceptably small, the
     323     * step is performed, success is returned and dt might be increased to make
     324     * the steps as large as possible. This method also updates t if a step is
     325     * performed.
     326     *
     327     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     328     *               Simple System concept.
     329     * \param in The state of the ODE which should be solved.
     330     * \param dxdt The derivative of state.
     331     * \param t The value of the time. Updated if the step is successfull.
     332     * \param out Used to store the result of the step.
     333     * \param dt The step size. Updated.
     334     * \return success if the step was accepted, fail otherwise.
    227335     */
    228336    template< class System , class StateIn , class DerivIn , class StateOut >
     
    252360            if( m_max_rel_error < 0.5 )
    253361            {
     362                // error should be > 0
     363                m_max_rel_error = max BOOST_PREVENT_MACRO_SUBSTITUTION ( pow( 5.0 , -m_stepper.stepper_order() ) , m_max_rel_error );
    254364                //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
    255365                t += dt;
    256                 dt *= min BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( m_max_rel_error ,
    257                                                                static_cast<value_type>(-1) / m_stepper.stepper_order() ) ,
    258                            static_cast<value_type>(5) );
     366                dt *= static_cast<value_type>(9)/static_cast<value_type>(10) * pow( m_max_rel_error ,
     367                                                                                    static_cast<value_type>(-1) / m_stepper.stepper_order() );
    259368                return success;
    260369            }
     
    267376    }
    268377
     378    /**
     379     * \brief Returns the error of the last step.
     380     *
     381     * returns The last error of the step.
     382     */
    269383    value_type last_error( void ) const
    270384    {
     
    273387
    274388
    275 
     389    /**
     390     * \brief Adjust the size of all temporaries in the stepper manually.
     391     * \param x A state from which the size of the temporaries to be resized is deduced.
     392     */
    276393    template< class StateType >
    277394    void adjust_size( const StateType &x )
     
    283400    }
    284401
     402    /**
     403     * \brief Returns the instance of the underlying stepper.
     404     * \returns The instance of the underlying stepper.
     405     */
    285406    stepper_type& stepper( void )
    286407    {
     
    288409    }
    289410
     411    /**
     412     * \brief Returns the instance of the underlying stepper.
     413     * \returns The instance of the underlying stepper.
     414     */
    290415    const stepper_type& stepper( void ) const
    291416    {
     
    355480    * try_step( sys , x , dxdt , t , dt )
    356481    * try_step( sys , in , dxdt_in , t , out , dxdt_out , dt )
     482 */
     483/**
     484 * \brief Implements step size control for Runge-Kutta FSAL steppers with
     485 * error estimation.
     486 *
     487 * This class implements the step size control for FSAL Runge-Kutta
     488 * steppers with error estimation.
     489 *
     490 * \tparam ErrorStepper The stepper type with error estimation, has to fullfill the ErrorStepper concept.
     491 * \tparam ErrorChecker The error checker
     492 * \tparam Resizer The resizer policy type.
    357493 */
    358494template<
     
    376512    typedef ErrorChecker error_checker_type;
    377513    typedef explicit_controlled_stepper_fsal_tag stepper_category;
     514
     515#ifndef DOXYGEN_SKIP
    378516    typedef typename stepper_type::wrapped_state_type wrapped_state_type;
    379517    typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
    380518
    381519    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
    382 
     520#endif // DOXYGEN_SKIP
     521
     522    /**
     523     * \brief Constructs the controlled Runge-Kutta stepper.
     524     * \param error_checker An instance of the error checker.
     525     * \param stepper An instance of the underlying stepper.
     526     */
    383527    controlled_runge_kutta(
    384528            const error_checker_type &error_checker = error_checker_type() ,
     
    394538     * The two overloads are needed in order to solve the forwarding problem
    395539     */
     540    /**
     541     * \brief Tries to perform one step.
     542     *
     543     * This method tries to do one step with step size dt. If the error estimate
     544     * is to large, the step is rejected and the method returns fail and the
     545     * step size dt is reduced. If the error estimate is acceptably small, the
     546     * step is performed, success is returned and dt might be increased to make
     547     * the steps as large as possible. This method also updates t if a step is
     548     * performed.
     549     *
     550     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     551     *               Simple System concept.
     552     * \param x The state of the ODE which should be solved. Overwritten if
     553     * the step is successfull.
     554     * \param t The value of the time. Updated if the step is successfull.
     555     * \param dt The step size. Updated.
     556     * \return success if the step was accepted, fail otherwise.
     557     */
    396558    template< class System , class StateInOut >
    397559    controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
     
    400562    }
    401563
     564
     565    /**
     566     * \brief Tries to perform one step. Solves the forwarding problem and
     567     * allows for using boost range as state_type.
     568     *
     569     * This method tries to do one step with step size dt. If the error estimate
     570     * is to large, the step is rejected and the method returns fail and the
     571     * step size dt is reduced. If the error estimate is acceptably small, the
     572     * step is performed, success is returned and dt might be increased to make
     573     * the steps as large as possible. This method also updates t if a step is
     574     * performed.
     575     *
     576     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     577     *               Simple System concept.
     578     * \param x The state of the ODE which should be solved. Overwritten if
     579     * the step is successfull. Can be a boost range.
     580     * \param t The value of the time. Updated if the step is successfull.
     581     * \param dt The step size. Updated.
     582     * \return success if the step was accepted, fail otherwise.
     583     */
    402584    template< class System , class StateInOut >
    403585    controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
     
    414596     *
    415597     * The disabler is needed to solve ambigous overloads
     598     */
     599    /**
     600     * \brief Tries to perform one step.
     601     *
     602     * \note This method is disabled if state_type=time_type to avoid ambiguity.
     603     *
     604     * This method tries to do one step with step size dt. If the error estimate
     605     * is to large, the step is rejected and the method returns fail and the
     606     * step size dt is reduced. If the error estimate is acceptably small, the
     607     * step is performed, success is returned and dt might be increased to make
     608     * the steps as large as possible. This method also updates t if a step is
     609     * performed.
     610     *
     611     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     612     *               Simple System concept.
     613     * \param in The state of the ODE which should be solved.
     614     * \param t The value of the time. Updated if the step is successfull.
     615     * \param out Used to store the result of the step.
     616     * \param dt The step size. Updated.
     617     * \return success if the step was accepted, fail otherwise.
    416618     */
    417619    template< class System , class StateIn , class StateOut >
     
    431633     *
    432634     * This version does not solve the forwarding problem, boost::range can not be used.
     635     */
     636    /**
     637     * \brief Tries to perform one step.
     638     *
     639     * This method tries to do one step with step size dt. If the error estimate
     640     * is to large, the step is rejected and the method returns fail and the
     641     * step size dt is reduced. If the error estimate is acceptably small, the
     642     * step is performed, success is returned and dt might be increased to make
     643     * the steps as large as possible. This method also updates t if a step is
     644     * performed.
     645     *
     646     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     647     *               Simple System concept.
     648     * \param x The state of the ODE which should be solved. Overwritten if
     649     * the step is successfull.
     650     * \param dxdt The derivative of state.
     651     * \param t The value of the time. Updated if the step is successfull.
     652     * \param dt The step size. Updated.
     653     * \return success if the step was accepted, fail otherwise.
    433654     */
    434655    template< class System , class StateInOut , class DerivInOut >
     
    451672     *
    452673     * This version does not solve the forwarding problem, boost::range can not be used.
     674     */
     675    /**
     676     * \brief Tries to perform one step.
     677     *
     678     * This method tries to do one step with step size dt. If the error estimate
     679     * is to large, the step is rejected and the method returns fail and the
     680     * step size dt is reduced. If the error estimate is acceptably small, the
     681     * step is performed, success is returned and dt might be increased to make
     682     * the steps as large as possible. This method also updates t if a step is
     683     * performed.
     684     *
     685     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     686     *               Simple System concept.
     687     * \param in The state of the ODE which should be solved.
     688     * \param dxdt The derivative of state.
     689     * \param t The value of the time. Updated if the step is successfull.
     690     * \param out Used to store the result of the step.
     691     * \param dt The step size. Updated.
     692     * \return success if the step was accepted, fail otherwise.
    453693     */
    454694    template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
     
    480720            if( max_rel_err < 0.5 )
    481721            {                //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
     722                // error should be > 0
     723                max_rel_err = max BOOST_PREVENT_MACRO_SUBSTITUTION ( pow( 5.0 , -m_stepper.stepper_order() ) , max_rel_err );
    482724                t += dt;
    483                 dt *= min BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) ) , static_cast<value_type>(5) );
     725                dt *= static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) );
    484726                return success;
    485727            }
     
    493735
    494736
    495 
     737    /**
     738     * \brief Resets the internal state of the underlying FSAL stepper.
     739     */
    496740    void reset( void )
    497741    {
     
    499743    }
    500744
     745    /**
     746     * \brief Initializes the internal state storing an internal copy of the derivative.
     747     *
     748     * \param deriv The initial derivative of the ODE.
     749     */
    501750    template< class DerivIn >
    502751    void initialize( const DerivIn &deriv )
     
    506755    }
    507756
     757    /**
     758     * \brief Initializes the internal state storing an internal copy of the derivative.
     759     *
     760     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     761     *               Simple System concept.
     762     * \param x The initial state of the ODE which should be solved.
     763     * \param t The initial time.
     764     */
    508765    template< class System , class StateIn >
    509766    void initialize( System system , const StateIn &x , time_type t )
     
    514771    }
    515772
     773    /**
     774     * \brief Returns true if the stepper has been initialized, false otherwise.
     775     *
     776     * \return true, if the stepper has been initialized, false otherwise.
     777     */
    516778    bool is_initialized( void ) const
    517779    {
     
    520782
    521783
    522 
     784    /**
     785     * \brief Adjust the size of all temporaries in the stepper manually.
     786     * \param x A state from which the size of the temporaries to be resized is deduced.
     787     */
    523788    template< class StateType >
    524789    void adjust_size( const StateType &x )
     
    531796
    532797
     798    /**
     799     * \brief Returns the instance of the underlying stepper.
     800     * \returns The instance of the underlying stepper.
     801     */
    533802    stepper_type& stepper( void )
    534803    {
     
    536805    }
    537806
     807    /**
     808     * \brief Returns the instance of the underlying stepper.
     809     * \returns The instance of the underlying stepper.
     810     */
    538811    const stepper_type& stepper( void ) const
    539812    {
     
    598871
    599872
     873/********** DOXYGEN **********/
     874
     875/**** DEFAULT ERROR CHECKER ****/
     876
     877/**
     878 * \class default_error_checker
     879 * \brief The default error checker to be used with Runge-Kutta error steppers
     880 *
     881 * This class provides the default mechanism to compare the error estimates
     882 * reported by Runge-Kutta error steppers with user defined error bounds.
     883 * It is used by the controlled_runge_kutta steppers.
     884 *
     885 * \tparam Value The value type.
     886 * \tparam Algebra The algebra type.
     887 * \tparam Operations The operations type.
     888 */
     889
     890    /**
     891     * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt )
     892     * \brief Constructs the error checker.
     893     *
     894     * The error is calculated as follows: ????
     895     *
     896     * \param eps_abs Absolute tolerance level.
     897     * \param eps_rel Relative tolerance level.
     898     * \param a_x Factor for the weight of the state.
     899     * \param a_dxdt Factor for the weight of the derivative.
     900     */
     901   
     902    /**
     903     * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
     904     * \brief Calculates the error level.
     905     *
     906     * If the returned error level is greater than 1, the estimated error was
     907     * larger than the permitted error bounds and the step should be repeated
     908     * with a smaller step size.
     909     *
     910     * \param x_old State at the beginning of the step.
     911     * \param dxdt_old Derivative at the beginning of the step.
     912     * \param x_err Error estimate.
     913     * \param dt Time step.
     914     * \return error
     915     */
     916
     917    /**
     918     * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
     919     * \brief Calculates the error level using a given algebra.
     920     *
     921     * If the returned error level is greater than 1, the estimated error was
     922     * larger than the permitted error bounds and the step should be repeated
     923     * with a smaller step size.
     924     *
     925     * \param algebra The algebra used for calculation of the error.
     926     * \param x_old State at the beginning of the step.
     927     * \param dxdt_old Derivative at the beginning of the step.
     928     * \param x_err Error estimate.
     929     * \param dt Time step.
     930     * \return error
     931     */
    600932
    601933
  • trunk/boost/numeric/odeint/stepper/controlled_step_result.hpp

    r81491 r81754  
    2525
    2626/**
    27  * \enum Enum representing the return values of the controlled steppers.
     27 * \enum controlled_step_result
     28 *
     29 * \brief Enum representing the return values of the controlled steppers.
    2830 */
    2931typedef enum
    3032{
    31     success , /** <The trial step war successful, hence the state and the time have been advanced. */
    32     fail      /** <The step was not succesful and might possibly be repeated with a small step size. */
     33    success , /**< The trial step was successful, hence the state and the time have been advanced. */
     34    fail      /**< The step was not succesful and might possibly be repeated with a small step size. */
    3335} controlled_step_result;
    3436
  • trunk/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp

    r81491 r81754  
    4343
    4444
    45 
    46 
     45/**
     46 * \brief The class representing dense-output Runge-Kutta steppers.
     47 * \note In this stepper, the initialize method has to be called before using
     48 * the do_step method.
     49 *
     50 * The dense-output functionality allows to interpolate the solution between
     51 * subsequent integration points using intermediate results obtained during the
     52 * computation. This version works based on a normal stepper without step-size
     53 * control.
     54 *
     55 *
     56 * \tparam Stepper The stepper type of the underlying algorithm.
     57 */
    4758template< class Stepper >
    4859class dense_output_runge_kutta< Stepper , stepper_tag >
     
    6879
    6980
    70 
     81    /**
     82     * \brief Constructs the dense_output_runge_kutta class. An instance of the
     83     * underlying stepper can be provided.
     84     * \param stepper An instance of the underlying stepper.
     85     */
    7186    dense_output_runge_kutta( const stepper_type &stepper = stepper_type() )
    7287    : m_stepper( stepper ) , m_resizer() ,
     
    7691
    7792
     93    /**
     94     * \brief Initializes the stepper. Has to be called before do_step can be
     95     * used to set the initial conditions and the step size.
     96     * \param x0 The inital state of the ODE which should be solved.
     97     * \param t0 The initial time, at which the step should be performed.
     98     * \param dt0 The step size.
     99     */
    78100    template< class StateType >
    79101    void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
     
    85107    }
    86108
     109    /**
     110     * \brief Does one time step.
     111     * \note initialize has to be called before using this method to set the
     112     * inital conditions x,t and the stepsize.
     113     * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fullfil the
     114     *               Simple System concept.
     115     * \return Pair with start and end time of the integration step.
     116     */
    87117    template< class System >
    88118    std::pair< time_type , time_type > do_step( System system )
     
    98128     * The next two overloads are needed to solve the forwarding problem
    99129     */
     130   
     131    /**
     132     * \brief Calculates the solution at an intermediate point.
     133     * \param t The time at which the solution should be calculated, has to be
     134     * in the current time interval.
     135     * \param x The output variable where the result is written into.
     136     */
    100137    template< class StateOut >
    101138    void calc_state( time_type t , StateOut &x ) const
     
    104141    }
    105142
     143    /**
     144     * \brief Calculates the solution at an intermediate point. Solves the forwarding problem
     145     * \param t The time at which the solution should be calculated, has to be
     146     * in the current time interval.
     147     * \param x The output variable where the result is written into, can be a boost range.
     148     */
    106149    template< class StateOut >
    107150    void calc_state( time_type t , const StateOut &x ) const
     
    110153    }
    111154
     155    /**
     156     * \brief Adjust the size of all temporaries in the stepper manually.
     157     * \param x A state from which the size of the temporaries to be resized is deduced.
     158     */
    112159    template< class StateType >
    113160    void adjust_size( const StateType &x )
     
    117164    }
    118165
     166    /**
     167     * \brief Returns the current state of the solution.
     168     * \return The current state of the solution x(t).
     169     */
    119170    const state_type& current_state( void ) const
    120171    {
     
    122173    }
    123174
     175    /**
     176     * \brief Returns the current time of the solution.
     177     * \return The current time of the solution t.
     178     */
    124179    time_type current_time( void ) const
    125180    {
     
    127182    }
    128183
     184    /**
     185     * \brief Returns the last state of the solution.
     186     * \return The last state of the solution x(t-dt).
     187     */
    129188    const state_type& previous_state( void ) const
    130189    {
     
    132191    }
    133192
     193    /**
     194     * \brief Returns the last time of the solution.
     195     * \return The last time of the solution t-dt.
     196     */
    134197    time_type previous_time( void ) const
    135198    {
     
    188251
    189252
    190 
     253/**
     254 * \brief The class representing dense-output Runge-Kutta steppers with FSAL property.
     255 *
     256 * The interface is the same as for dense_output_runge_kutta< Stepper , stepper_tag >.
     257 * This class provides dense output functionality based on methods with step size controlled
     258 *
     259 *
     260 * \tparam Stepper The stepper type of the underlying algorithm.
     261 */
    191262template< class Stepper >
    192263class dense_output_runge_kutta< Stepper , explicit_controlled_stepper_fsal_tag >
  • trunk/boost/numeric/odeint/stepper/euler.hpp

    r81491 r81754  
    3131
    3232
    33 
    34 /**
    35  * \class euler
    36  * \brief An implementation of the Euler method.
    37  *
    38  * The Euler method is a very simply solver for ordinary differential equations. This method should not be used
    39  * for real applications. It is only useful for demonstration purposes. Step size control is not provided but
    40  * trivial continous output is available.
    41  *
    42  * This class derives from explicit_stepper_base and inherits its interface via CRTP (current recurring template pattern),
    43  * see explicit_stepper_base
    44  *
    45  * \tparam State The state type.
    46  * \tparam Value The value type.
    47  * \tparam Deriv The type representing the time derivative of the state.
    48  * \tparam Time The time representing the independent variable - the time.
    49  * \tparam Algebra The algebra type.
    50  * \tparam Operations The operations type.
    51  * \tparam Resizer The resizer policy type.
    52  */
    5333template<
    5434class State ,
     
    6646  1 , State , Value , Deriv , Time , Algebra , Operations , Resizer >
    6747#else
    68 class euler : public explicit_stepper_base< euler< ... > , ... >
     48class euler : public explicit_stepper_base
    6949#endif
    7050{
     
    9171
    9272
     73    euler( const algebra_type &algebra = algebra_type() ) : stepper_base_type( algebra )
     74    { }
     75
     76    template< class System , class StateIn , class DerivIn , class StateOut >
     77    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
     78    {
     79        stepper_base_type::m_algebra.for_each3( out , in , dxdt ,
     80                typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) );
     81
     82    }
     83
     84    template< class StateOut , class StateIn1 , class StateIn2 >
     85    void calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 &current_state , time_type t_new ) const
     86    {
     87        const time_type delta = t - t_old;
     88        stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v ,
     89                typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , delta ) );
     90    }
     91
     92    template< class StateType >
     93    void adjust_size( const StateType &x )
     94    {
     95        stepper_base_type::adjust_size( x );
     96    }
     97};
     98
     99
     100
     101/********** DOXYGEN ***********/
     102
     103/**
     104 * \class euler
     105 * \brief An implementation of the Euler method.
     106 *
     107 * The Euler method is a very simply solver for ordinary differential equations. This method should not be used
     108 * for real applications. It is only useful for demonstration purposes. Step size control is not provided but
     109 * trivial continous output is available.
     110 *
     111 * This class derives from explicit_stepper_base and inherits its interface via CRTP (current recurring template pattern),
     112 * see explicit_stepper_base
     113 *
     114 * \tparam State The state type.
     115 * \tparam Value The value type.
     116 * \tparam Deriv The type representing the time derivative of the state.
     117 * \tparam Time The time representing the independent variable - the time.
     118 * \tparam Algebra The algebra type.
     119 * \tparam Operations The operations type.
     120 * \tparam Resizer The resizer policy type.
     121 */
     122
    93123    /**
     124     * \fn euler::euler( const algebra_type &algebra )
    94125     * \brief Constructs the euler class. This constructor can be used as a default
    95126     * constructor of the algebra has a default constructor.
    96127     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    97128     */
    98     euler( const algebra_type &algebra = algebra_type() ) : stepper_base_type( algebra )
    99     { }
    100 
     129   
    101130    /**
     131     * \fn euler::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    102132     * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method.
    103      * The result is updated out of place, hence the input is in `in` and the output in `out`. `do_step_impl` is
    104      * used by explicit_stepper_base.
     133     * The result is updated out of place, hence the input is in `in` and the output in `out`.
     134     * Access to this step functionality is provided by explicit_stepper_base and
     135     * `do_step_impl` should not be called directly.
    105136     *
    106137     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     
    112143     * \param dt The step size.
    113144     */
    114     template< class System , class StateIn , class DerivIn , class StateOut >
    115     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    116     {
    117         stepper_base_type::m_algebra.for_each3( out , in , dxdt ,
    118                 typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) );
    119 
    120     }
    121145
    122146
    123147    /**
     148     * \fn euler::calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 &current_state , time_type t_new ) const
    124149     * \brief This method is used for continous output and it calculates the state `x` at a time `t` from the
    125150     * knowledge of two states `old_state` and `current_state` at time points `t_old` and `t_new`.
    126151     */
    127     template< class StateOut , class StateIn1 , class StateIn2 >
    128     void calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 &current_state , time_type t_new ) const
    129     {
    130         const time_type delta = t - t_old;
    131         stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v ,
    132                 typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , delta ) );
    133     }
    134152
    135153    /**
     154     * \fn euler::adjust_size( const StateType &x )
    136155     * \brief Adjust the size of all temporaries in the stepper manually.
    137156     * \param x A state from which the size of the temporaries to be resized is deduced.
    138157     */
    139     template< class StateType >
    140     void adjust_size( const StateType &x )
    141     {
    142         stepper_base_type::adjust_size( x );
    143     }
    144 };
    145 
    146 
    147 
    148158
    149159} // odeint
  • trunk/boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp

    r81491 r81754  
    3737namespace odeint {
    3838
     39
     40template<
     41size_t StageCount,
     42size_t Order,
     43size_t StepperOrder ,
     44size_t ErrorOrder ,
     45class State ,
     46class Value = double ,
     47class Deriv = State ,
     48class Time = Value ,
     49class Algebra = range_algebra ,
     50class Operations = default_operations ,
     51class Resizer = initially_resizer
     52>
     53#ifndef DOXYGEN_SKIP
     54class explicit_error_generic_rk
     55: public explicit_error_stepper_base<
     56  explicit_error_generic_rk< StageCount , Order , StepperOrder , ErrorOrder , State ,
     57  Value , Deriv , Time , Algebra , Operations , Resizer > ,
     58  Order , StepperOrder , ErrorOrder , State , Value , Deriv , Time , Algebra ,
     59  Operations , Resizer >
     60#else
     61class explicit_error_generic_rk : public explicit_error_stepper_base
     62#endif
     63{
     64
     65public:
     66#ifndef DOXYGEN_SKIP
     67    typedef explicit_error_stepper_base<
     68            explicit_error_generic_rk< StageCount , Order , StepperOrder , ErrorOrder , State ,
     69            Value , Deriv , Time , Algebra , Operations , Resizer > ,
     70            Order , StepperOrder , ErrorOrder , State , Value , Deriv , Time , Algebra ,
     71            Operations , Resizer > stepper_base_type;
     72#else
     73    typedef explicit_stepper_base< ... > stepper_base_type;
     74#endif
     75    typedef typename stepper_base_type::state_type state_type;
     76    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
     77    typedef typename stepper_base_type::value_type value_type;
     78    typedef typename stepper_base_type::deriv_type deriv_type;
     79    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
     80    typedef typename stepper_base_type::time_type time_type;
     81    typedef typename stepper_base_type::algebra_type algebra_type;
     82    typedef typename stepper_base_type::operations_type operations_type;
     83    typedef typename stepper_base_type::resizer_type resizer_type;
     84#ifndef DOXYGEN_SKIP
     85    typedef explicit_error_generic_rk< StageCount , Order , StepperOrder , ErrorOrder , State ,
     86            Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type;
     87#endif
     88    typedef detail::generic_rk_algorithm< StageCount , Value , Algebra , Operations > rk_algorithm_type;
     89
     90    typedef typename rk_algorithm_type::coef_a_type coef_a_type;
     91    typedef typename rk_algorithm_type::coef_b_type coef_b_type;
     92    typedef typename rk_algorithm_type::coef_c_type coef_c_type;
     93
     94    static const size_t stage_count = StageCount;
     95
     96private:
     97
     98
     99public:
     100
     101    // we use an explicit_generic_rk to do the normal rk step
     102    // and add a separate calculation of the error estimate afterwards
     103    explicit_error_generic_rk( const coef_a_type &a ,
     104            const coef_b_type &b ,
     105            const coef_b_type &b2 ,
     106            const coef_c_type &c ,
     107            const algebra_type &algebra = algebra_type() )
     108    : stepper_base_type( algebra ) , m_rk_algorithm( a , b , c ) , m_b2( b2 )
     109    { }
     110
     111
     112    template< class System , class StateIn , class DerivIn , class StateOut , class Err >
     113    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
     114            time_type t , StateOut &out , time_type dt , Err &xerr )
     115    {
     116        // normal step
     117        do_step_impl( system , in , dxdt , t , out , dt );
     118
     119        // additionally, perform the error calculation
     120        detail::template generic_rk_call_algebra< StageCount , algebra_type >()( stepper_base_type::m_algebra ,
     121                xerr , dxdt , m_F , detail::generic_rk_scale_sum_err< StageCount , operations_type , value_type , time_type >( m_b2 , dt) );
     122    }
     123
     124
     125    template< class System , class StateIn , class DerivIn , class StateOut >
     126    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
     127            time_type t , StateOut &out , time_type dt )
     128    {
     129        m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
     130
     131        // actual calculation done in generic_rk.hpp
     132        m_rk_algorithm.do_step( stepper_base_type::m_algebra , system , in , dxdt , t , out , dt , m_x_tmp.m_v , m_F );
     133    }
     134
     135
     136    template< class StateIn >
     137    void adjust_size( const StateIn &x )
     138    {
     139        resize_impl( x );
     140        stepper_base_type::adjust_size( x );
     141    }
     142
     143
     144private:
     145
     146    template< class StateIn >
     147    bool resize_impl( const StateIn &x )
     148    {
     149        bool resized( false );
     150        resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() );
     151        for( size_t i = 0 ; i < StageCount-1 ; ++i )
     152        {
     153            resized |= adjust_size_by_resizeability( m_F[i] , x , typename is_resizeable<deriv_type>::type() );
     154        }
     155        return resized;
     156    }
     157
     158
     159    rk_algorithm_type m_rk_algorithm;
     160    coef_b_type m_b2;
     161
     162    resizer_type m_resizer;
     163
     164    wrapped_state_type m_x_tmp;
     165    wrapped_deriv_type m_F[StageCount-1];
     166
     167};
     168
     169
     170/********* DOXYGEN *********/
    39171
    40172/**
     
    61193 * \tparam Resizer The resizer policy type.
    62194 */
    63 template<
    64 size_t StageCount,
    65 size_t Order,
    66 size_t StepperOrder ,
    67 size_t ErrorOrder ,
    68 class State ,
    69 class Value = double ,
    70 class Deriv = State ,
    71 class Time = Value ,
    72 class Algebra = range_algebra ,
    73 class Operations = default_operations ,
    74 class Resizer = initially_resizer
    75 >
    76 class explicit_error_generic_rk
    77 : public explicit_error_stepper_base<
    78   explicit_error_generic_rk< StageCount , Order , StepperOrder , ErrorOrder , State ,
    79   Value , Deriv , Time , Algebra , Operations , Resizer > ,
    80   Order , StepperOrder , ErrorOrder , State , Value , Deriv , Time , Algebra ,
    81   Operations , Resizer >
    82 {
    83 
    84 public:
    85 
    86     typedef explicit_error_stepper_base<
    87             explicit_error_generic_rk< StageCount , Order , StepperOrder , ErrorOrder , State ,
    88             Value , Deriv , Time , Algebra , Operations , Resizer > ,
    89             Order , StepperOrder , ErrorOrder , State , Value , Deriv , Time , Algebra ,
    90             Operations , Resizer > stepper_base_type;
    91 
    92     typedef typename stepper_base_type::state_type state_type;
    93     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
    94     typedef typename stepper_base_type::value_type value_type;
    95     typedef typename stepper_base_type::deriv_type deriv_type;
    96     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
    97     typedef typename stepper_base_type::time_type time_type;
    98     typedef typename stepper_base_type::algebra_type algebra_type;
    99     typedef typename stepper_base_type::operations_type operations_type;
    100     typedef typename stepper_base_type::resizer_type resizer_type;
    101     //typedef typename stepper_base_type::stepper_type stepper_type;
    102     typedef explicit_error_generic_rk< StageCount , Order , StepperOrder , ErrorOrder , State ,
    103             Value , Deriv , Time , Algebra , Operations , Resizer > stepper_type;
    104 
    105     typedef detail::generic_rk_algorithm< StageCount , Value , Algebra , Operations > rk_algorithm_type;
    106 
    107     typedef typename rk_algorithm_type::coef_a_type coef_a_type;
    108     typedef typename rk_algorithm_type::coef_b_type coef_b_type;
    109     typedef typename rk_algorithm_type::coef_c_type coef_c_type;
    110 
    111     static const size_t stage_count = StageCount;
    112 
    113 private:
    114 
    115 
    116 public:
    117 
    118     // we use an explicit_generic_rk to do the normal rk step
    119     // and add a separate calculation of the error estimate afterwards
    120     explicit_error_generic_rk( const coef_a_type &a ,
    121             const coef_b_type &b ,
    122             const coef_b_type &b2 ,
    123             const coef_c_type &c ,
    124             const algebra_type &algebra = algebra_type() )
    125     : stepper_base_type( algebra ) , m_rk_algorithm( a , b , c ) , m_b2( b2 )
    126     { }
    127 
    128 
    129     /**
     195
     196
     197    /**
     198     * \fn explicit_error_generic_rk::explicit_error_generic_rk( const coef_a_type &a , const coef_b_type &b , const coef_b_type &b2 , const coef_c_type &c , const algebra_type &algebra )
     199     * \brief Constructs the explicit_error_generik_rk class with the given parameters a, b, b2 and c. See examples section for details on the coefficients.
     200     *
     201     * \param a Triangular matrix of parameters b in the Butcher tableau.
     202     * \param b Last row of the butcher tableau.
     203     * \param b2 Parameters for lower-order evaluation to estimate the error.
     204     * \param c Parameters to calculate the time points in the Butcher tableau.
     205     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
     206     */
     207
     208
     209    /**
     210     * \fn explicit_error_generic_rk::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt , Err &xerr )
    130211     * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method.
    131212     * The result is updated out-of-place, hence the input is in `in` and the output in `out`. Futhermore, an
     
    142223     */
    143224
    144     template< class System , class StateIn , class DerivIn , class StateOut , class Err >
    145     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
    146             time_type t , StateOut &out , time_type dt , Err &xerr )
    147     {
    148         // normal step
    149         do_step_impl( system , in , dxdt , t , out , dt );
    150 
    151         // additionally, perform the error calculation
    152         detail::template generic_rk_call_algebra< StageCount , algebra_type >()( stepper_base_type::m_algebra ,
    153                 xerr , dxdt , m_F , detail::generic_rk_scale_sum_err< StageCount , operations_type , value_type , time_type >( m_b2 , dt) );
    154     }
    155 
    156 
    157     /**
     225    /**
     226     * \fn explicit_error_generic_rk::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    158227     * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method.
    159      * The result is updated out-of-place, hence the input is in `in` and the output in `out`. `do_step_impl` is
    160      * used by explicit_error_stepper_base.
     228     * The result is updated out-of-place, hence the input is in `in` and the output in `out`.
     229     * Access to this step functionality is provided by explicit_stepper_base and
     230     * `do_step_impl` should not be called directly.
    161231     *
    162232     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     
    168238     * \param dt The step size.
    169239     */
    170     template< class System , class StateIn , class DerivIn , class StateOut >
    171     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
    172             time_type t , StateOut &out , time_type dt )
    173     {
    174         m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
    175 
    176         // actual calculation done in generic_rk.hpp
    177         m_rk_algorithm.do_step( stepper_base_type::m_algebra , system , in , dxdt , t , out , dt , m_x_tmp.m_v , m_F );
    178     }
    179 
    180 
    181 
    182     /**
     240
     241    /**
     242     * \fn explicit_error_generic_rk::adjust_size( const StateIn &x )
    183243     * \brief Adjust the size of all temporaries in the stepper manually.
    184244     * \param x A state from which the size of the temporaries to be resized is deduced.
    185245     */
    186     template< class StateType >
    187     void adjust_size( const StateType &x )
    188     {
    189         resize_impl( x );
    190         stepper_base_type::adjust_size( x );
    191     }
    192 
    193 
    194 private:
    195 
    196     template< class StateIn >
    197     bool resize_impl( const StateIn &x )
    198     {
    199         bool resized( false );
    200         resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() );
    201         for( size_t i = 0 ; i < StageCount-1 ; ++i )
    202         {
    203             resized |= adjust_size_by_resizeability( m_F[i] , x , typename is_resizeable<deriv_type>::type() );
    204         }
    205         return resized;
    206     }
    207 
    208 
    209     rk_algorithm_type m_rk_algorithm;
    210     coef_b_type m_b2;
    211 
    212     resizer_type m_resizer;
    213 
    214     wrapped_state_type m_x_tmp;
    215     wrapped_deriv_type m_F[StageCount-1];
    216 
    217 
    218 };
    219246
    220247}
  • trunk/boost/numeric/odeint/stepper/explicit_generic_rk.hpp

    r81491 r81754  
    7676#endif
    7777
     78
     79template<
     80size_t StageCount,
     81size_t Order,
     82class State ,
     83class Value ,
     84class Deriv ,
     85class Time ,
     86class Algebra ,
     87class Operations ,
     88class Resizer
     89>
     90#ifndef DOXYGEN_SKIP
     91class explicit_generic_rk : public explicit_stepper_base<
     92explicit_generic_rk< StageCount , Order , State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
     93Order , State , Value , Deriv , Time , Algebra , Operations , Resizer >
     94#else
     95class explicit_generic_rk : public explicit_stepper_base
     96#endif
     97{
     98
     99public:
     100
     101    #ifndef DOXYGEN_SKIP
     102    typedef explicit_stepper_base<
     103            explicit_generic_rk< StageCount , Order , State , Value , Deriv ,Time , Algebra , Operations , Resizer > ,
     104            Order , State , Value , Deriv , Time , Algebra ,
     105            Operations , Resizer > stepper_base_type;
     106    #else
     107    typedef explicit_stepper_base< ... > stepper_base_type;
     108    #endif
     109
     110    typedef typename stepper_base_type::state_type state_type;
     111    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
     112    typedef typename stepper_base_type::value_type value_type;
     113    typedef typename stepper_base_type::deriv_type deriv_type;
     114    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
     115    typedef typename stepper_base_type::time_type time_type;
     116    typedef typename stepper_base_type::algebra_type algebra_type;
     117    typedef typename stepper_base_type::operations_type operations_type;
     118    typedef typename stepper_base_type::resizer_type resizer_type;
     119
     120    #ifndef DOXYGEN_SKIP
     121    typedef explicit_generic_rk< StageCount , Order , State , Value , Deriv ,Time , Algebra , Operations , Resizer > stepper_type;
     122    #endif
     123
     124    typedef detail::generic_rk_algorithm< StageCount , Value , Algebra , Operations > rk_algorithm_type;
     125
     126    typedef typename rk_algorithm_type::coef_a_type coef_a_type;
     127    typedef typename rk_algorithm_type::coef_b_type coef_b_type;
     128    typedef typename rk_algorithm_type::coef_c_type coef_c_type;
     129
     130    #ifndef DOXYGEN_SKIP
     131    static const size_t stage_count = StageCount;
     132    #endif
     133
     134public:
     135
     136    explicit_generic_rk( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c ,
     137            const algebra_type &algebra = algebra_type() )
     138    : stepper_base_type( algebra ) , m_rk_algorithm( a , b , c )
     139    { }
     140
     141
     142    template< class System , class StateIn , class DerivIn , class StateOut >
     143    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
     144            time_type t , StateOut &out , time_type dt )
     145    {
     146        m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
     147
     148        // actual calculation done in generic_rk.hpp
     149        m_rk_algorithm.do_step( stepper_base_type::m_algebra , system , in , dxdt , t , out , dt , m_x_tmp.m_v , m_F );
     150    }
     151
     152    template< class StateIn >
     153    void adjust_size( const StateIn &x )
     154    {
     155        resize_impl( x );
     156        stepper_base_type::adjust_size( x );
     157    }
     158
     159private:
     160
     161    template< class StateIn >
     162    bool resize_impl( const StateIn &x )
     163    {
     164        bool resized( false );
     165        resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() );
     166        for( size_t i = 0 ; i < StageCount-1 ; ++i )
     167        {
     168            resized |= adjust_size_by_resizeability( m_F[i] , x , typename is_resizeable<deriv_type>::type() );
     169        }
     170        return resized;
     171    }
     172
     173
     174    rk_algorithm_type m_rk_algorithm;
     175
     176    resizer_type m_resizer;
     177
     178    wrapped_state_type m_x_tmp;
     179    wrapped_deriv_type m_F[StageCount-1];
     180
     181};
     182
     183
     184
     185/*********** DOXYGEN *************/
    78186
    79187/**
     
    97205 * \tparam Resizer The resizer policy type.
    98206 */
    99 template<
    100 size_t StageCount,
    101 size_t Order,
    102 class State ,
    103 class Value ,
    104 class Deriv ,
    105 class Time ,
    106 class Algebra ,
    107 class Operations ,
    108 class Resizer
    109 >
    110 class explicit_generic_rk : public explicit_stepper_base<
    111 explicit_generic_rk< StageCount , Order , State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
    112 Order , State , Value , Deriv , Time , Algebra , Operations , Resizer >
    113 {
    114 
    115 public:
    116 
    117     #ifndef DOXYGEN_SKIP
    118     typedef explicit_stepper_base<
    119             explicit_generic_rk< StageCount , Order , State , Value , Deriv ,Time , Algebra , Operations , Resizer > ,
    120             Order , State , Value , Deriv , Time , Algebra ,
    121             Operations , Resizer > stepper_base_type;
    122     #else
    123     typedef explicit_stepper_base< ... > stepper_base_type;
    124     #endif
    125 
    126     typedef typename stepper_base_type::state_type state_type;
    127     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
    128     typedef typename stepper_base_type::value_type value_type;
    129     typedef typename stepper_base_type::deriv_type deriv_type;
    130     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
    131     typedef typename stepper_base_type::time_type time_type;
    132     typedef typename stepper_base_type::algebra_type algebra_type;
    133     typedef typename stepper_base_type::operations_type operations_type;
    134     typedef typename stepper_base_type::resizer_type resizer_type;
    135 
    136     #ifndef DOXYGEN_SKIP
    137     typedef explicit_generic_rk< StageCount , Order , State , Value , Deriv ,Time , Algebra , Operations , Resizer > stepper_type;
    138     #endif
    139 
    140     typedef detail::generic_rk_algorithm< StageCount , Value , Algebra , Operations > rk_algorithm_type;
    141 
    142     typedef typename rk_algorithm_type::coef_a_type coef_a_type;
    143     typedef typename rk_algorithm_type::coef_b_type coef_b_type;
    144     typedef typename rk_algorithm_type::coef_c_type coef_c_type;
    145 
    146     #ifndef DOXYGEN_SKIP
    147     static const size_t stage_count = StageCount;
    148     #endif
    149 
    150 private:
    151 
    152 public:
    153207
    154208    /**
    155      * \brief Constructs the explicit_generic_rk class.
    156      * \param a The coefficients a in the butcher tableau, see the details section for an example.
    157      * \param b The coefficients b in the butcher tableau, see the details section for an example.
    158      * \param c The coefficients c in the butcher tableau, see the details section for an example.
     209     * \fn explicit_generic_rk::explicit_generic_rk( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c , const algebra_type &algebra )
     210     * \brief Constructs the explicit_generic_rk class. See examples section for details on the coefficients.
     211     * \param a Triangular matrix of parameters b in the Butcher tableau.
     212     * \param b Last row of the butcher tableau.
     213     * \param c Parameters to calculate the time points in the Butcher tableau.
    159214     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    160215     */
    161     explicit_generic_rk( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c ,
    162             const algebra_type &algebra = algebra_type() )
    163     : stepper_base_type( algebra ) , m_rk_algorithm( a , b , c )
    164     { }
    165 
    166 
     216   
    167217    /**
     218     * \fn explicit_generic_rk::do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
    168219     * \brief This method performs one step. The derivative `dxdt` of `in` at the time `t` is passed to the method.
    169      * The result is updated out of place, hence the input is in `in` and the output in `out`. `do_step_impl` is
    170      * used by explicit_stepper_base.
     220     * The result is updated out of place, hence the input is in `in` and the output in `out`.
     221     * Access to this step functionality is provided by explicit_stepper_base and
     222     * `do_step_impl` should not be called directly.
    171223     *
    172224     * \param system The system function to solve, hence the r.h.s. of the ODE. It must fullfil the
     
    178230     * \param dt The step size.
    179231     */
    180     template< class System , class StateIn , class DerivIn , class StateOut >
    181     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
    182             time_type t , StateOut &out , time_type dt )
    183     {
    184         m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
    185 
    186         // actual calculation done in generic_rk.hpp
    187         m_rk_algorithm.do_step( stepper_base_type::m_algebra , system , in , dxdt , t , out , dt , m_x_tmp.m_v , m_F );
    188     }
     232
    189233
    190234    /**
     235     * \fn explicit_generic_rk::adjust_size( const StateIn &x )
    191236     * \brief Adjust the size of all temporaries in the stepper manually.
    192237     * \param x A state from which the size of the temporaries to be resized is deduced.
    193238     */
    194     template< class StateType >
    195     void adjust_size( const StateType &x )
    196     {
    197         resize_impl( x );
    198         stepper_base_type::adjust_size( x );
    199     }
    200 
    201 private:
    202 
    203     template< class StateIn >
    204     bool resize_impl( const StateIn &x )
    205     {
    206         bool resized( false );
    207         resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() );
    208         for( size_t i = 0 ; i < StageCount-1 ; ++i )
    209         {
    210             resized |= adjust_size_by_resizeability( m_F[i] , x , typename is_resizeable<deriv_type>::type() );
    211         }
    212         return resized;
    213     }
    214 
    215 
    216     rk_algorithm_type m_rk_algorithm;
    217 
    218     resizer_type m_resizer;
    219 
    220     wrapped_state_type m_x_tmp;
    221     wrapped_deriv_type m_F[StageCount-1];
    222 
    223 };
    224239
    225240}
  • trunk/boost/numeric/odeint/stepper/modified_midpoint.hpp

    r81491 r81754  
    4141class Resizer = initially_resizer
    4242>
     43#ifndef DOXYGEN_SKIP
    4344class modified_midpoint
    4445: public explicit_stepper_base<
    4546  modified_midpoint< State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
    4647  2 , State , Value , Deriv , Time , Algebra , Operations , Resizer >
     48#else
     49class modified_midpoint : public explicit_stepper_base
     50#endif
    4751{
    4852
     
    153157 * This Stepper is for use in Bulirsch Stoer only. It DOES NOT meet any stepper concept.
    154158 */
    155 
    156159template<
    157160class State ,
     
    279282};
    280283
     284
     285
     286/********** DOXYGEN ***********/
     287
     288/**
     289 * \class modified_midpoint
     290 *
     291 * Implementation of the modified midpoint method with a configurable
     292 * number of intermediate steps. This class is used by the Bulirsch-Stoer
     293 * algorithm and is not meant for direct usage.
     294 */
     295
     296
     297/**
     298 * \class modified_midpoint_dense_out
     299 *
     300 * Implementation of the modified midpoint method with a configurable
     301 * number of intermediate steps. This class is used by the dense ouput
     302 * Bulirsch-Stoer algorithm and is not meant for direct usage.
     303 * \note This stepper is for internal use only and does not meet
     304 * any stepper concept.
     305 */
     306
     307
    281308}
    282309}
  • trunk/boost/numeric/odeint/stepper/runge_kutta4.hpp

    r81491 r81754  
    9898
    9999
     100template<
     101class State ,
     102class Value = double ,
     103class Deriv = State ,
     104class Time = Value ,
     105class Algebra = range_algebra ,
     106class Operations = default_operations ,
     107class Resizer = initially_resizer
     108>
     109#ifndef DOXYGEN_SKIP
     110class runge_kutta4 : public explicit_generic_rk< 4 , 4 , State , Value , Deriv , Time ,
     111Algebra , Operations , Resizer >
     112#else
     113class runge_kutta4 : public explicit_generic_rk
     114#endif
     115{
     116
     117public:
     118
     119#ifndef DOXYGEN_SKIP
     120    typedef explicit_generic_rk< 4 , 4 , State , Value , Deriv , Time ,
     121            Algebra , Operations , Resizer > stepper_base_type;
     122#endif
     123    typedef typename stepper_base_type::state_type state_type;
     124    typedef typename stepper_base_type::value_type value_type;
     125    typedef typename stepper_base_type::deriv_type deriv_type;
     126    typedef typename stepper_base_type::time_type time_type;
     127    typedef typename stepper_base_type::algebra_type algebra_type;
     128    typedef typename stepper_base_type::operations_type operations_type;
     129    typedef typename stepper_base_type::resizer_type resizer_type;
     130
     131    #ifndef DOXYGEN_SKIP
     132    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
     133    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
     134    typedef typename stepper_base_type::stepper_type stepper_type;
     135    #endif
     136
     137    runge_kutta4( const algebra_type &algebra = algebra_type() ) : stepper_base_type(
     138            boost::fusion::make_vector( rk4_coefficients_a1<Value>() , rk4_coefficients_a2<Value>() , rk4_coefficients_a3<Value>() ) ,
     139            rk4_coefficients_b<Value>() , rk4_coefficients_c<Value>() , algebra )
     140    { }
     141
     142};
     143
    100144/**
    101145 * \class runge_kutta4
     
    120164 * \tparam Resizer The resizer policy type.
    121165 */
    122 template<
    123 class State ,
    124 class Value = double ,
    125 class Deriv = State ,
    126 class Time = Value ,
    127 class Algebra = range_algebra ,
    128 class Operations = default_operations ,
    129 class Resizer = initially_resizer
    130 >
    131 class runge_kutta4 : public explicit_generic_rk< 4 , 4 , State , Value , Deriv , Time ,
    132 Algebra , Operations , Resizer >
    133 {
    134166
    135 public:
    136 
    137     typedef explicit_generic_rk< 4 , 4 , State , Value , Deriv , Time ,
    138             Algebra , Operations , Resizer > stepper_base_type;
    139     typedef typename stepper_base_type::state_type state_type;
    140     typedef typename stepper_base_type::value_type value_type;
    141     typedef typename stepper_base_type::deriv_type deriv_type;
    142     typedef typename stepper_base_type::time_type time_type;
    143     typedef typename stepper_base_type::algebra_type algebra_type;
    144     typedef typename stepper_base_type::operations_type operations_type;
    145     typedef typename stepper_base_type::resizer_type resizer_type;
    146 
    147     #ifndef DOXYGEN_SKIP
    148     typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
    149     typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
    150     typedef typename stepper_base_type::stepper_type stepper_type;
    151     #endif
    152 
    153 
    154     /**
    155      * \brief Constructs the runge_kutta4 class. This constructor can be used as a default
    156      * constructor if the algebra has a default constructor.
    157      * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
    158      */
    159     runge_kutta4( const algebra_type &algebra = algebra_type() ) : stepper_base_type(
    160             boost::fusion::make_vector( rk4_coefficients_a1<Value>() , rk4_coefficients_a2<Value>() , rk4_coefficients_a3<Value>() ) ,
    161             rk4_coefficients_b<Value>() , rk4_coefficients_c<Value>() , algebra )
    162     { }
    163 
    164 };
     167/**
     168 * \fn runge_kutta4::runge_kutta4( const algebra_type &algebra = algebra_type() )
     169 * \brief Constructs the runge_kutta4 class. This constructor can be used as a default
     170 * constructor if the algebra has a default constructor.
     171 * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
     172 */
    165173
    166174}
  • trunk/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp

    </
    r81491 r81754  
    3232namespace numeric {
    3333namespace odeint {
     34
     35template<
     36class State ,
     37class Value = double ,
     38class Deriv = State ,
     39class Time = Value ,
     40class Algebra = range_algebra ,
     41class Operations = default_operations ,
     42class Resizer = initially_resizer
     43>
     44#ifndef DOXYGEN_SKIP
     45class runge_kutta4_classic
     46: public explicit_stepper_base<
     47  runge_kutta4_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
     48  4 , State , Value , Deriv , Time , Algebra , Operations , Resizer >
     49#else
     50class runge_kutta4_classic : public explicit_stepper_base
     51#endif
     52{
     53
     54public :
     55
     56    #ifndef DOXYGEN_SKIP
     57    typedef explicit_stepper_base<
     58    runge_kutta4_classic< State , Value , Deriv , Time , Algebra , Operations , Resizer > ,
     59    4 , State , Value , Deriv , Time , Algebra , Operations , Resizer > stepper_base_type;
     60    #else
     61    typedef explicit_stepper_base< runge_kutta4_classic< ... > , ... > stepper_base_type;
     62    #endif
     63
     64    typedef typename stepper_base_type::state_type state_type;
     65    typedef typename stepper_base_type::value_type value_type;
     66    typedef typename stepper_base_type::deriv_type deriv_type;
     67    typedef typename stepper_base_type::time_type time_type;
     68    typedef typename stepper_base_type::algebra_type algebra_type;
     69    typedef typename stepper_base_type::operations_type operations_type;
     70    typedef typename stepper_base_type::resizer_type resizer_type;
     71
     72    #ifndef DOXYGEN_SKIP
     73    typedef typename stepper_base_type::stepper_type stepper_type;
     74    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
     75    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
     76    #endif // DOXYGEN_SKIP
     77
     78
     79
     80    runge_kutta4_classic( const algebra_type &algebra = algebra_type() ) : stepper_base_type( algebra )
     81    { }
     82
     83
     84    template< class System , class StateIn , class DerivIn , class StateOut >
     85    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
     86    {
     87        // ToDo : check if size of in,dxdt,out are equal?
     88
     89        static const value_type val1 = static_cast< value_type >( 1 );
     90
     91        m_resizer.adjust_size( in , detail::bind( &stepper_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
     92
     93        typename odeint::unwrap_reference< System >::type &sys = system;
     94
     95        const time_type dh = dt / static_cast< value_type >( 2 );
     96        const time_type th = t + dh;
     97
     98        // dt * dxdt = k1
     99        // m_x_tmp = x + dh*dxdt
     100        stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , dxdt ,
     101                typename operations_type::template scale_sum2< value_type , time_type >( val1 , dh ) );
     102
     103
     104        // dt * m_dxt = k2
     105        sys( m_x_tmp.m_v , m_dxt.m_v , th );
     106
     107        // m_x_tmp = x + dh*m_dxt
     108        stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , m_dxt.m_v ,
     109                typename operations_type::template scale_sum2< value_type , time_type >( val1 , dh ) );
     110
     111
     112        // dt * m_dxm = k3
     113        sys( m_x_tmp.m_v , m_dxm.m_v , th );
     114        //m_x_tmp = x + dt*m_dxm
     115        stepper_base_type::m_algebra.for_each3( m_x_tmp.m_v , in , m_dxm.m_v ,
     116                typename operations_type::template scale_sum2< value_type , time_type >( val1 , dt ) );
     117
     118
     119        // dt * m_dxh = k4
     120        sys( m_x_tmp.m_v , m_dxh.m_v , t + dt );
     121        //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
     122        time_type dt6 = dt / static_cast< value_type >( 6 );
     123        time_type dt3 = dt / static_cast< value_type >( 3 );
     124        stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v ,
     125                typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
     126    }
     127
     128    template< class StateType >
     129    void adjust_size( const StateType &x )
     130    {
     131        resize_impl( x );
     132        stepper_base_type::adjust_size( x );
     133    }
     134
     135private:
     136
     137    template< class StateIn >
     138    bool resize_impl( const StateIn &x )
     139    {
     140        bool resized = false;
     141        resized |= adjust_size_by_resizeability( m_x_tmp , x , typename is_resizeable<state_type>::type() );
     142        resized |= adjust_size_by_resizeability( m_dxm , x , typename is_resizeable<deriv_type>::type() );
     143        resized |= adjust_size_by_resizeability( m_dxt , x , typename is_resizeable<deriv_type>::type() );
     144        resized |= adjust_size_by_resizeability( m_dxh , x , typename is_resizeable<deriv_type>::type() );
     145        return resized;
     146    }
     147
     148
     149    resizer_type m_resizer;
     150
     151    wrapped_deriv_type m_dxt;
     152    wrapped_deriv_type m_dxm;
     153    wrapped_deriv_type m_dxh;
     154    wrapped_state_type m_x_tmp;
     155
     156};
     157
     158
     159/********* DOXYGEN *********/
    34160
    35161/**
     
    56182 * \tparam Resizer The resizer policy type.
    57183 */
    58 template<
    59 class State ,
    60 class Value = double ,
    61 class Deriv = State ,
    62 class Time = Value ,
    63 class Algebra = range_algebra ,