Subject: [Boost-commit] svn:boost r57733 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-17 17:12:03
Author: mariomulansky
Date: 2009-11-17 17:12:02 EST (Tue, 17 Nov 2009)
New Revision: 57733
URL: http://svn.boost.org/trac/boost/changeset/57733
Log:
typecast with static_const,
typedefs in generic rk
Binary files modified:
sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
Text files modified:
sandbox/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp | 121 ++++++++++++++++++---------------------
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 11 +-
sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp | 21 +++---
sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp | 33 +++++-----
sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp | 25 +++++--
5 files changed, 106 insertions(+), 105 deletions(-)
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 2009-11-17 17:12:02 EST (Tue, 17 Nov 2009)
@@ -1,4 +1,4 @@
-/* Boost odeint/integrator.hpp header file
+/* Boost odeint/integrator_adaptive.hpp header file
Copyright 2009 Karsten Ahnert
Copyright 2009 Mario Mulansky
@@ -25,20 +25,20 @@
template<
- class Stepper,
- class DynamicalSystem,
- class StepController,
- class Observer
- >
+ class Stepper,
+ class DynamicalSystem,
+ class StepController,
+ class Observer
+ >
size_t integrate_adaptive(
- Stepper &stepper,
- DynamicalSystem &system,
- StepController &controller,
- typename Stepper::time_type start_time,
- typename Stepper::time_type dt,
- typename Stepper::container_type &state,
- typename Stepper::time_type end_time,
- Observer &observer )
+ Stepper &stepper,
+ DynamicalSystem &system,
+ StepController &controller,
+ typename Stepper::time_type start_time,
+ typename Stepper::time_type dt,
+ typename Stepper::container_type &state,
+ typename Stepper::time_type end_time,
+ Observer &observer )
{
controlled_step_result result;
size_t iterations = 0;
@@ -58,18 +58,8 @@
iterations++;
}
-/* while( result != SUCCESS ) // as long as dt is too large/small
- {
- // do the controlled step
- result = controller.controlled_step( stepper, system, state, t, dt_ );
- if( result != STEP_SIZE_DECREASED )
- { // we did a step
- observer(t, state, system);
- iterations++;
- }*/
- if( !( t+dt_ > t) )
- throw; // we've reached machine precision with dt - no advancing in t
- //}
+ if( !( t+dt_ > t) )
+ throw; // we've reached machine precision with dt - no advancing in t
}
return iterations;
}
@@ -80,16 +70,15 @@
class StepController
>
size_t integrate_adaptive(
- Stepper &stepper,
- DynamicalSystem &system,
- StepController &controller,
- typename Stepper::time_type start_time,
- typename Stepper::time_type dt,
- typename Stepper::container_type &state,
- typename Stepper::time_type end_time
- )
+ Stepper &stepper,
+ DynamicalSystem &system,
+ StepController &controller,
+ typename Stepper::time_type start_time,
+ typename Stepper::time_type dt,
+ typename Stepper::container_type &state,
+ typename Stepper::time_type end_time )
{
- return integrate_adaptive(
+ return integrate_adaptive(
stepper , system , controller ,
start_time , dt , state , end_time ,
do_nothing_observer<
@@ -122,21 +111,22 @@
class TimeSequence,
class InsertIterator
>
- size_t integrate(Stepper &stepper,
- DynamicalSystem &system,
- StepController &controller,
- typename Stepper::container_type &state,
- TimeSequence ×,
- InsertIterator state_inserter,
- typename Stepper::time_type &dt)
+ size_t integrate(
+ Stepper &stepper,
+ DynamicalSystem &system,
+ StepController &controller,
+ typename Stepper::container_type &state,
+ TimeSequence ×,
+ InsertIterator state_inserter,
+ typename Stepper::time_type &dt)
{
- if( times.empty() ) return 0;
- else
- {
- state_copy_observer<InsertIterator, TimeSequence> observer(times, state_inserter);
- return integrate_adaptive(stepper, system, controller, times.front() ,
- dt, state, times.back() , observer);
- }
+ if( times.empty() ) return 0;
+ else
+ {
+ state_copy_observer<InsertIterator, TimeSequence> observer(times, state_inserter);
+ return integrate_adaptive(stepper, system, controller, times.front() ,
+ dt, state, times.back() , observer);
+ }
}
@@ -172,21 +162,24 @@
To avoid extensive chages in dt, the decrease factor is limited to 0.2 and
the increase factor to 5.0.
*/
- template< class Stepper,
- class DynamicalSystem,
- class InsertIterator ,
- class TimeSequence
- >
- size_t integrate(Stepper &stepper,
- DynamicalSystem &system,
- typename Stepper::container_type &x,
- TimeSequence ×,
- InsertIterator state_inserter,
- typename Stepper::time_type dt = 1E-4,
- typename Stepper::time_type eps_abs = 1E-6,
- typename Stepper::time_type eps_rel = 1E-7,
- typename Stepper::time_type a_x = 1.0 ,
- typename Stepper::time_type a_dxdt = 1.0)
+ template<
+ class Stepper,
+ class DynamicalSystem,
+ class InsertIterator ,
+ class TimeSequence
+ >
+ size_t integrate(
+ Stepper &stepper,
+ DynamicalSystem &system,
+ typename Stepper::container_type &x,
+ TimeSequence ×,
+ InsertIterator state_inserter,
+ typename Stepper::time_type dt = 1E-4,
+ typename Stepper::time_type eps_abs = 1E-6,
+ typename Stepper::time_type eps_rel = 1E-7,
+ typename Stepper::time_type a_x = 1.0 ,
+ typename Stepper::time_type a_dxdt = 1.0
+ )
{
// we use the standard controller for this adaptive integrator
step_controller_standard< typename Stepper::container_type,
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 2009-11-17 17:12:02 EST (Tue, 17 Nov 2009)
@@ -99,17 +99,16 @@
m_resizer.adjust_size( x , xerr );
m_xtemp = x;
- time_type dt2 = 0.5 * dt;
+ time_type dt2 = static_cast<time_type>(0.5) * dt;
next_step( system , m_xtemp , dxdt , t , dt );
next_step( system , x , dxdt , t , dt2 );
next_step( system , x , t+dt2 , dt2 );
- detail::it_algebra::assign_diff(
- xerr.begin() ,
- xerr.end() ,
- m_xtemp.begin() ,
- x.begin() );
+ detail::it_algebra::assign_diff( xerr.begin() ,
+ xerr.end() ,
+ m_xtemp.begin() ,
+ x.begin() );
}
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 2009-11-17 17:12:02 EST (Tue, 17 Nov 2009)
@@ -4,8 +4,8 @@
Copyright 2009 Mario Mulansky
Copyright 2009 Andre Bergner
- This file includes the explicit runge kutta solver for
- ordinary differential equations.
+ This file includes the explicit 4th order runge kutta
+ solver for ordinary differential equations.
It solves any ODE dx/dt = f(x,t).
@@ -75,10 +75,9 @@
// public interface
public:
- stepper_rk4( void )
- : current_size(0)
+ stepper_rk4( void )
{
- }
+ }
order_type order() const { return 4; }
@@ -91,14 +90,14 @@
{
using namespace detail::it_algebra;
- const time_type val1 = time_type( 1.0 );
+ const time_type val1 = static_cast<time_type>( 1.0 );
m_resizer.adjust_size( x , m_dxt );
m_resizer.adjust_size( x , m_dxm );
m_resizer.adjust_size( x , m_xt );
m_resizer.adjust_size( x , m_dxh );
- time_type dh = time_type( 0.5 ) * dt;
+ time_type dh = static_cast<time_type>( 0.5 ) * dt;
time_type th = t + dh;
// dt * dxdt = k1
@@ -126,10 +125,10 @@
//x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
scale_sum( x.begin(), x.end(),
val1, x.begin(),
- dt / time_type( 6.0 ), dxdt.begin(),
- dt / time_type( 3.0 ), m_dxt.begin(),
- dt / time_type( 3.0 ), m_dxm.begin(),
- dt / time_type( 6.0 ), m_dxh.begin() );
+ dt / static_cast<time_type>( 6.0 ), dxdt.begin(),
+ dt / static_cast<time_type>( 3.0 ), m_dxt.begin(),
+ dt / static_cast<time_type>( 3.0 ), m_dxm.begin(),
+ dt / static_cast<time_type>( 6.0 ), m_dxh.begin() );
}
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
==============================================================================
Binary files. No diff available.
Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp 2009-11-17 17:12:02 EST (Tue, 17 Nov 2009)
@@ -26,7 +26,6 @@
#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/resizer.hpp>
-#include <iostream>
namespace boost {
namespace numeric {
@@ -85,13 +84,15 @@
typedef std::vector< container_type > container_vector;
typedef std::vector< iterator > container_iterator_vector;
+ typedef std::vector< time_type > parameter_vector;
+ typedef std::vector< parameter_vector > parameter_matrix;
container_vector m_xvec;
container_iterator_vector m_xiter_vec;
container_type m_xtmp;
- const std::vector< time_type > m_a;
- const std::vector< std::vector<time_type> > m_b;
- const std::vector< time_type > m_c;
+ const parameter_vector m_a;
+ const parameter_matrix m_b;
+ const parameter_vector m_c;
order_type m_q;
@@ -111,10 +112,10 @@
void check_consitency()
{
- typename std::vector< time_type >::const_iterator a_iter = m_a.begin();
- typename std::vector< time_type >::const_iterator c_iter = m_c.begin();
- typename std::vector< std::vector<time_type> >::const_iterator b_iter = m_b.begin();
- typename std::vector<time_type>::const_iterator b_iter_iter;
+ typename parameter_vector::const_iterator a_iter = m_a.begin();
+ typename parameter_vector::const_iterator c_iter = m_c.begin();
+ typename parameter_matrix::const_iterator b_iter = m_b.begin();
+ typename parameter_vector::const_iterator b_iter_iter;
// check 1: a_i = sum_j b_ij
while( a_iter != m_a.end() ) {
@@ -184,9 +185,9 @@
Note, that a_0 = 1 (implicitely) and 0^0 = 1
so this means sum_i c_i = 1 at k=1
*/
- stepper_rk_generic( std::vector<time_type> &a,
- std::vector< std::vector<time_type> > &b,
- std::vector< time_type > &c)
+ stepper_rk_generic( parameter_vector &a,
+ parameter_matrix &b,
+ parameter_vector &c)
: m_a(a), m_b(b), m_c(c), m_q(c.size())
{
m_xvec.resize(m_q);
@@ -213,7 +214,7 @@
(*xiter_iter++) = (*x_iter++).begin();
while( x_iter != m_xvec.end() )
- {
+ {
m_resizer.adjust_size(x, (*x_iter));
(*xiter_iter++) = (*x_iter++).begin();
}
@@ -221,10 +222,10 @@
x_iter = m_xvec.begin()+1;
- typename std::vector< time_type >::const_iterator a_iter = m_a.begin();
- typename std::vector< std::vector<time_type> >::const_iterator b_iter = m_b.begin();
+ typename parameter_vector::const_iterator a_iter = m_a.begin();
+ typename parameter_matrix::const_iterator b_iter = m_b.begin();
while( x_iter != m_xvec.end() )
- {
+ {
reset_iter(m_xiter_vec.begin());
scale_sum_generic( m_xtmp.begin(), m_xtmp.end(),
(*b_iter).begin(), (*b_iter).end(), dt,
@@ -234,7 +235,7 @@
}
reset_iter(m_xiter_vec.begin());
- typename std::vector< time_type >::const_iterator c_iter = m_c.begin();
+ typename parameter_vector::const_iterator c_iter = m_c.begin();
scale_sum_generic( x.begin(), x.end(),
m_c.begin(), m_c.end(), dt,
x.begin(), m_xiter_vec.begin() );
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 2009-11-17 17:12:02 EST (Tue, 17 Nov 2009)
@@ -51,21 +51,29 @@
int main( int argc , char **argv )
{
- state_type x;
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 20.0;
+ state_type x1;
+ x1[0] = 1.0;
+ x1[1] = 0.0;
+ x1[2] = 20.0;
+ state_type x2(x1);
- vector<state_type> x_t_vec;
+
+ vector<state_type> x1_t_vec;
+ vector<state_type> x2_t_vec;
vector<double> times(time_points);
for( size_t i=0; i<time_points; i++ ) {
times[i] = 0.1*i;
}
stepper_half_step< stepper_euler< state_type > > euler;
- size_t steps = integrate( euler, lorenz, x, times, back_inserter(x_t_vec));
+ size_t steps = integrate( euler, lorenz, x1, times, back_inserter(x1_t_vec));
+
+ clog << "Euler Steps: " << steps << endl;
- clog << "Steps: " << steps << endl;
+ stepper_rk5_ck< state_type > rk5;
+ steps = integrate( rk5, lorenz, x2, times, back_inserter(x2_t_vec));
+
+ clog << "RK5 Steps: " << steps << endl;
cout.precision(5);
cout.setf(ios::fixed,ios::floatfield);
@@ -74,7 +82,8 @@
for( size_t i=0; i<time_points; i++ ) {
//cout << "current state: " ;
cout << times[i] << tab;
- cout << x_t_vec[i][0] << tab << x_t_vec[i][1] << tab << x_t_vec[i][2] << endl;
+ cout << x1_t_vec[i][0] << tab << x1_t_vec[i][1] << tab << x1_t_vec[i][2] << tab;
+ cout << x2_t_vec[i][0] << tab << x2_t_vec[i][1] << tab << x2_t_vec[i][2] << endl;
}
return 0;