boost/numeric/odeint/stepper/implicit_euler.hpp

   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
  68     implicit_euler( value_type epsilon = 1E-6 )
  69     : m_epsilon( epsilon )
  70     { }
  71
  72
  73     template< class System >
  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 >
 123     void adjust_size( const StateType &x )
 124     {
 125         resize_impl( x );
 126     }
 127
 128
 129 private:
 130
 131     template< class StateIn >
 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
 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 ) 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