Subject: [Boost-commit] svn:boost r64786 - in sandbox/odeint/branches/karsten: . boost/numeric/odeint/stepper libs/numeric/odeint/test
From: mario.mulansky_at_[hidden]
Date: 2010-08-13 10:35:12
Author: mariomulansky
Date: 2010-08-13 10:35:09 EDT (Fri, 13 Aug 2010)
New Revision: 64786
URL: http://svn.boost.org/trac/boost/changeset/64786
Log:
added implicit euler (not working, yet)
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/.project | 1 +
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 5 +++++
2 files changed, 6 insertions(+), 0 deletions(-)
Modified: sandbox/odeint/branches/karsten/.project
==============================================================================
--- sandbox/odeint/branches/karsten/.project (original)
+++ sandbox/odeint/branches/karsten/.project 2010-08-13 10:35:09 EDT (Fri, 13 Aug 2010)
@@ -3,6 +3,7 @@
<name>karsten</name>
<comment></comment>
<projects>
+ <project>boost</project>
</projects>
<buildSpec>
<buildCommand>
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp 2010-08-13 10:35:09 EDT (Fri, 13 Aug 2010)
@@ -0,0 +1,106 @@
+/*
+ boost header: BOOST_NUMERIC_ODEINT/implicit_euler.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 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_BOOST_NUMERIC_ODEINT_IMPLICIT_EULER_HPP_INCLUDED
+#define BOOST_BOOST_NUMERIC_ODEINT_IMPLICIT_EULER_HPP_INCLUDED
+
+#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 implicit_euler
+{
+
+public:
+
+ typedef ValueType value_type;
+ typedef boost::numeric::ublas::vector< value_type > state_type;
+ typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+ typedef boost::numeric::ublas::permutation_matrix< std::size_t > pmatrix_type;
+
+ implicit_euler( const value_type epsilon = 1E-6 ) : m_epsilon( epsilon )
+ { }
+
+ template< class System , class Jacobi >
+ void do_step( System &system , Jacobi &jacobi , state_type &x , value_type t , value_type dt )
+ {
+ m_dxdt.resize( x.size() );
+ m_x.resize( x.size() );
+ m_b.resize( x.size() );
+ m_jacobi.resize( x.size() , x.size() );
+ m_pm.resize( x.size() );
+
+ t += dt;
+
+ // apply first Newton step
+ system( x , m_dxdt , t );
+
+ m_b = dt * m_dxdt;
+
+ jacobi( x , m_jacobi , t );
+ m_jacobi *= dt;
+ m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
+
+ solve( m_b , m_jacobi );
+
+ m_x = x - m_b;
+
+ // iterate Newton until some precision is reached
+ while( boost::numeric::ublas::norm_2( m_b ) > m_epsilon )
+ {
+ system( m_x , m_dxdt , t );
+ m_b = x - m_x - dt*m_dxdt;
+
+ jacobi( m_x , m_jacobi , t );
+ m_jacobi *= dt;
+ m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
+
+ solve( m_b , m_jacobi );
+
+ m_x -= m_b;
+ }
+
+ }
+
+private:
+
+ void solve( state_type &x , matrix_type &m )
+ {
+ int res = boost::numeric::ublas::lu_factorize( m , m_pm );
+ if( res != 0 ) exit(0);
+
+ boost::numeric::ublas::lu_substitute( m , m_pm , x );
+ }
+
+private:
+ state_type m_dxdt;
+ state_type m_x;
+ state_type m_b;
+ matrix_type m_jacobi;
+ pmatrix_type m_pm;
+
+ const value_type m_epsilon;
+
+};
+
+
+} // odeint
+} // numeric
+} // boost
+
+
+#endif //BOOST_BOOST_NUMERIC_ODEINT_IMPLICIT_EULER_HPP_INCLUDED
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile 2010-08-13 10:35:09 EDT (Fri, 13 Aug 2010)
@@ -40,6 +40,11 @@
:
: <link>static
]
+ [ run check_implicit_euler.cpp
+ :
+ :
+ : <link>static
+ ]
;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp 2010-08-13 10:35:09 EDT (Fri, 13 Aug 2010)
@@ -0,0 +1,57 @@
+/* Boost check_implicit_euler.cpp test file
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ This file tests the use of the euler stepper
+
+ 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)
+*/
+
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/implicit_euler.hpp>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef double value_type;
+typedef boost::numeric::ublas::vector< value_type > state_type;
+typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+
+
+void system( state_type &x , state_type &dxdt , const value_type t )
+{
+ dxdt( 0 ) = x( 0 ) + 2 * x( 1 );
+ dxdt( 1 ) = x( 1 );
+}
+
+void jacobi( state_type &x , matrix_type &jacobi , const value_type t )
+{
+ jacobi( 0 , 0 ) = 1;
+ jacobi( 0 , 1 ) = 2;
+ jacobi( 1 , 0 ) = 0;
+ jacobi( 1 , 1 ) = 1;
+}
+
+
+BOOST_AUTO_TEST_CASE( test_euler )
+{
+ implicit_euler< value_type > stepper;
+ state_type x( 2 );
+ x(0) = 0.0; x(1) = 1.0;
+
+ value_type eps = 1E-12;
+
+ stepper.do_step( system , jacobi , x , 0.0 , 0.1 );
+
+ BOOST_CHECK_MESSAGE( abs( x(0) - 0.1 ) < eps , x[0] - 0.1 );
+
+}