|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r71826 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: lorenz_special taylor_v4
From: karsten.ahnert_at_[hidden]
Date: 2011-05-08 11:56:12
Author: karsten
Date: 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
New Revision: 71826
URL: http://svn.boost.org/trac/boost/changeset/71826
Log:
continued taylor
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_original.f (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/Jamfile
- copied, changed from r71691, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/Jamfile
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp
- copied, changed from r71691, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp
- copied, changed from r71691, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_main.cpp
Removed:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz.f
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/Jamfile | 2 -
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp | 80 +++++++++++++++------------------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp | 46 ++++++++---------------
3 files changed, 47 insertions(+), 81 deletions(-)
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz.f
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz.f 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
+++ (empty file)
@@ -1,140 +0,0 @@
- IMPLICIT REAL*8 (A-H,O-Z)
- call sections
- stop
- end
-
- Subroutine SECTIONS
- IMPLICIT REAL*8 (A-H,O-Z)
- real*8 AX(77777),AY(77777),AZ(77777),AT(77777)
- DIMENSION DX(0:40),DY(0:40),DZ(0:40)
-C This subroutine provides the values of consecutive
-c intersection coordinates of the solution of the Lorenz equations
-c with the plane Z=const
-c (keeping them in the arrays AX and AY). Array DT contains the
-c time values at the moments of intersections.
-c These data are written into the file 'lorsec.dat'
-
- COMMON /ST/ DX,DY,DZ,Q,NO
- open(9,file='lorsec.dat',status='unknown',access='sequential',
- : form='formatted')
-C Introduce the parameters:
- P=10.d0
- B=8.d0/3.d0
- R=28.d0
-C Insert the initial point:
- write(*,"(' Initial X ? '$)")
- read*,x
- write(*,"(' Initial Y ? '$)")
- read*,y
- write(*,"(' Initial Z ? '$)")
- read*,z
- t=0
-C Insert the desired number of the Poincare sections
- write(*,"(' How many section points to find ? '$)")
- read*,JSECT
- isect=0
-C Insert the coordinate Z of the section plane
- write(*,"(' Coordinate Z of the secant plane = ? '$)")
- read*,zsect
- 10 xp=x
- yp=y
- zp=z
- call step (P,R,B,X,Y,Z,DT)
- jjj=jjj+1
- T=T+DT
- IF(z.lt.zsect.and.zp.ge.zsect)then
-C computation of the section coordinates
- td=0
- tu=DT/q
- 12 tt=td+0.5d0*(tu-td)
- QZ=0.
- DO I=0,NO-1
- QZ=(QZ+DZ(NO-I))*tt
- ENDDO
- if(zp+QZ.gt.zsect)td=tt
- if(zp+QZ.le.zsect)tU=tt
- if(tu-td.Gt.1.d-10*DT/Q)goto 12
- QX=0.
- QY=0.
- DO I=0,NO-1
- QX=(QX+DX(NO-I))*tt
- QY=(QY+DY(NO-I))*tt
- ENDDO
- isect=isect+1
- ax(isect)=xp+qx
- ay(isect)=yp+qy
- at(isect)=T-DT+tt*Q
- write(9,'(i10,3f20.11)') isect,ax(isect),ay(isect),at(isect)
- if(jsect/20*(20*isect/jsect).eq.isect)
- : write(*,"(i3,'%'$)")100*isect/jsect
- ENDIF
- if(isect.lt.jsect)goto 10
- print*,'ready'
- close(9)
- return
- end
-
- subroutine step (P,R,B,X,Y,Z,DT)
-c This subroutine performs one integration step for the Lorenz
-c equations by integrating them with the method of order NO,
-c and relative error TO. Here P,R and B are the usual parameters
-c of the Lorenz equations. The input values of X,Y,Z are the current
-c values of the respective variables; as the output values they
-c contain the new values. The neccesary stepsize is computed
-c automatically and can be recovered from the output value DT.
- IMPLICIT REAL*8 (A-H,O-Z)
- DIMENSION DX(0:40),DY(0:40),DZ(0:40)
- COMMON /ST/ DX,DY,DZ,Q,NO
- DATA Q/1.d0/
- TO=1.D-17
- NO=25
- DX(0)=X
- DY(0)=Y
- DZ(0)=Z
- Q1=1.
- DO 19 ND=1,NO
- DXY=0
- DXZ=0
- ND1=ND-1
- QQ=Q/ND
- DO I=0,ND1
- DXY=DXY+DX(I)*DY(ND1-I)
- DXZ=DXZ+DX(I)*DZ(ND1-I)
- enddo
- DX(ND)=P*(DY(ND1)-DX(ND1))*QQ
- DY(ND)=(R*DX(ND1)-DY(ND1)-DXZ)*QQ
- DZ(ND)=(DXY-B*DZ(ND1))*QQ
- 13 Q1=1.
- Q2=DABS(DX(ND))+DABS(DY(ND))+DABS(DZ(ND))
- IF(Q2.LT.1.D-19)Q1=1.5
- IF(Q2.GT.1.D19)Q1=0.6
- IF(Q1.LT.0.95.OR.Q1.GT.1.05)THEN
- Q2=Q1
- DO M=1,ND
- DX(M)=DX(M)*Q2
- DY(M)=DY(M)*Q2
- DZ(M)=DZ(M)*Q2
- Q2=Q2*Q1
- enddo
- Q=Q*Q1
- GOTO 13
- ENDIF
- 19 CONTINUE
- EX=DABS(DX(NO))/(DABS(X)+1.D-35)
- EY=DABS(DY(NO))/(DABS(Y)+1.D-35)
- EZ=DABS(DZ(NO))/(DABS(Z)+1.D-35)
- DT=(TO/DMAX1(EX,EY,EZ))**(1.d0/NO)
- QZ=0.
- QY=0.
- QX=0.
- DO I=0,NO-1
- QX=(QX+DX(NO-I))*DT
- QY=(QY+DY(NO-I))*DT
- QZ=(QZ+DZ(NO-I))*DT
- ENDDO
- X=X+QX
- Y=Y+QY
- Z=Z+QZ
- DT=DT*Q
- return
- END
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_original.f
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_original.f 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
@@ -0,0 +1,140 @@
+ IMPLICIT REAL*8 (A-H,O-Z)
+ call sections
+ stop
+ end
+
+ Subroutine SECTIONS
+ IMPLICIT REAL*8 (A-H,O-Z)
+ real*8 AX(77777),AY(77777),AZ(77777),AT(77777)
+ DIMENSION DX(0:40),DY(0:40),DZ(0:40)
+C This subroutine provides the values of consecutive
+c intersection coordinates of the solution of the Lorenz equations
+c with the plane Z=const
+c (keeping them in the arrays AX and AY). Array DT contains the
+c time values at the moments of intersections.
+c These data are written into the file 'lorsec.dat'
+
+ COMMON /ST/ DX,DY,DZ,Q,NO
+ open(9,file='lorsec.dat',status='unknown',access='sequential',
+ : form='formatted')
+C Introduce the parameters:
+ P=10.d0
+ B=8.d0/3.d0
+ R=28.d0
+C Insert the initial point:
+ write(*,"(' Initial X ? '$)")
+ read*,x
+ write(*,"(' Initial Y ? '$)")
+ read*,y
+ write(*,"(' Initial Z ? '$)")
+ read*,z
+ t=0
+C Insert the desired number of the Poincare sections
+ write(*,"(' How many section points to find ? '$)")
+ read*,JSECT
+ isect=0
+C Insert the coordinate Z of the section plane
+ write(*,"(' Coordinate Z of the secant plane = ? '$)")
+ read*,zsect
+ 10 xp=x
+ yp=y
+ zp=z
+ call step (P,R,B,X,Y,Z,DT)
+ jjj=jjj+1
+ T=T+DT
+ IF(z.lt.zsect.and.zp.ge.zsect)then
+C computation of the section coordinates
+ td=0
+ tu=DT/q
+ 12 tt=td+0.5d0*(tu-td)
+ QZ=0.
+ DO I=0,NO-1
+ QZ=(QZ+DZ(NO-I))*tt
+ ENDDO
+ if(zp+QZ.gt.zsect)td=tt
+ if(zp+QZ.le.zsect)tU=tt
+ if(tu-td.Gt.1.d-10*DT/Q)goto 12
+ QX=0.
+ QY=0.
+ DO I=0,NO-1
+ QX=(QX+DX(NO-I))*tt
+ QY=(QY+DY(NO-I))*tt
+ ENDDO
+ isect=isect+1
+ ax(isect)=xp+qx
+ ay(isect)=yp+qy
+ at(isect)=T-DT+tt*Q
+ write(9,'(i10,3f20.11)') isect,ax(isect),ay(isect),at(isect)
+ if(jsect/20*(20*isect/jsect).eq.isect)
+ : write(*,"(i3,'%'$)")100*isect/jsect
+ ENDIF
+ if(isect.lt.jsect)goto 10
+ print*,'ready'
+ close(9)
+ return
+ end
+
+ subroutine step (P,R,B,X,Y,Z,DT)
+c This subroutine performs one integration step for the Lorenz
+c equations by integrating them with the method of order NO,
+c and relative error TO. Here P,R and B are the usual parameters
+c of the Lorenz equations. The input values of X,Y,Z are the current
+c values of the respective variables; as the output values they
+c contain the new values. The neccesary stepsize is computed
+c automatically and can be recovered from the output value DT.
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION DX(0:40),DY(0:40),DZ(0:40)
+ COMMON /ST/ DX,DY,DZ,Q,NO
+ DATA Q/1.d0/
+ TO=1.D-17
+ NO=25
+ DX(0)=X
+ DY(0)=Y
+ DZ(0)=Z
+ Q1=1.
+ DO 19 ND=1,NO
+ DXY=0
+ DXZ=0
+ ND1=ND-1
+ QQ=Q/ND
+ DO I=0,ND1
+ DXY=DXY+DX(I)*DY(ND1-I)
+ DXZ=DXZ+DX(I)*DZ(ND1-I)
+ enddo
+ DX(ND)=P*(DY(ND1)-DX(ND1))*QQ
+ DY(ND)=(R*DX(ND1)-DY(ND1)-DXZ)*QQ
+ DZ(ND)=(DXY-B*DZ(ND1))*QQ
+ 13 Q1=1.
+ Q2=DABS(DX(ND))+DABS(DY(ND))+DABS(DZ(ND))
+ IF(Q2.LT.1.D-19)Q1=1.5
+ IF(Q2.GT.1.D19)Q1=0.6
+ IF(Q1.LT.0.95.OR.Q1.GT.1.05)THEN
+ Q2=Q1
+ DO M=1,ND
+ DX(M)=DX(M)*Q2
+ DY(M)=DY(M)*Q2
+ DZ(M)=DZ(M)*Q2
+ Q2=Q2*Q1
+ enddo
+ Q=Q*Q1
+ GOTO 13
+ ENDIF
+ 19 CONTINUE
+ EX=DABS(DX(NO))/(DABS(X)+1.D-35)
+ EY=DABS(DY(NO))/(DABS(Y)+1.D-35)
+ EZ=DABS(DZ(NO))/(DABS(Z)+1.D-35)
+ DT=(TO/DMAX1(EX,EY,EZ))**(1.d0/NO)
+ QZ=0.
+ QY=0.
+ QX=0.
+ DO I=0,NO-1
+ QX=(QX+DX(NO-I))*DT
+ QY=(QY+DY(NO-I))*DT
+ QZ=(QZ+DZ(NO-I))*DT
+ ENDDO
+ X=X+QX
+ Y=Y+QY
+ Z=Z+QZ
+ DT=DT*Q
+ return
+ END
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/lorenz_special/lorenz_special.f 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
@@ -0,0 +1,104 @@
+ IMPLICIT REAL*8 (A-H,O-Z)
+ call integration
+ stop
+ end
+
+ Subroutine INTEGRATION
+ IMPLICIT REAL*8 (A-H,O-Z)
+ real*8 AX(77777),AY(77777),AZ(77777),AT(77777)
+ DIMENSION DX(0:40),DY(0:40),DZ(0:40)
+C This subroutine provides the values of consecutive
+c intersection coordinates of the solution of the Lorenz equations
+c with the plane Z=const
+c (keeping them in the arrays AX and AY). Array DT contains the
+c time values at the moments of intersections.
+c These data are written into the file 'lorsec.dat'
+
+ COMMON /ST/ DX,DY,DZ,Q,NO
+ INTEGER counter
+C Introduce the parameters:
+ P=10.d0
+ B=8.d0/3.d0
+ R=28.d0
+C Insert the initial point:
+ x=10.d0
+ y=10.d0
+ z=10.d0
+ t=0.d0
+ tend=5000.d0
+ counter=0
+ 10 xp=x
+ yp=y
+ zp=z
+ call step (P,R,B,X,Y,Z,DT)
+ T=T+DT
+ counter = counter + 1
+ if(t.lt.tend) goto 10
+ write(*,"(i8)")counter
+ return
+ end
+
+ subroutine step (P,R,B,X,Y,Z,DT)
+c This subroutine performs one integration step for the Lorenz
+c equations by integrating them with the method of order NO,
+c and relative error TO. Here P,R and B are the usual parameters
+c of the Lorenz equations. The input values of X,Y,Z are the current
+c values of the respective variables; as the output values they
+c contain the new values. The neccesary stepsize is computed
+c automatically and can be recovered from the output value DT.
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION DX(0:40),DY(0:40),DZ(0:40)
+ COMMON /ST/ DX,DY,DZ,Q,NO
+ DATA Q/1.d0/
+ TO=1.D-17
+ NO=25
+ DX(0)=X
+ DY(0)=Y
+ DZ(0)=Z
+ Q1=1.
+ DO 19 ND=1,NO
+ DXY=0
+ DXZ=0
+ ND1=ND-1
+ QQ=Q/ND
+ DO I=0,ND1
+ DXY=DXY+DX(I)*DY(ND1-I)
+ DXZ=DXZ+DX(I)*DZ(ND1-I)
+ enddo
+ DX(ND)=P*(DY(ND1)-DX(ND1))*QQ
+ DY(ND)=(R*DX(ND1)-DY(ND1)-DXZ)*QQ
+ DZ(ND)=(DXY-B*DZ(ND1))*QQ
+ 13 Q1=1.
+ Q2=DABS(DX(ND))+DABS(DY(ND))+DABS(DZ(ND))
+ IF(Q2.LT.1.D-19)Q1=1.5
+ IF(Q2.GT.1.D19)Q1=0.6
+ IF(Q1.LT.0.95.OR.Q1.GT.1.05)THEN
+ Q2=Q1
+ DO M=1,ND
+ DX(M)=DX(M)*Q2
+ DY(M)=DY(M)*Q2
+ DZ(M)=DZ(M)*Q2
+ Q2=Q2*Q1
+ enddo
+ Q=Q*Q1
+ GOTO 13
+ ENDIF
+ 19 CONTINUE
+ EX=DABS(DX(NO))/(DABS(X)+1.D-35)
+ EY=DABS(DY(NO))/(DABS(Y)+1.D-35)
+ EZ=DABS(DZ(NO))/(DABS(Z)+1.D-35)
+ DT=(TO/DMAX1(EX,EY,EZ))**(1.d0/NO)
+ QZ=0.
+ QY=0.
+ QX=0.
+ DO I=0,NO-1
+ QX=(QX+DX(NO-I))*DT
+ QY=(QY+DY(NO-I))*DT
+ QZ=(QZ+DZ(NO-I))*DT
+ ENDDO
+ X=X+QX
+ Y=Y+QY
+ Z=Z+QZ
+ DT=DT*Q
+ return
+ END
Copied: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/Jamfile (from r71691, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/Jamfile)
==============================================================================
--- /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/Jamfile 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
@@ -11,5 +11,3 @@
;
exe taylor_main : taylor_main.cpp ;
-exe taylor_lorenz_evol_direct : taylor_lorenz_evol_direct.cpp ;
-exe taylor_lorenz_evol_adaptive_statistics : taylor_lorenz_evol_adaptive_statistics.cpp ;
Copied: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp (from r71691, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor.hpp 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
@@ -31,7 +31,15 @@
// proto includes
#include <boost/proto/proto.hpp>
-
+/*
+ * OK # dt entfernen aus context
+ * OK # do_step zu try_step
+ * # Skalierung
+ * # Fehler als Member parameter
+ * # Fehler anhand der hoechsten Ableitung berechnen
+ * # Schrittweite bestimmen
+ * # Schritt ausfuehren
+ */
namespace boost {
@@ -70,12 +78,11 @@
typedef std::tr1::array< state_type , MaxOrder > deriv_type;
size_t which;
- double dt;
const state_type &x;
deriv_type &derivs;
- taylor_context( size_t _which , double _dt , const state_type &_x , deriv_type &_derivs )
- : which( _which ) , dt( _dt ) , x( _x ) , derivs( _derivs ) { }
+ taylor_context( size_t _which , const state_type &_x , deriv_type &_derivs )
+ : which( _which ) , x( _x ) , derivs( _derivs ) { }
};
@@ -141,8 +148,8 @@
double tmp = 0.0;
for( size_t k=0 ; k<=data.which ; ++k )
{
- data_type data1( k , data.dt , data.x , data.derivs );
- data_type data2( data.which - k , data.dt , data.x , data.derivs );
+ data_type data1( k , data.x , data.derivs );
+ data_type data2( data.which - k , data.x , data.derivs );
tmp += boost::math::binomial_coefficient< double >( data.which , k )
* Grammar()( proto::left( expr ) , state , data1 )
@@ -337,8 +344,8 @@
System m_sys;
taylor_context_type m_data;
- eval_derivs( System sys , const State &x , deriv_type &derivs , size_t which , double dt )
- : m_sys( sys ) , m_data( which , dt , x , derivs ) { }
+ eval_derivs( System sys , const State &x , deriv_type &derivs , size_t which )
+ : m_sys( sys ) , m_data( which , x , derivs ) { }
template< class Index >
void operator()( Index )
@@ -347,18 +354,15 @@
const expr_type &expr = boost::fusion::at< Index >( m_sys );
double deriv = taylor_transform()( expr , 0.0 , m_data );
- m_data.derivs[ m_data.which ][ Index::value ] = m_data.dt / double( m_data.which + 1 ) * deriv;
+ m_data.derivs[ m_data.which ][ Index::value ] = 1.0 / double( m_data.which + 1 ) * deriv;
}
};
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 )
+ 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 , dt );
+ return eval_derivs< System , State , Order >( sys , x , derivs , i );
}
-
-
- struct optimize_tree : proto::or_< proto::_ > { };
}
template< size_t N , size_t Order >
@@ -395,71 +399,49 @@
template< class System >
- void do_step( System sys , state_type &x , time_type t , time_type dt )
+ void try_step( System sys , state_type &x , time_type &t , time_type &dt )
{
BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
- do_step( sys , x , t , x , dt );
+ try_step( sys , x , t , x , dt );
}
template< class System >
- void do_step( System sys , const state_type &in , time_type t , state_type &out , time_type dt )
+ void try_step( System sys , const state_type &in , time_type &t , state_type &out , time_type &dt )
{
BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
- typedef typename boost::fusion::result_of::transform< System , taylor_adf::optimize_tree >::type optimized_system_type;
-
- optimized_system_type optimized_system = boost::fusion::transform( sys , taylor_adf::optimize_tree() );
-
for( size_t i=0 ; i<Order ; ++i )
{
- boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( optimized_system , in , m_derivs , i , dt ) );
+ boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , m_derivs , i ) );
}
for( size_t i=0 ; i<N ; ++i )
{
- out[i] = in[i];
+ value_type tmp = 0.0;
for( size_t k=0 ; k<order_value ; ++k )
- out[i] += m_derivs[k][i];
+ tmp += dt * ( tmp + m_derivs[k][order_value-1-i] );
+ out[i] = in[i] + tmp;
}
+
+ t += dt;
}
- template< class System >
- void do_step( System sys , state_type &x , time_type t , time_type dt , state_type &xerr )
+ const derivs_type& get_last_derivs( void ) const
{
- BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
- BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
- BOOST_STATIC_ASSERT(( order_value > 1 ));
-
- do_step( sys , x , t , x , dt , xerr );
+ return m_derivs;
}
template< class System >
- void do_step( System sys , const state_type &in , time_type t , state_type &out , time_type dt , state_type &xerr )
+ void eval_derivs( System sys , const state_type &in , derivs_type &der ) const
{
- BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
- BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
- BOOST_STATIC_ASSERT(( order_value > 1 ));
-
for( size_t i=0 ; i<Order ; ++i )
{
- boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , m_derivs , i , dt ) );
- }
-
- for( size_t i=0 ; i<N ; ++i )
- {
- out[i] = in[i];
- for( size_t k=0 ; k<order_value ; ++k )
- out[i] += m_derivs[k][i];
- xerr[i] = m_derivs[order_value-1][i];
+ boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , der , i ) );
}
- }
- const derivs_type& get_last_derivs( void ) const
- {
- return m_derivs;
}
Copied: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp (from r71691, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_main.cpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_main.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v4/taylor_main.cpp 2011-05-08 11:56:10 EDT (Sun, 08 May 2011)
@@ -27,12 +27,6 @@
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];
-}
@@ -51,33 +45,25 @@
{
cout.precision( 14 );
- taylor_type t;
-
- state_type in = {{ 10.0 , 10.0 , 10.0 }} , dxdt = {{ 0.0 , 0.0 , 0.0 }} , xerr = {{ 0.0 , 0.0 , 0.0 }} , out = {{ 0.0 ,0.0 , 0.0 }};
-
- lorenz( in , dxdt , 0.0 );
-
- cout << in << endl;
- cout << dxdt << endl << endl;
-
- t.do_step(
- fusion::make_vector
- (
- sigma * ( arg2 - arg1 ) ,
- R * arg1 - arg2 - arg1 * arg3 ,
- arg1 * arg2 - b * arg3
- ) ,
- in , 0.0 , out , 0.1 , xerr );
+ taylor_type stepper;
+ state_type x = {{ 10.0 , 10.0 , 10.0 }} ;
+ double t = 0.0;
+ double dt = 0.01;
+ while( t < 10.0 )
+ {
+ stepper.try_step(
+ fusion::make_vector
+ (
+ sigma * ( arg2 - arg1 ) ,
+ R * arg1 - arg2 - arg1 * arg3 ,
+ arg1 * arg2 - b * arg3
+ ) ,
+ x , t , dt );
- cout << endl;
- cout << in << endl;
- cout << xerr << endl;
- cout << out << endl << endl;
- const derivs_type &derivs = t.get_last_derivs();
- for( size_t i=0 ; i<derivs.size() ; ++i )
- cout << derivs[i] << endl;
+ cout << t << "\t" << x << 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