mirror of
https://github.com/boostorg/odeint.git
synced 2025-05-11 13:34:09 +00:00
first implementation of error stepper rk_ck
This commit is contained in:
mariomulansky
parent
02b6c43d77
commit
044376a849
@ -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
|
||||
|
||||
@ -25,7 +25,7 @@ struct standard_algebra
|
||||
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 @@ struct standard_algebra
|
||||
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 @@ struct standard_algebra
|
||||
|
||||
|
||||
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 @@ struct standard_algebra
|
||||
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 @@ struct standard_algebra
|
||||
|
||||
// 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
|
||||
|
||||
@ -52,6 +52,63 @@ struct standard_operations
|
||||
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;
|
||||
}
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -63,7 +63,7 @@ public :
|
||||
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 ) );
|
||||
}
|
||||
|
||||
|
||||
|
||||
159
boost/numeric/odeint/stepper/runge_kutta_error_ck.hpp
Normal file
159
boost/numeric/odeint/stepper/runge_kutta_error_ck.hpp
Normal file
@ -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
|
||||
@ -70,7 +70,19 @@ void check_stepper_concept( Stepper &stepper , System system , typename Stepper:
|
||||
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 @@ void test_euler_with_array( void )
|
||||
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");
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user