Subject: [Boost-commit] svn:boost r71848 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: . lorenz_special taylor_v4
From: karsten.ahnert_at_[hidden]
Date: 2011-05-09 14:44:00
Author: karsten
Date: 2011-05-09 14:43:57 EDT (Mon, 09 May 2011)
New Revision: 71848
URL: http://svn.boost.org/trac/boost/changeset/71848
Log:
continuing taylor, cpp code is as fast as the fortran code, hehe
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb | 253 +++++++++++++++++++++++++++++----------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f | 10
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp | 98 +++++++++------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp | 42 +++---
4 files changed, 275 insertions(+), 128 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz.nb 2011-05-09 14:43:57 EDT (Mon, 09 May 2011)
@@ -10,10 +10,10 @@
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[ 145, 7]
-NotebookDataLength[ 10200, 344]
-NotebookOptionsPosition[ 9018, 300]
-NotebookOutlinePosition[ 9356, 315]
-CellTagsIndexPosition[ 9313, 312]
+NotebookDataLength[ 14490, 465]
+NotebookOptionsPosition[ 13070, 414]
+NotebookOutlinePosition[ 13407, 429]
+CellTagsIndexPosition[ 13364, 426]
WindowFrame->Normal*)
(* Beginning of Notebook Content *)
@@ -82,7 +82,8 @@
RowBox[{"{",
RowBox[{"0", ",", "170", ",",
FractionBox["220", "3"]}], "}"}]], "Output",
- CellChangeTimes->{3.510824811744413*^9, 3.510824844343238*^9}]
+ CellChangeTimes->{3.510824811744413*^9, 3.510824844343238*^9,
+ 3.513935237148198*^9}]
}, Open ]],
Cell[CellGroupData[{
@@ -121,7 +122,8 @@
RowBox[{"-",
FractionBox["2710", "3"]}], ",",
FractionBox["13540", "9"]}], "}"}]], "Output",
- CellChangeTimes->{{3.510824862798667*^9, 3.510824892253048*^9}}]
+ CellChangeTimes->{{3.510824862798667*^9, 3.510824892253048*^9},
+ 3.513935237968691*^9}]
}, Open ]],
Cell[CellGroupData[{
@@ -178,7 +180,80 @@
FractionBox["78100", "3"]}], ",",
FractionBox["148130", "9"], ",",
FractionBox["106780", "27"]}], "}"}]], "Output",
- CellChangeTimes->{3.510824928988983*^9}]
+ CellChangeTimes->{3.510824928988983*^9, 3.513935238940962*^9}]
+}, Open ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"ddddx", "=",
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{
+ RowBox[{"D", "[",
+ RowBox[{"vec", ",",
+ RowBox[{"{",
+ RowBox[{"t", ",", "3"}], "}"}]}], "]"}], "/.",
+ 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"}]}]}], "/.",
+ RowBox[{
+ RowBox[{
+ RowBox[{"x", "'''"}], "[", "t", "]"}], "\[Rule]",
+ RowBox[{
+ RowBox[{"-", "78100"}], "/", "3"}]}]}], "/.",
+ RowBox[{
+ RowBox[{
+ RowBox[{"y", "'''"}], "[", "t", "]"}], "\[Rule]",
+ RowBox[{"148130", "/", "9"}]}]}], "/.",
+ RowBox[{
+ RowBox[{
+ RowBox[{"z", "'''"}], "[", "t", "]"}], "\[Rule]",
+ RowBox[{"106780", "/", "27"}]}]}]}]], "Input",
+ CellChangeTimes->{{3.513935187550949*^9, 3.513935230954701*^9}}],
+
+Cell[BoxData[
+ RowBox[{"{",
+ RowBox[{
+ FractionBox["3824300", "9"], ",",
+ RowBox[{"-",
+ FractionBox["24262390", "27"]}], ",",
+ FractionBox["61617460", "81"]}], "}"}]], "Output",
+ CellChangeTimes->{{3.513935233233439*^9, 3.513935239844738*^9}}]
}, Open ]],
Cell[CellGroupData[{
@@ -186,75 +261,102 @@
Cell[BoxData[{
RowBox[{"N", "[", "dx", "]"}], "\[IndentingNewLine]",
RowBox[{"N", "[", "ddx", "]"}], "\[IndentingNewLine]",
- RowBox[{"N", "[", "dddx", "]"}]}], "Input",
- CellChangeTimes->{{3.510824932150373*^9, 3.510824938416258*^9}}],
+ RowBox[{"N", "[", "dddx", "]"}], "\[IndentingNewLine]",
+ RowBox[{"N", "[", "ddddx", "]"}]}], "Input",
+ CellChangeTimes->{{3.510824932150373*^9, 3.510824938416258*^9}, {
+ 3.513935242881462*^9, 3.513935243956193*^9}}],
Cell[BoxData[
RowBox[{"{",
RowBox[{"0.`", ",", "170.`", ",", "73.33333333333333`"}], "}"}]], "Output",
- CellChangeTimes->{3.510824939061275*^9}],
+ CellChangeTimes->{3.510824939061275*^9, 3.5139352446122847`*^9}],
Cell[BoxData[
RowBox[{"{",
RowBox[{"1700.`", ",",
RowBox[{"-", "903.3333333333334`"}], ",", "1504.4444444444443`"}],
"}"}]], "Output",
- CellChangeTimes->{3.5108249391281633`*^9}],
+ CellChangeTimes->{3.510824939061275*^9, 3.513935244616696*^9}],
Cell[BoxData[
RowBox[{"{",
RowBox[{
RowBox[{"-", "26033.333333333332`"}], ",", "16458.88888888889`", ",",
"3954.814814814815`"}], "}"}]], "Output",
- CellChangeTimes->{3.510824939152525*^9}]
+ CellChangeTimes->{3.510824939061275*^9, 3.513935244620942*^9}],
+
+Cell[BoxData[
+ RowBox[{"{",
+ RowBox[{"424922.22222222225`", ",",
+ RowBox[{"-", "898607.0370370371`"}], ",", "760709.3827160494`"}],
+ "}"}]], "Output",
+ CellChangeTimes->{3.510824939061275*^9, 3.513935244625566*^9}]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[{
RowBox[{
- RowBox[{"dt", "=", "0.1"}], ";"}], "\[IndentingNewLine]",
+ RowBox[{"dt", "=", "1"}], ";"}], "\[IndentingNewLine]",
RowBox[{"dx1", "=",
- RowBox[{
- RowBox[{"1", "/", "10"}], " ", "dx"}]}], "\[IndentingNewLine]",
+ RowBox[{"dt", " ", "dx"}]}], "\[IndentingNewLine]",
RowBox[{"ddx1", "=",
RowBox[{
RowBox[{
- RowBox[{"1", "/", "100"}], "/", "2"}], "ddx"}]}], "\[IndentingNewLine]",
+ RowBox[{"dt", "^", "2"}], "/", "2"}], "ddx"}]}], "\[IndentingNewLine]",
RowBox[{"dddx1", "=",
RowBox[{
RowBox[{
- RowBox[{"1", "/", "1000"}], "/", "6"}], "dddx"}]}]}], "Input",
+ RowBox[{"dt", "^", "3"}], "/",
+ RowBox[{"3", "!"}]}], "dddx"}]}], "\[IndentingNewLine]",
+ RowBox[{"ddddx1", "=",
+ RowBox[{
+ RowBox[{
+ RowBox[{"dt", "^", "4"}], "/",
+ RowBox[{"4", "!"}]}], "ddddx"}]}]}], "Input",
CellChangeTimes->{{3.5108251832490597`*^9, 3.510825322790757*^9}, {
3.510825737885933*^9, 3.510825758949971*^9}, {3.510825802567717*^9,
- 3.510825805399754*^9}}],
+ 3.510825805399754*^9}, {3.5139352586584597`*^9, 3.513935304517981*^9}}],
Cell[BoxData[
RowBox[{"{",
- RowBox[{"0", ",", "17", ",",
- FractionBox["22", "3"]}], "}"}]], "Output",
+ RowBox[{"0", ",", "170", ",",
+ FractionBox["220", "3"]}], "}"}]], "Output",
CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9},
- 3.510825323504698*^9, 3.510825760025854*^9, 3.510825805994877*^9}],
+ 3.510825323504698*^9, 3.510825760025854*^9, 3.510825805994877*^9, {
+ 3.5139352825257673`*^9, 3.5139353051164513`*^9}}],
+
+Cell[BoxData[
+ RowBox[{"{",
+ RowBox[{"850", ",",
+ RowBox[{"-",
+ FractionBox["1355", "3"]}], ",",
+ FractionBox["6770", "9"]}], "}"}]], "Output",
+ CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9},
+ 3.510825323504698*^9, 3.510825760025854*^9, 3.510825805994877*^9, {
+ 3.5139352825257673`*^9, 3.513935305122323*^9}}],
Cell[BoxData[
RowBox[{"{",
RowBox[{
- FractionBox["17", "2"], ",",
RowBox[{"-",
- FractionBox["271", "60"]}], ",",
- FractionBox["677", "90"]}], "}"}]], "Output",
+ FractionBox["39050", "9"]}], ",",
+ FractionBox["74065", "27"], ",",
+ FractionBox["53390", "81"]}], "}"}]], "Output",
CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9},
- 3.510825323504698*^9, 3.510825760025854*^9, 3.51082580600082*^9}],
+ 3.510825323504698*^9, 3.510825760025854*^9, 3.510825805994877*^9, {
+ 3.5139352825257673`*^9, 3.5139353051274033`*^9}}],
Cell[BoxData[
RowBox[{"{",
RowBox[{
+ FractionBox["956075", "54"], ",",
RowBox[{"-",
- FractionBox["781", "180"]}], ",",
- FractionBox["14813", "5400"], ",",
- FractionBox["5339", "8100"]}], "}"}]], "Output",
+ FractionBox["12131195", "324"]}], ",",
+ FractionBox["15404365", "486"]}], "}"}]], "Output",
CellChangeTimes->{{3.510825209963582*^9, 3.5108252877049*^9},
- 3.510825323504698*^9, 3.510825760025854*^9, 3.5108258060061197`*^9}]
+ 3.510825323504698*^9, 3.510825760025854*^9, 3.510825805994877*^9, {
+ 3.5139352825257673`*^9, 3.5139353051323833`*^9}}]
}, Open ]],
Cell[CellGroupData[{
@@ -265,41 +367,53 @@
RowBox[{"N", "[",
RowBox[{"ddx1", ",", "20"}], "]"}], "\[IndentingNewLine]",
RowBox[{"N", "[",
- RowBox[{"dddx1", ",", "20"}], "]"}]}], "Input",
+ RowBox[{"dddx1", ",", "20"}], "]"}], "\[IndentingNewLine]",
+ RowBox[{"N", "[",
+ RowBox[{"ddddx1", ",", "20"}], "]"}]}], "Input",
CellChangeTimes->{{3.5108252915121098`*^9, 3.510825294371749*^9}, {
- 3.510825326753573*^9, 3.5108253314518747`*^9}}],
+ 3.510825326753573*^9, 3.5108253314518747`*^9}, {3.513935286019421*^9,
+ 3.513935289521921*^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}],
+ "0", ",", "170.`20.", ",",
+ "73.33333333333333333333333333333333333333`20."}], "}"}]], "Output",
+ CellChangeTimes->{
+ 3.510825295082189*^9, 3.510825332149527*^9, 3.510825761579811*^9,
+ 3.510825806943115*^9, {3.513935290210574*^9, 3.513935306844468*^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}],
+ RowBox[{"850.`20.", ",",
+ RowBox[{"-", "451.66666666666666666666666666666666666667`20."}], ",",
+ "752.22222222222222222222222222222222222222`20."}], "}"}]], "Output",
+ CellChangeTimes->{
+ 3.510825295082189*^9, 3.510825332149527*^9, 3.510825761579811*^9,
+ 3.510825806943115*^9, {3.513935290210574*^9, 3.513935306849414*^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}]
+ RowBox[{"-", "4338.88888888888888888888888888888888888889`20."}], ",",
+ "2743.14814814814814814814814814814814814815`20.", ",",
+ "659.13580246913580246913580246913580246914`20."}], "}"}]], "Output",
+ CellChangeTimes->{
+ 3.510825295082189*^9, 3.510825332149527*^9, 3.510825761579811*^9,
+ 3.510825806943115*^9, {3.513935290210574*^9, 3.513935306853644*^9}}],
+
+Cell[BoxData[
+ RowBox[{"{",
+ RowBox[{"17705.09259259259259259259259259259259259259`20.", ",",
+ RowBox[{"-", "37441.95987654320987654320987654320987654321`20."}], ",",
+ "31696.22427983539094650205761316872427983539`20."}], "}"}]], "Output",
+ CellChangeTimes->{
+ 3.510825295082189*^9, 3.510825332149527*^9, 3.510825761579811*^9,
+ 3.510825806943115*^9, {3.513935290210574*^9, 3.513935306892515*^9}}]
}, Open ]]
},
WindowSize->{640, 750},
-WindowMargins->{{148, Automatic}, {Automatic, 27}},
+WindowMargins->{{146, Automatic}, {Automatic, 2}},
FrontEndVersion->"7.0 for Linux x86 (64-bit) (November 11, 2008)",
StyleDefinitions->"Default.nb"
]
@@ -318,33 +432,40 @@
Cell[1638, 55, 292, 8, 32, "Input"],
Cell[CellGroupData[{
Cell[1955, 67, 384, 11, 32, "Input"],
-Cell[2342, 80, 173, 4, 46, "Output"]
+Cell[2342, 80, 198, 5, 46, "Output"]
+}, Open ]],
+Cell[CellGroupData[{
+Cell[2577, 90, 830, 26, 32, "Input"],
+Cell[3410, 118, 249, 7, 46, "Output"]
}, Open ]],
Cell[CellGroupData[{
-Cell[2552, 89, 830, 26, 55, "Input"],
-Cell[3385, 117, 223, 6, 46, "Output"]
+Cell[3696, 130, 1387, 43, 32, "Input"],
+Cell[5086, 175, 248, 7, 46, "Output"]
}, Open ]],
Cell[CellGroupData[{
-Cell[3645, 128, 1387, 43, 99, "Input"],
-Cell[5035, 173, 226, 7, 46, "Output"]
+Cell[5371, 187, 1958, 59, 143, "Input"],
+Cell[7332, 248, 257, 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"]
+Cell[7626, 260, 346, 6, 99, "Input"],
+Cell[7975, 268, 172, 3, 31, "Output"],
+Cell[8150, 273, 209, 5, 31, "Output"],
+Cell[8362, 280, 222, 5, 31, "Output"],
+Cell[8587, 287, 221, 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"]
+Cell[8845, 297, 777, 21, 121, "Input"],
+Cell[9625, 320, 297, 6, 46, "Output"],
+Cell[9925, 328, 341, 8, 46, "Output"],
+Cell[10269, 338, 371, 9, 46, "Output"],
+Cell[10643, 349, 381, 9, 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"]
+Cell[11061, 363, 499, 11, 99, "Input"],
+Cell[11563, 376, 302, 7, 31, "Output"],
+Cell[11868, 385, 364, 7, 52, "Output"],
+Cell[12235, 394, 408, 8, 31, "Output"],
+Cell[12646, 404, 408, 7, 31, "Output"]
}, Open ]]
}
]
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f 2011-05-09 14:43:57 EDT (Mon, 09 May 2011)
@@ -3,6 +3,8 @@
stop
end
+
+
Subroutine INTEGRATION
IMPLICIT REAL*8 (A-H,O-Z)
real*8 AX(77777),AY(77777),AZ(77777),AT(77777)
@@ -16,6 +18,8 @@
COMMON /ST/ DX,DY,DZ,Q,NO
INTEGER counter
+ open(9,file='lorenz.dat',status='unknown',access='sequential',
+ : form='formatted')
C Introduce the parameters:
P=10.d0
B=8.d0/3.d0
@@ -27,11 +31,9 @@
t=0.d0
tend=5000.d0
counter=0
- 10 xp=x
- yp=y
- zp=z
- call step (P,R,B,X,Y,Z,DT)
+10 call step (P,R,B,X,Y,Z,DT)
T=T+DT
+ write(9,'(4f20.11)') t,x,y,z
counter = counter + 1
if(t.lt.tend) goto 10
write(*,"(i8)")counter
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp 2011-05-09 14:43:57 EDT (Mon, 09 May 2011)
@@ -159,9 +159,13 @@
data_type data1( k , data.x , data.derivs , data.dt_fac );
data_type data2( data.which - k , data.x , data.derivs , data.dt_fac );
- tmp += boost::math::binomial_coefficient< double >( data.which , k )
- * Grammar()( proto::left( expr ) , state , data1 )
- * Grammar()( proto::right( expr ) , state , data2 );
+// tmp += boost::math::binomial_coefficient< double >( data.which , k )
+// * Grammar()( proto::left( expr ) , state , data1 )
+// * Grammar()( proto::right( expr ) , state , data2 );
+
+ tmp += Grammar()( proto::left( expr ) , state , data1 )
+ * Grammar()( proto::right( expr ) , state , data2 );
+
}
return tmp;
@@ -388,6 +392,8 @@
typedef state_type deriv_type;
typedef std::tr1::array< state_type , order_value > derivs_type;
+ taylor( void ) : m_derivs() , m_dt_fac( 1.0 ) { }
+
order_type order( void ) const
{
return order_value;
@@ -424,38 +430,41 @@
BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
- double dt_fac = dt;
- eval_derivs( sys , in , m_derivs , dt_fac );
-
-// double max_error = 0.0;
-// for( size_t i=0 ; i<N ; ++i )
-// {
-// double error = std::abs( m_derivs[order_value-1][i] ) / ( std::abs( in[i] ) + 1.0e-35 );
-// max_error = std::max( error , max_error );
-// }
+ eval_derivs( sys , in , m_derivs , m_dt_fac );
-// dt = pow( dt0 / max_error , 1.0 / double( order_value ) );
+ double max_error = 0.0;
+ for( size_t i=0 ; i<dim ; ++i )
+ {
+ double error = std::abs( m_derivs[order_value-1][i] ) / ( std::abs( in[i] ) + 1.0e-35 );
+ max_error = std::max( error , max_error );
+ }
-// clog << dt << tab << dt0 << tab << max_error << tab << dt_fac << endl;
+ dt = pow( dt0 / max_error , 1.0 / double( order_value ) );
+// clog << dt << tab << max_error << tab << in[0] << tab << in[1] << tab << in[2] << tab;
+// clog << m_derivs[0][0] << tab << m_derivs[0][1] << tab << m_derivs[0][2] << tab << m_dt_fac << endl;
- for( size_t i=0 ; i<N ; ++i )
+ for( size_t i=0 ; i<dim ; ++i )
{
value_type tmp = 0.0;
for( size_t k=0 ; k<order_value ; ++k )
- tmp += 1.0 * ( tmp + m_derivs[order_value-1-k][i] );
+ tmp = dt * ( tmp + m_derivs[order_value-1-k][i] );
out[i] = in[i] + tmp;
}
+// clog << dt << tab << dt0 << tab << max_error << tab << m_dt_fac;
+// for( size_t j=0 ; j<dim ; ++j ) clog << tab << out[j];
+// clog << endl;
+
// clog << endl;
// for( size_t i=0 ; i<4 ; ++i )
// {
-// for( size_t j=0 ; j<N ; ++j )
+// for( size_t j=0 ; j<dim ; ++j )
// clog << m_derivs[i][j] << "\t";
// clog << endl;
// }
// clog << endl;
-// dt *= dt_fac;
+ dt *= m_dt_fac;
t += dt;
}
@@ -472,10 +481,15 @@
const double min_fac = 1.5;
const double max_fac = 0.6;
- for( size_t i=0 ; i<Order ; ++i )
+ for( size_t i=0 ; i<order_value ; ++i )
{
boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , der , dt_fac , i ) );
+// clog << i << tab << "Deriv : ";
+// for( size_t j=0 ; j<dim ; ++j )
+// clog << tab << der[i][j];
+// clog << endl;
+
/*
* OK # Fehler bestimmen
* OK # Fehler grenzen pruefen
@@ -483,25 +497,30 @@
* OK # dt_fac aendern
* OK # schon berechnete Ableitungen skalieren
*/
-// while( true )
-// {
-// double err = 0.0;
-// for( size_t j=0 ; j<N ; ++j ) err += std::abs( der[i][j] );
-//
-// if( err < min_error )
-// {
-// scale_derivs( der , i , min_fac );
-// dt_fac *= min_fac;
-// continue;
-// }
-// if( err > max_error )
-// {
-// scale_derivs( der , i , max_fac );
-// dt_fac *= max_fac;
-// continue;
-// }
-// break;
-// }
+ while( true )
+ {
+ double err = 0.0;
+ for( size_t j=0 ; j<dim ; ++j ) err += std::abs( der[i][j] );
+// clog << i;
+// for( size_t j=0 ; j<dim ; ++j ) clog << tab << der[i][j];
+// clog << tab << err << endl;
+
+ if( err < min_error )
+ {
+// clog << i << tab << "min_error : " << err << tab << dt_fac << endl;
+ scale_derivs( der , i , min_fac );
+ dt_fac *= min_fac;
+ continue;
+ }
+ if( err > max_error )
+ {
+// clog << i << tab << "max_error : " << err << tab << dt_fac << endl;
+ scale_derivs( der , i , max_fac );
+ dt_fac *= max_fac;
+ continue;
+ }
+ break;
+ }
}
@@ -513,7 +532,7 @@
for( size_t i=0 ; i<=order ; ++i )
{
scale *= fac;
- for( size_t j=0 ; j<N ; ++j )
+ for( size_t j=0 ; j<dim ; ++j )
der[i][j] *= scale;
}
}
@@ -523,6 +542,7 @@
private:
derivs_type m_derivs;
+ double m_dt_fac;
};
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp 2011-05-09 14:43:57 EDT (Mon, 09 May 2011)
@@ -6,6 +6,7 @@
*/
#include <iostream>
+#include <fstream>
#include "taylor.hpp"
@@ -19,7 +20,7 @@
return out;
}
-typedef boost::numeric::odeint::taylor< 3 , 10 > taylor_type;
+typedef boost::numeric::odeint::taylor< 3 , 25 > taylor_type;
typedef taylor_type::state_type state_type;
typedef taylor_type::derivs_type derivs_type;
@@ -50,22 +51,8 @@
state_type x = {{ 10.0 , 10.0 , 10.0 }} ;
double t = 0.0;
- double dt = 0.001;
- for( size_t i=0 ; i<10000 ; ++i )
- {
- stepper.try_step(
- fusion::make_vector
- (
- sigma * ( arg2 - arg1 ) ,
- R * arg1 - arg2 - arg1 * arg3 ,
- arg1 * arg2 - b * arg3
- ) ,
- x , t , dt );
- cout << i << "\t" << t << "\t" << x << "\n";
- }
-
-
-// while( t < 10.0 )
+ double dt = 1.0;
+// for( size_t i=0 ; i<10000 ; ++i )
// {
// stepper.try_step(
// fusion::make_vector
@@ -75,9 +62,26 @@
// arg1 * arg2 - b * arg3
// ) ,
// x , t , dt );
-//
-// cout << t << "\t" << x << endl;
+// cout << i << "\t" << t << "\t" << x << "\n";
// }
+ ofstream fout( "lorenz.dat" );
+ size_t count = 0;
+ while( t < 5000.0 )
+ {
+ stepper.try_step(
+ fusion::make_vector
+ (
+ sigma * ( arg2 - arg1 ) ,
+ R * arg1 - arg2 - arg1 * arg3 ,
+ arg1 * arg2 - b * arg3
+ ) ,
+ x , t , dt );
+
+ fout << t << "\t" << x << endl;
+ ++count;
+ }
+ clog << count << endl;
+
return 0;
}