Subject: [Boost-commit] svn:boost r71115 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: . taylor_v1
From: karsten.ahnert_at_[hidden]
Date: 2011-04-08 09:21:14
Author: karsten
Date: 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
New Revision: 71115
URL: http://svn.boost.org/trac/boost/changeset/71115
Log:
rearranging ideas/taylor
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/Jamfile
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor.hpp
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_controller.hpp
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_controller.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_lorenz_evol_adaptive.cpp
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_lorenz_evol_adaptive_statistics.cpp
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive_statistics.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_lorenz_evol_direct.cpp
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_direct.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_main.cpp
- copied unchanged from r71114, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp
Removed:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_controller.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive_statistics.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_direct.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,16 +0,0 @@
-# (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 ;
-exe taylor_lorenz_evol_direct : taylor_lorenz_evol_direct.cpp ;
-exe taylor_lorenz_evol_adaptive : taylor_lorenz_evol_adaptive.cpp ;
-exe taylor_lorenz_evol_adaptive_statistics : taylor_lorenz_evol_adaptive_statistics.cpp ;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
@@ -0,0 +1,20 @@
+Version 2:
+
+* Requirements:
+ * Optimize a * x
+ * Optimize a + x
+ * Optimize x * y
+ * Optimize g(x)
+
+* Ansatz :
+ * mult< terminal< double >( a ) , expr > to mult_scalar< expr )( a )
+ * plus< terminal< double >( a ) , expr > to shift_scalar< expr )( a )
+ * mult< expr1 , expr2 > to mult_leibnitz< expr1 , expr2 >
+ * unary_op< expr > to unary_op_faa_di_bruno_bruno< expr >
+
+* Possible solutions :
+ * Proto transform to create
+ deriv[which] = deriv_transform()( expr , x , derivs , which );
+ * custom eval context to create
+ deriv[which] = proto::eval( expr , my_context( x , derivs , which ) );
+ * start with 1D Odes
\ No newline at end of file
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,372 +0,0 @@
-/*
- * 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_ */
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_controller.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_controller.hpp 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,148 +0,0 @@
-/*
- boost header: NUMERIC_ODEINT_STEPPER/controlled_error_stepper.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_STEPPER_TAYLOR_CONTROLLER_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_CONTROLLER_HPP_INCLUDED
-
-#include <cmath>
-
-#include <boost/numeric/odeint/algebra/range_algebra.hpp>
-#include <boost/numeric/odeint/algebra/default_operations.hpp>
-#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
-#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
-
-
-
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
-/*
- * explicit stepper version
- */
-template
-<
- class TaylorStepper
->
-class taylor_controller
-{
-
-public:
-
- typedef TaylorStepper stepper_type;
- typedef typename stepper_type::state_type state_type;
- typedef typename stepper_type::value_type value_type;
- typedef typename stepper_type::deriv_type deriv_type;
- typedef typename stepper_type::time_type time_type;
- typedef typename stepper_type::order_type order_type;
- typedef default_error_checker< value_type > error_checker_type;
- typedef controlled_stepper_tag stepper_category;
-
-
- taylor_controller(
- value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
- value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
- value_type a_x = static_cast< value_type >( 1.0 ) ,
- value_type a_dxdt = static_cast< value_type >( 1.0 ) ,
- const stepper_type &stepper = stepper_type()
- )
- : m_stepper( stepper ) ,
- m_error_checker( eps_abs , eps_rel , a_x , a_dxdt ) ,
- m_xerr() , m_xnew()
- {
- }
-
-
- template< class System >
- controlled_step_result try_step( System system , state_type &x , time_type &t , time_type &dt )
- {
- controlled_step_result res = try_step( system , x , t , m_xnew , dt );
- if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged ) )
- {
- x = m_xnew;
- }
- return res;
- }
-
-
- template< class System >
- controlled_step_result try_step( System system , state_type &in , time_type &t , state_type &out , time_type &dt )
- {
- using std::max;
- using std::min;
- using std::pow;
-
- // do one step with error calculation
- m_stepper.do_step( system , in , t , out , dt , m_xerr );
-
- value_type max_rel_error = m_error_checker.error( in , m_stepper.get_last_derivs()[0] , m_xerr , dt );
-
- if( max_rel_error > 1.1 )
- {
- // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
- dt *= max( 0.9 * pow( max_rel_error , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
- return step_size_decreased;
- }
- else
- {
- if( max_rel_error < 0.5 )
- {
- //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
- t += dt;
- dt *= min( 0.9 * pow( max_rel_error , -1.0 / m_stepper.stepper_order() ) , 5.0 );
- return success_step_size_increased;
- }
- else
- {
- t += dt;
- return success_step_size_unchanged;
- }
- }
- }
-
-
- stepper_type& stepper( void )
- {
- return m_stepper;
- }
-
- const stepper_type& stepper( void ) const
- {
- return m_stepper;
- }
-
-
-private:
-
-
-
- stepper_type m_stepper;
- error_checker_type m_error_checker;
- state_type m_xerr;
- state_type m_xnew;
-};
-
-
-
-
-
-
-
-
-} // odeint
-} // numeric
-} // boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_CONTROLLER_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive.cpp 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,106 +0,0 @@
-/*
- * main.cpp
- *
- * Created on: Apr 2, 2011
- * Author: karsten
- */
-
-#include <iostream>
-#include <fstream>
-
-#include "taylor.hpp"
-#include "taylor_controller.hpp"
-
-#include <boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp>
-#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
-#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
-
-#include <boost/fusion/include/make_vector.hpp>
-#include <boost/spirit/include/phoenix.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 , 20 > taylor_type;
-typedef taylor_type::state_type state_type;
-typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type > rk54_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;
-namespace phoenix = boost::phoenix;
-
-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;
-
-struct streaming_observer
-{
- std::ostream& m_out;
-
- streaming_observer( std::ostream &out ) : m_out( out ) { }
-
- template< class State >
- void operator()( const State &x , double t ) const
- {
- m_out << t;
- for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i];
- m_out << "\n";
- }
-};
-
-
-int main( int argc , char **argv )
-{
- double eps_abs = 1.0e-6;
- double eps_rel = 1.0e-6;
- if( argc > 2 )
- {
- eps_abs = atof( argv[1] );
- eps_rel = atof( argv[2] );
- }
-
-
- rk54_type rk54_plain;
- controlled_error_stepper< rk54_type > rk54( rk54_plain , default_error_checker< double >( eps_abs , eps_rel ) );
- taylor_controller< taylor_type > taylor_controller( eps_abs , eps_rel );
-
- state_type x1 = {{ 10.0 , 10.0 , 10.0 }} , x2 = x1;
-
- ofstream fout( "lorenz_rk54.dat" );
- size_t steps_rk54 = integrate_adaptive( rk54 , lorenz , x1 , 0.0 , 50.0 , 0.1 , streaming_observer( fout ) );
- clog << "Steps RK 54 : " << steps_rk54 << endl;;
-
- ofstream fout2( "lorenz_taylor.dat" );
- size_t steps_taylor = integrate_adaptive( taylor_controller ,
- fusion::make_vector
- (
- sigma * ( arg2 - arg1 ) ,
- R * arg1 - arg2 - arg1 * arg3 ,
- arg1 * arg2 - b * arg3
- ) , x2 , 0.0 , 50.0 , 0.1 , streaming_observer( fout2 ) );
- clog << "Steps Taylor : " << steps_taylor << endl;
-
- return 0;
-}
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive_statistics.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_adaptive_statistics.cpp 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,166 +0,0 @@
-/*
- * main.cpp
- *
- * Created on: Apr 2, 2011
- * Author: karsten
- */
-
-#include <iostream>
-#include <fstream>
-
-#include "taylor.hpp"
-#include "taylor_controller.hpp"
-
-#include <boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp>
-#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
-#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
-
-#include <boost/fusion/include/make_vector.hpp>
-#include <boost/spirit/include/phoenix.hpp>
-#include <boost/assign.hpp>
-#include <boost/timer.hpp>
-
-#include <boost/mpl/range_c.hpp>
-#include <boost/mpl/for_each.hpp>
-
-#define tab "\t"
-
-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 std::tr1::array< double , 3 > state_type;
-
-typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type > rk54_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;
-namespace phoenix = boost::phoenix;
-namespace mpl = boost::mpl;
-
-using namespace std;
-using namespace boost::numeric::odeint;
-using namespace boost::assign;
-
-using boost::numeric::odeint::taylor_adf::arg1;
-using boost::numeric::odeint::taylor_adf::arg2;
-using boost::numeric::odeint::taylor_adf::arg3;
-
-struct run
-{
- vector< double > m_eps_abs_values;
- vector< double > m_eps_rel_values;
- double m_t_end;
-
- run( const vector< double > &eps_abs_values , const vector< double > &eps_rel_values , double t_end )
- : m_eps_abs_values( eps_abs_values ) , m_eps_rel_values( eps_rel_values ) , m_t_end( t_end ){ }
-
- template< class Order >
- void operator()( Order ) const
- {
- static const size_t order = Order::value;
- typedef boost::numeric::odeint::taylor< 3 , order > taylor_type;
-
- boost::timer timer;
-
- char str[512] = "";
- sprintf( str , "dat/lorenz_taylor_%02d.dat" , int( order ) );
- ofstream fout( str );
-
- for( size_t i=0 ; i<m_eps_abs_values.size() ; ++i )
- {
- for( size_t j=0 ; j<m_eps_rel_values.size() ; ++j )
- {
- double eps_abs = m_eps_abs_values[i];
- double eps_rel = m_eps_rel_values[j];
-
- taylor_controller< taylor_type > taylor_controller( eps_abs , eps_rel );
-
- state_type x = {{ 10.0 , 10.0 , 10.0 }};
-
- timer.restart();
- size_t steps_taylor = integrate_adaptive( taylor_controller ,
- fusion::make_vector
- (
- sigma * ( arg2 - arg1 ) ,
- R * arg1 - arg2 - arg1 * arg3 ,
- arg1 * arg2 - b * arg3
- ) , x , 0.0 , m_t_end , 0.1 );
- double time_taylor = timer.elapsed();
-
- fout << i << tab << j << tab << eps_abs << tab << eps_rel << tab;
- fout << steps_taylor << tab << time_taylor;
- fout << endl;
- }
- fout << endl;
- }
-
-
- }
-};
-
-
-
-
-int main( int argc , char **argv )
-{
- if( argc != 2 )
- {
- cerr << "usage taylor_lorenz_eval_adaptive_statistics t_end" << endl;
- return -1;
- }
- double t_end = atof( argv[1] );
-
- vector< double > eps_abs_values;
- vector< double > eps_rel_values;
-
- eps_abs_values += 1.0e1 , 1.0 , 1.0e-1 , 1.0e-2 , 1.0e-3 , 1.0e-4 , 1.0e-5 , 1.0e-6 , 1.0e-7 , 1.0e-8 , 1.0e-9 , 1.0e-10 , 1.0e-11 , 1.0e-12 , 1.0e-13 , 1.0e-14;
- eps_rel_values += 1.0e-1 , 1.0e-2 , 1.0e-3 , 1.0e-4 , 1.0e-5 , 1.0e-6 , 1.0e-7 , 1.0e-8 , 1.0e-9 , 1.0e-10 , 1.0e-11 , 1.0e-12 , 1.0e-13 , 1.0e-14;
-
-
- typedef mpl::range_c< size_t , 5 , 30 > order_values;
- mpl::for_each< order_values >( run( eps_abs_values , eps_rel_values , t_end ) );
-
-
- boost::timer timer;
- ofstream fout( "dat/lorenz_rk54.dat" );
- for( size_t i=0 ; i<eps_abs_values.size() ; ++i )
- {
- for( size_t j=0 ; j<eps_rel_values.size() ; ++j )
- {
- double eps_abs = eps_abs_values[i];
- double eps_rel = eps_rel_values[j];
-
- rk54_type rk54_plain;
- controlled_error_stepper< rk54_type > rk54( rk54_plain , default_error_checker< double >( eps_abs , eps_rel ) );
-
- state_type x = {{ 10.0 , 10.0 , 10.0 }};
-
- timer.restart();
- size_t steps_rk54 = integrate_adaptive( rk54 , lorenz , x , 0.0 , t_end , 0.1 );
- double time_rk54 = timer.elapsed();
-
- fout << i << tab << j << tab << eps_abs << tab << eps_rel << tab;
- fout << steps_rk54 << tab << time_rk54 << tab;
- fout << endl;
- }
- fout << endl;
- }
-
- return 0;
-}
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_direct.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol_direct.cpp 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,80 +0,0 @@
-/*
- * main.cpp
- *
- * Created on: Apr 2, 2011
- * Author: karsten
- */
-
-#include <iostream>
-
-#include "taylor.hpp"
-
-#include <boost/numeric/odeint/stepper/explicit_rk4.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 , 11 > taylor_type;
-typedef taylor_type::state_type state_type;
-typedef boost::numeric::odeint::explicit_rk4< state_type > rk4_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 )
-{
- taylor_type tay;
- rk4_type rk4;
-
-
- state_type x1 = {{ 10.0 , 10.0 , 10.0 }} , x2 = x1;
-
- const double dt = 0.1;
- double t = 0.0;
- for( size_t i=0 ; i<10000 ; ++i , t += dt )
- {
- tay.do_step(
- fusion::make_vector
- (
- sigma * ( arg2 - arg1 ) ,
- R * arg1 - arg2 - arg1 * arg3 ,
- arg1 * arg2 - b * arg3
- ) ,
- x1 , t , dt );
-
- rk4.do_step( lorenz , x2 , t , dt );
-
- cout << t << "\t" << x1 << "\t" << x2 << "\n";
- }
-
-
- return 0;
-}
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp 2011-04-08 09:21:12 EDT (Fri, 08 Apr 2011)
+++ (empty file)
@@ -1,83 +0,0 @@
-/*
- * 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;
-}