|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r68531 - in sandbox/odeint/branches/karsten: . boost/numeric/odeint/stepper boost/numeric/odeint/stepper/base
From: karsten.ahnert_at_[hidden]
Date: 2011-01-28 08:39:45
Author: karsten
Date: 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
New Revision: 68531
URL: http://svn.boost.org/trac/boost/changeset/68531
Log:
* dense output dopri5
Text files modified:
sandbox/odeint/branches/karsten/Jamroot | 2
sandbox/odeint/branches/karsten/TODO | 2
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp | 6 ---
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 2 -
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_controlled_explicit_fsal.hpp | 17 ++-------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit.hpp | 12 +-----
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp | 71 ++++++++++++++++++++++++++++++++++++++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_euler.hpp | 6 +-
8 files changed, 82 insertions(+), 36 deletions(-)
Modified: sandbox/odeint/branches/karsten/Jamroot
==============================================================================
--- sandbox/odeint/branches/karsten/Jamroot (original)
+++ sandbox/odeint/branches/karsten/Jamroot 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -34,7 +34,7 @@
# docs:
-build-project libs/numeric/odeint/doc ;
+# build-project libs/numeric/odeint/doc ;
Modified: sandbox/odeint/branches/karsten/TODO
==============================================================================
--- sandbox/odeint/branches/karsten/TODO (original)
+++ sandbox/odeint/branches/karsten/TODO 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -7,6 +7,7 @@
* test, if copy construct of stepper_base is called when explicit_euler is used
* test gsl
* test explicit stepper with ranges
+ * test units with dense output
OK * change standard_operations::rel_error in order to word with units and test it
OK * include implicit euler
OK * call via std::pair< deriv , jacobi >
@@ -36,7 +37,6 @@
* check header guards
* check copyright note
* documente every file in the preamble
-* include adjust_size_by_policy in all steppers and call this routine inside do_step()
* check once more, if all contructor, destructors and assign-operators are present
* Integrate functions
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -134,12 +134,6 @@
m_size_adjuster.adjust_size( x );
}
- template< class StateType >
- void adjust_size_by_policy( const StateType &x )
- {
- m_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- }
-
protected:
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -262,8 +262,6 @@
}
-
-
stepper_type& stepper( void )
{
return m_stepper;
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_controlled_explicit_fsal.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_controlled_explicit_fsal.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_controlled_explicit_fsal.hpp 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -126,7 +126,8 @@
template< class StateType >
void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
{
- adjust_size_by_policy( x0 );
+ m_state_adjuster.adjust_size_by_policy( x0 , adjust_size_policy() );
+ m_deriv_adjuster.adjust_size_by_policy( x0 , adjust_size_policy() );
boost::numeric::odeint::copy( x0 , *m_current_state );
m_t = t0;
m_dt = dt0;
@@ -152,7 +153,7 @@
{
res = m_stepper.try_step( system , *m_current_state , *m_current_deriv , m_t , *m_old_state , *m_old_deriv , m_dt );
if( count++ == max_count )
- throw std::overflow_error( "dense_output_dopri5 : too much iterations!");
+ throw std::overflow_error( "dense_output_controlled_explicit_fsal : too much iterations!");
}
while( res == step_size_decreased );
std::swap( m_current_state , m_old_state );
@@ -163,7 +164,7 @@
template< class StateOut >
void calc_state( const time_type &t , StateOut &x )
{
-// m_stepper.calc_state( x , t , *m_old_state , m_t_old );
+ m_stepper.stepper().calc_state( t , x , *m_old_state , *m_old_deriv , m_t_old , *m_current_state , *m_current_deriv , m_t );
}
template< class StateType >
@@ -174,16 +175,6 @@
m_stepper.adjust_size( x );
}
- template< class StateType >
- void adjust_size_by_policy( const StateType &x )
- {
- m_state_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- m_deriv_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- // ToDo : implement this in all stepper
-// m_stepper.adjust_size_by_policy( x );
- }
-
-
const state_type& current_state( void ) const
{
return *m_current_state;
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_explicit.hpp 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -100,7 +100,7 @@
template< class StateType >
void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 )
{
- adjust_size_by_policy( x0 );
+ m_size_adjuster.adjust_size_by_policy( x0 , adjust_size_policy() );
boost::numeric::odeint::copy( x0 , *m_current_state );
m_t = t0;
m_dt = dt0;
@@ -119,7 +119,7 @@
template< class StateOut >
void calc_state( const time_type &t , StateOut &x )
{
- m_stepper.calc_state( x , t , *m_old_state , m_t_old );
+ m_stepper.calc_state( x , t , *m_old_state , m_t_old , *m_current_state , m_t );
}
template< class StateType >
@@ -129,14 +129,6 @@
m_stepper.adjust_size( x );
}
- template< class StateType >
- void adjust_size_by_policy( const StateType &x )
- {
- m_size_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- m_stepper.adjust_size_by_policy( x );
- }
-
-
const state_type& current_state( void ) const
{
return *m_current_state;
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -212,6 +212,77 @@
}
+ /*
+ * Calculates Dense-Output for Dopri5
+ *
+ * See Hairer, Norsett, Wanner: Solving Ordinary Differential Equations, Nonstiff Problems. I, p.191/192
+ *
+ * y(t+theta) = y(t) + h * sum_i^7 b_i(theta) * k_i
+ *
+ * A = theta^2 * ( 3 - 2 theta )
+ * B = theta^2 * ( theta - 1 )
+ * C = theta^2 * ( theta - 1 )^2
+ * D = theta * ( theta - 1 )^2
+ *
+ * b_1( theta ) = A * b_1 - C * X1( theta ) + D
+ * b_2( theta ) = 0
+ * b_3( theta ) = A * b_3 + C * X3( theta )
+ * b_4( theta ) = A * b_4 - C * X4( theta )
+ * b_5( theta ) = A * b_5 + C * X5( theta )
+ * b_6( theta ) = A * b_6 - C * X6( theta )
+ * b_7( theta ) = B + C * X7( theta )
+ *
+ * An alternative Method is described in:
+ *
+ * www-m2.ma.tum.de/homepages/simeon/numerik3/kap3.ps
+ */
+ template< class StateOut , class StateIn1 , class DerivIn1 , class StateIn2 , class DerivIn2 >
+ void calc_state( const time_type &t , StateOut &x ,
+ const StateIn1 &x_old , const DerivIn1 &deriv_old , const time_type &t_old ,
+ const StateIn2 &x_new , const DerivIn2 &deriv_new , const time_type &t_new )
+ {
+ const value_type b1 = static_cast<value_type> ( 35.0 ) / static_cast<value_type>( 384.0 );
+ const value_type b3 = static_cast<value_type> ( 500.0 ) / static_cast<value_type>( 1113.0 );
+ const value_type b4 = static_cast<value_type> ( 125.0 ) / static_cast<value_type>( 192.0 );
+ const value_type b5 = static_cast<value_type> ( -2187.0 ) / static_cast<value_type>( 6784.0 );
+ const value_type b6 = static_cast<value_type> ( 11.0 ) / static_cast<value_type>( 84.0 );
+
+ time_type dt = ( t_new - t_old );
+ value_type theta = ( t - t_old ) / dt;
+ value_type X1 = static_cast< value_type >( 5.0 ) * ( static_cast< value_type >( 2558722523.0 ) - static_cast< value_type >( 31403016.0 ) * theta ) / static_cast< value_type >( 11282082432.0 );
+ value_type X3 = static_cast< value_type >( 100.0 ) * ( static_cast< value_type >( 882725551.0 ) - static_cast< value_type >( 15701508.0 ) * theta ) / static_cast< value_type >( 32700410799.0 );
+ value_type X4 = static_cast< value_type >( 25.0 ) * ( static_cast< value_type >( 443332067.0 ) - static_cast< value_type >( 31403016.0 ) * theta ) / static_cast< value_type >( 1880347072.0 ) ;
+ value_type X5 = static_cast< value_type >( 32805.0 ) * ( static_cast< value_type >( 23143187.0 ) - static_cast< value_type >( 3489224.0 ) * theta ) / static_cast< value_type >( 199316789632.0 );
+ value_type X6 = static_cast< value_type >( 55.0 ) * ( static_cast< value_type >( 29972135.0 ) - static_cast< value_type >( 7076736.0 ) * theta ) / static_cast< value_type >( 822651844.0 );
+ value_type X7 = static_cast< value_type >( 10.0 ) * ( static_cast< value_type >( 7414447.0 ) - static_cast< value_type >( 829305.0 ) * theta ) / static_cast< value_type >( 29380423.0 );
+
+ value_type theta_m_1 = theta - static_cast< value_type >( 1.0 );
+ value_type theta_sq = theta * theta;
+ value_type A = theta_sq * ( static_cast< value_type >( 3.0 ) - static_cast< value_type >( 2.0 ) * theta );
+ value_type B = theta_sq * theta_m_1;
+ value_type C = theta_sq * theta_m_1 * theta_m_1;
+ value_type D = theta * theta_m_1 * theta_m_1;
+
+ value_type b1_theta = A * b1 - C * X1 + D;
+ value_type b3_theta = A * b3 + C * X3;
+ value_type b4_theta = A * b4 - C * X4;
+ value_type b5_theta = A * b5 + C * X5;
+ value_type b6_theta = A * b6 - C * X6;
+ value_type b7_theta = B + C * X7;
+
+// const state_type &k1 = *m_old_deriv;
+// const state_type &k3 = dopri5().m_k3;
+// const state_type &k4 = dopri5().m_k4;
+// const state_type &k5 = dopri5().m_k5;
+// const state_type &k6 = dopri5().m_k6;
+// const state_type &k7 = *m_current_deriv;
+
+ typename algebra_type::for_each8()( x , *x_old , *deriv_old , m_k3 , m_k4 , m_k5 , m_k6 , deriv_new ,
+ typename operations_type::template scale_sum7< value_type , time_type , time_type , time_type , time_type , time_type , time_type >( 1.0 , dt * b1_theta , dt * b3_theta , dt * b4_theta , dt * b5_theta , dt * b6_theta , dt * b7_theta ) );
+ }
+
+
+
void adjust_size( const state_type &x )
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 2011-01-28 08:39:41 EST (Fri, 28 Jan 2011)
@@ -53,11 +53,11 @@
}
- template< class StateOut , class StateIn >
- void calc_state( StateOut &x , const time_type &t , const StateIn &old_state , const time_type &t_old )
+ template< class StateOut , class StateIn1 , class StateIn2 >
+ void calc_state( StateOut &x , const time_type &t , const StateIn1 &old_state , const time_type &t_old , const StateIn2 ¤t_state , const time_type &t_new )
{
time_type delta = t - t_old;
- typename algebra_type::for_each3()( x , old_state , stepper_base_type::m_dxdt , typename operations_type::template scale_sum2< time_type , time_type >( 1.0 , delta ) );
+ typename algebra_type::for_each3()( x , old_state , stepper_base_type::m_dxdt , typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , delta ) );
}
};
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