Chris@16: /*
Chris@16:  [auto_generated]
Chris@16:  boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp
Chris@16: 
Chris@16:  [begin_description]
Chris@16:  Default Integrate adaptive implementation.
Chris@16:  [end_description]
Chris@16: 
Chris@16:  Copyright 2009-2011 Karsten Ahnert
Chris@16:  Copyright 2009-2011 Mario Mulansky
Chris@16: 
Chris@16:  Distributed under the Boost Software License, Version 1.0.
Chris@16:  (See accompanying file LICENSE_1_0.txt or
Chris@16:  copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  */
Chris@16: 
Chris@16: 
Chris@16: #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
Chris@16: #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED
Chris@16: 
Chris@16: #include <stdexcept>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
Chris@16: #include <boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp>
Chris@16: #include <boost/numeric/odeint/util/bind.hpp>
Chris@16: #include <boost/numeric/odeint/util/unwrap_reference.hpp>
Chris@16: #include <boost/numeric/odeint/util/copy.hpp>
Chris@16: 
Chris@16: #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
Chris@16: 
Chris@16: 
Chris@16: #include <iostream>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace numeric {
Chris@16: namespace odeint {
Chris@16: namespace detail {
Chris@16: 
Chris@16: // forward declaration
Chris@16: template< class Stepper , class System , class State , class Time , class Observer>
Chris@16: Time integrate_n_steps(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time start_time , Time dt , size_t num_of_steps ,
Chris@16:         Observer observer , stepper_tag );
Chris@16: 
Chris@16: /*
Chris@16:  * integrate_adaptive for simple stepper is basically an integrate_const + some last step
Chris@16:  */
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_adaptive(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time start_time , Time end_time , Time dt ,
Chris@16:         Observer observer , stepper_tag
Chris@16: )
Chris@16: {
Chris@16:     size_t steps = static_cast< size_t >( (end_time-start_time)/dt );
Chris@16:     Time end = detail::integrate_n_steps( stepper , system , start_state , start_time ,
Chris@16:                                           dt , steps , observer , stepper_tag() );
Chris@16:     if( less_with_sign( end , end_time , dt ) )
Chris@16:     {   //make a last step to end exactly at end_time
Chris@16:         stepper.do_step( system , start_state , end , end_time - end );
Chris@16:         steps++;
Chris@16:         typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@16:         obs( start_state , end_time );
Chris@16:     }
Chris@16:     return steps;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: /*
Chris@16:  * classical integrate adaptive
Chris@16:  */
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_adaptive(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time &start_time , Time end_time , Time &dt ,
Chris@16:         Observer observer , controlled_stepper_tag
Chris@16: )
Chris@16: {
Chris@16:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@16: 
Chris@16:     const size_t max_attempts = 1000;
Chris@16:     const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
Chris@16:     size_t count = 0;
Chris@16:     while( less_with_sign( start_time , end_time , dt ) )
Chris@16:     {
Chris@16:         obs( start_state , start_time );
Chris@16:         if( less_with_sign( end_time , start_time + dt , dt ) )
Chris@16:         {
Chris@16:             dt = end_time - start_time;
Chris@16:         }
Chris@16: 
Chris@16:         size_t trials = 0;
Chris@16:         controlled_step_result res = success;
Chris@16:         do
Chris@16:         {
Chris@16:             res = stepper.try_step( system , start_state , start_time , dt );
Chris@16:             ++trials;
Chris@16:         }
Chris@16:         while( ( res == fail ) && ( trials < max_attempts ) );
Chris@16:         if( trials == max_attempts ) throw std::overflow_error( error_string );
Chris@16: 
Chris@16:         ++count;
Chris@16:     }
Chris@16:     obs( start_state , start_time );
Chris@16:     return count;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: /*
Chris@16:  * integrate adaptive for dense output steppers
Chris@16:  *
Chris@16:  * step size control is used if the stepper supports it
Chris@16:  */
Chris@16: template< class Stepper , class System , class State , class Time , class Observer >
Chris@16: size_t integrate_adaptive(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@16:         Time start_time , Time end_time , Time dt ,
Chris@16:         Observer observer , dense_output_stepper_tag )
Chris@16: {
Chris@16:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@16: 
Chris@16:     size_t count = 0;
Chris@16:     stepper.initialize( start_state , start_time , dt );
Chris@16: 
Chris@16:     while( less_with_sign( stepper.current_time() , end_time , stepper.current_time_step() ) )
Chris@16:     {
Chris@16:         while( less_eq_with_sign( stepper.current_time() + stepper.current_time_step() ,
Chris@16:                end_time ,
Chris@16:                stepper.current_time_step() ) )
Chris@16:         {   //make sure we don't go beyond the end_time
Chris@16:             obs( stepper.current_state() , stepper.current_time() );
Chris@16:             stepper.do_step( system );
Chris@16:             ++count;
Chris@16:         }
Chris@16:         stepper.initialize( stepper.current_state() , stepper.current_time() , end_time - stepper.current_time() );
Chris@16:     }
Chris@16:     obs( stepper.current_state() , stepper.current_time() );
Chris@16:     // overwrite start_state with the final point
Chris@16:     boost::numeric::odeint::copy( stepper.current_state() , start_state );
Chris@16:     return count;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: } // namespace detail
Chris@16: } // namespace odeint
Chris@16: } // namespace numeric
Chris@16: } // namespace boost
Chris@16: 
Chris@16: 
Chris@16: #endif // BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_ADAPTIVE_HPP_INCLUDED