|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r63928 - in sandbox/odeint/branches/karsten: boost/numeric boost/numeric/odeint/algebra boost/numeric/odeint/stepper libs/numeric/odeint/test
From: mario.mulansky_at_[hidden]
Date: 2010-07-12 11:54:17
Author: mariomulansky
Date: 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
New Revision: 63928
URL: http://svn.boost.org/trac/boost/changeset/63928
Log:
first implementation of error stepper rk_ck
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/runge_kutta_error_ck.hpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp | 1
sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_algebra.hpp | 142 ++++++++++++++++++++++++++++++++++++++-
sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_operations.hpp | 57 ++++++++++++++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_euler.hpp | 2
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp | 20 +++++
5 files changed, 215 insertions(+), 7 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
@@ -18,5 +18,6 @@
#include <boost/config.hpp>
#include <boost/numeric/odeint/stepper/explicit_euler.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_error_ck.hpp>
#endif // BOOST_NUMERIC_ODEINT_HPP
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_algebra.hpp 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
@@ -25,7 +25,7 @@
typedef Container container_type;
template< class StateType1 , class StateType2 , class Operation >
- static void transform2( StateType1 &s1 , StateType2 &s2 , Operation op )
+ static void for_each2( StateType1 &s1 , StateType2 &s2 , Operation op )
{
// ToDo : check that number of arguments of the operation is equal 2
@@ -34,13 +34,13 @@
BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType2 >::type , container_type >::value ));
// ToDo : pack into detail namespace
- transform2( boost::begin( s1 ) , boost::end( s1 ) ,
+ for_each2( boost::begin( s1 ) , boost::end( s1 ) ,
boost::begin( s2 ) , op );
}
// ToDo : pack into namespace detail
template< class Iterator1 , class Iterator2 , class Operation >
- static void transform2( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Operation op )
+ static void for_each2( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Operation op )
{
for( ; first1 != last1 ; )
op( *first1++ , *first2++ );
@@ -48,7 +48,7 @@
template< class StateType1 , class StateType2 , class StateType3 , class Operation >
- static void transform3( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , Operation op )
+ static void for_each3( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , Operation op )
{
// ToDo : check that number of arguments of the operation is equal 3
@@ -58,7 +58,7 @@
BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType3 >::type , container_type >::value ));
// ToDo : pack into detail namespace
- transform2( boost::begin( s1 ) , boost::end( s1 ) ,
+ for_each3( boost::begin( s1 ) , boost::end( s1 ) ,
boost::begin( s2 ) ,
boost::begin( s3 ) ,
op );
@@ -66,12 +66,142 @@
// ToDo : pack into namespace detail
template< class Iterator1 , class Iterator2 , class Iterator3 , class Operation >
- static void transform3( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3, Operation op )
+ static void for_each3( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3, Operation op )
{
for( ; first1 != last1 ; )
op( *first1++ , *first2++ , *first3++ );
}
+
+
+ template< class StateType1 , class StateType2 , class StateType3 , class StateType4 , class Operation >
+ static void for_each4( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , StateType4 &s4 , Operation op )
+ {
+ // ToDo : check that number of arguments of the operation is equal 3
+
+ // ToDo : generate macro
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType1 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType2 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType3 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType4 >::type , container_type >::value ));
+
+ // ToDo : pack into detail namespace
+ for_each4( boost::begin( s1 ) , boost::end( s1 ) ,
+ boost::begin( s2 ) ,
+ boost::begin( s3 ) ,
+ boost::begin( s4 ) ,
+ op );
+ }
+
+ // ToDo : pack into namespace detail
+ template< class Iterator1 , class Iterator2 , class Iterator3 , class Iterator4 , class Operation >
+ static void for_each4( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3, Iterator4 first4, Operation op )
+ {
+ for( ; first1 != last1 ; )
+ op( *first1++ , *first2++ , *first3++ , *first4++ );
+ }
+
+
+
+ template< class StateType1 , class StateType2 , class StateType3 , class StateType4 , class StateType5 , class Operation >
+ static void for_each5( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , StateType4 &s4 , StateType5 &s5 , Operation op )
+ {
+ // ToDo : check that number of arguments of the operation is equal 3
+
+ // ToDo : generate macro
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType1 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType2 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType3 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType4 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType5 >::type , container_type >::value ));
+
+ // ToDo : pack into detail namespace
+ for_each5( boost::begin( s1 ) , boost::end( s1 ) ,
+ boost::begin( s2 ) ,
+ boost::begin( s3 ) ,
+ boost::begin( s4 ) ,
+ boost::begin( s5 ) ,
+ op );
+ }
+
+ // ToDo : pack into namespace detail
+ template< class Iterator1 , class Iterator2 , class Iterator3 , class Iterator4 , class Iterator5 , class Operation >
+ static void for_each5( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3,
+ Iterator4 first4, Iterator5 first5, Operation op )
+ {
+ for( ; first1 != last1 ; )
+ op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ );
+ }
+
+
+
+ template< class StateType1 , class StateType2 , class StateType3 , class StateType4 , class StateType5 , class StateType6 , class Operation >
+ static void for_each6( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , StateType4 &s4 , StateType5 &s5 , StateType6 &s6 , Operation op )
+ {
+ // ToDo : check that number of arguments of the operation is equal 6
+
+ // ToDo : generate macro
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType1 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType2 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType3 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType4 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType5 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType6 >::type , container_type >::value ));
+
+ // ToDo : pack into detail namespace
+ for_each6( boost::begin( s1 ) , boost::end( s1 ) ,
+ boost::begin( s2 ) ,
+ boost::begin( s3 ) ,
+ boost::begin( s4 ) ,
+ boost::begin( s5 ) ,
+ boost::begin( s6 ) ,
+ op );
+ }
+
+ // ToDo : pack into namespace detail
+ template< class Iterator1 , class Iterator2 , class Iterator3 , class Iterator4 , class Iterator5 , class Iterator6 , class Operation >
+ static void for_each6( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3,
+ Iterator4 first4, Iterator5 first5, Iterator6 first6 , Operation op )
+ {
+ for( ; first1 != last1 ; )
+ op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ , *first6++ );
+ }
+
+
+ template< class StateType1 , class StateType2 , class StateType3 , class StateType4 , class StateType5 , class StateType6 ,class StateType7 , class Operation >
+ static void for_each7( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , StateType4 &s4 , StateType5 &s5 , StateType6 &s6 , StateType7 &s7 , Operation op )
+ {
+ // ToDo : check that number of arguments of the operation is equal 6
+
+ // ToDo : generate macro
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType1 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType2 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType3 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType4 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType5 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType6 >::type , container_type >::value ));
+ BOOST_STATIC_ASSERT(( boost::is_same< typename boost::remove_const< StateType7 >::type , container_type >::value ));
+
+ // ToDo : pack into detail namespace
+ for_each7( boost::begin( s1 ) , boost::end( s1 ) ,
+ boost::begin( s2 ) ,
+ boost::begin( s3 ) ,
+ boost::begin( s4 ) ,
+ boost::begin( s5 ) ,
+ boost::begin( s6 ) ,
+ boost::begin( s7 ) ,
+ op );
+ }
+
+ // ToDo : pack into namespace detail
+ template< class Iterator1 , class Iterator2 , class Iterator3 , class Iterator4 , class Iterator5 , class Iterator6 , class Iterator7 , class Operation >
+ static void for_each7( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3,
+ Iterator4 first4, Iterator5 first5, Iterator6 first6 , Iterator7 first7 , Operation op )
+ {
+ for( ; first1 != last1 ; )
+ op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ , *first6++ , *first7++ );
+ }
+
};
} // odeint
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_operations.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_operations.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_operations.hpp 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
@@ -52,6 +52,63 @@
t1 = m_alpha1 * t2 + m_alpha2 * t3;
}
};
+
+ struct scale_sum3
+ {
+ time_type m_alpha1 , m_alpha2 , m_alpha3;
+
+ scale_sum3( time_type alpha1 , time_type alpha2 , time_type alpha3 )
+ : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) { }
+
+ template< class T1 , class T2 , class T3 , class T4 >
+ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 ) const
+ {
+ t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4;
+ }
+ };
+
+ struct scale_sum4
+ {
+ time_type m_alpha1 , m_alpha2 , m_alpha3 , m_alpha4;
+
+ scale_sum4( time_type alpha1 , time_type alpha2 , time_type alpha3 , time_type alpha4)
+ : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) { }
+
+ template< class T1 , class T2 , class T3 , class T4 , class T5 >
+ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5) const
+ {
+ t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5;
+ }
+ };
+
+ struct scale_sum5
+ {
+ time_type m_alpha1 , m_alpha2 , m_alpha3 , m_alpha4 , m_alpha5;
+
+ scale_sum5( time_type alpha1 , time_type alpha2 , time_type alpha3 , time_type alpha4 , time_type alpha5)
+ : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) { }
+
+ template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 >
+ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6) const
+ {
+ t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6;
+ }
+ };
+
+ struct scale_sum6
+ {
+ time_type m_alpha1 , m_alpha2 , m_alpha3 , m_alpha4 , m_alpha5 , m_alpha6;
+
+ scale_sum6( time_type alpha1 , time_type alpha2 , time_type alpha3 , time_type alpha4 , time_type alpha5 , time_type alpha6 )
+ : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ){ }
+
+ template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 >
+ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 ,const T7 &t7) const
+ {
+ t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7;
+ }
+ };
+
};
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_euler.hpp 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
@@ -63,7 +63,7 @@
template< class System >
void do_step( System system , state_type &x , const state_type &dxdt , time_type t , time_type dt )
{
- algebra_type::transform2( x , dxdt , typename operations_type::increment( dt ) );
+ algebra_type::for_each2( x , dxdt , typename operations_type::increment( dt ) );
}
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/runge_kutta_error_ck.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/runge_kutta_error_ck.hpp 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
@@ -0,0 +1,159 @@
+/*
+ boost header: BOOST_NUMERIC_ODEINT/runge_kutta_error.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_RUNGE_KUTTA_ERROR_HPP_INCLUDED
+#define BOOST_BOOST_NUMERIC_ODEINT_RUNGE_KUTTA_ERROR_HPP_INCLUDED
+
+#include <boost/numeric/odeint/algebra/standard_algebra.hpp>
+#include <boost/numeric/odeint/algebra/standard_operations.hpp>
+#include <boost/numeric/odeint/algebra/standard_resize.hpp>
+
+#include <boost/numeric/odeint/stepper/explicit_stepper_base.hpp>
+#include <boost/numeric/odeint/stepper/detail/macros.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+/*
+ * ToDO: write error_stepper_base
+ */
+
+template<
+ class State ,
+ class Time = double ,
+ class Algebra = standard_algebra< State > ,
+ class Operations = standard_operations< Time > ,
+ class AdjustSizePolicy = adjust_size_initially_tag
+ >
+class runge_kutta_error_ck
+: public explicit_stepper_base<
+ runge_kutta_error_ck< State , Time , Algebra , Operations , AdjustSizePolicy > ,
+ 5 , State , Time , Algebra , Operations , AdjustSizePolicy >
+{
+
+public :
+
+ BOOST_ODEINT_EXPLICIT_STEPPERS_TYPEDEFS( runge_kutta_error_ck , 5 );
+
+ friend class explicit_stepper_base< runge_kutta_error_ck< State , Time , Algebra , Operations , AdjustSizePolicy > ,
+ 5 , State , Time , Algebra , Operations , AdjustSizePolicy >;
+
+ runge_kutta_error_ck( void ) : m_dxdt() , m_x1() , m_x2() , m_x3() , m_x4() , m_x5() , m_x6()
+ { }
+
+ template< class System >
+ void do_step( System system , state_type &x , time_type t , time_type dt , state_type &xerr)
+ {
+ this->adjust_size_by_policy( x , adjust_size_policy() );
+ system( x , m_dxdt ,t );
+ do_step( system , x , m_dxdt , t , dt , xerr);
+ }
+
+
+ template< class System >
+ void do_step( System system , state_type &x , const state_type &dxdt , time_type t , time_type dt , state_type &xerr )
+ {
+ /* ToDo: separate resize m_dxdt and m_x1..6 */
+ this->adjust_size_by_policy( x , adjust_size_policy() );
+
+ const time_type a2 = static_cast<time_type> ( 0.2 );
+ const time_type a3 = static_cast<time_type> ( 0.3 );
+ const time_type a4 = static_cast<time_type> ( 0.6 );
+ const time_type a5 = static_cast<time_type> ( 1.0 );
+ const time_type a6 = static_cast<time_type> ( 0.875 );
+
+ const time_type b21 = static_cast<time_type> ( 0.2 );
+ const time_type b31 = static_cast<time_type> ( 3.0 ) / static_cast<time_type>( 40.0 );
+ const time_type b32 = static_cast<time_type> ( 9.0 ) / static_cast<time_type>( 40.0 );
+ const time_type b41 = static_cast<time_type> ( 0.3 );
+ const time_type b42 = static_cast<time_type> ( -0.9 );
+ const time_type b43 = static_cast<time_type> ( 1.2 );
+ const time_type b51 = static_cast<time_type> ( -11.0 ) / static_cast<time_type>( 54.0 );
+ const time_type b52 = static_cast<time_type> ( 2.5 );
+ const time_type b53 = static_cast<time_type> ( -70.0 ) / static_cast<time_type>( 27.0 );
+ const time_type b54 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 27.0 );
+ const time_type b61 = static_cast<time_type> ( 1631.0 ) / static_cast<time_type>( 55296.0 );
+ const time_type b62 = static_cast<time_type> ( 175.0 ) / static_cast<time_type>( 512.0 );
+ const time_type b63 = static_cast<time_type> ( 575.0 ) / static_cast<time_type>( 13824.0 );
+ const time_type b64 = static_cast<time_type> ( 44275.0 ) / static_cast<time_type>( 110592.0 );
+ const time_type b65 = static_cast<time_type> ( 253.0 ) / static_cast<time_type>( 4096.0 );
+
+ const time_type c1 = static_cast<time_type> ( 37.0 ) / static_cast<time_type>( 378.0 );
+ const time_type c3 = static_cast<time_type> ( 250.0 ) / static_cast<time_type>( 621.0 );
+ const time_type c4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 594.0 );
+ const time_type c6 = static_cast<time_type> ( 512.0 ) / static_cast<time_type>( 1771.0 );
+
+ const time_type dc1 = c1 - static_cast<time_type> ( 2825.0 ) / static_cast<time_type>( 27648 );
+ const time_type dc3 = c3 - static_cast<time_type> ( 18575.0 ) / static_cast<time_type>( 48384.0 );
+ const time_type dc4 = c4 - static_cast<time_type> ( 13525.0 ) / static_cast<time_type>( 55296.0 );
+ const time_type dc5 = static_cast<time_type> ( -277.0 ) / static_cast<time_type>( 14336.0 );
+ const time_type dc6 = c6 - static_cast<time_type> ( 0.25 );
+
+ //m_x1 = x + dt*b21*dxdt
+ algebra_type::for_each3( m_x1 , x , dxdt ,
+ typename operations_type::scale_sum2( 1.0 , dt*b21 ) );
+
+ system( m_x1 , m_x2 , t + dt*a2 );
+ // m_x1 = x + dt*b31*dxdt + dt*b32*m_x2
+ algebra_type::for_each4( m_x1 , x , dxdt , m_x2 ,
+ typename operations_type::scale_sum3( 1.0 , dt*b31 , dt*b32 ));
+
+ system( m_x1 , m_x3 , t + dt*a3 );
+ // m_x1 = x + dt * (b41*dxdt + b42*m_x2 + b43*m_x3)
+ algebra_type::for_each5( m_x1 , x , dxdt , m_x2 , m_x3 ,
+ typename operations_type::scale_sum4( 1.0 , dt*b41 , dt*b42 , dt*b43 ));
+
+ system( m_x1, m_x4 , t + dt*a4 );
+ algebra_type::for_each6( m_x1 , x , dxdt , m_x2 , m_x3 , m_x4 ,
+ typename operations_type::scale_sum5( 1.0 , dt*b51 , dt*b52 , dt*b53 , dt*b54 ));
+
+ system( m_x1 , m_x5 , t + dt*a5 );
+ algebra_type::for_each7( m_x1 , x , dxdt , m_x2 , m_x3 , m_x4 , m_x5 ,
+ typename operations_type::scale_sum6( 1.0 , dt*b61 , dt*b62 , dt*b63 , dt*b64 , dt*b65 ));
+
+ system( m_x1 , m_x6 , t + dt*a6 );
+ algebra_type::for_each6( x , x , dxdt , m_x3 , m_x4 , m_x6 ,
+ typename operations_type::scale_sum5( 1.0 , dt*c1 , dt*c3 , dt*c4 , dt*c6 ));
+
+ //error estimate
+ algebra_type::for_each6( xerr , dxdt , m_x3 , m_x4 , m_x5 , m_x6 ,
+ typename operations_type::scale_sum5( dt*dc1 , dt*dc3 , dt*dc4 , dt*dc5 , dt*dc6 ));
+
+ }
+
+
+
+private:
+
+ void adjust_size_impl( const state_type &x )
+ {
+ adjust_size( x , m_dxdt );
+ adjust_size( x , m_x1 );
+ adjust_size( x , m_x2 );
+ adjust_size( x , m_x3 );
+ adjust_size( x , m_x4 );
+ adjust_size( x , m_x5 );
+ adjust_size( x , m_x6 );
+ }
+
+ state_type m_dxdt;
+ state_type m_x1, m_x2, m_x3, m_x4, m_x5, m_x6;
+
+
+};
+
+
+} // odeint
+} // numeric
+} // boost
+
+#endif //BOOST_BOOST_NUMERIC_ODEINT_RUNGE_KUTTA_ERROR_HPP_INCLUDED
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp 2010-07-12 11:54:16 EDT (Mon, 12 Jul 2010)
@@ -70,7 +70,19 @@
BOOST_CHECK_SMALL( fabs( xval - 0.1 ) , eps );
}
+template< class Stepper , class System >
+void check_error_stepper_concept( Stepper &stepper , System system ,
+ typename Stepper::state_type &x , typename Stepper::state_type &xerr )
+{
+ typedef Stepper stepper_type;
+ typedef typename stepper_type::state_type container_type;
+ typedef typename stepper_type::order_type order_type;
+ typedef typename stepper_type::time_type time_type;
+ stepper.do_step( system , x , 0.0 , 0.1 , xerr);
+ double xval = * boost::begin( x );
+ BOOST_CHECK_SMALL( fabs( xval - 0.1 ) , eps );
+}
//template< class ErrorStepper >
@@ -210,6 +222,14 @@
check_stepper_concept( euler , constant_system4 , x );
}
+void test_runge_kutta_error_ck_with_vector( void )
+{
+ state_type1 x( 1 , 0.0 );
+ state_type1 xerr( 1 , 0.0 );
+ runge_kutta_error_ck< state_type1 > rk_ck;
+ check_error_stepper_concept( rk_ck , constant_system1 , x , xerr );
+}
+
test_suite* init_unit_test_suite( int argc, char* argv[] )
{
test_suite *test = BOOST_TEST_SUITE("check stepper concepts");
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk