1 /* 2 [auto_generated] 3 boost/numeric/odeint/stepper/implicit_euler.hpp 4 5 [begin_description] 6 Impementation of the implicit Euler method. Works with ublas::vector as state type. 7 [end_description] 8 9 Copyright 2010-2012 Mario Mulansky 10 Copyright 2010-2012 Karsten Ahnert 11 Copyright 2012 Christoph Koke 12 13 Distributed under the Boost Software License, Version 1.0. 14 (See accompanying file LICENSE_1_0.txt or 15 copy at http://www.boost.org/LICENSE_1_0.txt) 16 */ 17 18 19 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED 20 #define BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED 21 22 23 #include <utility> 24 25 #include <boost/numeric/odeint/util/bind.hpp> 26 #include <boost/numeric/odeint/util/unwrap_reference.hpp> 27 #include <boost/numeric/odeint/stepper/stepper_categories.hpp> 28 29 #include <boost/numeric/odeint/util/ublas_wrapper.hpp> 30 #include <boost/numeric/odeint/util/is_resizeable.hpp> 31 #include <boost/numeric/odeint/util/resizer.hpp> 32 33 #include <boost/numeric/ublas/vector.hpp> 34 #include <boost/numeric/ublas/matrix.hpp> 35 #include <boost/numeric/ublas/lu.hpp> 36 37 namespace boost { 38 namespace numeric { 39 namespace odeint { 40 41 42 43 44 45 46 47 48 template< class ValueType , class Resizer = initially_resizer > 49 class implicit_euler 50 { 51 52 public: 53 54 typedef ValueType value_type; 55 typedef value_type time_type; 56 typedef boost::numeric::ublas::vector< value_type > state_type; 57 typedef state_wrapper< state_type > wrapped_state_type; 58 typedef state_type deriv_type; 59 typedef state_wrapper< deriv_type > wrapped_deriv_type; 60 typedef boost::numeric::ublas::matrix< value_type > matrix_type; 61 typedef state_wrapper< matrix_type > wrapped_matrix_type; 62 typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type; 63 typedef state_wrapper< pmatrix_type > wrapped_pmatrix_type; 64 typedef Resizer resizer_type; 65 typedef stepper_tag stepper_category; 66 typedef implicit_euler< ValueType , Resizer > stepper_type; 67 implicit_euler(value_type epsilon=1E-6)68 implicit_euler( value_type epsilon = 1E-6 ) 69 : m_epsilon( epsilon ) 70 { } 71 72 73 template< class System > do_step(System system,state_type & x,time_type t,time_type dt)74 void do_step( System system , state_type &x , time_type t , time_type dt ) 75 { 76 typedef typename odeint::unwrap_reference< System >::type system_type; 77 typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type; 78 typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type; 79 system_type &sys = system; 80 deriv_func_type &deriv_func = sys.first; 81 jacobi_func_type &jacobi_func = sys.second; 82 83 m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<state_type> , detail::ref( *this ) , detail::_1 ) ); 84 85 for( size_t i=0 ; i<x.size() ; ++i ) 86 m_pm.m_v[i] = i; 87 88 t += dt; 89 90 // apply first Newton step 91 deriv_func( x , m_dxdt.m_v , t ); 92 93 m_b.m_v = dt * m_dxdt.m_v; 94 95 jacobi_func( x , m_jacobi.m_v , t ); 96 m_jacobi.m_v *= dt; 97 m_jacobi.m_v -= boost::numeric::ublas::identity_matrix< value_type >( x.size() ); 98 99 solve( m_b.m_v , m_jacobi.m_v ); 100 101 m_x.m_v = x - m_b.m_v; 102 103 // iterate Newton until some precision is reached 104 // ToDo: maybe we should apply only one Newton step -> linear implicit one-step scheme 105 while( boost::numeric::ublas::norm_2( m_b.m_v ) > m_epsilon ) 106 { 107 deriv_func( m_x.m_v , m_dxdt.m_v , t ); 108 m_b.m_v = x - m_x.m_v + dt*m_dxdt.m_v; 109 110 // simplified version, only the first Jacobian is used 111 // jacobi( m_x , m_jacobi , t ); 112 // m_jacobi *= dt; 113 // m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() ); 114 115 solve( m_b.m_v , m_jacobi.m_v ); 116 117 m_x.m_v -= m_b.m_v; 118 } 119 x = m_x.m_v; 120 } 121 122 template< class StateType > adjust_size(const StateType & x)123 void adjust_size( const StateType &x ) 124 { 125 resize_impl( x ); 126 } 127 128 129 private: 130 131 template< class StateIn > resize_impl(const StateIn & x)132 bool resize_impl( const StateIn &x ) 133 { 134 bool resized = false; 135 resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() ); 136 resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() ); 137 resized |= adjust_size_by_resizeability( m_b , x , typename is_resizeable<deriv_type>::type() ); 138 resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() ); 139 resized |= adjust_size_by_resizeability( m_pm , x , typename is_resizeable<pmatrix_type>::type() ); 140 return resized; 141 } 142 143 solve(state_type & x,matrix_type & m)144 void solve( state_type &x , matrix_type &m ) 145 { 146 int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v ); 147 if( res != 0 ) std::exit(0); 148 boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x ); 149 } 150 151 private: 152 153 value_type m_epsilon; 154 resizer_type m_resizer; 155 wrapped_deriv_type m_dxdt; 156 wrapped_state_type m_x; 157 wrapped_deriv_type m_b; 158 wrapped_matrix_type m_jacobi; 159 wrapped_pmatrix_type m_pm; 160 161 162 }; 163 164 165 } // odeint 166 } // numeric 167 } // boost 168 169 170 #endif // BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED 171