]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 ) 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 |