Subject: [Boost-commit] svn:boost r60390 - in sandbox/odeint: . boost/numeric/odeint libs/numeric/odeint/examples libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-03-09 17:36:15
Author: karsten
Date: 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
New Revision: 60390
URL: http://svn.boost.org/trac/boost/changeset/60390
Log:
bug fixes in the stepper, check of the integrate functions
Removed:
sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp
Binary files modified:
sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
Text files modified:
sandbox/odeint/ToDo | 17 +++++---
sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp | 18 +++++---
sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp | 77 +++++++++++++++++++++++----------------
sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp | 28 ++++++++-----
sandbox/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp | 15 ++++---
sandbox/odeint/boost/numeric/odeint/integrator_constant_stepsize.hpp | 2 +
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 2
sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp | 4 +-
sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp | 3 +
sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp | 2
sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp | 31 +++++++--------
sandbox/odeint/libs/numeric/odeint/examples/Jamfile | 1
sandbox/odeint/libs/numeric/odeint/examples/doc_integrate.cpp | 4 +-
sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp | 3 +
sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp | 5 +-
sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp | 6 +-
sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp | 7 +--
sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp | 2
sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp | 2
19 files changed, 133 insertions(+), 96 deletions(-)
Modified: sandbox/odeint/ToDo
==============================================================================
--- sandbox/odeint/ToDo (original)
+++ sandbox/odeint/ToDo 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -12,22 +12,27 @@
* typedef typename traits_type::container_type container_type statt typedef Container container_type - DONE
* Concept checks, für vernünftige Fehlermeldungen
* auch für checks ob abs() vorhanden ist
-* check in all steppers the order
-* rk78 die fehler implementieren
-* references entfernen
-* check: orders and adjust_size() in :
+* check: orders and adjust_size() and something more in :
* stepper_euler.hpp - DONE
* stepper_half_step.hpp - DONE
* stepper_midpoint.hpp - DONE, changed stepcount -> step_number
* stepper_rk4_classical.hpp - DONE
* stepper_rk4.hpp - DONE
* stepper_rk5_ck.hpp - DONE
- * stepper_rk78_fehlberg.hpp
- * stepper_rk_generic.hpp
+ * check the error orders
+ * stepper_rk78_fehlberg.hpp - DONE
+ * rk78 die fehler implementieren - DONE
+ * check the error orders
+ * controlled_stepper_standard.hpp - DONE
+ * contruction of the internal stepper from its standard constructor - DONE
+ * controlled_stepper_bs.hpp
+
+
* die auskommentierten Iteratoren typedefs entfernen
+
Mario:
* code cleanup:
* error_checker in controlled_stepper_standard - DONE
Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -20,13 +20,11 @@
#include <limits>
#include <vector>
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/stepper_midpoint.hpp>
#include <boost/numeric/odeint/error_checker_standard.hpp>
#include <boost/numeric/odeint/container_traits.hpp>
-#include <iostream>
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
namespace boost {
namespace numeric {
@@ -126,6 +124,16 @@
}
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size(x, m_xerr);
+ traits_type::adjust_size(x, m_x_mp);
+ traits_type::adjust_size(x, m_x_scale);
+ traits_type::adjust_size(x, m_dxdt);
+ m_stepper_mp.adjust_size( x );
+ }
+
+
template< class DynamicalSystem >
controlled_step_result try_step(
DynamicalSystem &system ,
@@ -134,9 +142,6 @@
time_type &t ,
time_type &dt )
{
- traits_type::adjust_size(x, m_xerr);
- traits_type::adjust_size(x, m_x_mp);
- traits_type::adjust_size(x, m_x_scale);
// get error scale
m_error_checker.fill_scale(x, dxdt, dt, m_x_scale);
@@ -268,7 +273,6 @@
time_type &t,
time_type &dt )
{
- traits_type::adjust_size(x, m_dxdt);
system(x, m_dxdt, t);
return try_step(system, x, m_dxdt, t, dt );
}
Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -18,18 +18,19 @@
#include <algorithm>
#include <complex>
-#include <boost/concept_check.hpp>
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/error_checker_standard.hpp>
#include <boost/numeric/odeint/container_traits.hpp>
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
+
+
+
namespace boost {
namespace numeric {
namespace odeint {
- typedef enum{
+ typedef enum
+ {
success ,
step_size_decreased ,
step_size_increased
@@ -74,17 +75,22 @@
// typedef typename stepper_type::iterator iterator;
// typedef typename stepper_type::const_iterator const_iterator;
+ typedef error_checker_standard<
+ container_type, time_type , traits_type
+ > error_checker_type;
+
// private members
private:
- stepper_type &m_stepper;
- error_checker_standard< container_type, time_type , traits_type > m_error_checker;
+ stepper_type m_stepper;
+ error_checker_type m_error_checker;
time_type m_eps_abs;
time_type m_eps_rel;
time_type m_a_x;
time_type m_a_dxdt;
+
container_type m_dxdt;
container_type m_x_tmp;
container_type m_x_err;
@@ -93,20 +99,43 @@
// public functions
public:
+
+ order_type order_error_step() { return m_stepper.order_error_step(); }
controlled_stepper_standard(
- ErrorStepper &stepper,
- time_type abs_err, time_type rel_err,
- time_type factor_x, time_type factor_dxdt )
- : m_stepper(stepper),
- m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
+ time_type abs_err, time_type rel_err,
+ time_type factor_x, time_type factor_dxdt )
+ : m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
m_eps_abs(abs_err),
m_eps_rel(rel_err),
m_a_x(factor_x),
m_a_dxdt(factor_dxdt)
- { }
+ {
+ }
+
+ controlled_stepper_standard(
+ const container_type &x ,
+ time_type abs_err, time_type rel_err,
+ time_type factor_x, time_type factor_dxdt )
+ : m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
+ m_eps_abs(abs_err),
+ m_eps_rel(rel_err),
+ m_a_x(factor_x),
+ m_a_dxdt(factor_dxdt)
+ {
+ adjust_size( x );
+ }
+
+ void adjust_size( const container_type &x )
+ {
+ traits_type::adjust_size( x , m_x_err );
+ traits_type::adjust_size( x , m_x_scale );
+ traits_type::adjust_size( x , m_dxdt );
+ m_stepper.adjust_size( x );
+ }
+
+
- order_type order() { return m_stepper.order(); }
/* Tries a controlled step with the given stepsize dt. If dt is too large,
x remains unchanged, an appropriate stepsize is assigned to dt and
@@ -119,13 +148,10 @@
controlled_step_result try_step(
DynamicalSystem &system,
container_type &x,
- container_type &dxdt,
+ const container_type &dxdt,
time_type &t,
time_type &dt )
{
- traits_type::adjust_size( x , m_x_err );
- traits_type::adjust_size( x , m_x_scale );
-
m_error_checker.fill_scale( x , dxdt , dt , m_x_scale );
m_x_tmp = x;
@@ -137,7 +163,7 @@
{
// error too large - decrease dt
// limit scaling factor to 0.2
- dt *= std::max( 0.9*pow(max_rel_err , -1.0/(m_stepper.order_error()-1.0)),
+ dt *= std::max( 0.9 * pow(max_rel_err , -1.0/(m_stepper.order_error()-1.0)),
0.2 );
// reset state
@@ -151,7 +177,7 @@
//error too small - increase dt and keep the evolution
t += dt;
// limit scaling factor to 5.0
- dt *= std::min( 0.9*pow(max_rel_err , -1.0/m_stepper.order()), 5.0 );
+ dt *= std::min( 0.9*pow(max_rel_err , -1.0/m_stepper.order_error_step()), 5.0 );
return step_size_increased;
}
else
@@ -169,23 +195,12 @@
time_type &t,
time_type &dt )
{
- traits_type::adjust_size( x , m_dxdt );
system( x , m_dxdt , t );
return try_step( system , x , m_dxdt , t , dt );
}
};
- template< class ErrorStepper >
- controlled_stepper_standard< ErrorStepper > make_controlled_stepper_standard(
- ErrorStepper &stepper,
- typename ErrorStepper::time_type abs_err, typename ErrorStepper::time_type rel_err,
- typename ErrorStepper::time_type factor_x, typename ErrorStepper::time_type factor_dxdt )
- {
- controlled_stepper_standard< ErrorStepper > controlled_stepper(
- stepper , abs_err , rel_err , factor_x , factor_dxdt );
- return controlled_stepper;
- };
Modified: sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -26,50 +26,58 @@
template< class Container,
class Time,
class Traits = container_traits< Container > >
- class error_checker_standard {
+ class error_checker_standard
+ {
public:
+
typedef Time time_type;
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
+// typedef typename traits_type::iterator iterator;
+// typedef typename traits_type::const_iterator const_iterator;
private:
+
time_type m_eps_abs;
time_type m_eps_rel;
time_type m_a_x;
time_type m_a_dxdt;
public:
+
// constructor
error_checker_standard(
time_type abs_err, time_type rel_err,
time_type factor_x, time_type factor_dxdt )
: m_eps_abs( abs_err ), m_eps_rel( rel_err ),
m_a_x( factor_x ), m_a_dxdt( factor_dxdt )
- { }
+ {
+ }
+
void fill_scale(
- container_type &x,
- container_type &dxdt,
+ const container_type &x,
+ const container_type &dxdt,
time_type dt,
- container_type &scale )
+ container_type &scale ) const
{
detail::it_algebra::weighted_scale( traits_type::begin(scale),
traits_type::end(scale),
traits_type::begin(x),
- traits_type::end(dxdt),
+ traits_type::begin(dxdt),
m_eps_abs,
m_eps_rel,
m_a_x,
m_a_x*dt );
}
- time_type get_max_error_ratio( container_type &x_err, container_type &scale)
+ time_type get_max_error_ratio(
+ const container_type &x_err,
+ const container_type &scale) const
{
return detail::it_algebra::max_ratio( traits_type::begin(x_err),
traits_type::end(x_err),
@@ -78,8 +86,6 @@
}
const time_type get_epsilon() { return std::max(m_eps_rel, m_eps_abs); }
-
-
};
}
Modified: sandbox/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -37,10 +37,13 @@
typename ControlledStepper::time_type dt,
Observer &observer )
{
+ typedef typename ControlledStepper::time_type time_type;
+
+ stepper.adjust_size( state );
+
controlled_step_result result;
size_t iterations = 0;
- typename ControlledStepper::time_type t = start_time;
- typename ControlledStepper::time_type dt_ = dt;
+ time_type t = start_time , dt_ = dt;
observer(t, state, system);
@@ -150,12 +153,12 @@
T a_dxdt = 1.0
)
{
- typedef stepper_rk5_ck< ContainerType , T > stepper_type;
// we use cash karp stepper as base stepper
- stepper_type stepper_cash_karp;
+ typedef stepper_rk5_ck< ContainerType , T > stepper_type;
+
// we use the standard controller for this adaptive integrator
- controlled_stepper_standard< stepper_type >
- controlled_stepper(stepper_cash_karp, eps_abs, eps_rel, a_x, a_dxdt );
+ controlled_stepper_standard< stepper_type >
+ controlled_stepper( eps_abs, eps_rel, a_x, a_dxdt );
// initialized with values from above
// call the normal integrator
Modified: sandbox/odeint/boost/numeric/odeint/integrator_constant_stepsize.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/integrator_constant_stepsize.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/integrator_constant_stepsize.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -35,6 +35,7 @@
Observer &observer
)
{
+ stepper.adjust_size( state );
size_t iteration = 0;
while( start_time < end_time )
{
@@ -89,6 +90,7 @@
Observer &observer
)
{
+ stepper.adjust_size( state );
size_t iteration = 0;
while( iteration < num_of_steps )
{
Modified: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -132,7 +132,7 @@
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
- container_type &dxdt ,
+ const container_type &dxdt ,
time_type t ,
time_type dt ,
container_type &xerr )
Modified: sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -106,7 +106,7 @@
void midpoint_step(
DynamicalSystem &system ,
container_type &x ,
- container_type &dxdt ,
+ const container_type &dxdt ,
time_type t ,
time_type dt ,
container_type &x_out )
@@ -161,7 +161,7 @@
void do_step(
DynamicalSystem &system ,
container_type &x ,
- container_type &dxdt ,
+ const container_type &dxdt ,
time_type t ,
time_type dt )
{
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -75,6 +75,7 @@
// constructor, which adjusts the internal containers
stepper_rk4( const container_type &x )
{
+ adjust_size( x );
}
void adjust_size( const container_type &x )
@@ -91,7 +92,7 @@
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
- container_type &dxdt ,
+ const container_type &dxdt ,
time_type t ,
time_type dt )
{
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -86,7 +86,7 @@
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
- container_type &dxdt ,
+ const container_type &dxdt ,
time_type t ,
time_type dt )
{
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
==============================================================================
Binary files. No diff available.
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -99,7 +99,7 @@
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
- container_type &dxdt ,
+ const container_type &dxdt ,
time_type t ,
time_type dt )
{
@@ -359,8 +359,6 @@
}
-/*
-
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
container_type &x ,
@@ -372,26 +370,29 @@
const time_type cc01 = static_cast<time_type>( 41.0 / 840.0 );
const time_type cc06 = static_cast<time_type>( 34.0 / 105.0 );
const time_type cc07 = static_cast<time_type>( 9.0 / 35.0 );
- const time_type cc08 = static_cast<time_type>( cc08 );
+ const time_type cc08 = static_cast<time_type>( cc07 );
const time_type cc09 = static_cast<time_type>( 9.0 / 280.0 );
const time_type cc10 = static_cast<time_type>( cc09 );
const time_type cc11 = static_cast<time_type>( cc01 );
+ using namespace detail::it_algebra;
+
xerr = x;
do_step( system , xerr , dxdt , t , dt );
// now, k1-k13 are calculated and stored in m_k01 - m_k13
- scale_sum( x.begin() , x.end() ,
- static_cast<time_type>(1.0) , x.begin(),
- dt * cc01 , dxdt ,
- dt * cc06 , m_k06 ,
- dt * cc07 , m_k07 ,
- dt * cc08 , m_k08 ,
- dt * cc09 , m_k09 ,
- dt * cc10 , m_k10 ,
- dt * cc11 , m_k11 );
+ scale_sum( traits_type::begin(x) , traits_type::end(x) ,
+ static_cast<time_type>(1.0) , traits_type::begin(x) ,
+ dt * cc01 , traits_type::begin(dxdt) ,
+ dt * cc06 , traits_type::begin(m_k06) ,
+ dt * cc07 , traits_type::begin(m_k07) ,
+ dt * cc08 , traits_type::begin(m_k08) ,
+ dt * cc09 , traits_type::begin(m_k09) ,
+ dt * cc10 , traits_type::begin(m_k10) ,
+ dt * cc11 , traits_type::begin(m_k11) );
- increment( xerr.begin() , xerr.end() , x.begin() , static_cast<time_type>(-1.0) );
+ increment( traits_type::begin(xerr) , traits_type::end(xerr) ,
+ traits_type::begin(x) , static_cast<time_type>(-1.0) );
}
template< class DynamicalSystem >
@@ -401,11 +402,9 @@
time_type dt ,
container_type &xerr )
{
- traits_type::adjust_size( x , m_dxdt );
system( x , m_dxdt , t );
do_step( system , x , m_dxdt , t , dt , xerr );
}
-*/
};
} // namespace odeint
Modified: sandbox/odeint/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/Jamfile 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -15,7 +15,6 @@
;
exe harmonic_oscillator : harmonic_oscillator.cpp ;
-# exe test_container_and_stepper : test_container_and_stepper.cpp ;
exe hamiltonian_test : hamiltonian_test.cpp ;
exe lorenz_cmp_rk4_rk_generic : lorenz_cmp_rk4_rk_generic.cpp ;
Modified: sandbox/odeint/libs/numeric/odeint/examples/doc_integrate.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/doc_integrate.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/doc_integrate.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -47,11 +47,11 @@
//[ define_conntrolled_stepper
controlled_stepper_standard< stepper_rk5_ck< state_type > >
- controlled_rk5( rk5 , 1E-6 , 1E-7 , 1.0 , 1.0 );
+ controlled_rk5( 1E-6 , 1E-7 , 1.0 , 1.0 );
//]
//[ integrate_adapt
- integrate_adaptive( make_controlled_stepper_standard( rk5 , 1E-6 , 1E-7 , 1.0 , 1.0 ),
+ integrate_adaptive( controlled_rk5 ,
harmonic_oscillator, x, 0.0, 10.0, 0.01 );
//]
Modified: sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -127,6 +127,9 @@
odeint::stepper_euler<harm_osc::state_type> euler;
odeint::stepper_rk4<harm_osc::state_type> rk4;
+ euler.adjust_size( x_euler );
+ rk4.adjust_size( x_rk4 );
+
// Write the output header. The '#' symbol creates a comment in gnuplot
cout << "# Output from the simple harmonic oscillator example." << endl;
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -24,7 +24,8 @@
#include <tr1/array>
#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/controlled_stepper_bs.hpp>
+#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
+ // #include <boost/numeric/odeint/controlled_stepper_bs.hpp>
#define tab "\t"
@@ -86,7 +87,7 @@
x[2] = 20.0;
stepper_rk5_ck< state_type > rk5;
- controlled_stepper_standard< stepper_rk5_ck< state_type > > controlled_rk5( rk5, eps_abs , eps_rel, 1.0, 1.0 );
+ controlled_stepper_standard< stepper_rk5_ck< state_type > > controlled_rk5( eps_abs , eps_rel, 1.0, 1.0 );
output_observer rk5_obs("lorenz_rk5.dat");
size_t steps = integrate_adaptive( controlled_rk5, lorenz, x, 0.0, end_time, 1E-2, rk5_obs );
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -68,7 +68,7 @@
stepper_half_step< stepper_euler< state_type > > euler;
controlled_stepper_standard< stepper_half_step< stepper_euler< state_type > > >
- euler_controlled( euler , eps_abs, eps_rel, 1.0, 1.0);
+ euler_controlled( eps_abs, eps_rel, 1.0, 1.0);
size_t steps = integrate( euler_controlled, lorenz, x1,
0.0, 10.0, 1E-4,
back_inserter(t1_vec),
@@ -78,7 +78,7 @@
stepper_half_step< stepper_rk4< state_type > > rk4;
controlled_stepper_standard< stepper_half_step< stepper_rk4< state_type > > >
- rk4_controlled( rk4 , eps_abs, eps_rel, 1.0, 1.0);
+ rk4_controlled( eps_abs, eps_rel, 1.0, 1.0);
steps = integrate( rk4_controlled, lorenz, x2, 0.0, 10.0, 1E-4,
back_inserter(t2_vec),
back_inserter(x2_t_vec));
@@ -88,7 +88,7 @@
stepper_rk5_ck< state_type > rk5;
controlled_stepper_standard< stepper_rk5_ck< state_type > >
- rk5_controlled( rk5 , eps_abs, eps_rel, 1.0, 1.0);
+ rk5_controlled( eps_abs, eps_rel, 1.0, 1.0);
steps = integrate( rk5_controlled, lorenz, x3, 0.0, 10.0, 1E-4,
back_inserter(t3_vec),
back_inserter(x3_t_vec));
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -69,8 +69,8 @@
x1[2] = 0.0;
state_type2 x2 = {{ 1.0 , 0.0 , 0.0 }};
- stepper_rk4< state_type1 > stepper1;
- stepper_rk4< state_type2 > stepper2;
+ stepper_rk4< state_type1 > stepper1( x1 );
+ stepper_rk4< state_type2 > stepper2( x2 );
clock_t start , end;
double t;
@@ -78,8 +78,7 @@
start= clock();
t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- stepper1.do_step
-( lorenz1 , x1 , t , dt );
+ stepper1.do_step( lorenz1 , x1 , t , dt );
end = clock();
cout << "vector : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
cout << "x: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
Modified: sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -56,7 +56,7 @@
stepper_half_step< stepper_rk4< state_type > > stepper;
controlled_stepper_standard< stepper_half_step< stepper_rk4< state_type > > >
- controlled_stepper( stepper, 1E-6 , 1E-7 , 1.0 , 1.0 );
+ controlled_stepper( 1E-6 , 1E-7 , 1.0 , 1.0 );
size_t steps = integrate( controlled_stepper, my_system, x,
0.0, 100.0,1E-4,
Deleted: sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
+++ (empty file)
@@ -1,268 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_stepper.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- A simple compile test program
-
- Furthermore, the usage of std::tr1::array and std::vector in odeint is
- shown and the performance of both containers is compared.
-
- 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 <vector>
-#include <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/container_traits_blitz_array.hpp>
-#include <boost/numeric/odeint/container_traits_mtl4_dense2d.hpp>
-#include <boost/numeric/odeint/container_traits_ublas_matrix.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::numeric::odeint;
-
-const size_t n1 = 3;
-const size_t n2 = 5;
-const double dt = 0.01;
-
-typedef std::vector< double > state_type1;
-typedef std::tr1::array< double , n1 > state_type2;
-typedef blitz::Array< double , 2 > state_type3;
-typedef mtl::dense2D< double > state_type4;
-typedef boost::numeric::ublas::matrix< double > state_type5;
-
-void deriv1( state_type1 &x , state_type1 &dxdt , double t )
-{
- const double sigma = 10.0 , R = 28.0 , b = 8.0 / 3.0;
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-void deriv2( state_type2 &x , state_type2 &dxdt , double t )
-{
- const double sigma = 10.0 , R = 28.0 , b = 8.0 / 3.0;
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-void deriv3( state_type3 &x , state_type3 &dxdt , double t )
-{
- for( int i=0 ; i<int(n1) ; ++i )
- for( int j=0 ; j<int(n2) ; ++j )
- dxdt( i , j ) = drand48();
-}
-
-void deriv4( state_type4 &x , state_type4 &dxdt , double t )
-{
- for( size_t i=0 ; i<n1 ; ++i )
- for( size_t j=0 ; j<n2 ; ++j )
- dxdt( i , j ) = drand48();
-}
-
-void deriv5( state_type5 &x , state_type5 &dxdt , double t )
-{
- for( size_t i=0 ; i<n1 ; ++i )
- for( size_t j=0 ; j<n2 ; ++j )
- dxdt( i , j ) = drand48();
-}
-
-template<
- class Stepper1 ,
- class Stepper2 ,
- class Stepper3 ,
- class Stepper4 ,
- class Stepper5
- >
-void test_steppers(
- Stepper1 &stepper1 ,
- Stepper2 &stepper2 ,
- Stepper3 &stepper3 ,
- Stepper4 &stepper4 ,
- Stepper5 &stepper5
- )
-{
- state_type1 x1( n1 , 0.0 ) , dxdt1( x1 );
- state_type2 x2 = {{ 1.0 , 0.0 , 0.0 }} , dxdt2( x2 );
- state_type3 x3( n1 , n2 ) , dxdt3( x3 );
- state_type4 x4( n1 , n2 ) , dxdt4( x4 );
- state_type5 x5( n1 , n2 ) , dxdt5( x5 );
-
- stepper1.do_step( deriv1 , x1 , 0.0 , dt );
- stepper2.do_step( deriv2 , x2 , 0.0 , dt );
- stepper3.do_step( deriv3 , x3 , 0.0 , dt );
- stepper4.do_step( deriv4 , x4 , 0.0 , dt );
- stepper5.do_step( deriv5 , x5 , 0.0 , dt );
-
- stepper1.do_step( deriv1 , x1 , dxdt1 , 0.0 , dt );
- stepper2.do_step( deriv2 , x2 , dxdt2 , 0.0 , dt );
- stepper3.do_step( deriv3 , x3 , dxdt3 , 0.0 , dt );
- stepper4.do_step( deriv4 , x4 , dxdt4 , 0.0 , dt );
- stepper5.do_step( deriv5 , x5 , dxdt5 , 0.0 , dt );
-}
-
-
-
-template<
- class Stepper1 ,
- class Stepper2 ,
- class Stepper3 ,
- class Stepper4 ,
- class Stepper5
- >
-void test_error_steppers(
- Stepper1 &stepper1 ,
- Stepper2 &stepper2 ,
- Stepper3 &stepper3 ,
- Stepper4 &stepper4 ,
- Stepper5 &stepper5
- )
-{
- state_type1 x1( n1 , 0.0 ) , dxdt1( x1 ) , error1( x1 );
- state_type2 x2 = {{ 1.0 , 0.0 , 0.0 }} , dxdt2( x2 ) , error2( x2 );
- state_type3 x3( n1 , n2 ) , dxdt3( x3 ) , error3( x3 );
- state_type4 x4( n1 , n2 ) , dxdt4( x4 ) , error4( x4 );
- state_type5 x5( n1 , n2 ) , dxdt5( x5 ) , error5( x5 );
-
- stepper1.do_step( deriv1 , x1 , 0.0 , dt , error1 );
- stepper2.do_step( deriv2 , x2 , 0.0 , dt , error2 );
- stepper3.do_step( deriv3 , x3 , 0.0 , dt , error3 );
- stepper4.do_step( deriv4 , x4 , 0.0 , dt , error4 );
- stepper5.do_step( deriv5 , x5 , 0.0 , dt , error5 );
-
- stepper1.do_step( deriv1 , x1 , dxdt1 , 0.0 , dt , error1 );
- stepper2.do_step( deriv2 , x2 , dxdt2 , 0.0 , dt , error2 );
- stepper3.do_step( deriv3 , x3 , dxdt3 , 0.0 , dt , error3 );
- stepper4.do_step( deriv4 , x4 , dxdt4 , 0.0 , dt , error4 );
- stepper5.do_step( deriv5 , x5 , dxdt5 , 0.0 , dt , error5 );
-}
-
-template< class Stepper , class Deriv , class Container >
-void test_stepper( Stepper &stepper , Deriv &deriv , Container &x , Container &dxdt )
-{
- typename Stepper::container_type xx;
- typename Stepper::order_type order;
- typename Stepper::time_type time;
- typename Stepper::traits_type traits;
- typename Stepper::value_type value;
- typename Stepper::iterator iterator;
- typename Stepper::const_iterator const_iteratio;
-
- stepper.do_step( deriv , x , 0.0 , dt );
- stepper.do_step( deriv , x , dxdt , 0.0 , dt );
-}
-
-
-
-int main( int argc , char **argv )
-{
- cout << "Hello world!" << endl;
-
- state_type1 x1( n1 , 0.0 ) , dxdt1( x1 ) , error1( x1 );
- state_type2 x2 = {{ 1.0 , 0.0 , 0.0 }} , dxdt2( x2 ) , error2( x2 );
- state_type3 x3( n1 , n2 ) , dxdt3( x3 ) , error3( x3 );
- state_type4 x4( n1 , n2 ) , dxdt4( x4 ) , error4( x4 );
- state_type5 x5( n1 , n2 ) , dxdt5( x5 ) , error5( x5 );
-
- stepper_euler< state_type1 > s1;
- stepper_euler< state_type2 > s2;
- stepper_euler< state_type3 > s3;
- stepper_euler< state_type4 > s4;
- stepper_euler< state_type5 > s5;
-
- test_stepper( s1 , deriv1 , x1 , dxdt1 );
-
-
- {
- stepper_euler< state_type1 > s1;
- stepper_euler< state_type2 > s2;
- stepper_euler< state_type3 > s3;
- stepper_euler< state_type4 > s4;
- stepper_euler< state_type5 > s5;
- test_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "Euler Stepper OK!" << endl;
-
-
- {
- stepper_half_step< stepper_euler< state_type1 > > s1;
- stepper_half_step< stepper_euler< state_type2 > > s2;
- stepper_half_step< stepper_euler< state_type3 > > s3;
- stepper_half_step< stepper_euler< state_type4 > > s4;
- stepper_half_step< stepper_euler< state_type5 > > s5;
- test_steppers( s1 , s2 , s3 , s4 , s5 );
- test_error_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "Half Step Stepper OK!" << endl;
-
-
- {
- stepper_rk4< state_type1 > s1;
- stepper_rk4< state_type2 > s2;
- stepper_rk4< state_type3 > s3;
- stepper_rk4< state_type4 > s4;
- stepper_rk4< state_type5 > s5;
- test_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "RK4 Stepper OK!" << endl;
-
-
-
-
- {
- stepper_rk4_classical< state_type1 > s1;
- stepper_rk4_classical< state_type2 > s2;
- stepper_rk4_classical< state_type3 > s3;
- stepper_rk4_classical< state_type4 > s4;
- stepper_rk4_classical< state_type5 > s5;
- test_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "RK4 Classical Stepper OK!" << endl;
-
-
-
- {
- stepper_rk5_ck< state_type1 > s1;
- stepper_rk5_ck< state_type2 > s2;
- stepper_rk5_ck< state_type3 > s3;
- stepper_rk5_ck< state_type4 > s4;
- stepper_rk5_ck< state_type5 > s5;
- test_error_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "RK5 CK Stepper OK!" << endl;
-
-
- {
- stepper_midpoint< state_type1 > s1;
- stepper_midpoint< state_type2 > s2;
- stepper_midpoint< state_type3 > s3;
- stepper_midpoint< state_type4 > s4;
- stepper_midpoint< state_type5 > s5;
- test_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "Midpoint Stepper OK!" << endl;
-
-
- {
- stepper_rk78_fehlberg< state_type1 > s1;
- stepper_rk78_fehlberg< state_type2 > s2;
- stepper_rk78_fehlberg< state_type3 > s3;
- stepper_rk78_fehlberg< state_type4 > s4;
- stepper_rk78_fehlberg< state_type5 > s5;
- test_steppers( s1 , s2 , s3 , s4 , s5 );
- }
- cout << "RK78 Stepper OK!" << endl;
-
-
- cout << "Everything is allright!" << endl;
-
- return 0;
-}
Modified: sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/test/check_stepper_concepts.cpp 2010-03-09 17:36:13 EST (Tue, 09 Mar 2010)
@@ -156,7 +156,7 @@
{
stepper_rk78_fehlberg< std::vector<double> > stepper;
check_stepper_concept( stepper , 8 );
-// check_error_stepper_concept( stepper , 7 , 8 );
+ check_error_stepper_concept( stepper , 7 , 8 );
}