Subject: [Boost-commit] svn:boost r66426 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher
From: karsten.ahnert_at_[hidden]
Date: 2010-11-07 08:48:20
Author: karsten
Date: 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
New Revision: 66426
URL: http://svn.boost.org/trac/boost/changeset/66426
Log:
* butcher tableau steppers introduced in ideas
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Makefile (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/convert_value.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/diagnostic.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/predefined_steppers.hpp (contents, props changed)
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Makefile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Makefile 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,22 @@
+CC = g++
+
+RATIO_ROOT = $(HOME)/boost/chrono
+
+CXXFLAGS = -O3 -I$(BOOST_ROOT) -I$(ODEINT_ROOT) -I$(RATIO_ROOT)
+
+HEADERS = algebra.hpp convert_value.hpp diagnostic.hpp explicit_runge_kutta.hpp predefined_steppers.hpp
+
+all : butcher_performance butcher_lorenz butcher_vdp
+
+butcher_performance : butcher_performance.o
+butcher_performance.o : butcher_performance.cpp $(HEADERS)
+
+butcher_vdp : butcher_vdp.o
+butcher_vdp.o : butcher_vdp.cpp $(HEADERS)
+
+butcher_lorenz : butcher_lorenz.o
+butcher_lorenz.o : butcher_lorenz.cpp $(HEADERS)
+
+
+clean :
+ -rm *.o *~ butcher_performance butcher_vdp butcher_lorenz
\ No newline at end of file
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,103 @@
+/*
+ * algebra.hpp
+ *
+ * Created on: Nov 7, 2010
+ * Author: karsten
+ */
+
+#ifndef ALGEBRA_HPP_
+#define ALGEBRA_HPP_
+
+#include <boost/mpl/int.hpp>
+#include <boost/mpl/at.hpp>
+
+#include "convert_value.hpp"
+
+namespace mpl = boost::mpl;
+
+template< class state_type , class coef_tye , class index >
+struct algebra
+{
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ {}
+};
+
+template< class state_type , class coef_type >
+struct algebra< state_type , coef_type , mpl::int_< 0 > >
+{
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ {
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ x_tmp[i] = x[i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i];
+ }
+ }
+};
+
+template< class state_type , class coef_type >
+struct algebra< state_type , coef_type , mpl::int_< 1 > >
+{
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ {
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ x_tmp[i] = x[i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i];
+ }
+ }
+};
+
+template< class state_type , class coef_type >
+struct algebra< state_type , coef_type , mpl::int_< 2 > >
+{
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ {
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ x_tmp[i] = x[i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i];
+ }
+
+ }
+};
+
+template< class state_type , class coef_type >
+struct algebra< state_type , coef_type , mpl::int_< 3 > >
+{
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ {
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ x_tmp[i] = x[i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i];
+ }
+ }
+};
+
+template< class state_type , class coef_type >
+struct algebra< state_type , coef_type , mpl::int_< 4 > >
+{
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ {
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ x_tmp[i] = x[i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i] +
+ dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value() * k_vector[4][i];
+ }
+ }
+};
+
+
+
+#endif /* ALGEBRA_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,66 @@
+/*
+ * butcher_test.cpp
+ *
+ * Created on: Nov 5, 2010
+ * Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#include "predefined_steppers.hpp"
+
+#define tab "\t"
+
+using namespace std;
+
+typedef std::tr1::array< double , 3 > state_type;
+typedef mpl_rk4_stepper< state_type > stepper_type;
+typedef boost::numeric::odeint::stepper_rk4< state_type > stepper_std_type;
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const 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];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+ stepper_type stepper;
+ stepper_std_type stepper_std;
+ state_type x = {{ 1.0 , 1.0 , 2.0 }};
+ state_type x_std = x;
+
+ cout << "Tableau : " << endl;
+ stepper.print_tableau();
+ cout << endl;
+
+ double t = 0.0 , dt = 0.01;
+ ofstream fout( "lorenz.dat" );
+ for( size_t i=0 ; i<10000 ; ++i )
+ {
+ fout << t << tab;
+ fout << x[0] << tab << x[1] << tab << x[2] << tab;
+ fout << x_std[0] << tab << x_std[1] << tab << x_std[2] << "\n";
+
+ stepper.do_step( lorenz , x , t , dt );
+ stepper_std.do_step( lorenz , x_std , t , dt );
+ t += dt;
+ }
+
+
+ return 0;
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,100 @@
+/*
+ * butcher_test.cpp
+ *
+ * Created on: Nov 5, 2010
+ * Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#include "predefined_steppers.hpp"
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::accumulators;
+
+typedef accumulator_set<
+ double , stats< tag::mean , tag::variance >
+ > accumulator_type;
+
+ostream& operator<<( ostream& out , accumulator_type &acc )
+{
+ out << boost::accumulators::mean( acc ) << tab;
+// out << boost::accumulators::variance( acc ) << tab;
+ return out;
+}
+
+typedef boost::timer timer_type;
+
+
+typedef std::tr1::array< double , 3 > state_type;
+typedef mpl_euler_stepper< state_type > stepper_type;
+typedef boost::numeric::odeint::stepper_euler< state_type > stepper_std_type;
+
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double 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];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+ stepper_type stepper;
+ stepper_std_type stepper_std;
+
+ const size_t num_of_steps = 10000000;
+ const size_t dt = 0.01;
+
+ accumulator_type acc1 , acc2;
+ timer_type timer;
+
+ srand48( 12312354 );
+
+ while( true )
+ {
+ state_type x = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
+ state_type x_std = x;
+ double t = 0.0 ;
+ double t_std = 0.0;
+
+ timer.restart();
+ for( size_t i=0 ; i<num_of_steps ; ++i,t+=dt )
+ stepper.do_step( lorenz , x , t , dt );
+ acc1( timer.elapsed() );
+
+ timer.restart();
+ for( size_t i=0 ; i<num_of_steps ; ++i,t_std+=dt )
+ stepper_std.do_step( lorenz , x_std , t_std , dt );
+ acc2( timer.elapsed() );
+
+ clog.precision( 3 );
+ clog.width( 5 );
+ clog << acc1 << " " << acc2 << " " << x[0] << " " << x_std[0] << endl;
+ }
+
+
+
+ return 0;
+}
+
+/*
+ * Compile with :
+ * g++ -O3 -I$BOOST_ROOT -I$HOME/boost/chrono -I$ODEINT_ROOT butcher_performance.cpp
+ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,69 @@
+/*
+ * butcher_test.cpp
+ *
+ * Created on: Nov 5, 2010
+ * Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#include "predefined_steppers.hpp"
+
+#define tab "\t"
+
+using namespace std;
+
+typedef std::tr1::array< double , 2 > state_type;
+typedef mpl_rk4_stepper< state_type > stepper_type;
+typedef boost::numeric::odeint::stepper_rk4< state_type > stepper_std_type;
+
+void vdp( const state_type &x , state_type &dxdt , double t )
+{
+ const double mu = 0.2;
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] + mu * ( 1.0 - x[0]*x[0] ) * x[1];
+}
+
+typedef mpl_rk4_stepper< state_type > stepper_type;
+typedef boost::numeric::odeint::stepper_rk4< state_type > stepper_std_type;
+
+
+
+
+int main( int argc , char **argv )
+{
+ stepper_type midpoint;
+ stepper_std_type midpoint_std;
+ state_type x = {{ 1.0 , 1.0 }};
+ state_type x_std = x;
+
+ cout << "Tableau : " << endl;
+ midpoint.print_tableau();
+ cout << endl;
+
+ double t = 0.0 , dt = 0.01;
+ ofstream fout( "vdp.dat" );
+ for( size_t i=0 ; i<10000 ; ++i )
+ {
+ fout << t << tab;
+ fout << x[0] << tab << x[1] << tab ;
+ fout << x_std[0] << tab << x_std[1] << "\n";
+
+ midpoint.do_step( vdp , x , t , dt );
+ midpoint_std.do_step( vdp , x_std , t , dt );
+ t += dt;
+ }
+
+
+
+// cout << double( one_half::num ) / double( one_half::den ) << endl;
+// cout << conv< one_half >::get_value() << endl;
+// cout << conv< one >::get_value() << endl;
+
+ return 0;
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/convert_value.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/convert_value.hpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,36 @@
+/*
+ * convert_values.hpp
+ *
+ * Created on: Nov 7, 2010
+ * Author: karsten
+ */
+
+#ifndef CONVERT_VALUE_HPP_
+#define CONVERT_VALUE_HPP_
+
+
+#include <boost/ratio.hpp>
+
+
+template< class T >
+struct convert_value
+{
+ static double get_value( void )
+ {
+ return double( T::value );
+ }
+};
+
+
+template< long N , long D >
+struct convert_value< boost::ratio< N , D > >
+{
+ typedef typename boost::ratio< N , D > number;
+ static double get_value( void )
+ {
+ return double( number::num ) / double( number::den );
+ }
+};
+
+
+#endif /* CONVERT_VALUE_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/diagnostic.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/diagnostic.hpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,59 @@
+/*
+ * diagnostic.hpp
+ *
+ * Created on: Nov 7, 2010
+ * Author: karsten
+ */
+
+#ifndef DIAGNOSTIC_HPP_
+#define DIAGNOSTIC_HPP_
+
+#include <iostream>
+
+#include <boost/mpl/for_each.hpp>
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/int.hpp>
+#include <boost/mpl/size.hpp>
+
+#include "convert_value.hpp"
+
+namespace mpl = boost::mpl;
+
+struct print_value
+{
+ template< class T >
+ void operator()( T )
+ {
+ std::cout << convert_value< T >::get_value() << " ";
+ }
+};
+
+struct print_vector
+{
+ template< class Vec >
+ void operator()( Vec )
+ {
+ mpl::for_each< Vec >( print_value() );
+ std::cout << "\n";
+ }
+};
+
+struct print_tableau_vector
+{
+ template< class Vec >
+ void operator()( Vec )
+ {
+ typedef typename mpl::at< Vec , mpl::int_< 0 > >::type index;
+ typedef typename mpl::at< Vec , mpl::int_< 1 > >::type c;
+ typedef typename mpl::at< Vec , mpl::int_< 2 > >::type coef;
+
+ std::cout << "Size : " << mpl::size< Vec >::value << "\t";
+ std::cout << convert_value< index >::get_value() << " ";
+ std::cout << convert_value< c >::get_value() << " ";
+ mpl::for_each< coef >( print_value() );
+ std::cout << "\n";
+ }
+};
+
+
+#endif /* DIAGNOSTIC_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,159 @@
+/*
+ * explicit_runge_kutta.hpp
+ *
+ * Created on: Nov 5, 2010
+ * Author: karsten
+ */
+
+#ifndef EXPLICIT_RUNGE_KUTTA_HPP_
+#define EXPLICIT_RUNGE_KUTTA_HPP_
+
+#include <boost/static_assert.hpp>
+#include <boost/ratio.hpp>
+#include <boost/type_traits/is_same.hpp>
+
+#include <boost/mpl/vector.hpp>
+#include <boost/mpl/size.hpp>
+#include <boost/mpl/push_back.hpp>
+#include <boost/mpl/for_each.hpp>
+#include <boost/mpl/zip_view.hpp>
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/int.hpp>
+#include <boost/mpl/at.hpp>
+
+#include "convert_value.hpp"
+#include "diagnostic.hpp"
+#include "algebra.hpp"
+
+namespace mpl = boost::mpl;
+
+
+/*
+ * c_1 |
+ * c_2 | a_{2,1}
+ * c_3 | a_{3,1} a_{3,2}
+ * ...
+ * c_s | a_{s,1} a_{s,2} ... a_{s,s-1}
+ * -----------------------------------
+ * | b_1 b_2 b_{s-1} b_s
+ *
+ * ButcherA : the coefficients a_ij
+ *
+ * ButcherC : the time coefficients c_i
+ *
+ * ButcherC : the sum coefficients b_j
+ */
+template< class State , class ButcherA , class ButcherC , class ButcherB >
+class explicit_runge_kutta
+{
+public:
+
+ typedef ButcherA butcher_a;
+ typedef ButcherB butcher_b;
+ typedef ButcherC butcher_c;
+
+ static const size_t dim = mpl::size< ButcherC >::value;
+
+ typedef mpl::range_c< size_t , 0 , dim > indices;
+ typedef typename mpl::push_back< butcher_a , butcher_b >::type butcher_right;
+ typedef mpl::zip_view< mpl::vector< indices , butcher_c , butcher_right > > butcher_tableau;
+
+
+ typedef double time_type;
+ typedef State state_type;
+
+
+ /*
+ * t , dt , system, x, k[]
+ */
+ template< class System >
+ struct calculate_stage
+ {
+ System &system;
+ state_type &x , &x_tmp;
+ state_type *k_vector;
+ double t;
+ double dt;
+
+ calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , double _t , double _dt )
+ : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
+ {}
+
+ /*
+ * for( size_t j=0 ; j<dim ; ++j )
+ * {
+ * system( x[j] , k[j] , c[j] )
+ * x[j+1] = sum( i , a[j][i] , k[j] )
+ * }
+ */
+
+ template< class Stage > void operator()( Stage v )
+ {
+ typedef typename mpl::at< Stage , mpl::int_< 0 > >::type index_type;
+ typedef typename mpl::at< Stage , mpl::int_< 1 > >::type time_value_type;
+ typedef typename mpl::at< Stage , mpl::int_< 2 > >::type coef_type;
+
+ const static size_t index = index_type::value;
+ const static double time_value = convert_value< time_value_type >::get_value();
+
+ BOOST_STATIC_ASSERT(( ( mpl::size< coef_type >::value - 1 ) == index ));
+
+// cout << index << " " << time_value << " " << endl;
+
+ if( index == 0 ) system( x , k_vector[index] , t + time_value * dt );
+ else system( x_tmp , k_vector[index] , t + time_value * dt );
+
+ if( index == ( dim - 1 ) )
+ algebra< state_type , coef_type , mpl::int_< index > >::do_step( x , x , k_vector , dt );
+ else
+ algebra< state_type , coef_type , mpl::int_< index > >::do_step( x_tmp , x , k_vector , dt );
+ };
+ };
+
+
+ explicit_runge_kutta( void )
+ {
+ }
+
+ template< class System >
+ void do_step( System &system , state_type &x , double t , double dt )
+ {
+ mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
+ }
+
+
+
+ void print_butcher_a( void )
+ {
+ mpl::for_each< butcher_a >( print_vector() );
+ }
+
+ void print_butcher_b( void )
+ {
+ mpl::for_each< butcher_b >( print_value() );
+ }
+
+ void print_butcher_c( void )
+ {
+ mpl::for_each< butcher_c >( print_value() );
+ }
+
+ void print_tableau( void )
+ {
+ mpl::for_each< butcher_tableau >( print_tableau_vector() );
+ }
+
+
+
+private:
+
+ state_type m_k_vector[dim];
+ state_type m_x_tmp;
+
+ BOOST_STATIC_ASSERT(( dim > 0 ));
+ BOOST_STATIC_ASSERT(( dim == size_t( mpl::size< butcher_b >::value ) ));
+ BOOST_STATIC_ASSERT(( ( dim - 1 ) == mpl::size< butcher_a >::value ));
+};
+
+
+#endif /* EXPLICIT_RUNGE_KUTTA_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/predefined_steppers.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/predefined_steppers.hpp 2010-11-07 08:48:11 EST (Sun, 07 Nov 2010)
@@ -0,0 +1,77 @@
+/*
+ * predefined_steppers.hpp
+ *
+ * Created on: Nov 7, 2010
+ * Author: karsten
+ */
+
+#ifndef PREDEFINED_STEPPERS_HPP_
+#define PREDEFINED_STEPPERS_HPP_
+
+#include <boost/mpl/int.hpp>
+#include <boost/ratio.hpp>
+
+#include "explicit_runge_kutta.hpp"
+
+typedef mpl::int_< 0 > null;
+typedef mpl::int_< 1 > one;
+typedef boost::ratio< 1 , 2 > one_half;
+typedef boost::ratio< 1 , 3 > one_third;
+typedef boost::ratio< 1 , 6 > one_sixth;
+
+
+
+/*
+ * euler :
+ * 0 |
+ * | 1
+ */
+
+typedef mpl::vector< null > euler_c;
+typedef mpl::vector<> euler_a;
+typedef mpl::vector< one > euler_b;
+
+template< class state_type >
+struct mpl_euler_stepper : public explicit_runge_kutta< state_type , euler_a , euler_c , euler_b > { };
+
+
+/*
+ * midpoint :
+ * 0 |
+ * 1/2 | 1/2
+ * -----------
+ * | 0 1
+ */
+
+typedef mpl::vector< null , one_half > midpoint_c;
+typedef mpl::vector< mpl::vector< one_half > > midpoint_a;
+typedef mpl::vector< null , one > midpoint_b;
+
+template< class state_type >
+struct mpl_midpoint_stepper : public explicit_runge_kutta< state_type , midpoint_a , midpoint_c , midpoint_b > { };
+
+
+
+
+/*
+ * classical Runge Kutta
+ * 0 |
+ * 1/2 | 1/2
+ * 1/2 | 0 1/2
+ * 1 | 0 0 1
+ * ------------------
+ * | 1/6 1/3 1/6 1/6
+ *
+ */
+typedef mpl::vector< null , one_half , one_half , one > rk4_c;
+typedef mpl::vector<
+ mpl::vector< one_half > ,
+ mpl::vector< null , one_half > ,
+ mpl::vector< null , null , one >
+ > rk4_a;
+typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;
+
+template< class state_type >
+struct mpl_rk4_stepper : public explicit_runge_kutta< state_type , rk4_a , rk4_c , rk4_b > { };
+
+#endif /* PREDEFINED_STEPPERS_HPP_ */