|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r66378 - in sandbox/odeint/branches/karsten: boost/numeric/odeint/algebra boost/numeric/odeint/algebra/detail boost/numeric/odeint/stepper libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-11-03 08:41:30
Author: karsten
Date: 2010-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
New Revision: 66378
URL: http://svn.boost.org/trac/boost/changeset/66378
Log:
* dense output dopri5 finalized
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp | 9 +++
sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_algebra.hpp | 24 +++++++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/standard_operations.hpp | 16 ++++++
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp | 100 ++++++++++++++++++++++++++-------------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_dopri5.hpp | 3 +
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp | 10 ----
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp | 8 ++
7 files changed, 124 insertions(+), 46 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp 2010-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -70,6 +70,15 @@
op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ , *first6++ , *first7++ );
}
+template< class Iterator1 , class Iterator2 , class Iterator3 , class Iterator4 , class Iterator5 , class Iterator6 , class Iterator7 , class Iterator8 , class Operation >
+inline void for_each8( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3,
+ Iterator4 first4, Iterator5 first5, Iterator6 first6 , Iterator7 first7 , Iterator8 first8 , Operation op )
+{
+ for( ; first1 != last1 ; )
+ op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ , *first6++ , *first7++ , *first8++ );
+}
+
+
} // detail
} // odeint
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-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -143,6 +143,30 @@
op );
}
+ template< class StateType1 , class StateType2 , class StateType3 , class StateType4 , class StateType5 , class StateType6 ,class StateType7 , class StateType8 , class Operation >
+ static void for_each8( StateType1 &s1 , StateType2 &s2 , StateType3 &s3 , StateType4 &s4 , StateType5 &s5 , StateType6 &s6 , StateType7 &s7 , StateType8 &s8 , Operation op )
+ {
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType1 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType2 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType3 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType4 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType5 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType6 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType7 , container_type );
+ BOOST_ODEINT_CHECK_CONTAINER_TYPE( StateType8 , container_type );
+
+ detail::for_each8( boost::begin( s1 ) , boost::end( s1 ) ,
+ boost::begin( s2 ) ,
+ boost::begin( s3 ) ,
+ boost::begin( s4 ) ,
+ boost::begin( s5 ) ,
+ boost::begin( s6 ) ,
+ boost::begin( s7 ) ,
+ boost::begin( s8 ) ,
+ op );
+ }
+
+
template< class ValueType , class StateType , class Reduction >
static ValueType reduce( StateType &s , Reduction red , ValueType init)
{
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-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -98,6 +98,22 @@
}
};
+ struct scale_sum7
+ {
+ const time_type m_alpha1 , m_alpha2 , m_alpha3 , m_alpha4 , m_alpha5 , m_alpha6 , m_alpha7;
+
+ scale_sum7( const time_type alpha1 , const time_type alpha2 , const time_type alpha3 , const time_type alpha4 ,
+ const time_type alpha5 , const time_type alpha6 , const time_type alpha7 )
+ : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) { }
+
+ template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 >
+ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 ) const
+ {
+ t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8;
+ }
+ };
+
+
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/dense_output_dopri5.hpp 2010-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -127,56 +127,88 @@
* C = theta^2 * ( theta - 1 )^2
* D = theta * ( theta - 1 )^2
*
- * b_1( theta ) = A * b_1 + C * X1( theta ) + D
+ * 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_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_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
*/
void calc_state( time_type t , state_type &x )
{
+ const time_type b1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
+ const time_type b3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
+ const time_type b4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
+ const time_type b5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
+ const time_type b6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
+
+ time_type dt = ( m_t - m_t_old );
+ time_type theta = ( t - m_t_old ) / dt;
+ time_type X1 = static_cast< time_type >( 5.0 ) * ( static_cast< time_type >( 2558722523.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 11282082432.0 );
+ time_type X3 = static_cast< time_type >( 100.0 ) * ( static_cast< time_type >( 882725551.0 ) - static_cast< time_type >( 15701508.0 ) * theta ) / static_cast< time_type >( 32700410799.0 );
+ time_type X4 = static_cast< time_type >( 25.0 ) * ( static_cast< time_type >( 443332067.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 1880347072.0 ) ;
+ time_type X5 = static_cast< time_type >( 32805.0 ) * ( static_cast< time_type >( 23143187.0 ) - static_cast< time_type >( 3489224.0 ) * theta ) / static_cast< time_type >( 199316789632.0 );
+ time_type X6 = static_cast< time_type >( 55.0 ) * ( static_cast< time_type >( 29972135.0 ) - static_cast< time_type >( 7076736.0 ) * theta ) / static_cast< time_type >( 822651844.0 );
+ time_type X7 = static_cast< time_type >( 10.0 ) * ( static_cast< time_type >( 7414447.0 ) - static_cast< time_type >( 829305.0 ) * theta ) / static_cast< time_type >( 29380423.0 );
+
+ time_type theta_m_1 = theta - static_cast< time_type >( 1.0 );
+ time_type theta_sq = theta * theta;
+ time_type A = theta_sq * ( static_cast< time_type >( 3.0 ) - static_cast< time_type >( 2.0 ) * theta );
+ time_type B = theta_sq * theta_m_1;
+ time_type C = theta_sq * theta_m_1 * theta_m_1;
+ time_type D = theta * theta_m_1 * theta_m_1;
+
+ time_type b1_theta = A * b1 - C * X1 + D;
+ time_type b3_theta = A * b3 + C * X3;
+ time_type b4_theta = A * b4 - C * X4;
+ time_type b5_theta = A * b5 + C * X5;
+ time_type b6_theta = A * b6 - C * X6;
+ time_type b7_theta = B + C * X7;
+
+ const state_type &k1 = *m_old_deriv;
+ const state_type &k3 = dopri5().m_x3;
+ const state_type &k4 = dopri5().m_x4;
+ const state_type &k5 = dopri5().m_x5;
+ const state_type &k6 = dopri5().m_x6;
+ const state_type &k7 = *m_current_deriv;
+
+ algebra_type::for_each8( x , *m_old_state , k1 , k3 , k4 , k5 , k6 , k7 ,
+ typename operations_type::scale_sum7( 1.0 , dt * b1_theta , dt * b3_theta , dt * b4_theta , dt * b5_theta , dt * b6_theta , dt * b7_theta ) );
+// algebra_type::for_each3( x , *m_old_state , k1 ,
+// typename operations_type::scale_sum2( 1.0 , dt * theta ) );
+ }
-// const time_type b1 = static_cast<time_type> ( 35.0 ) / static_cast<time_type>( 384.0 );
-// const time_type b3 = static_cast<time_type> ( 500.0 ) / static_cast<time_type>( 1113.0 );
-// const time_type b4 = static_cast<time_type> ( 125.0 ) / static_cast<time_type>( 192.0 );
-// const time_type b5 = static_cast<time_type> ( -2187.0 ) / static_cast<time_type>( 6784.0 );
-// const time_type b6 = static_cast<time_type> ( 11.0 ) / static_cast<time_type>( 84.0 );
-//
-//
-// time_type theta = t - m_t_old;
-// time_type X1 = static_cast< time_type >( 5.0 ) * ( static_cast< time_type >( 2558722523.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 11282082432.0 );
-// time_type X3 = static_cast< time_type >( 100.0 ) * ( static_cast< time_type >( 882725551.0 ) - static_cast< time_type >( 15701508.0 ) * theta ) / static_cast< time_type >( 32700410799.0 );
-// time_type X4 = static_cast< time_type >( 25.0 ) * ( static_cast< time_type >( 443332067.0 ) - static_cast< time_type >( 31403016.0 ) * theta ) / static_cast< time_type >( 1880347072.0 ) ;
-// time_type X5 = static_cast< time_type >( 32805.0 ) * ( static_cast< time_type >( 23143187.0 ) - static_cast< time_type >( 3489224.0 ) * theta ) / static_cast< time_type >( 199316789632.0 );
-// time_type X6 = static_cast< time_type >( 55.0 ) * ( static_cast< time_type >( 29972135.0 ) - static_cast< time_type >( 7076736.0 ) * theta ) / static_cast< time_type >( 822651844.0 );
-// time_type X7 = static_cast< time_type >( 10.0 ) * ( static_cast< time_type >( 7414447.0 ) - static_cast< time_type >( 829305.0 ) * theta ) / static_cast< time_type >( 29380423.0 );
-//
-// time_type theta_m_1 = theta - static_cast< time_type >( 1.0 );
-// time_type theta_sq = theta * theta;
-// time_type A = theta_sq * ( static_cast< time_type >( 3.0 ) - static_cast< time_type >( 2.0 ) * theta );
-// time_type B = theta_sq * theta_m_1;
-// time_type C = theta_sq * theta_m_1 * theta_m_1;
-// time_type D = theta * theta_m_1 * theta_m_1;
-
-// time_type b1_theta = A *
-
-
-
+ const state_type& current_state( void ) const
+ {
+ return *m_current_state;
+ }
+ const time_type& current_time( void ) const
+ {
+ return m_t;
+ }
-// algebra_type::for_each3( x , *m_old_state , m_euler.m_dxdt , typename operations_type::scale_sum2( 1.0 , delta ) );
+ const time_type& previous_state( void ) const
+ {
+ return *m_old_state;
}
- const state_type& current_state( void ) const { return *m_current_state; }
- const time_type& current_time( void ) const { return m_t; }
- const time_type& previous_state( void ) const { return *m_old_state; }
- const time_type& previous_time( void ) const { return m_t_old; }
+ const time_type& previous_time( void ) const
+ {
+ return m_t_old;
+ }
private:
+ dopri5_type& dopri5( void ) { return m_stepper.stepper(); }
+ const dopri5_type& dopri5( void ) const { return m_stepper.stepper(); }
+
stepper_type &m_stepper;
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 2010-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -42,6 +42,9 @@
public :
+ template< class ControlledStepper >
+ friend class dense_output_dopri5;
+
BOOST_ODEINT_EXPLICIT_STEPPERS_AND_ERROR_STEPPERS_TYPEDEFS( explicit_error_dopri5 , 5 , 5 , 4 );
typedef explicit_error_stepper_fsal_tag stepper_category;
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-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -16,9 +16,6 @@
#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 {
@@ -58,15 +55,10 @@
m_jacobi *= dt;
m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
- 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
@@ -85,8 +77,6 @@
solve( m_b , jacobi_tmp );
- std::clog << m_b << std::endl;
-
m_x -= m_b;
}
x = m_x;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/dense_output_stepper_evolution.cpp 2010-11-03 08:41:19 EDT (Wed, 03 Nov 2010)
@@ -60,12 +60,16 @@
state_type orig = sys_solution( t , x0 );
state_type diff = {{ orig[0] - x[0] , orig[1] - x[1] }};
double diff_abs = sqrt( diff[0] * diff[0] + diff[1] * diff[1] );
- state_out << t << tab << x << tab << orig << tab << diff << diff_abs << endl;
+ state_out << t << tab << x << tab << orig << tab << diff << tab << diff_abs << endl;
}
else
{
stepper.do_step( sys );
- stepper_out << stepper.current_time() << "\t" << stepper.current_state() << std::endl;
+ state_type orig = sys_solution( stepper.current_time() , x0 );
+ const state_type &xx = stepper.current_state();
+ state_type diff = {{ orig[0] - xx[0] , orig[1] - xx[1] }};
+ double diff_abs = sqrt( diff[0] * diff[0] + diff[1] * diff[1] );
+ stepper_out << stepper.current_time() << "\t" << xx << tab << orig << tab << diff << tab << diff_abs << std::endl;
continue;
}
t += dt;
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