Subject: [Boost-commit] svn:boost r58036 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/examples libs/numeric/odeint/stuff/test
From: mario.mulansky_at_[hidden]
Date: 2009-11-29 15:36:33
Author: mariomulansky
Date: 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
New Revision: 58036
URL: http://svn.boost.org/trac/boost/changeset/58036
Log:
bulirsch stoer stepper with example - yeah
Added:
sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp (contents, props changed)
sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp (contents, props changed)
sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_bs.cpp (contents, props changed)
Text files modified:
sandbox/odeint/boost/numeric/odeint.hpp | 3 +
sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp | 14 +++++--
sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp | 3 -
sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp | 68 ++++++++++++++++++++++++++++++---------
sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp | 13 ++++--
5 files changed, 73 insertions(+), 28 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -25,6 +25,9 @@
#include <boost/numeric/odeint/stepper_rk78_fehlberg.hpp>
#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
+#include <boost/numeric/odeint/controlled_stepper_bs.hpp>
+
+#include <boost/numeric/odeint/error_checker_standard.hpp>
#include <boost/numeric/odeint/integrator_adaptive_stepsize.hpp>
#include <boost/numeric/odeint/integrator_constant_stepsize.hpp>
Added: sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -0,0 +1,371 @@
+/* Boost odeint/controlled_stepper_bs.hpp header file
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ This file includes the explicit Burlisch Stoer
+ solver for ordinary differential equations.
+
+ It solves any ODE dx/dt = f(x,t) via
+ x(t+dt) = x(t) + dt*f(x,t)
+
+ 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_NUMERIC_ODEINT_CONTROLLED_STEPPER_BS_HPP
+#define BOOST_NUMERIC_ODEINT_CONTROLLED_STEPPER_BS_HPP
+
+#include <limits>
+
+#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/resizer.hpp>
+#include <boost/numeric/odeint/stepper_midpoint.hpp>
+#include <boost/numeric/odeint/error_checker_standard.hpp>
+
+#include <iostream>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+ // typedef enum{success, step_size_decreased, step_size_increased} controlled_step_result;
+
+ template<
+ class Container ,
+ class Time = double ,
+ class Resizer = resizer< Container >
+ >
+ class controlled_stepper_bs
+ {
+ // provide basic typedefs
+ public:
+
+ typedef Container container_type;
+ typedef Resizer resizer_type;
+ typedef Time time_type;
+ typedef const unsigned short order_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::iterator iterator;
+
+
+ // check the concept of the ContainerType
+ private:
+
+ BOOST_CLASS_REQUIRE( container_type ,
+ boost::numeric::odeint, Container );
+
+
+ // private memebers
+ private:
+ stepper_midpoint< Container, Time, Resizer > m_stepper_mp;
+ resizer_type m_resizer;
+ error_checker_standard<container_type, time_type> m_error_checker;
+
+ const unsigned short m_k_max;
+
+ const time_type m_safety1;
+ const time_type m_safety2;
+ const time_type m_max_dt_factor;
+ const time_type m_min_step_scale;
+ const time_type m_max_step_scale;
+
+ bool m_continuous_calls;
+ bool m_decreased_step_during_last_call;
+
+ time_type m_dt_last;
+ time_type m_t_last;
+ time_type m_current_eps;
+
+ unsigned short m_current_k_max;
+ unsigned short m_current_k_opt;
+
+ container_type m_x0;
+ container_type m_xerr;
+ container_type m_x_mp;
+ container_type m_x_scale;
+ container_type m_dxdt;
+
+ typedef std::vector< time_type > value_vector;
+ typedef std::vector< std::vector< time_type > > value_matrix;
+ typedef std::vector< unsigned short > us_vector;
+
+ value_vector m_error;
+ value_vector m_a;
+ value_matrix m_alpha;
+ us_vector m_interval_sequence;
+
+ value_vector m_times;
+ std::vector< container_type > m_d;
+ container_type m_c;
+
+ // public functions
+ public:
+
+ // constructor
+ controlled_stepper_bs(
+ 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_k_max(8),
+ m_safety1(0.25), m_safety2(0.7),
+ m_max_dt_factor( 0.1 ), m_min_step_scale(5E-5), m_max_step_scale(0.7),
+ m_continuous_calls(false), m_decreased_step_during_last_call( false ),
+ m_dt_last( 1.0E30),
+ m_current_eps( -1.0 )
+ {
+ m_error.resize(m_k_max);
+ m_a.resize(m_k_max+1);
+ m_alpha.resize(m_k_max); // k_max * k_max matrix
+ typename value_matrix::iterator it = m_alpha.begin();
+ while( it != m_alpha.end() )
+ (*it++).resize(m_k_max);
+ m_interval_sequence.resize(m_k_max+1);
+ for( unsigned short i = 1; i <= m_k_max+1; i++ )
+ m_interval_sequence[i-1] = (2*i);
+
+ m_times.resize(m_k_max);
+ m_d.resize(m_k_max);
+ }
+
+
+ template< class DynamicalSystem >
+ controlled_step_result try_step(
+ DynamicalSystem &system ,
+ container_type &x ,
+ container_type &dxdt ,
+ time_type &t ,
+ time_type &dt )
+ {
+ m_resizer.adjust_size(x, m_xerr);
+ m_resizer.adjust_size(x, m_x_mp);
+ m_resizer.adjust_size(x, m_x_scale);
+
+ // get error scale
+ m_error_checker.fill_scale(x, dxdt, dt, m_x_scale);
+
+ //std::clog << " x scale: " << m_x_scale[0] << '\t' << m_x_scale[1] << std::endl;
+
+ m_x0 = x; // save starting state
+ time_type max_eps = m_error_checker.get_epsilon();
+ if( m_current_eps != max_eps )
+ { // error changed -> recalculate tableau
+ initialize_step_adjust_tableau( max_eps );
+ m_current_eps = max_eps;
+ m_dt_last = m_t_last = std::numeric_limits< value_type >::max(); // unrealistic
+ }
+ // if t and dt are exactly the parameters from last step we are called continuously
+ bool continuous_call = ((t == m_t_last) || (dt == m_dt_last ));
+ if( !continuous_call )
+ m_current_k_opt = m_current_k_max;
+
+ bool converged = false;
+ value_type step_scale = 1.0;
+ unsigned short k_conv = 0;
+
+ for( unsigned short k=0; k<=m_current_k_max; k++ )
+ {
+ unsigned short stepcount = m_interval_sequence[k];
+ //out-of-place midpoint step
+ m_stepper_mp.set_stepcount(stepcount);
+ m_stepper_mp.do_step(system, m_x0, dxdt, t, dt, m_x_mp);
+ //std::clog << "x_mp: " << k << '\t' << m_x_mp[0] << '\t' << m_x_mp[1] << std::endl;
+ time_type t_est = (dt/stepcount)*(dt/stepcount);
+ extrapolate(k, t_est, m_x_mp, x, m_xerr);
+ //std::clog << "Error: " << k << '\t' << m_xerr[0] << '\t' << m_xerr[1] << std::endl;
+ if( k != 0 )
+ {
+ value_type max_err = m_error_checker.get_max_error_ratio(m_xerr, m_x_scale);
+ m_error[k-1] = std::pow( max_err/m_safety1, 1.0/(2*k+1) );
+ if( (k >= m_current_k_opt-1) || !continuous_call )
+ { //we're in the order window where convergence is expected
+ if( max_err < 1.0 )
+ {
+ k_conv = k;
+ converged = true;
+ //std::clog << "Converged at: " << k << '\t' << max_err << std::endl;
+ break;
+ } else {
+ converged = false;
+ if( (k == m_current_k_max) || (k == m_current_k_opt+1) )
+ {
+ step_scale = m_safety2/m_error[k-1];
+ break;
+ } else if( (k == m_current_k_opt) &&
+ (m_alpha[m_current_k_opt-1][m_current_k_opt] < m_error[k-1] ) )
+ {
+ step_scale = static_cast<value_type>(1.0)/m_error[k-1];
+ break;
+ } else if( (m_current_k_opt == m_current_k_max) &&
+ (m_alpha[k-1][m_current_k_max-1] < m_error[k-1]) )
+ {
+ step_scale = m_alpha[k-1][m_current_k_max-1]*m_safety2/m_error[k-1];
+ //std::clog << " Case 3 " << m_alpha[k-1][m_current_k_max-1] << '\t' << m_error[k-1] << '\t' << step_scale << std::endl;
+ break;
+ } else if( m_alpha[k-1][m_current_k_opt] < m_error[k-1] )
+ {
+ step_scale = m_alpha[k-1][m_current_k_opt-1]/m_error[k-1];
+ break;
+ }
+ }
+ }
+ }
+ }
+ if( !converged ) { // dt was too large - no converge up to order k_max
+
+ // step_scale is always > 1
+ step_scale = std::max(step_scale, m_min_step_scale); // at least m_min ...
+ step_scale = std::min(step_scale, m_max_step_scale); // at most m_max ...
+ dt *= step_scale;
+ m_dt_last = dt;
+ m_t_last = t;
+ m_decreased_step_during_last_call = true;
+ x = m_x0; // copy the state back
+ return step_size_decreased;
+
+ } else { //converged
+
+ t += dt; // step accomplished
+
+ time_type work_min = std::numeric_limits< time_type >::max();
+ for( unsigned short k=0; k < k_conv ; k++ )
+ { // compute optimal convergence order and corresponding stepsize
+ time_type factor = std::max(m_error[k], m_max_dt_factor);
+ time_type work = factor * m_a[k+1];
+ if( work < work_min )
+ {
+ step_scale = factor;
+ work_min = work;
+ m_current_k_opt = k+1;
+ }
+ }
+ m_dt_last = dt/step_scale;
+
+ //std::clog << "Internal: " << dt << '\t' << k_conv-1 << '\t' << m_current_k_opt << '\t' << m_current_k_max << '\t' << m_decreased_step_during_last_call << std::endl;
+
+ if( (m_current_k_opt >= k_conv) &&
+ (m_current_k_opt != m_current_k_max) &&
+ !m_decreased_step_during_last_call )
+ { // check possible order increase, if step was not decreased before
+ time_type factor = std::max(
+ step_scale/m_alpha[m_current_k_opt-1][m_current_k_opt],
+ m_max_dt_factor );
+ if( m_a[m_current_k_opt+1]*factor <= work_min )
+ {
+ m_dt_last = dt/factor;
+ m_current_k_opt++;
+ }
+ }
+ dt = m_dt_last;
+ m_t_last = t;
+ m_decreased_step_during_last_call = false;
+ //std::clog << "Internal: " << dt << '\t' << m_current_k_opt << std::endl;
+ return step_size_increased;
+ }
+ }
+
+ template< class System >
+ controlled_step_result try_step(
+ System &system,
+ container_type &x,
+ time_type &t,
+ time_type &dt )
+ {
+ m_resizer.adjust_size(x, m_dxdt);
+ system(x, m_dxdt, t);
+ return try_step(system, x, m_dxdt, t, dt );
+ }
+
+
+ private: // private functions
+
+ void initialize_step_adjust_tableau( time_type eps )
+ {
+ m_a[0] = m_interval_sequence[0]+1;
+ for( unsigned short k=0; k<m_k_max; k++ )
+ m_a[k+1] = m_a[k] + m_interval_sequence[k+1];
+ for( unsigned short i=1; i<m_k_max; i++ )
+ {
+ for( unsigned short k=0; k<i; k++ )
+ {
+ m_alpha[k][i] = std::pow(
+ m_safety1*eps,
+ (m_a[k+1]-m_a[i+1])/
+ ((m_a[i+1]-m_a[0]+1.0)*(2*k+3))
+ );
+ }
+ }
+ m_current_k_opt = m_k_max-1;
+ for( unsigned short k=1; k<m_k_max-1; k++ )
+ {
+ if( m_a[k+1] > m_a[k]*m_alpha[k-1][k] )
+ {
+ m_current_k_opt = k;
+ break;
+ }
+ }
+ m_current_k_max = m_current_k_opt;
+ }
+
+ void extrapolate(
+ unsigned short k_est,
+ time_type t_est,
+ container_type &x_est,
+ container_type &x,
+ container_type &x_err )
+ {
+ //m_resizer.adjust_size(x, m_c);
+ //std::vector< container_type > m_d_iter = m_d.begin();
+ //while( m_d_iter != m_d.end() )
+ // m_resizer.adjust_size(x, (*m_d_iter++));
+
+ m_times[k_est] = t_est;
+ x_err = x = x_est;
+
+ const iterator x_end = x.end();
+
+ if( k_est == 0 )
+ {
+ m_d[0] = x_est;
+ } else
+ {
+ m_c = x_est;
+ for( unsigned short k=0; k<k_est; k++ )
+ {
+ value_type delta = static_cast<value_type>(1.0) /
+ (m_times[k_est-k-1]-static_cast<value_type>(t_est));
+ value_type val1 = static_cast<value_type>(t_est)*delta;
+ value_type val2 = m_times[k_est-k-1]*delta;
+
+ //std::clog << " values: " << delta << '\t' << val1 << '\t' << val2 << std::endl;
+
+ iterator x_iter = x.begin();
+ iterator x_err_iter = x_err.begin();
+ iterator d_k_iter = m_d[k].begin();
+ iterator c_iter = m_c.begin();
+ while( x_iter != x_end )
+ {
+ //std::clog << " extrapolate: " << '\t' << *x_iter << '\t' << *x_err_iter << '\t' << *d_k_iter << std::endl;
+ value_type q = *d_k_iter;
+ *d_k_iter++ = *x_err_iter;
+ delta = *c_iter - q;
+ *x_err_iter = val1 * delta;
+ *c_iter++ = val2 * delta;
+ *x_iter++ += *x_err_iter++;
+ }
+ //std::clog << " extrapolate: " << '\t' << x[1] << '\t' << x_err[1] << '\t' << m_d[k][1] << std::endl;
+ m_d[k_est] = x_err;
+ }
+ }
+ }
+
+ };
+}
+}
+}
+
+#endif
Added: sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -0,0 +1,79 @@
+
+
+#ifndef BOOST_NUMERIC_ODEINT_ERROR_CHECKER_STANDARD_HPP
+#define BOOST_NUMERIC_ODEINT_ERROR_CHECKER_STANDARD_HPP
+
+#include <cmath>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+ template< class Container, class Time >
+ class error_checker_standard {
+
+
+ public:
+ typedef Container container_type;
+ typedef Time time_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::iterator 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,
+ time_type dt,
+ container_type &scale )
+ {
+ iterator x_iter = x.begin();
+ const iterator x_end = x.end();
+ iterator dxdt_iter = dxdt.begin();
+ iterator scale_iter = scale.begin();
+ while( x_iter != x_end )
+ {
+ *scale_iter++ = m_eps_abs +
+ m_eps_rel * (m_a_x * std::abs(*x_iter++) +
+ m_a_dxdt * dt * std::abs(*dxdt_iter++));
+ }
+ }
+
+ time_type get_max_error_ratio( container_type &x_err, container_type &scale)
+ {
+ iterator x_err_iter = x_err.begin();
+ const iterator x_err_end = x_err.end();
+ iterator scale_iter = scale.begin();
+ time_type max_rel_error = static_cast<time_type>(0.0);
+ while( x_err_iter != x_err_end )
+ {
+ max_rel_error = std::max(
+ static_cast<time_type>(std::abs(*x_err_iter++)/(*scale_iter++)) ,
+ max_rel_error);
+ }
+ return max_rel_error;
+ }
+
+ const time_type get_epsilon() { return std::max(m_eps_rel, m_eps_abs); }
+
+
+ };
+
+}
+}
+}
+
+#endif
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 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -23,7 +23,7 @@
#include <boost/numeric/odeint/concepts/state_concept.hpp>
#include <boost/numeric/odeint/resizer.hpp>
-
+#include <iostream>
namespace boost {
namespace numeric {
@@ -105,12 +105,14 @@
using namespace detail::it_algebra;
// m_x0 = x + h*dxdt
- scale_sum( m_x0.begin(), m_x0.end(),
+ scale_sum( m_x1.begin(), m_x1.end(),
t_1, x.begin(),
h, dxdt.begin() );
- system( m_x0, m_dxdt, th );
+ system( m_x1, m_dxdt, th );
- m_x1 = x;
+ m_x0 = x;
+
+ //std::clog<< "mp: " << 0 << '\t' << m_x1[0] << '\t' << m_x1[1] << '\t' << m_dxdt[0] << '\t' << m_dxdt[1] << std::endl;
unsigned short i = 1;
while( i != m_stepcount )
@@ -123,9 +125,11 @@
system( m_x1, m_dxdt, th);
i++;
}
+
+ //std::clog<< "mp: " << i << '\t' << m_x1[0] << '\t' << m_x1[1] << '\t' << m_dxdt[0] << '\t' << m_dxdt[1] << std::endl;
// last step
- // x = 0.5*( m_x0 + m_x1 + h*m_dxdt
+ // x = 0.5*( m_x0 + m_x1 + h*m_dxdt )
scale_sum( x_out.begin(), x_out.end(),
t_05, m_x0.begin(),
t_05, m_x1.begin(),
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 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -58,8 +58,7 @@
// public interface
public:
- order_type order() const { return 8; }
-
+ order_type order() const { return 7; }
template< class DynamicalSystem >
void do_step( DynamicalSystem &system ,
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 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -18,11 +18,13 @@
*/
#include <iostream>
+#include <fstream>
#include <vector>
#include <list>
#include <tr1/array>
#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/controlled_stepper_bs.hpp>
#define tab "\t"
@@ -34,43 +36,77 @@
const double R = 28.0;
const double b = 8.0 / 3.0;
-const size_t olen = 10000;
-
-const double eps_abs = 1E-2;
-const double eps_rel = 1E-3;;
-
typedef array<double, 3> state_type;
+size_t function_calls = 0;
+
void lorenz( state_type &x , state_type &dxdt , double t )
{
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];
+ function_calls++;
}
-void print_state( double t, state_type &x, ... )
+
+class output_observer
{
- cout << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
-}
+
+ ofstream m_file;
+
+public:
+ output_observer( const char* file_name )
+ : m_file( file_name )
+ {
+ m_file.precision(10);
+ //m_file.setf(ios::fixed,ios::floatfield);
+ }
+
+ ~output_observer()
+ {
+ m_file.close();
+ }
+
+ void operator()( double t, state_type &x, ... )
+ {
+ m_file << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
+ }
+};
int main( int argc , char **argv )
{
+ const double end_time = 25.0;
+
+ const double eps_abs = 1E-10;
+ const double eps_rel = 1E-10;;
+
state_type x;
x[0] = 1.0;
x[1] = 0.0;
x[2] = 20.0;
- stepper_half_step< stepper_euler< state_type > > euler;
- controlled_stepper_standard< stepper_half_step< stepper_euler< state_type > > >
- controlled_stepper( euler, eps_abs , eps_rel, 1.0, 1.0);
-
- cout.precision(5);
- cout.setf(ios::fixed,ios::floatfield);
+ 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 );
+ output_observer rk5_obs("lorenz_rk5.dat");
+ size_t steps = integrate_adaptive( controlled_rk5, lorenz, x, 0.0, end_time, 1E-2, rk5_obs );
+
+ clog << "RK5: " << steps << " steps. (" << function_calls << " function calls)" << endl;
+
+ x[0] = 1.0;
+ x[1] = 0.0;
+ x[2] = 20.0;
+
+ function_calls = 0;
+
+ controlled_stepper_bs< state_type > controlled_bs(eps_abs, eps_rel, 1.0, 1.0);
- size_t steps = integrate_adaptive( controlled_stepper, lorenz, x, 0.0, 10.0, 1E-2, print_state );
+ output_observer bs_obs("lorenz_bs.dat");
+ steps = integrate_adaptive( controlled_bs, lorenz, x, 0.0, end_time, 1E-2, bs_obs );
- clog << "Number of steps: " << steps << endl;
+ clog << "BS: " << steps << " steps. (" << function_calls << " function calls)" << endl;
+
+
return 0;
}
Added: sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_bs.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_bs.cpp 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -0,0 +1,51 @@
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+#include <cmath>
+
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+typedef std::tr1::array< double , 2 > state_type;
+
+
+const double gam = .15;
+
+void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t)
+{
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - gam*x[1];
+}
+
+
+int main( int argc , char **argv )
+{
+ state_type x1 = {{ 1.0, 0.0 }};
+
+ const size_t step_count = 100;
+ controlled_stepper_bs< state_type > stepper_bs(0.0, 1E-11, 1.0, 1.0);
+ double t = 0.0;
+ double dt = 0.01;
+ cout << " Everything initialized - starting time evolution " << endl;
+ cout.precision(16);
+ clog.precision(16);
+ for( size_t step=0; step < step_count; step++ )
+ {
+ stepper_bs.try_step(harmonic_oscillator, x1, t, dt);
+ //clog << " ####################################################### " << endl;
+ cout << t << tab << dt << tab << x1[0] << tab << x1[1] << endl;
+ //clog << " ####################################################### " << endl;
+ }
+}
+
+/*
+ Compile with
+ g++ -Wall -O3 -I$BOOST_ROOT -I../../../../../ stepper_bs.cpp
+*/
+
Modified: sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp 2009-11-29 15:36:32 EST (Sun, 29 Nov 2009)
@@ -41,11 +41,14 @@
double t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
{
- stepper_euler.next_step( harmonic_oscillator , x1 , t , dt );
- stepper_mp.next_step( harmonic_oscillator , x2 , t , dt , 2 );
- stepper_mp.next_step( harmonic_oscillator , x3 , t , dt , 4 );
- stepper_mp.next_step( harmonic_oscillator , x4 , t , dt , 8 );
- stepper_rk4.next_step( harmonic_oscillator , x5 , t , dt );
+ stepper_euler.do_step( harmonic_oscillator , x1 , t , dt );
+ stepper_mp.set_stepcount(2);
+ stepper_mp.do_step( harmonic_oscillator , x2 , t , dt );
+ stepper_mp.set_stepcount(4);
+ stepper_mp.do_step( harmonic_oscillator , x3 , t , dt );
+ stepper_mp.set_stepcount(8);
+ stepper_mp.do_step( harmonic_oscillator , x4 , t , dt );
+ stepper_rk4.do_step( harmonic_oscillator , x5 , t , dt );
cout<< t << tab << x1[0]*x1[0] + x1[1]*x1[1];
cout<< tab << x2[0]*x2[0] + x2[1]*x2[1];
cout<< tab << x3[0]*x3[0] + x3[1]*x3[1];