boost/odeint
2
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2026-01-26 06:42:23 +00:00
Code Issues Projects Releases Wiki Activity

+ fusion stepper

Browse Source
This commit is contained in:
mariomulansky
2010-11-21 17:29:28 +00:00
parent fb1f32778a
commit f12fd3c00a
3 changed files with 298 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
libs/numeric/odeint/ideas/generic_stepper/algebra.hpp
View File

@@ -42,4 +42,19 @@ struct algebra< state_type , 1 >
};
template< typename state_type>
struct algebra< state_type , 2 >
{
typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
typedef typename boost::numeric::odeint::standard_operations< double > std_op;
static void foreach( state_type &x_tmp , const state_type &x , const std::vector< double > &a , const state_type *k_vector , const double dt )
{
std_algebra::for_each4( x_tmp , x , k_vector[0] , k_vector[1] , std_op::scale_sum3( 1.0 , a[0]*dt , a[1]*dt ) );
}
};
#endif /* ALGEBRA_HPP_ */

273
libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp Normal file
View File

@@ -0,0 +1,273 @@
/*
* fusion_stepper.hpp
*
* Created on: Nov 21, 2010
* Author: mario
*/
#ifndef FUSION_STEPPER_HPP_
#define FUSION_STEPPER_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 <boost/mpl/copy.hpp>
#include <boost/mpl/insert_range.hpp>
#include <vector>
#include <boost/fusion/container.hpp>
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/include/mpl.hpp>
#include <typeinfo>
#include "algebra.hpp"
namespace mpl = boost::mpl;
namespace fusion = boost::fusion;
using namespace std;
template< typename State , size_t stage >
class stepper_stage
{
public:
typedef double value_type;
typedef State state_type;
typedef vector< vector< double > > parameter_array_type2D;
typedef vector< double > parameter_array_type1D;
void init( const parameter_array_type2D &a , const parameter_array_type1D &c )
{
m_a_row = a[stage-1];
m_c = c[stage-1];
}
template< typename System >
void operator() ( System &system , state_type &x_tmp , const state_type &x , state_type *k_vector ,
double t , const double dt ) const
{
system( x_tmp , k_vector[stage-1] , t + m_c * dt );
algebra<state_type , stage>::foreach< stage >( x_tmp , x , m_a_row , k_vector , dt);
}
private:
vector< value_type > m_a_row;
value_type m_c;
};
template< typename State >
class stepper_stage< State , 1 >
{
public:
typedef double value_type;
typedef State state_type;
typedef vector< vector< double > > parameter_array_type2D;
typedef vector< double > parameter_array_type1D;
void init( const parameter_array_type2D &a , const parameter_array_type1D &c )
{
m_a_row = a[0];
m_c = c[0];
}
template< typename System >
void operator() ( System &system , state_type &x_tmp , const state_type &x , state_type *k_vector ,
double t , const double dt ) const
{
system( x , k_vector[0] , t + m_c * dt );
algebra<state_type , 1>::foreach( x_tmp , x , m_a_row , k_vector , dt);
}
private:
vector< value_type > m_a_row;
value_type m_c;
};
template< typename State , size_t stage >
class stepper_last_stage
{
public:
typedef double value_type;
typedef State state_type;
typedef vector< vector< double > > parameter_array_type2D;
typedef vector< double > parameter_array_type1D;
void init( const parameter_array_type1D &b , const parameter_array_type1D &c )
{
m_b = b;
m_c = c[stage-1];
}
template< typename System >
void operator() ( System &system , state_type &x_tmp , state_type &x , state_type *k_vector , double t , const double dt )
{
system( x_tmp , k_vector[stage-1] , t + m_c * dt );
algebra<state_type , stage>::foreach( x , x , m_b , k_vector , dt);
}
private:
vector< value_type > m_b;
value_type m_c;
};
template< typename State >
class stepper_last_stage< State , 1 >
{
public:
typedef double value_type;
typedef State state_type;
typedef vector< vector< double > > parameter_array_type2D;
typedef vector< double > parameter_array_type1D;
void init( const parameter_array_type1D &b , const parameter_array_type1D &c )
{
m_b = b;
m_c = c[0];
}
template< typename System >
void operator() ( System &system , state_type &x_tmp , state_type &x , state_type k_vector[1] , double t , const double dt )
{
system( x , k_vector[0] , t + m_c * dt );
algebra<state_type , 1>::foreach( x , x , m_b , k_vector , dt);
}
private:
vector< value_type > m_b;
value_type m_c;
};
template< class T , class Constant >
struct stepper_stage_wrapper
{
typedef stepper_stage< T , Constant::value > type;
};
/* Runge Kutta Stepper - consisting of several stages */
template< typename State , size_t stage_count >
class runge_kutta_stepper
{
typedef State state_type;
typedef vector< vector< double > > parameter_array_type2D;
typedef vector< double > parameter_array_type1D;
template< class System >
struct calculate_stage
{
System &system;
state_type &x , &x_tmp;
state_type *k_vector;
const double t;
const double dt;
calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
const double _t , const double _dt )
: system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
{}
template< class Stage >
void operator()( Stage const &stage ) const
{
stage( system , x_tmp , x , k_vector , t , dt );
}
};
struct init_stage {
const parameter_array_type2D &m_a;
const parameter_array_type1D &m_c;
init_stage( const parameter_array_type2D &a , const parameter_array_type1D &c )
: m_a( a ) , m_c( c )
{ }
template< typename Stage >
void operator() ( Stage &stage ) const
{
stage.init( m_a , m_c );
}
};
public:
typedef mpl::range_c< size_t , 1 , stage_count > indices;
typedef typename fusion::result_of::as_vector <
typename mpl::copy // create mpl::vector< stepper_stage< 0 >, stepper_stage< 1 > , .. stepper_stage< stage_count-1 > >
<
indices ,
mpl::inserter
<
mpl::vector0<> ,
mpl::push_back< mpl::_1 , stepper_stage_wrapper< State , mpl::_2 > >
>
>::type
>::type stage_vector;
typedef stepper_last_stage< State , stage_count > last_stage_type;
runge_kutta_stepper( const parameter_array_type2D &a ,
const parameter_array_type1D &b ,
const parameter_array_type1D &c )
{
fusion::for_each( m_stages , init_stage( a , c ) );
m_last_stage.init( b , c );
}
template< class System >
void do_step( System &system , state_type &x , double t , const double dt )
{
fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
m_last_stage( system , m_x_tmp , x , m_k_vector , t , dt );
}
private:
state_type m_k_vector[stage_count];
state_type m_x_tmp;
stage_vector m_stages;
last_stage_type m_last_stage;
};
#endif /* FUSION_STEPPER_HPP_ */

