Subject: [Boost-commit] svn:boost r64836 - in sandbox/odeint/branches/karsten: . boost/numeric/odeint/algebra boost/numeric/odeint/stepper libs/numeric/odeint/examples libs/numeric/odeint/test
From: mario.mulansky_at_[hidden]
Date: 2010-08-15 16:54:25
Author: mariomulansky
Date: 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
New Revision: 64836
URL: http://svn.boost.org/trac/boost/changeset/64836
Log:
+stiff example
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/stiff_system.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/.project | 1 +
sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_resize.hpp | 12 +++++++++++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp | 34 +++++++++++++++++++++++++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile | 3 ++-
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 2 +-
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp | 7 ++++---
6 files changed, 44 insertions(+), 15 deletions(-)
Modified: sandbox/odeint/branches/karsten/.project
==============================================================================
--- sandbox/odeint/branches/karsten/.project (original)
+++ sandbox/odeint/branches/karsten/.project 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -4,6 +4,7 @@
<comment></comment>
<projects>
<project>boost</project>
+ <project>boost_1_41_0</project>
</projects>
<buildSpec>
<buildCommand>
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_resize.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_resize.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_resize.hpp 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -15,6 +15,8 @@
#include <vector>
#include <list>
+/* ToDo: boost ublas dependency here? */
+#include <boost/numeric/ublas/vector.hpp>
namespace boost {
namespace numeric {
@@ -53,7 +55,15 @@
const static bool value = type::value;
};
-
+/*
+ * specialization for boost::numeric::ublas::vector
+ */
+template< class T >
+struct is_resizeable< boost::numeric::ublas::vector< T > >
+{
+ struct type : public boost::true_type { };
+ const static bool value = type::value;
+};
template< class Container >
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -16,6 +16,9 @@
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
+#include <iostream>
+#include <boost/numeric/ublas/io.hpp>
+
namespace boost {
namespace numeric {
namespace odeint {
@@ -30,9 +33,9 @@
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;
+ typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
- implicit_euler( const value_type epsilon = 1E-6 ) : m_epsilon( epsilon )
+ implicit_euler( const value_type epsilon = 1E-6 ) : m_epsilon( epsilon ) , m_pm( 1 )
{ }
template< class System , class Jacobi >
@@ -42,7 +45,7 @@
m_x.resize( x.size() );
m_b.resize( x.size() );
m_jacobi.resize( x.size() , x.size() );
- m_pm.resize( x.size() );
+ m_pm = pmatrix_type( x.size() ); // no resize because we also need default filling
t += dt;
@@ -55,25 +58,38 @@
m_jacobi *= dt;
m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
- solve( m_b , m_jacobi );
+ std::clog << m_jacobi << std::endl;
+ std::clog << m_b << std::endl;
+
+ matrix_type jacobi_tmp( m_jacobi );
+
+ solve( m_b , jacobi_tmp );
+
+ std::clog << m_b << std::endl;
m_x = x - m_b;
// 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_epsilon )
{
system( m_x , m_dxdt , t );
- m_b = x - m_x - dt*m_dxdt;
+ m_b = x - m_x + dt*m_dxdt;
+ /* we use simplified newton where the jacobi matrix is evaluated only once
jacobi( m_x , m_jacobi , t );
m_jacobi *= dt;
m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
+ */
+ jacobi_tmp = m_jacobi;
+
+ solve( m_b , jacobi_tmp );
- solve( m_b , m_jacobi );
+ std::clog << m_b << std::endl;
m_x -= m_b;
}
-
+ x = m_x;
}
private:
@@ -87,14 +103,14 @@
}
private:
+ const value_type m_epsilon;
+
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;
-
};
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -7,7 +7,7 @@
project
: requirements
<include>../../../..
- <include>$BOOST_ROOT
+ <include>$(BOOST_ROOT)
<define>BOOST_ALL_NO_LIB=1
: build-dir .
;
@@ -15,3 +15,4 @@
# exe harmonic_oscillator : harmonic_oscillator.cpp ;
# exe solar_system : solar_system.cpp point_type.hpp ;
+exe stiff_system : stiff_system.cpp ;
\ No newline at end of file
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/stiff_system.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/stiff_system.cpp 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -0,0 +1,65 @@
+/* Boost implicit_euler.cpp example file
+
+ Copyright 2010 Karsten Ahnert
+ Copyright 2010 Mario Mulansky
+
+ This file shows the use of the implicit 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 <iostream>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/implicit_euler.hpp>
+#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+
+#define tab '\t'
+
+typedef double value_type;
+typedef boost::numeric::ublas::vector< value_type > state_type;
+typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+
+void stiff_system( state_type &x , state_type &dxdt , value_type t )
+{
+ dxdt( 0 ) = -101.0 * x( 0 ) - 100.0 * x( 1 );
+ dxdt( 1 ) = x( 0 );
+}
+
+void jacobi( state_type &x , matrix_type &jacobi , value_type t )
+{
+ jacobi( 0 , 0 ) = -101.0;
+ jacobi( 0 , 1 ) = -100.0;
+ jacobi( 1 , 0 ) = 1.0;
+ jacobi( 1 , 1 ) = 0.0;
+}
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+int main( void )
+{
+ explicit_euler< state_type , value_type /*, vector_space_algebra */> expl_euler;
+ implicit_euler< value_type > impl_euler;
+
+ state_type x1( 2 );
+ x1( 0 ) = 1.0; x1( 1 ) = 0.0;
+ state_type x2( x1 );
+
+ const double dt = 0.01; // for dt >= 0.01 the euler method gets unstable
+ const size_t steps = 1000;
+
+ for( size_t step = 0 ; step < steps ; ++step )
+ {
+ clog << step << " of " << 100 << endl;
+ cout << step*dt << tab << x1( 0 ) << tab << x1( 1 ) << tab << x2( 0 ) << tab << x2( 1 ) << endl;
+ expl_euler.do_step( stiff_system , x1 , step*dt , dt );
+ impl_euler.do_step( stiff_system , jacobi , x2 , step*dt , dt );
+ }
+ cout << steps*dt << tab << x1( 0 ) << tab << x1( 1 ) << tab << x2( 0 ) << tab << x2( 1 ) << endl;
+}
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-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -21,7 +21,7 @@
# <link>static
# <link>shared:<define>BOOST_TEST_DYN_LINK=1
<include>../../../..
- <include>$BOOST_ROOT
+ <include>$(BOOST_ROOT)
;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_implicit_euler.cpp 2010-08-15 16:54:23 EDT (Sun, 15 Aug 2010)
@@ -10,6 +10,7 @@
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
+#define BOOST_TEST_MODULE test_implicit_euler
#include <boost/test/unit_test.hpp>
@@ -27,7 +28,7 @@
typedef boost::numeric::ublas::matrix< value_type > matrix_type;
-void system( state_type &x , state_type &dxdt , const value_type t )
+void sys( state_type &x , state_type &dxdt , const value_type t )
{
dxdt( 0 ) = x( 0 ) + 2 * x( 1 );
dxdt( 1 ) = x( 1 );
@@ -38,7 +39,7 @@
jacobi( 0 , 0 ) = 1;
jacobi( 0 , 1 ) = 2;
jacobi( 1 , 0 ) = 0;
- jacobi( 1 , 1 ) = 1;
+ jacobi( 1 , 1 ) = 3;
}
@@ -50,7 +51,7 @@
value_type eps = 1E-12;
- stepper.do_step( system , jacobi , x , 0.0 , 0.1 );
+ stepper.do_step( sys , jacobi , x , 0.0 , 0.1 );
BOOST_CHECK_MESSAGE( abs( x(0) - 0.1 ) < eps , x[0] - 0.1 );