Subject: [Boost-commit] svn:boost r71116 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: taylor_v1 taylor_v2
From: karsten.ahnert_at_[hidden]
Date: 2011-04-08 09:23:42
Author: karsten
Date: 2011-04-08 09:23:41 EDT (Fri, 08 Apr 2011)
New Revision: 71116
URL: http://svn.boost.org/trac/boost/changeset/71116
Log:
preparing taylor v2
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/Jamfile | 2 +-
1 files changed, 1 insertions(+), 1 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/Jamfile 2011-04-08 09:23:41 EDT (Fri, 08 Apr 2011)
@@ -7,7 +7,7 @@
project
: requirements
<define>BOOST_ALL_NO_LIB=1
- <include>../../../../..
+ <include>../../../../../..
;
exe taylor_main : taylor_main.cpp ;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile 2011-04-08 09:23:41 EDT (Fri, 08 Apr 2011)
@@ -0,0 +1,13 @@
+# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
+# 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)
+
+
+
+project
+ : requirements
+ <define>BOOST_ALL_NO_LIB=1
+ <include>../../../../../..
+ ;
+
+exe taylor_main : taylor_main.cpp ;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp 2011-04-08 09:23:41 EDT (Fri, 08 Apr 2011)
@@ -0,0 +1,372 @@
+/*
+ * taylor.hpp
+ *
+ * Created on: Apr 2, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_
+#define BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_
+
+#include <tr1/array>
+
+// general boost includes
+#include <boost/static_assert.hpp>
+#include <boost/math/special_functions/binomial.hpp>
+#include <boost/math/special_functions/factorials.hpp>
+
+// fusion includes
+#include <boost/fusion/include/is_sequence.hpp>
+#include <boost/fusion/include/size.hpp>
+#include <boost/fusion/include/at.hpp>
+
+// mpl includes
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/for_each.hpp>
+
+// proto includes
+#include <boost/proto/proto.hpp>
+
+
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+namespace taylor_adf
+{
+
+ namespace proto = boost::proto;
+
+ template<int I> struct placeholder {};
+
+ proto::terminal< placeholder< 0 > >::type const arg1 = {{}};
+ proto::terminal< placeholder< 1 > >::type const arg2 = {{}};
+ proto::terminal< placeholder< 2 > >::type const arg3 = {{}};
+
+
+
+ template< class State , size_t Order >
+ struct context_derivs : proto::callable_context< context_derivs< State , Order > const >
+ {
+ typedef double result_type;
+
+ const State &m_x;
+ std::tr1::array< State , Order > &m_derivs;
+ size_t m_which;
+
+ context_derivs( const State &x , std::tr1::array< State , Order > &derivs , size_t which )
+ : m_x( x ) , m_derivs( derivs ) , m_which( which ) { }
+
+
+
+ template< int I >
+ double operator()( proto::tag::terminal , placeholder<I> ) const
+ {
+ return ( m_which == 0 ) ? m_x[I] : m_derivs[m_which-1][I];
+ }
+
+ double operator()( proto::tag::terminal , double x ) const
+ {
+ return ( m_which == 0 ) ? x : 0.0;
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::plus , L const &l , R const &r ) const
+ {
+ return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) +
+ proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::minus , L const &l , R const &r ) const
+ {
+ return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) -
+ proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::multiplies , L const &l , R const &r ) const
+ {
+ double tmp = 0.0;
+ for( size_t k=0 ; k<=m_which ; ++k )
+ {
+ tmp += boost::math::binomial_coefficient< double >( m_which , k )
+ * proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , k ) )
+ * proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which - k ) );
+ }
+ return tmp;
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::divides , L const &l , R const &r ) const
+ {
+ return 0.0;
+ }
+
+ };
+
+ template< class System , class State , size_t Order >
+ struct eval_derivs
+ {
+ typedef std::tr1::array< State , Order > derivs_type;
+
+ System m_sys;
+ const State &m_x;
+ derivs_type &m_derivs;
+ size_t m_which;
+ double m_dt;
+
+ eval_derivs( System sys , const State &x , derivs_type &derivs , size_t which , double dt )
+ : m_sys( sys ) , m_x( x ) , m_derivs( derivs ) , m_which( which ) , m_dt( dt ) { }
+
+ template< class Index >
+ void operator()( Index )
+ {
+ m_derivs[m_which][ Index::value ] = m_dt / double( m_which + 1 ) * boost::proto::eval(
+ boost::fusion::at< Index >( m_sys ) ,
+ taylor_adf::context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ }
+ };
+
+ template< class System , class State , size_t Order >
+ eval_derivs< System , State , Order > make_eval_derivs( System sys , const State &x , std::tr1::array< State , Order > &derivs , size_t i , double dt )
+ {
+ return eval_derivs< System , State , Order >( sys , x , derivs , i , dt );
+ }
+
+
+
+}
+
+template< size_t N , size_t Order >
+class taylor
+{
+public:
+
+ static const size_t dim = N;
+ static const size_t order_value = Order;
+ typedef double value_type;
+ typedef value_type time_type;
+ typedef unsigned short order_type;
+
+ typedef std::tr1::array< value_type , dim > state_type;
+ typedef state_type deriv_type;
+ typedef std::tr1::array< state_type , order_value > derivs_type;
+
+ order_type order( void ) const
+ {
+ return order_value;
+ }
+
+ order_type stepper_order( void ) const
+ {
+ return order_value;
+ }
+
+ order_type error_order( void ) const
+ {
+ return order_value - 1;
+ }
+
+
+
+
+ template< class System >
+ void do_step( System sys , state_type &x , time_type t , time_type dt )
+ {
+ BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+ BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+
+ do_step( sys , x , t , x , dt );
+ }
+
+ template< class System >
+ void do_step( System sys , const state_type &in , time_type t , state_type &out , time_type dt )
+ {
+ BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+ BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+
+ for( size_t i=0 ; i<Order ; ++i )
+ {
+ boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , m_derivs , i , dt ) );
+ }
+
+ for( size_t i=0 ; i<N ; ++i )
+ {
+ out[i] = in[i];
+ for( size_t k=0 ; k<order_value ; ++k )
+ out[i] += m_derivs[k][i];
+ }
+ }
+
+ template< class System >
+ void do_step( System sys , state_type &x , time_type t , time_type dt , state_type &xerr )
+ {
+ BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+ BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+ BOOST_STATIC_ASSERT(( order_value > 1 ));
+
+ do_step( sys , x , t , x , dt , xerr );
+ }
+
+ template< class System >
+ void do_step( System sys , const state_type &in , time_type t , state_type &out , time_type dt , state_type &xerr )
+ {
+ BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+ BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+ BOOST_STATIC_ASSERT(( order_value > 1 ));
+
+ for( size_t i=0 ; i<Order ; ++i )
+ {
+ boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , m_derivs , i , dt ) );
+ }
+
+ for( size_t i=0 ; i<N ; ++i )
+ {
+ out[i] = in[i];
+ for( size_t k=0 ; k<order_value ; ++k )
+ out[i] += m_derivs[k][i];
+ xerr[i] = m_derivs[order_value-1][i];
+ }
+ }
+
+ const derivs_type& get_last_derivs( void ) const
+ {
+ return m_derivs;
+ }
+
+
+
+private:
+
+ derivs_type m_derivs;
+
+};
+
+
+
+namespace taylor_adf2
+{
+
+ namespace proto = boost::proto;
+
+ template<int I> struct placeholder {};
+
+ proto::terminal< placeholder< 0 > >::type const arg1 = {{}};
+ proto::terminal< placeholder< 1 > >::type const arg2 = {{}};
+ proto::terminal< placeholder< 2 > >::type const arg3 = {{}};
+
+
+
+ template< class State , size_t Order >
+ struct context_derivs : proto::callable_context< context_derivs< State , Order > const >
+ {
+ typedef double result_type;
+
+ const State &m_x;
+ std::tr1::array< State , Order > &m_derivs;
+ size_t m_which;
+
+ context_derivs( const State &x , std::tr1::array< State , Order > &derivs , size_t which )
+ : m_x( x ) , m_derivs( derivs ) , m_which( which ) { }
+
+
+
+ template< int I >
+ double operator()( proto::tag::terminal , placeholder<I> ) const
+ {
+ return ( m_which == 0 ) ? m_x[I] : m_derivs[m_which-1][I];
+ }
+
+ double operator()( proto::tag::terminal , double x ) const
+ {
+ return ( m_which == 0 ) ? x : 0.0;
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::plus , L const &l , R const &r ) const
+ {
+ return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) +
+ proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::minus , L const &l , R const &r ) const
+ {
+ return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) -
+ proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::multiplies , L const &l , R const &r ) const
+ {
+ double tmp = 0.0;
+ for( size_t k=0 ; k<=m_which ; ++k )
+ {
+ tmp += boost::math::binomial_coefficient< double >( m_which , k )
+ * proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , k ) )
+ * proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which - k ) );
+ }
+ return tmp;
+ }
+
+
+ template< typename L , typename R >
+ double operator ()( proto::tag::divides , L const &l , R const &r ) const
+ {
+ return 0.0;
+ }
+
+ };
+
+ template< class System , class State , size_t Order >
+ struct eval_derivs
+ {
+ typedef std::tr1::array< State , Order > derivs_type;
+
+ System m_sys;
+ const State &m_x;
+ derivs_type &m_derivs;
+ size_t m_which;
+
+ eval_derivs( System sys , const State &x , derivs_type &derivs , size_t which )
+ : m_sys( sys ) , m_x( x ) , m_derivs( derivs ) , m_which( which ) { }
+
+ template< class Index >
+ void operator()( Index )
+ {
+ m_derivs[m_which][ Index::value ] = boost::proto::eval(
+ boost::fusion::at< Index >( m_sys ) ,
+ taylor_adf::context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ }
+ };
+
+ template< class System , class State , size_t Order >
+ eval_derivs< System , State , Order > make_eval_derivs( System sys , const State &x , std::tr1::array< State , Order > &derivs , size_t i )
+ {
+ return eval_derivs< System , State , Order >( sys , x , derivs , i );
+ }
+
+
+
+}
+
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif /* BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp 2011-04-08 09:23:41 EDT (Fri, 08 Apr 2011)
@@ -0,0 +1,83 @@
+/*
+ * main.cpp
+ *
+ * Created on: Apr 2, 2011
+ * Author: karsten
+ */
+
+#include <iostream>
+
+#include "taylor.hpp"
+
+#include <boost/fusion/include/make_vector.hpp>
+
+template< typename T , size_t N >
+std::ostream& operator<<( std::ostream& out , const std::tr1::array< T , N > &x )
+{
+ if( !x.empty() ) out << x[0];
+ for( size_t i=1 ; i<x.size() ; ++i ) out << "\t" << x[i];
+ return out;
+}
+
+typedef boost::numeric::odeint::taylor< 3 , 10 > taylor_type;
+typedef taylor_type::state_type state_type;
+typedef taylor_type::derivs_type derivs_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];
+}
+
+
+
+
+
+namespace fusion = boost::fusion;
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+using boost::numeric::odeint::taylor_adf::arg1;
+using boost::numeric::odeint::taylor_adf::arg2;
+using boost::numeric::odeint::taylor_adf::arg3;
+
+int main( int argc , char **argv )
+{
+ cout.precision( 14 );
+
+ taylor_type t;
+
+ state_type in = {{ 10.0 , 10.0 , 10.0 }} , dxdt = {{ 0.0 , 0.0 , 0.0 }} , xerr = {{ 0.0 , 0.0 , 0.0 }} , out = {{ 0.0 ,0.0 , 0.0 }};
+
+ lorenz( in , dxdt , 0.0 );
+
+ cout << in << endl;
+ cout << dxdt << endl << endl;
+
+ t.do_step(
+ fusion::make_vector
+ (
+ sigma * ( arg2 - arg1 ) ,
+ R * arg1 - arg2 - arg1 * arg3 ,
+ arg1 * arg2 - b * arg3
+ ) ,
+ in , 0.0 , out , 0.1 , xerr );
+
+
+
+ cout << in << endl;
+ cout << dxdt << endl;
+ cout << xerr << endl;
+ cout << out << endl << endl;
+ const derivs_type &derivs = t.get_last_derivs();
+ for( size_t i=0 ; i<derivs.size() ; ++i )
+ cout << derivs[i] << endl;
+
+ return 0;
+}