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@101:  Copyright 2011-2013 Karsten Ahnert
Chris@101:  Copyright 2011-2012 Mario Mulansky
Chris@101:  Copyright 2012 Christoph Koke
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@101: #include <boost/throw_exception.hpp>
Chris@101: 
Chris@16: #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
Chris@16: #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
Chris@101: #include <boost/numeric/odeint/integrate/detail/integrate_const.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@101: size_t integrate_const(
Chris@16:         Stepper stepper , System system , State &start_state ,
Chris@101:         Time start_time , Time end_time , Time dt ,
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@101:     size_t steps = detail::integrate_const( stepper , system , start_state , start_time ,
Chris@101:                                             end_time , dt , observer , stepper_tag() );
Chris@101:     typename odeint::unwrap_reference< Observer >::type &obs = observer;
Chris@101:     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
Chris@101: 
Chris@101:     Time end = start_time + dt*steps;
Chris@16:     if( less_with_sign( end , end_time , dt ) )
Chris@16:     {   //make a last step to end exactly at end_time
Chris@101:         st.do_step( system , start_state , end , end_time - end );
Chris@16:         steps++;
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@101:     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
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@101:         if( less_with_sign( end_time , static_cast<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@101:             res = st.try_step( system , start_state , start_time , dt );
Chris@16:             ++trials;
Chris@16:         }
Chris@16:         while( ( res == fail ) && ( trials < max_attempts ) );
Chris@101:         if( trials == max_attempts ) BOOST_THROW_EXCEPTION( 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@101:     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
Chris@16: 
Chris@16:     size_t count = 0;
Chris@101:     st.initialize( start_state , start_time , dt );
Chris@16: 
Chris@101:     while( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
Chris@16:     {
Chris@101:         while( less_eq_with_sign( static_cast<Time>(st.current_time() + st.current_time_step()) ,
Chris@16:                end_time ,
Chris@101:                st.current_time_step() ) )
Chris@16:         {   //make sure we don't go beyond the end_time
Chris@101:             obs( st.current_state() , st.current_time() );
Chris@101:             st.do_step( system );
Chris@16:             ++count;
Chris@16:         }
Chris@101:         // calculate time step to arrive exactly at end time
Chris@101:         st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
Chris@16:     }
Chris@101:     obs( st.current_state() , st.current_time() );
Chris@16:     // overwrite start_state with the final point
Chris@101:     boost::numeric::odeint::copy( st.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