ceph/src/boost/libs/numeric/odeint/include/boost/numeric/odeint/stepper/rosenbrock4.hpp

   1 /*
   2  [auto_generated]
   3  boost/numeric/odeint/stepper/rosenbrock4.hpp
   4
   5  [begin_description]
   6  Implementation of the Rosenbrock 4 method for solving stiff ODEs. Note, that a
   7  controller and a dense-output stepper exist for this method,
   8  [end_description]
   9
  10  Copyright 2011-2013 Karsten Ahnert
  11  Copyright 2011-2012 Mario Mulansky
  12  Copyright 2012 Christoph Koke
  13
  14  Distributed under the Boost Software License, Version 1.0.
  15  (See accompanying file LICENSE_1_0.txt or
  16  copy at http://www.boost.org/LICENSE_1_0.txt)
  17  */
  18
  19
  20 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_HPP_INCLUDED
  21 #define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_HPP_INCLUDED
  22
  23
  24 #include <boost/numeric/odeint/util/bind.hpp>
  25 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  26 #include <boost/numeric/ublas/vector.hpp>
  27 #include <boost/numeric/ublas/matrix.hpp>
  28 #include <boost/numeric/ublas/lu.hpp>
  29
  30 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  31
  32 #include <boost/numeric/odeint/util/ublas_wrapper.hpp>
  33 #include <boost/numeric/odeint/util/is_resizeable.hpp>
  34 #include <boost/numeric/odeint/util/resizer.hpp>
  35
  36 #include <boost/numeric/ublas/vector.hpp>
  37 #include <boost/numeric/ublas/matrix.hpp>
  38 #include <boost/numeric/ublas/lu.hpp>
  39
  40
  41 namespace boost {
  42 namespace numeric {
  43 namespace odeint {
  44
  45
  46 /*
  47  * ToDo:
  48  *
  49  * 2. Interfacing for odeint, check if controlled_error_stepper can be used
  50  * 3. dense output
  51  */
  52
  53
  54
  55 template< class Value >
  56 struct default_rosenbrock_coefficients
  57 {
  58     typedef Value value_type;
  59     typedef unsigned short order_type;
  60
  61     default_rosenbrock_coefficients( void )
  62     : gamma ( static_cast< value_type >( 0.25 ) ) ,
  63       d1 ( static_cast< value_type >( 0.25 ) ) ,
  64       d2 ( static_cast< value_type >( -0.1043 ) ) ,
  65       d3 ( static_cast< value_type >( 0.1035 ) ) ,
  66       d4 ( static_cast< value_type >( 0.3620000000000023e-01 ) ) ,
  67       c2 ( static_cast< value_type >( 0.386 ) ) ,
  68       c3 ( static_cast< value_type >( 0.21 ) ) ,
  69       c4 ( static_cast< value_type >( 0.63 ) ) ,
  70       c21 ( static_cast< value_type >( -0.5668800000000000e+01 ) ) ,
  71       a21 ( static_cast< value_type >( 0.1544000000000000e+01 ) ) ,
  72       c31 ( static_cast< value_type >( -0.2430093356833875e+01 ) ) ,
  73       c32 ( static_cast< value_type >( -0.2063599157091915e+00 ) ) ,
  74       a31 ( static_cast< value_type >( 0.9466785280815826e+00 ) ) ,
  75       a32 ( static_cast< value_type >( 0.2557011698983284e+00 ) ) ,
  76       c41 ( static_cast< value_type >( -0.1073529058151375e+00 ) ) ,
  77       c42 ( static_cast< value_type >( -0.9594562251023355e+01 ) ) ,
  78       c43 ( static_cast< value_type >( -0.2047028614809616e+02 ) ) ,
  79       a41 ( static_cast< value_type >( 0.3314825187068521e+01 ) ) ,
  80       a42 ( static_cast< value_type >( 0.2896124015972201e+01 ) ) ,
  81       a43 ( static_cast< value_type >( 0.9986419139977817e+00 ) ) ,
  82       c51 ( static_cast< value_type >( 0.7496443313967647e+01 ) ) ,
  83       c52 ( static_cast< value_type >( -0.1024680431464352e+02 ) ) ,
  84       c53 ( static_cast< value_type >( -0.3399990352819905e+02 ) ) ,
  85       c54 ( static_cast< value_type >(  0.1170890893206160e+02 ) ) ,
  86       a51 ( static_cast< value_type >( 0.1221224509226641e+01 ) ) ,
  87       a52 ( static_cast< value_type >( 0.6019134481288629e+01 ) ) ,
  88       a53 ( static_cast< value_type >( 0.1253708332932087e+02 ) ) ,
  89       a54 ( static_cast< value_type >( -0.6878860361058950e+00 ) ) ,
  90       c61 ( static_cast< value_type >( 0.8083246795921522e+01 ) ) ,
  91       c62 ( static_cast< value_type >( -0.7981132988064893e+01 ) ) ,
  92       c63 ( static_cast< value_type >( -0.3152159432874371e+02 ) ) ,
  93       c64 ( static_cast< value_type >( 0.1631930543123136e+02 ) ) ,
  94       c65 ( static_cast< value_type >( -0.6058818238834054e+01 ) ) ,
  95       d21 ( static_cast< value_type >( 0.1012623508344586e+02 ) ) ,
  96       d22 ( static_cast< value_type >( -0.7487995877610167e+01 ) ) ,
  97       d23 ( static_cast< value_type >( -0.3480091861555747e+02 ) ) ,
  98       d24 ( static_cast< value_type >( -0.7992771707568823e+01 ) ) ,
  99       d25 ( static_cast< value_type >( 0.1025137723295662e+01 ) ) ,
 100       d31 ( static_cast< value_type >( -0.6762803392801253e+00 ) ) ,
 101       d32 ( static_cast< value_type >( 0.6087714651680015e+01 ) ) ,
 102       d33 ( static_cast< value_type >( 0.1643084320892478e+02 ) ) ,
 103       d34 ( static_cast< value_type >( 0.2476722511418386e+02 ) ) ,
 104       d35 ( static_cast< value_type >( -0.6594389125716872e+01 ) )
 105     {}
 106
 107     const value_type gamma;
 108     const value_type d1 , d2 , d3 , d4;
 109     const value_type c2 , c3 , c4;
 110     const value_type c21 ;
 111     const value_type a21;
 112     const value_type c31 , c32;
 113     const value_type a31 , a32;
 114     const value_type c41 , c42 , c43;
 115     const value_type a41 , a42 , a43;
 116     const value_type c51 , c52 , c53 , c54;
 117     const value_type a51 , a52 , a53 , a54;
 118     const value_type c61 , c62 , c63 , c64 , c65;
 119     const value_type d21 , d22 , d23 , d24 , d25;
 120     const value_type d31 , d32 , d33 , d34 , d35;
 121
 122     static const order_type stepper_order = 4;
 123     static const order_type error_order = 3;
 124 };
 125
 126
 127
 128 template< class Value , class Coefficients = default_rosenbrock_coefficients< Value > , class Resizer = initially_resizer >
 129 class rosenbrock4
 130 {
 131 private:
 132
 133 public:
 134
 135     typedef Value value_type;
 136     typedef boost::numeric::ublas::vector< value_type > state_type;
 137     typedef state_type deriv_type;
 138     typedef value_type time_type;
 139     typedef boost::numeric::ublas::matrix< value_type > matrix_type;
 140     typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
 141     typedef Resizer resizer_type;
 142     typedef Coefficients rosenbrock_coefficients;
 143     typedef stepper_tag stepper_category;
 144     typedef unsigned short order_type;
 145
 146     typedef state_wrapper< state_type > wrapped_state_type;
 147     typedef state_wrapper< deriv_type > wrapped_deriv_type;
 148     typedef state_wrapper< matrix_type > wrapped_matrix_type;
 149     typedef state_wrapper< pmatrix_type > wrapped_pmatrix_type;
 150
 151     typedef rosenbrock4< Value , Coefficients , Resizer > stepper_type;
 152
 153     const static order_type stepper_order = rosenbrock_coefficients::stepper_order;
 154     const static order_type error_order = rosenbrock_coefficients::error_order;
 155
 156     rosenbrock4( void )
 157     : m_resizer() , m_x_err_resizer() ,
 158       m_jac() , m_pm() ,
 159       m_dfdt() , m_dxdt() , m_dxdtnew() ,
 160       m_g1() , m_g2() , m_g3() , m_g4() , m_g5() ,
 161       m_cont3() , m_cont4() , m_xtmp() , m_x_err() ,
 162       m_coef()
 163     { }
 164
 165
 166     order_type order() const { return stepper_order; }
 167
 168     template< class System >
 169     void do_step( System system , const state_type &x , time_type t , state_type &xout , time_type dt , state_type &xerr )
 170     {
 171         // get the system and jacobi function
 172         typedef typename odeint::unwrap_reference< System >::type system_type;
 173         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
 174         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
 175         system_type &sys = system;
 176         deriv_func_type &deriv_func = sys.first;
 177         jacobi_func_type &jacobi_func = sys.second;
 178
 179         const size_t n = x.size();
 180
 181         m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<state_type> , detail::ref( *this ) , detail::_1 ) );
 182
 183         for( size_t i=0 ; i<n ; ++i )
 184             m_pm.m_v( i ) = i;
 185
 186         deriv_func( x , m_dxdt.m_v , t );
 187         jacobi_func( x , m_jac.m_v , t , m_dfdt.m_v );
 188
 189         m_jac.m_v *= -1.0;
 190         m_jac.m_v += 1.0 / m_coef.gamma / dt * boost::numeric::ublas::identity_matrix< value_type >( n );
 191         boost::numeric::ublas::lu_factorize( m_jac.m_v , m_pm.m_v );
 192
 193         for( size_t i=0 ; i<n ; ++i )
 194             m_g1.m_v[i] = m_dxdt.m_v[i] + dt * m_coef.d1 * m_dfdt.m_v[i];
 195         boost::numeric::ublas::lu_substitute( m_jac.m_v , m_pm.m_v , m_g1.m_v );
 196
 197
 198         for( size_t i=0 ; i<n ; ++i )
 199             m_xtmp.m_v[i] = x[i] + m_coef.a21 * m_g1.m_v[i];
 200         deriv_func( m_xtmp.m_v , m_dxdtnew.m_v , t + m_coef.c2 * dt );
 201         for( size_t i=0 ; i<n ; ++i )
 202             m_g2.m_v[i] = m_dxdtnew.m_v[i] + dt * m_coef.d2 * m_dfdt.m_v[i] + m_coef.c21 * m_g1.m_v[i] / dt;
 203         boost::numeric::ublas::lu_substitute( m_jac.m_v , m_pm.m_v , m_g2.m_v );
 204
 205
 206         for( size_t i=0 ; i<n ; ++i )
 207             m_xtmp.m_v[i] = x[i] + m_coef.a31 * m_g1.m_v[i] + m_coef.a32 * m_g2.m_v[i];
 208         deriv_func( m_xtmp.m_v , m_dxdtnew.m_v , t + m_coef.c3 * dt );
 209         for( size_t i=0 ; i<n ; ++i )
 210             m_g3.m_v[i] = m_dxdtnew.m_v[i] + dt * m_coef.d3 * m_dfdt.m_v[i] + ( m_coef.c31 * m_g1.m_v[i] + m_coef.c32 * m_g2.m_v[i] ) / dt;
 211         boost::numeric::ublas::lu_substitute( m_jac.m_v , m_pm.m_v , m_g3.m_v );
 212
 213
 214         for( size_t i=0 ; i<n ; ++i )
 215             m_xtmp.m_v[i] = x[i] + m_coef.a41 * m_g1.m_v[i] + m_coef.a42 * m_g2.m_v[i] + m_coef.a43 * m_g3.m_v[i];
 216         deriv_func( m_xtmp.m_v , m_dxdtnew.m_v , t + m_coef.c4 * dt );
 217         for( size_t i=0 ; i<n ; ++i )
 218             m_g4.m_v[i] = m_dxdtnew.m_v[i] + dt * m_coef.d4 * m_dfdt.m_v[i] + ( m_coef.c41 * m_g1.m_v[i] + m_coef.c42 * m_g2.m_v[i] + m_coef.c43 * m_g3.m_v[i] ) / dt;
 219         boost::numeric::ublas::lu_substitute( m_jac.m_v , m_pm.m_v , m_g4.m_v );
 220
 221
 222         for( size_t i=0 ; i<n ; ++i )
 223             m_xtmp.m_v[i] = x[i] + m_coef.a51 * m_g1.m_v[i] + m_coef.a52 * m_g2.m_v[i] + m_coef.a53 * m_g3.m_v[i] + m_coef.a54 * m_g4.m_v[i];
 224         deriv_func( m_xtmp.m_v , m_dxdtnew.m_v , t + dt );
 225         for( size_t i=0 ; i<n ; ++i )
 226             m_g5.m_v[i] = m_dxdtnew.m_v[i] + ( m_coef.c51 * m_g1.m_v[i] + m_coef.c52 * m_g2.m_v[i] + m_coef.c53 * m_g3.m_v[i] + m_coef.c54 * m_g4.m_v[i] ) / dt;
 227         boost::numeric::ublas::lu_substitute( m_jac.m_v , m_pm.m_v , m_g5.m_v );
 228
 229         for( size_t i=0 ; i<n ; ++i )
 230             m_xtmp.m_v[i] += m_g5.m_v[i];
 231         deriv_func( m_xtmp.m_v , m_dxdtnew.m_v , t + dt );
 232         for( size_t i=0 ; i<n ; ++i )
 233             xerr[i] = m_dxdtnew.m_v[i] + ( m_coef.c61 * m_g1.m_v[i] + m_coef.c62 * m_g2.m_v[i] + m_coef.c63 * m_g3.m_v[i] + m_coef.c64 * m_g4.m_v[i] + m_coef.c65 * m_g5.m_v[i] ) / dt;
 234         boost::numeric::ublas::lu_substitute( m_jac.m_v , m_pm.m_v , xerr );
 235
 236         for( size_t i=0 ; i<n ; ++i )
 237             xout[i] = m_xtmp.m_v[i] + xerr[i];
 238     }
 239
 240     template< class System >
 241     void do_step( System system , state_type &x , time_type t , time_type dt , state_type &xerr )
 242     {
 243         do_step( system , x , t , x , dt , xerr );
 244     }
 245
 246     /*
 247      * do_step without error output - just calls above functions with and neglects the error estimate
 248      */
 249     template< class System >
 250     void do_step( System system , const state_type &x , time_type t , state_type &xout , time_type dt )
 251     {
 252         m_x_err_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_x_err<state_type> , detail::ref( *this ) , detail::_1 ) );
 253         do_step( system , x , t , xout , dt , m_x_err.m_v );
 254     }
 255
 256     template< class System >
 257     void do_step( System system , state_type &x , time_type t , time_type dt )
 258     {
 259         m_x_err_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_x_err<state_type> , detail::ref( *this ) , detail::_1 ) );
 260         do_step( system , x , t , dt , m_x_err.m_v );
 261     }
 262
 263     void prepare_dense_output()
 264     {
 265         const size_t n = m_g1.m_v.size();
 266         for( size_t i=0 ; i<n ; ++i )
 267         {
 268             m_cont3.m_v[i] = m_coef.d21 * m_g1.m_v[i] + m_coef.d22 * m_g2.m_v[i] + m_coef.d23 * m_g3.m_v[i] + m_coef.d24 * m_g4.m_v[i] + m_coef.d25 * m_g5.m_v[i];
 269             m_cont4.m_v[i] = m_coef.d31 * m_g1.m_v[i] + m_coef.d32 * m_g2.m_v[i] + m_coef.d33 * m_g3.m_v[i] + m_coef.d34 * m_g4.m_v[i] + m_coef.d35 * m_g5.m_v[i];
 270         }
 271     }
 272
 273
 274     void calc_state( time_type t , state_type &x ,
 275             const state_type &x_old , time_type t_old ,
 276             const state_type &x_new , time_type t_new )
 277     {
 278         const size_t n = m_g1.m_v.size();
 279         time_type dt = t_new - t_old;
 280         time_type s = ( t - t_old ) / dt;
 281         time_type s1 = 1.0 - s;
 282         for( size_t i=0 ; i<n ; ++i )
 283             x[i] = x_old[i] * s1 + s * ( x_new[i] + s1 * ( m_cont3.m_v[i] + s * m_cont4.m_v[i] ) );
 284     }
 285
 286
 287
 288     template< class StateType >
 289     void adjust_size( const StateType &x )
 290     {
 291         resize_impl( x );
 292         resize_x_err( x );
 293     }
 294
 295
 296 protected:
 297
 298     template< class StateIn >
 299     bool resize_impl( const StateIn &x )
 300     {
 301         bool resized = false;
 302         resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
 303         resized |= adjust_size_by_resizeability( m_dfdt , x , typename is_resizeable<deriv_type>::type() );
 304         resized |= adjust_size_by_resizeability( m_dxdtnew , x , typename is_resizeable<deriv_type>::type() );
 305         resized |= adjust_size_by_resizeability( m_xtmp , x , typename is_resizeable<state_type>::type() );
 306         resized |= adjust_size_by_resizeability( m_g1 , x , typename is_resizeable<state_type>::type() );
 307         resized |= adjust_size_by_resizeability( m_g2 , x , typename is_resizeable<state_type>::type() );
 308         resized |= adjust_size_by_resizeability( m_g3 , x , typename is_resizeable<state_type>::type() );
 309         resized |= adjust_size_by_resizeability( m_g4 , x , typename is_resizeable<state_type>::type() );
 310         resized |= adjust_size_by_resizeability( m_g5 , x , typename is_resizeable<state_type>::type() );
 311         resized |= adjust_size_by_resizeability( m_cont3 , x , typename is_resizeable<state_type>::type() );
 312         resized |= adjust_size_by_resizeability( m_cont4 , x , typename is_resizeable<state_type>::type() );
 313         resized |= adjust_size_by_resizeability( m_jac , x , typename is_resizeable<matrix_type>::type() );
 314         resized |= adjust_size_by_resizeability( m_pm , x , typename is_resizeable<pmatrix_type>::type() );
 315         return resized;
 316     }
 317
 318     template< class StateIn >
 319     bool resize_x_err( const StateIn &x )
 320     {
 321         return adjust_size_by_resizeability( m_x_err , x , typename is_resizeable<state_type>::type() );
 322     }
 323
 324 private:
 325
 326
 327     resizer_type m_resizer;
 328     resizer_type m_x_err_resizer;
 329
 330     wrapped_matrix_type m_jac;
 331     wrapped_pmatrix_type m_pm;
 332     wrapped_deriv_type m_dfdt , m_dxdt , m_dxdtnew;
 333     wrapped_state_type m_g1 , m_g2 , m_g3 , m_g4 , m_g5;
 334     wrapped_state_type m_cont3 , m_cont4;
 335     wrapped_state_type m_xtmp;
 336     wrapped_state_type m_x_err;
 337
 338     const rosenbrock_coefficients m_coef;
 339 };
 340
 341
 342 } // namespace odeint
 343 } // namespace numeric
 344 } // namespace boost
 345
 346 #endif // BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_HPP_INCLUDED