mirror of
https://github.com/boostorg/odeint.git
synced 2026-01-26 06:42:23 +00:00
Ideas taylor added, a ode solver using the taylor expansion of the rhs of the ode
This commit is contained in:
karsten
14
libs/numeric/odeint/ideas/taylor/Jamfile
Normal file
14
libs/numeric/odeint/ideas/taylor/Jamfile
Normal file
@@ -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 ;
|
||||
353
libs/numeric/odeint/ideas/taylor/lorenz.nb
Normal file
353
libs/numeric/odeint/ideas/taylor/lorenz.nb
Normal file
@@ -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 *)
|
||||
314
libs/numeric/odeint/ideas/taylor/taylor.hpp
Normal file
314
libs/numeric/odeint/ideas/taylor/taylor.hpp
Normal file
@@ -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_ */
|
||||
79
libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp
Normal file
79
libs/numeric/odeint/ideas/taylor/taylor_lorenz_evol.cpp
Normal file
@@ -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;
|
||||
}
|
||||
76
libs/numeric/odeint/ideas/taylor/taylor_main.cpp
Normal file
76
libs/numeric/odeint/ideas/taylor/taylor_main.cpp
Normal file
@@ -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;
|
||||
}
|
||||
Reference in New Issue
Block a user