integrate/detail/integrate_n_steps.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp
0004
0005  [begin_description]
0006  integrate steps implementation
0007  [end_description]
0008
0009  Copyright 2012-2015 Mario Mulansky
0010  Copyright 2012 Christoph Koke
0011  Copyright 2012 Karsten Ahnert
0012
0013  Distributed under the Boost Software License, Version 1.0.
0014  (See accompanying file LICENSE_1_0.txt or
0015  copy at http://www.boost.org/LICENSE_1_0.txt)
0016  */
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_N_STEPS_HPP_INCLUDED
0019 #define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_N_STEPS_HPP_INCLUDED
0020
0021 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0022 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0023 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
0024 #include <boost/numeric/odeint/util/unit_helper.hpp>
0025
0026 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
0027
0028 namespace boost {
0029 namespace numeric {
0030 namespace odeint {
0031 namespace detail {
0032
0033 // forward declaration
0034 template< class Stepper , class System , class State , class Time , class Observer >
0035 size_t integrate_adaptive_checked(
0036         Stepper stepper , System system , State &start_state ,
0037         Time &start_time , Time end_time , Time &dt ,
0038         Observer observer, controlled_stepper_tag
0039 );
0040
0041
0042 /* basic version */
0043 template< class Stepper , class System , class State , class Time , class Observer>
0044 Time integrate_n_steps(
0045         Stepper stepper , System system , State &start_state ,
0046         Time start_time , Time dt , size_t num_of_steps ,
0047         Observer observer , stepper_tag )
0048 {
0049     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0050     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
0051
0052     Time time = start_time;
0053
0054     for( size_t step = 0; step < num_of_steps ; ++step )
0055     {
0056         obs( start_state , time );
0057         st.do_step( system , start_state , time , dt );
0058         // direct computation of the time avoids error propagation happening when using time += dt
0059         // we need clumsy type analysis to get boost units working here
0060         time = start_time + static_cast< typename unit_value_type<Time>::type >( step+1 ) * dt;
0061     }
0062     obs( start_state , time );
0063
0064     return time;
0065 }
0066
0067
0068 /* controlled version */
0069 template< class Stepper , class System , class State , class Time , class Observer >
0070 Time integrate_n_steps(
0071         Stepper stepper , System system , State &start_state ,
0072         Time start_time , Time dt , size_t num_of_steps ,
0073         Observer observer , controlled_stepper_tag )
0074 {
0075     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0076
0077     Time time = start_time;
0078     Time time_step = dt;
0079
0080     for( size_t step = 0; step < num_of_steps ; ++step )
0081     {
0082         obs( start_state , time );
0083         // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
0084         detail::integrate_adaptive(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt,
0085                                    null_observer(), controlled_stepper_tag());
0086         // direct computation of the time avoids error propagation happening when using time += dt
0087         // we need clumsy type analysis to get boost units working here
0088         time = start_time + static_cast< typename unit_value_type<Time>::type >(step+1) * time_step;
0089     }
0090     obs( start_state , time );
0091
0092     return time;
0093 }
0094
0095
0096 /* dense output version */
0097 template< class Stepper , class System , class State , class Time , class Observer >
0098 Time integrate_n_steps(
0099         Stepper stepper , System system , State &start_state ,
0100         Time start_time , Time dt , size_t num_of_steps ,
0101         Observer observer , dense_output_stepper_tag )
0102 {
0103     typename odeint::unwrap_reference< Observer >::type &obs = observer;
0104     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
0105
0106     Time time = start_time;
0107     const Time end_time = start_time + static_cast< typename unit_value_type<Time>::type >(num_of_steps) * dt;
0108
0109     st.initialize( start_state , time , dt );
0110
0111     size_t step = 0;
0112
0113     while( step < num_of_steps )
0114     {
0115         while( less_with_sign( time , st.current_time() , st.current_time_step() ) )
0116         {
0117             st.calc_state( time , start_state );
0118             obs( start_state , time );
0119             ++step;
0120             // direct computation of the time avoids error propagation happening when using time += dt
0121             // we need clumsy type analysis to get boost units working here
0122             time = start_time + static_cast< typename unit_value_type<Time>::type >(step) * dt;
0123         }
0124
0125         // we have not reached the end, do another real step
0126         if( less_with_sign( static_cast<Time>(st.current_time()+st.current_time_step()) ,
0127                             end_time ,
0128                             st.current_time_step() ) )
0129         {
0130             st.do_step( system );
0131         }
0132         else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
0133         { // do the last step ending exactly on the end point
0134             st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
0135             st.do_step( system );
0136         }
0137     }
0138
0139     // make sure we really end exactly where we should end
0140     while( st.current_time() < end_time )
0141     {
0142         if( less_with_sign( end_time ,
0143                             static_cast<Time>(st.current_time()+st.current_time_step()) ,
0144                             st.current_time_step() ) )
0145             st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
0146         st.do_step( system );
0147     }
0148
0149     // observation at final point
0150     obs( st.current_state() , end_time );
0151
0152     return time;
0153 }
0154
0155
0156 }
0157 }
0158 }
0159 }
0160
0161 #endif /* BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_N_STEPS_HPP_INCLUDED */