[Boost-commit] svn:boost r64786 - in sandbox/odeint/branches/karsten: . boost/numeric/odeint/stepper libs/numeric/odeint/test

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