163 lines
4.7 KiB
Plaintext
163 lines
4.7 KiB
Plaintext
/*
|
|
[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
|