1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
/*
[begin_description]
Modification of the implicit Euler method, works with the MTL4 matrix library only.
[end_description]
Copyright 2012-2013 Andreas Angelopoulos
Copyright 2012-2013 Karsten Ahnert
Copyright 2012-2013 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
#include <utility>
#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
namespace boost {
namespace numeric {
namespace odeint {
template< class ValueType , class Resizer = initially_resizer >
class implicit_euler_mtl4
{
public:
typedef ValueType value_type;
typedef value_type time_type;
typedef mtl::dense_vector<value_type> state_type;
typedef state_wrapper< state_type > wrapped_state_type;
typedef state_type deriv_type;
typedef state_wrapper< deriv_type > wrapped_deriv_type;
typedef mtl::compressed2D< value_type > matrix_type;
typedef state_wrapper< matrix_type > wrapped_matrix_type;
typedef Resizer resizer_type;
typedef stepper_tag stepper_category;
typedef implicit_euler_mtl4< ValueType , Resizer > stepper_type;
implicit_euler_mtl4( const value_type epsilon = 1E-6 )
: m_epsilon( epsilon ) , m_resizer() ,
m_dxdt() , m_x() ,
m_identity() , m_jacobi()
{ }
template< class System >
void do_step( System system , state_type &x , time_type t , time_type dt )
{
typedef typename odeint::unwrap_reference< System >::type system_type;
typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
system_type &sys = system;
deriv_func_type &deriv_func = sys.first;
jacobi_func_type &jacobi_func = sys.second;
m_resizer.adjust_size( x , detail::bind(
&stepper_type::template resize_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
m_identity.m_v = 1;
t += dt;
m_x.m_v = x;
deriv_func( x , m_dxdt.m_v , t );
jacobi_func( x , m_jacobi.m_v , t );
m_dxdt.m_v *= -dt;
m_jacobi.m_v *= dt;
m_jacobi.m_v -= m_identity.m_v ;
// using ilu_0 preconditioning -incomplete LU factorisation
// itl::pc::diagonal<matrix_type,double> L(m_jacobi.m_v);
itl::pc::ilu_0<matrix_type> L( m_jacobi.m_v );
solve( m_jacobi.m_v , m_x.m_v , m_dxdt.m_v , L );
x+= m_x.m_v;
}
template< class StateType >
void adjust_size( const StateType &x )
{
resize_impl( x );
}
private:
/*
Applying approximate iterative linear solvers
default solver is Biconjugate gradient stabilized method
itl::bicgstab(A, x, b, L, iter);
*/
template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner>
void solve(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& L, int max_iteractions =500)
{
// Termination criterion: r < 1e-6 * b or N iterations
itl::basic_iteration< double > iter( b , max_iteractions , 1e-6 );
itl::bicgstab( A , x , b , L , iter );
}
template< class StateIn >
bool resize_impl( const StateIn &x )
{
bool resized = false;
resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
resized |= adjust_size_by_resizeability( m_identity , x , typename is_resizeable<matrix_type>::type() );
resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
return resized;
}
private:
value_type m_epsilon;
resizer_type m_resizer;
wrapped_deriv_type m_dxdt;
wrapped_state_type m_x;
wrapped_matrix_type m_identity;
wrapped_matrix_type m_jacobi;
};
} // odeint
} // numeric
} // boost
#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
|
