odeint/stepper/implicit_euler.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/implicit_euler.hpp
0004
0005  [begin_description]
0006  Impementation of the implicit Euler method. Works with ublas::vector as state type.
0007  [end_description]
0008
0009  Copyright 2010-2012 Mario Mulansky
0010  Copyright 2010-2012 Karsten Ahnert
0011  Copyright 2012 Christoph Koke
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
0019 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED
0021
0022
0023 #include <utility>
0024
0025 #include <boost/numeric/odeint/util/bind.hpp>
0026 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0027 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0028
0029 #include <boost/numeric/odeint/util/ublas_wrapper.hpp>
0030 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0031 #include <boost/numeric/odeint/util/resizer.hpp>
0032
0033 #include <boost/numeric/ublas/vector.hpp>
0034 #include <boost/numeric/ublas/matrix.hpp>
0035 #include <boost/numeric/ublas/lu.hpp>
0036
0037 namespace boost {
0038 namespace numeric {
0039 namespace odeint {
0040
0041
0042
0043
0044
0045
0046
0047
0048 template< class ValueType , class Resizer = initially_resizer >
0049 class implicit_euler
0050 {
0051
0052 public:
0053
0054     typedef ValueType value_type;
0055     typedef value_type time_type;
0056     typedef boost::numeric::ublas::vector< value_type > state_type;
0057     typedef state_wrapper< state_type > wrapped_state_type;
0058     typedef state_type deriv_type;
0059     typedef state_wrapper< deriv_type > wrapped_deriv_type;
0060     typedef boost::numeric::ublas::matrix< value_type > matrix_type;
0061     typedef state_wrapper< matrix_type > wrapped_matrix_type;
0062     typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
0063     typedef state_wrapper< pmatrix_type > wrapped_pmatrix_type;
0064     typedef Resizer resizer_type;
0065     typedef stepper_tag stepper_category;
0066     typedef implicit_euler< ValueType , Resizer > stepper_type;
0067
0068     implicit_euler( value_type epsilon = 1E-6 )
0069     : m_epsilon( epsilon )
0070     { }
0071
0072
0073     template< class System >
0074     void do_step( System system , state_type &x , time_type t , time_type dt )
0075     {
0076         typedef typename odeint::unwrap_reference< System >::type system_type;
0077         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
0078         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
0079         system_type &sys = system;
0080         deriv_func_type &deriv_func = sys.first;
0081         jacobi_func_type &jacobi_func = sys.second;
0082
0083         m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<state_type> , detail::ref( *this ) , detail::_1 ) );
0084
0085         for( size_t i=0 ; i<x.size() ; ++i )
0086             m_pm.m_v[i] = i;
0087
0088         t += dt;
0089
0090         // apply first Newton step
0091         deriv_func( x , m_dxdt.m_v , t );
0092
0093         m_b.m_v = dt * m_dxdt.m_v;
0094
0095         jacobi_func( x , m_jacobi.m_v  , t );
0096         m_jacobi.m_v *= dt;
0097         m_jacobi.m_v -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
0098
0099         solve( m_b.m_v , m_jacobi.m_v );
0100
0101         m_x.m_v = x - m_b.m_v;
0102
0103         // iterate Newton until some precision is reached
0104         // ToDo: maybe we should apply only one Newton step -> linear implicit one-step scheme
0105         while( boost::numeric::ublas::norm_2( m_b.m_v ) > m_epsilon )
0106         {
0107             deriv_func( m_x.m_v , m_dxdt.m_v , t );
0108             m_b.m_v = x - m_x.m_v + dt*m_dxdt.m_v;
0109
0110             // simplified version, only the first Jacobian is used
0111             //            jacobi( m_x , m_jacobi , t );
0112             //            m_jacobi *= dt;
0113             //            m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
0114
0115             solve( m_b.m_v , m_jacobi.m_v );
0116
0117             m_x.m_v -= m_b.m_v;
0118         }
0119         x = m_x.m_v;
0120     }
0121
0122     template< class StateType >
0123     void adjust_size( const StateType &x )
0124     {
0125         resize_impl( x );
0126     }
0127
0128
0129 private:
0130
0131     template< class StateIn >
0132     bool resize_impl( const StateIn &x )
0133     {
0134         bool resized = false;
0135         resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
0136         resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
0137         resized |= adjust_size_by_resizeability( m_b , x , typename is_resizeable<deriv_type>::type() );
0138         resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
0139         resized |= adjust_size_by_resizeability( m_pm , x , typename is_resizeable<pmatrix_type>::type() );
0140         return resized;
0141     }
0142
0143
0144     void solve( state_type &x , matrix_type &m )
0145     {
0146         int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v );
0147         if( res != 0 ) std::exit(0);
0148         boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x );
0149     }
0150
0151 private:
0152
0153     value_type m_epsilon;
0154     resizer_type m_resizer;
0155     wrapped_deriv_type m_dxdt;
0156     wrapped_state_type m_x;
0157     wrapped_deriv_type m_b;
0158     wrapped_matrix_type m_jacobi;
0159     wrapped_pmatrix_type m_pm;
0160
0161
0162 };
0163
0164
0165 } // odeint
0166 } // numeric
0167 } // boost
0168
0169
0170 #endif // BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED