From 73d76a49ae4a430dfc3df875d5836f79a04d41bd Mon Sep 17 00:00:00 2001 From: karsten Date: Sun, 3 Apr 2011 15:13:01 +0000 Subject: [PATCH] Ideas taylor added, a ode solver using the taylor expansion of the rhs of the ode --- libs/numeric/odeint/ideas/taylor/Jamfile | 14 + libs/numeric/odeint/ideas/taylor/lorenz.nb | 353 ++++++++++++++++++ libs/numeric/odeint/ideas/taylor/taylor.hpp | 314 ++++++++++++++++ .../ideas/taylor/taylor_lorenz_evol.cpp | 79 ++++ .../odeint/ideas/taylor/taylor_main.cpp | 76 ++++ 5 files changed, 836 insertions(+) create mode 100644 libs/numeric/odeint/ideas/taylor/Jamfile create mode 100644 libs/numeric/odeint/ideas/taylor/lorenz.nb create mode 100644 libs/numeric/odeint/ideas/taylor/taylor.hpp create mode 100644 libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp create mode 100644 libs/numeric/odeint/ideas/taylor/taylor_main.cpp diff --git a/libs/numeric/odeint/ideas/taylor/Jamfile b/libs/numeric/odeint/ideas/taylor/Jamfile new file mode 100644 index 00000000..c8a326fe --- /dev/null +++ b/libs/numeric/odeint/ideas/taylor/Jamfile @@ -0,0 +1,14 @@ +# (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 + BOOST_ALL_NO_LIB=1 + ../../../../.. + ; + +exe taylor_main : taylor_main.cpp ; +exe taylor_lorenz_evol : taylor_lorenz_evol.cpp ; diff --git a/libs/numeric/odeint/ideas/taylor/lorenz.nb b/libs/numeric/odeint/ideas/taylor/lorenz.nb new file mode 100644 index 00000000..333cf1fc --- /dev/null +++ b/libs/numeric/odeint/ideas/taylor/lorenz.nb @@ -0,0 +1,353 @@ +(* Content-type: application/mathematica *) + +(*** Wolfram Notebook File ***) +(* http://www.wolfram.com/nb *) + +(* CreatedBy='Mathematica 7.0' *) + +(*CacheID: 234*) +(* Internal cache information: +NotebookFileLineBreakTest +NotebookFileLineBreakTest +NotebookDataPosition[ 145, 7] +NotebookDataLength[ 10200, 344] +NotebookOptionsPosition[ 9018, 300] +NotebookOutlinePosition[ 9356, 315] +CellTagsIndexPosition[ 9313, 312] +WindowFrame->Normal*) + +(* Beginning of Notebook Content *) +Notebook[{ +Cell[BoxData[{ + RowBox[{ + RowBox[{"\[Sigma]", "=", "10"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"R", "=", "28"}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"b", "=", + RowBox[{"8", "/", "3"}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"fx", "[", "t_", "]"}], ":=", + RowBox[{"\[Sigma]", + RowBox[{"(", + RowBox[{ + RowBox[{"y", "[", "t", "]"}], "-", + RowBox[{"x", "[", "t", "]"}]}], ")"}]}]}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"fy", "[", "t_", "]"}], ":=", + RowBox[{ + RowBox[{"R", " ", + RowBox[{"x", "[", "t", "]"}]}], "-", + RowBox[{"y", "[", "t", "]"}], "-", + RowBox[{ + RowBox[{"x", "[", "t", "]"}], " ", + RowBox[{"z", "[", "t", "]"}]}]}]}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{"fz", "[", "t_", "]"}], ":=", + RowBox[{ + RowBox[{ + RowBox[{"x", "[", "t", "]"}], + RowBox[{"y", "[", "t", "]"}]}], "-", + RowBox[{"b", " ", + RowBox[{"z", "[", "t", "]"}]}]}]}]}], "Input", + CellChangeTimes->{{3.510824673195559*^9, 3.5108246955668364`*^9}, { + 3.510824728655014*^9, 3.510824774844098*^9}}], + +Cell[BoxData[ + RowBox[{ + RowBox[{"vec", "=", + RowBox[{"{", + RowBox[{ + RowBox[{"fx", "[", "t", "]"}], ",", + RowBox[{"fy", "[", "t", "]"}], ",", + RowBox[{"fz", "[", "t", "]"}]}], "}"}]}], ";"}]], "Input", + CellChangeTimes->{{3.5108247955762157`*^9, 3.510824808339299*^9}}], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"dx", "=", + RowBox[{ + RowBox[{ + RowBox[{"vec", "/.", + RowBox[{ + RowBox[{"x", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{"y", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{"z", "[", "t", "]"}], "\[Rule]", "10"}]}]}]], "Input", + CellChangeTimes->{{3.5108248103120823`*^9, 3.510824843576672*^9}}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"0", ",", "170", ",", + FractionBox["220", "3"]}], "}"}]], "Output", + CellChangeTimes->{3.510824811744413*^9, 3.510824844343238*^9}] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"ddx", "=", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{"D", "[", + RowBox[{"vec", ",", "t"}], "]"}], "/.", + RowBox[{ + RowBox[{"x", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{"y", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{"z", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"x", "'"}], "[", "t", "]"}], "\[Rule]", "0"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"y", "'"}], "[", "t", "]"}], "\[Rule]", "170"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"z", "'"}], "[", "t", "]"}], "\[Rule]", + RowBox[{"220", "/", "3"}]}]}]}]], "Input", + CellChangeTimes->{{3.510824847695149*^9, 3.510824889430078*^9}}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"1700", ",", + RowBox[{"-", + FractionBox["2710", "3"]}], ",", + FractionBox["13540", "9"]}], "}"}]], "Output", + CellChangeTimes->{{3.510824862798667*^9, 3.510824892253048*^9}}] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"dddx", "=", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{"D", "[", + RowBox[{"vec", ",", + RowBox[{"{", + RowBox[{"t", ",", "2"}], "}"}]}], "]"}], "/.", + RowBox[{ + RowBox[{"x", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{"y", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{"z", "[", "t", "]"}], "\[Rule]", "10"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"x", "'"}], "[", "t", "]"}], "\[Rule]", "0"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"y", "'"}], "[", "t", "]"}], "\[Rule]", "170"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"z", "'"}], "[", "t", "]"}], "\[Rule]", + RowBox[{"220", "/", "3"}]}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"x", "''"}], "[", "t", "]"}], "\[Rule]", "1700"}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"y", "''"}], "[", "t", "]"}], "\[Rule]", + RowBox[{ + RowBox[{"-", "2710"}], "/", "3"}]}]}], "/.", + RowBox[{ + RowBox[{ + RowBox[{"z", "''"}], "[", "t", "]"}], "\[Rule]", + RowBox[{"13540", "/", "9"}]}]}]}]], "Input", + CellChangeTimes->{{3.5108248942284517`*^9, 3.510824926935175*^9}}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{"-", + FractionBox["78100", "3"]}], ",", + FractionBox["148130", "9"], ",", + FractionBox["106780", "27"]}], "}"}]], "Output", + CellChangeTimes->{3.510824928988983*^9}] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{"N", "[", "dx", "]"}], "\[IndentingNewLine]", + RowBox[{"N", "[", "ddx", "]"}], "\[IndentingNewLine]", + RowBox[{"N", "[", "dddx", "]"}]}], "Input", + CellChangeTimes->{{3.510824932150373*^9, 3.510824938416258*^9}}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"0.`", ",", "170.`", ",", "73.33333333333333`"}], "}"}]], "Output", + CellChangeTimes->{3.510824939061275*^9}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"1700.`", ",", + RowBox[{"-", "903.3333333333334`"}], ",", "1504.4444444444443`"}], + "}"}]], "Output", + CellChangeTimes->{3.5108249391281633`*^9}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{"-", "26033.333333333332`"}], ",", "16458.88888888889`", ",", + "3954.814814814815`"}], "}"}]], "Output", + CellChangeTimes->{3.510824939152525*^9}] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{ + RowBox[{"dt", "=", "0.1"}], ";"}], "\[IndentingNewLine]", + RowBox[{"dx1", "=", + RowBox[{ + RowBox[{"1", "/", "10"}], " ", "dx"}]}], "\[IndentingNewLine]", + RowBox[{"ddx1", "=", + RowBox[{ + RowBox[{ + RowBox[{"1", "/", "100"}], "/", "2"}], "ddx"}]}], "\[IndentingNewLine]", + RowBox[{"dddx1", "=", + RowBox[{ + RowBox[{ + RowBox[{"1", "/", "1000"}], "/", "6"}], "dddx"}]}]}], "Input", + CellChangeTimes->{{3.5108251832490597`*^9, 3.510825322790757*^9}, { + 3.510825737885933*^9, 3.510825758949971*^9}, {3.510825802567717*^9, + 3.510825805399754*^9}}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"0", ",", "17", ",", + FractionBox["22", "3"]}], "}"}]], "Output", + CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9}, + 3.510825323504698*^9, 3.510825760025854*^9, 3.510825805994877*^9}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + FractionBox["17", "2"], ",", + RowBox[{"-", + FractionBox["271", "60"]}], ",", + FractionBox["677", "90"]}], "}"}]], "Output", + CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9}, + 3.510825323504698*^9, 3.510825760025854*^9, 3.51082580600082*^9}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{"-", + FractionBox["781", "180"]}], ",", + FractionBox["14813", "5400"], ",", + FractionBox["5339", "8100"]}], "}"}]], "Output", + CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9}, + 3.510825323504698*^9, 3.510825760025854*^9, 3.5108258060061197`*^9}] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[{ + RowBox[{"N", "[", + RowBox[{"dx1", ",", "20"}], "]"}], "\[IndentingNewLine]", + RowBox[{"N", "[", + RowBox[{"ddx1", ",", "20"}], "]"}], "\[IndentingNewLine]", + RowBox[{"N", "[", + RowBox[{"dddx1", ",", "20"}], "]"}]}], "Input", + CellChangeTimes->{{3.5108252915121098`*^9, 3.510825294371749*^9}, { + 3.510825326753573*^9, 3.5108253314518747`*^9}}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + "0", ",", "17.`20.", ",", "7.33333333333333333333333333333333333333`20."}], + "}"}]], "Output", + CellChangeTimes->{3.510825295082189*^9, 3.510825332149527*^9, + 3.510825761579811*^9, 3.510825806943115*^9}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"8.5`20.", ",", + RowBox[{ + "-", "4.5166666666666666666666666666666666666658830037661184749947`20."}], + ",", "7.52222222222222222222222222222222222222`20."}], "}"}]], "Output", + CellChangeTimes->{3.510825295082189*^9, 3.510825332149527*^9, + 3.510825761579811*^9, 3.510825806948196*^9}], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{ + "-", "4.3388888888888888888888888888888888888896072465477247312549`20."}], + ",", "2.7431481481481481481481481481481481481481481481481481481482`20.", + ",", "0.6591358024691358024691358024691358024691358024691358024691`20."}], + "}"}]], "Output", + CellChangeTimes->{3.510825295082189*^9, 3.510825332149527*^9, + 3.510825761579811*^9, 3.510825806952417*^9}] +}, Open ]] +}, +WindowSize->{640, 750}, +WindowMargins->{{148, Automatic}, {Automatic, 27}}, +FrontEndVersion->"7.0 for Linux x86 (64-bit) (November 11, 2008)", +StyleDefinitions->"Default.nb" +] +(* End of Notebook Content *) + +(* Internal cache information *) +(*CellTagsOutline +CellTagsIndex->{} +*) +(*CellTagsIndex +CellTagsIndex->{} +*) +(*NotebookFileOutline +Notebook[{ +Cell[545, 20, 1090, 33, 143, "Input"], +Cell[1638, 55, 292, 8, 32, "Input"], +Cell[CellGroupData[{ +Cell[1955, 67, 384, 11, 32, "Input"], +Cell[2342, 80, 173, 4, 46, "Output"] +}, Open ]], +Cell[CellGroupData[{ +Cell[2552, 89, 830, 26, 55, "Input"], +Cell[3385, 117, 223, 6, 46, "Output"] +}, Open ]], +Cell[CellGroupData[{ +Cell[3645, 128, 1387, 43, 99, "Input"], +Cell[5035, 173, 226, 7, 46, "Output"] +}, Open ]], +Cell[CellGroupData[{ +Cell[5298, 185, 238, 4, 77, "Input"], +Cell[5539, 191, 148, 3, 31, "Output"], +Cell[5690, 196, 189, 5, 31, "Output"], +Cell[5882, 203, 200, 5, 31, "Output"] +}, Open ]], +Cell[CellGroupData[{ +Cell[6119, 213, 591, 16, 99, "Input"], +Cell[6713, 231, 241, 5, 46, "Output"], +Cell[6957, 238, 309, 8, 46, "Output"], +Cell[7269, 248, 322, 8, 46, "Output"] +}, Open ]], +Cell[CellGroupData[{ +Cell[7628, 261, 366, 8, 77, "Input"], +Cell[7997, 271, 249, 6, 31, "Output"], +Cell[8249, 279, 333, 7, 52, "Output"], +Cell[8585, 288, 417, 9, 52, "Output"] +}, Open ]] +} +] +*) + +(* End of internal cache information *) diff --git a/libs/numeric/odeint/ideas/taylor/taylor.hpp b/libs/numeric/odeint/ideas/taylor/taylor.hpp new file mode 100644 index 00000000..d69bf9fd --- /dev/null +++ b/libs/numeric/odeint/ideas/taylor/taylor.hpp @@ -0,0 +1,314 @@ +/* + * taylor.hpp + * + * Created on: Apr 2, 2011 + * Author: karsten + */ + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_ +#define BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_ + +#include + +// general boost includes +#include +#include +#include + +// fusion includes +#include +#include +#include + +// mpl includes +#include +#include + +// proto includes +#include + + + + +namespace boost { +namespace numeric { +namespace odeint { + + +namespace taylor_adf +{ + + namespace proto = boost::proto; + + template 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 ) 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 = Order; + typedef double value_type; + typedef value_type time_type; + + typedef std::tr1::array< value_type , dim > state_type; + typedef std::tr1::array< state_type , order > derivs_type; + + + + template< class System > + void do_step( System sys , state_type &x , value_type t , value_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 >( make_eval_derivs( sys , x , m_derivs , i , dt ) ); + } + + for( size_t i=0 ; 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 ) 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_ */ diff --git a/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp b/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp new file mode 100644 index 00000000..357cb6ed --- /dev/null +++ b/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp @@ -0,0 +1,79 @@ +/* + * main.cpp + * + * Created on: Apr 2, 2011 + * Author: karsten + */ + +#include + +#include "taylor.hpp" + +#include +#include + +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 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; +} diff --git a/libs/numeric/odeint/ideas/taylor/taylor_main.cpp b/libs/numeric/odeint/ideas/taylor/taylor_main.cpp new file mode 100644 index 00000000..f50dd80a --- /dev/null +++ b/libs/numeric/odeint/ideas/taylor/taylor_main.cpp @@ -0,0 +1,76 @@ +/* + * main.cpp + * + * Created on: Apr 2, 2011 + * Author: karsten + */ + +#include + +#include "taylor.hpp" + +#include + +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 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 ) +{ + taylor_type t; + + state_type x = {{ 10.0 , 10.0 , 10.0 }} , dxdt = {{ 0.0 , 0.0 , 0.0 }} , xerr = {{ 0.0 , 0.0 , 0.0 }}; + + lorenz( x , dxdt , 0.0 ); + + cout << x << endl; + cout << dxdt << endl << endl; + + t.do_step( + fusion::make_vector + ( + sigma * ( arg2 - arg1 ) , + R * arg1 - arg2 - arg1 * arg3 , + arg1 * arg2 - b * arg3 + ) , + x , 0.0 , 0.1 ); + + cout.precision( 14 ); + const derivs_type &derivs = t.get_last_derivs(); + for( size_t i=0 ; i