11
libs/numeric/odeint/ideas/generic_stepper/test.cpp
View File

@@ -9,12 +9,13 @@
#include <vector>
#include <tr1/array>
#include "generic_stepper.hpp"
#include "fusion_stepper.hpp"
using namespace std;
typedef tr1::array< double , 3 > state_type;
typedef runge_kutta_stepper< state_type , 1 > euler_stepper;
typedef runge_kutta_stepper< state_type , 2 > midpoint_stepper;
const double sigma = 10.0;
@@ -37,9 +38,17 @@ int main( void )
vector< double > c( 1 , 0.0 );
euler_stepper euler( a , b , c );
vector< vector< double > > a2( 1 , vector<double>( 1 , 0.5 ) );
vector< double > b2( 2 , 0.0 ); b2[1] = 0.5;
vector< double > c2( 1 , 0.0 ); c2[1] = 1.0;
midpoint_stepper midpoint( a2 , b2 , c2 );
state_type x = {{ 1.0 , 1.0 , 2.0 }};
state_type x2 = {{ 1.0 , 1.0 , 2.0 }};
double t = 0.0;
euler.do_step( lorenz , x , t , dt );
midpoint.do_step( lorenz , x2 , t , dt );
cout << x[0] << " , " << x[1] << " , " << x[2] << endl;
cout << x2[0] << " , " << x2[1] << " , " << x2[2] << endl;
}