FM / stepper

/usr/include/boost/numeric/odeint/stepper/

/usr/include/boost/numeric/odeint/stepper

Name	Size	Mode	Actions
base/	-	0755	rm
detail/	-	0755	rm
generation/	-	0755	rm
adams_bashforth.hpp	15264	0644	edit dl rm
adams_bashforth_moulton.hpp	12015	0644	edit dl rm
adams_moulton.hpp	6112	0644	edit dl rm
adaptive_adams_bashforth_moulton.hpp	8311	0644	edit dl rm
bulirsch_stoer.hpp	26189	0644	edit dl rm
bulirsch_stoer_dense_out.hpp	33859	0644	edit dl rm
controlled_adams_bashforth_moulton.hpp	10080	0644	edit dl rm
controlled_runge_kutta.hpp	38301	0644	edit dl rm
controlled_step_result.hpp	1092	0644	edit dl rm
dense_output_runge_kutta.hpp	14590	0644	edit dl rm
euler.hpp	6587	0644	edit dl rm
explicit_error_generic_rk.hpp	10708	0644	edit dl rm
explicit_generic_rk.hpp	8491	0644	edit dl rm
extrapolation_stepper.hpp	11155	0644	edit dl rm
generation.hpp	1477	0644	edit dl rm
implicit_euler.hpp	5236	0644	edit dl rm
modified_midpoint.hpp	10566	0644	edit dl rm
rosenbrock4.hpp	14602	0644	edit dl rm
rosenbrock4_controller.hpp	7290	0644	edit dl rm
rosenbrock4_dense_output.hpp	5595	0644	edit dl rm
runge_kutta4.hpp	5860	0644	edit dl rm
runge_kutta4_classic.hpp	9255	0644	edit dl rm
runge_kutta_cash_karp54.hpp	8676	0644	edit dl rm
runge_kutta_cash_karp54_classic.hpp	14318	0644	edit dl rm
runge_kutta_dopri5.hpp	21484	0644	edit dl rm
runge_kutta_fehlberg78.hpp	15005	0644	edit dl rm
stepper_categories.hpp	2162	0644	edit dl rm
symplectic_euler.hpp	3937	0644	edit dl rm
symplectic_rkn_sb3a_m4_mclachlan.hpp	5120	0644	edit dl rm
symplectic_rkn_sb3a_mclachlan.hpp	5389	0644	edit dl rm
velocity_verlet.hpp	13544	0644	edit dl rm

Edit: /usr/include/boost/numeric/odeint/stepper/implicit_euler.hpp (5236B)

/*
 [auto_generated]
 boost/numeric/odeint/stepper/implicit_euler.hpp

[begin_description]
 Impementation of the implicit Euler method. Works with ublas::vector as state type.
 [end_description]

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_STEPPER_IMPLICIT_EULER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_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/util/ublas_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>

namespace boost {
namespace numeric {
namespace odeint {

template< class ValueType , class Resizer = initially_resizer >
class implicit_euler
{

public:

typedef ValueType value_type;
    typedef value_type time_type;
    typedef boost::numeric::ublas::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 boost::numeric::ublas::matrix< value_type > matrix_type;
    typedef state_wrapper< matrix_type > wrapped_matrix_type;
    typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
    typedef state_wrapper< pmatrix_type > wrapped_pmatrix_type;
    typedef Resizer resizer_type;
    typedef stepper_tag stepper_category;
    typedef implicit_euler< ValueType , Resizer > stepper_type;

implicit_euler( value_type epsilon = 1E-6 )
    : m_epsilon( epsilon ) 
    { }

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 ) );

for( size_t i=0 ; i<x.size() ; ++i )
            m_pm.m_v[i] = i;

t += dt;

// apply first Newton step
        deriv_func( x , m_dxdt.m_v , t );

m_b.m_v = dt * m_dxdt.m_v;

jacobi_func( x , m_jacobi.m_v  , t );
        m_jacobi.m_v *= dt;
        m_jacobi.m_v -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );

solve( m_b.m_v , m_jacobi.m_v );

m_x.m_v = x - m_b.m_v;

// iterate Newton until some precision is reached
        // ToDo: maybe we should apply only one Newton step -> linear implicit one-step scheme
        while( boost::numeric::ublas::norm_2( m_b.m_v ) > m_epsilon )
        {
            deriv_func( m_x.m_v , m_dxdt.m_v , t );
            m_b.m_v = x - m_x.m_v + dt*m_dxdt.m_v;

// simplified version, only the first Jacobian is used
            //            jacobi( m_x , m_jacobi , t );
            //            m_jacobi *= dt;
            //            m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );

solve( m_b.m_v , m_jacobi.m_v );

m_x.m_v -= m_b.m_v;
        }
        x = m_x.m_v;
    }

template< class StateType >
    void adjust_size( const StateType &x )
    {
        resize_impl( x );
    }

private:

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_b , x , typename is_resizeable<deriv_type>::type() );
        resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
        resized |= adjust_size_by_resizeability( m_pm , x , typename is_resizeable<pmatrix_type>::type() );
        return resized;
    }

void solve( state_type &x , matrix_type &m )
    {
        int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v );
        if( res != 0 ) std::exit(0);
        boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x );
    }

private:

value_type m_epsilon;
    resizer_type m_resizer;
    wrapped_deriv_type m_dxdt;
    wrapped_state_type m_x;
    wrapped_deriv_type m_b;
    wrapped_matrix_type m_jacobi;
    wrapped_pmatrix_type m_pm;

};

} // odeint
} // numeric
} // boost

#endif // BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED

cmd: