|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r70937 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor
From: karsten.ahnert_at_[hidden]
Date: 2011-04-03 11:13:02
Author: karsten
Date: 2011-04-03 11:13:01 EDT (Sun, 03 Apr 2011)
New Revision: 70937
URL: http://svn.boost.org/trac/boost/changeset/70937
Log:
Ideas taylor added, a ode solver using the taylor expansion of the rhs of the ode
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp (contents, props changed)
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Jamfile 2011-04-03 11:13:01 EDT (Sun, 03 Apr 2011)
@@ -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
+ <define>BOOST_ALL_NO_LIB=1
+ <include>../../../../..
+ ;
+
+exe taylor_main : taylor_main.cpp ;
+exe taylor_lorenz_evol : taylor_lorenz_evol.cpp ;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb 2011-04-03 11:13:01 EDT (Sun, 03 Apr 2011)
@@ -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 *)
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor.hpp 2011-04-03 11:13:01 EDT (Sun, 03 Apr 2011)
@@ -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 <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 = 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<Order ; ++i )
+ {
+ boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , x , m_derivs , i , dt ) );
+ }
+
+ for( size_t i=0 ; i<N ; ++i )
+ {
+ for( size_t k=0 ; k<order ; ++k )
+ x[i] += m_derivs[k][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_lorenz_evol.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp 2011-04-03 11:13:01 EDT (Sun, 03 Apr 2011)
@@ -0,0 +1,79 @@
+/*
+ * 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;
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_main.cpp 2011-04-03 11:13:01 EDT (Sun, 03 Apr 2011)
@@ -0,0 +1,76 @@
+/*
+ * 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 )
+{
+ 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<derivs.size() ; ++i )
+ cout << derivs[i] << endl;
+
+ return 0;
+}
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk