|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r67959 - in sandbox/odeint/branches/karsten: boost/numeric/odeint/stepper libs/numeric/odeint/ideas/rosenbrock4
From: karsten.ahnert_at_[hidden]
Date: 2011-01-11 09:22:40
Author: karsten
Date: 2011-01-11 09:22:39 EST (Tue, 11 Jan 2011)
New Revision: 67959
URL: http://svn.boost.org/trac/boost/changeset/67959
Log:
rosenbrock continued
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4_stepper.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 18 +++++++---
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile | 4 ++
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp | 67 +++++++++++++++++++++++++++++++--------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp | 43 +++++++++++++++++-------
4 files changed, 98 insertions(+), 34 deletions(-)
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp 2011-01-11 09:22:39 EST (Tue, 11 Jan 2011)
@@ -77,7 +77,7 @@
)
: m_stepper( stepper ) , m_error_checker( error_checker ) ,
m_dxdt_size_adjuster() , m_xerr_size_adjuster() , m_xnew_size_adjuster() ,
- m_dxdt() , m_xerr() , m_xnew()
+ m_dxdt() , m_xerr() , m_xnew() , m_max_rel_error()
{
boost::numeric::odeint::construct( m_dxdt );
boost::numeric::odeint::construct( m_xerr );
@@ -142,21 +142,21 @@
// do one step with error calculation
m_stepper.do_step( sys , in , dxdt , t , out , dt , m_xerr );
- time_type max_rel_err = m_error_checker.error( in , dxdt , m_xerr , dt );
+ m_max_rel_error = m_error_checker.error( in , dxdt , m_xerr , dt );
- if( max_rel_err > 1.1 )
+ if( m_max_rel_error > 1.1 )
{
// error too large - decrease dt ,limit scaling factor to 0.2 and reset state
- dt *= max( 0.9 * pow( max_rel_err , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
+ dt *= max( 0.9 * pow( m_max_rel_error , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
return step_size_decreased;
}
else
{
- if( max_rel_err < 0.5 )
+ if( m_max_rel_error < 0.5 )
{
//error too small - increase dt and keep the evolution and limit scaling factor to 5.0
t += dt;
- dt *= min( 0.9 * pow( max_rel_err , -1.0 / m_stepper.stepper_order() ) , 5.0 );
+ dt *= min( 0.9 * pow( m_max_rel_error , -1.0 / m_stepper.stepper_order() ) , 5.0 );
return success_step_size_increased;
}
else
@@ -167,6 +167,11 @@
}
}
+ time_type last_error( void ) const
+ {
+ return m_max_rel_error;
+ }
+
@@ -204,6 +209,7 @@
state_type m_dxdt;
state_type m_xerr;
state_type m_xnew;
+ time_type m_max_rel_error;
};
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile 2011-01-11 09:22:39 EST (Tue, 11 Jan 2011)
@@ -18,4 +18,8 @@
exe rosenbrock4
: rosenbrock4.cpp
+ ;
+
+exe rosenbrock4_stepper
+ : rosenbrock4_stepper.cpp
;
\ No newline at end of file
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp 2011-01-11 09:22:39 EST (Tue, 11 Jan 2011)
@@ -22,25 +22,54 @@
typedef stepper_type::matrix_type matrix_type;
typedef rosenbrock4_controller< time_type > controlled_stepper_type;
+//template< class StateType >
+//void system( const StateType &x , StateType &dxdt , time_type t )
+//{
+// dxdt[0] = -0.013 * x[0] - 1000.0 * x[0] * x[2];
+// dxdt[1] = -2500.0 * x[1] * x[2];
+// dxdt[2] = -0.013 * x[0] - 1000.0 * x[0] * x[2] - 2500.0 * x[1] * x[2];
+//}
+//
+//void jacobi( const state_type &x , matrix_type &J , time_type t , state_type &dfdt )
+//{
+// J( 0 , 0 ) = -0.013 - 1000.0 * x[2];
+// J( 0 , 1 ) = 0.0;
+// J( 0 , 2 ) = -1000.0 * x[0];
+// J( 1 , 0 ) = 0.0;
+// J( 1 , 1 ) = -2500.0 * x[2];
+// J( 1 , 2 ) = -2500.0 * x[1];
+// J( 2 , 0 ) = -0.013 - 1000.0 * x[2];
+// J( 2 , 1 ) = -2500.0 * x[2];
+// J( 2 , 2 ) = -1000.0 * x[0] - 2500.0 * x[1];
+//
+// dfdt[0] = 0.0;
+// dfdt[1] = 0.0;
+// dfdt[2] = 0.0;
+//}
+
+const time_type sigma = 10.0;
+const time_type R = 28.0;
+const time_type b = 8.0 / 3.0;
+
template< class StateType >
void system( const StateType &x , StateType &dxdt , time_type t )
{
- dxdt[0] = -0.013 * x[0] - 1000.0 * x[0] * x[2];
- dxdt[1] = -2500.0 * x[1] * x[2];
- dxdt[2] = -0.013 * x[0] - 1000.0 * x[0] * x[2] - 2500.0 * x[1] * x[2];
+ 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];
}
void jacobi( const state_type &x , matrix_type &J , time_type t , state_type &dfdt )
{
- J( 0 , 0 ) = -0.013 - 1000.0 * x[2];
- J( 0 , 1 ) = 0.0;
- J( 0 , 2 ) = -1000.0 * x[0];
- J( 1 , 0 ) = 0.0;
- J( 1 , 1 ) = -2500.0 * x[2];
- J( 1 , 2 ) = -2500.0 * x[1];
- J( 2 , 0 ) = -0.013 - 1000.0 * x[2];
- J( 2 , 1 ) = -2500.0 * x[2];
- J( 2 , 2 ) = -1000.0 * x[0] - 2500.0 * x[1];
+ J( 0 , 0 ) = -sigma;
+ J( 0 , 1 ) = sigma;
+ J( 0 , 2 ) = 0.0;
+ J( 1 , 0 ) = R - x[2];
+ J( 1 , 1 ) = -1.0;
+ J( 1 , 2 ) = -x[0];
+ J( 2 , 0 ) = x[1];
+ J( 2 , 1 ) = x[0];
+ J( 2 , 2 ) = -b;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
@@ -48,6 +77,7 @@
}
+
int main( int argc , char **argv )
{
if( true )
@@ -65,8 +95,12 @@
size_t count = 0;
while( t < 50.0 )
{
- clog << t << "\t" << dt << "\n";
- fout << t << "\t" << dt << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << std::endl;
+// clog << t << "\t" << dt << "\t" << controlled_stepper.last_error() << "\n";
+ fout << t << "\t" << dt << "\t" << controlled_stepper.last_error() << "\t";
+ fout << x[0] << "\t" << x[1] << "\t" << x[2] << "\t";
+// fout << xerr[0] << "\t" << xerr[1] << "\t" << xerr[2] << "\t";
+ fout <<std::endl;
+
size_t trials = 0;
while( trials < 100 )
{
@@ -104,7 +138,10 @@
size_t count = 0;
while( t < 50.0 )
{
- fout << t << "\t" << dt << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << std::endl;
+ fout << t << "\t" << dt << "\t" << stepper.last_error() << "\t";
+ fout << x[0] << "\t" << x[1] << "\t" << x[2] << "\t";
+ fout <<std::endl;
+
size_t trials = 0;
while( trials < 100 )
{
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp 2011-01-11 09:22:39 EST (Tue, 11 Jan 2011)
@@ -20,6 +20,15 @@
namespace odeint {
+/*
+ * ToDo:
+ *
+ * 1. roll out parameters
+ * 2. how to call, with jacobi and system?
+ * 3. Introduce adjust size
+ * 4. Interfacing for odeint, check if controlled_error_stepper can be used
+ * 5. dense output
+ */
template< class Value >
class rosenbrock4
@@ -35,23 +44,26 @@
{
}
+// void do_step( System &system , Jacobi &jacobi , ... );
+// void do_step( pair< System& , Jacobi& > system , ... );
+// void do_step( pair< System , Jacobi > & , ... );
+
template< class System , class Jacobi >
void do_step( System &system , Jacobi &jacobi , const state_type &x , time_type t , state_type &xout , time_type dt , state_type &xerr )
{
- const double gamma = 1.0;
- const double d1 = 1.0 , d2 = 1.0 , d3 = 1.0 , d4 = 1.0;
- const double c2 = 1.0 , c3 = 1.0 , c4 = 1.0;
- const double c21 = 1.0;
- const double a21 = 1.0;
- const double c31 = 1.0 , c32 = 1.0;
- const double a31 = 1.0 , a32 = 1.0;
- const double c41 = 1.0 , c42 = 1.0 , c43 = 1.0;
- const double a41 = 1.0 , a42 = 1.0 , a43 = 1.0;
- const double c51 = 1.0 , c52 = 1.0 , c53 = 1.0 , c54 = 1.0 ;
- const double a51 = 1.0 , a52 = 1.0 , a53 = 1.0 , a54 = 1.0 ;
- const double c61 = 1.0 , c62 = 1.0 , c63 = 1.0 , c64 = 1.0 , c65 = 1.0;
-
+ const double gamma = 0.25;
+ const double d1 = 0.25 , d2 = -0.1043 , d3 = 0.1035 , d4 = 0.3620000000000023e-01;
+ const double c2 = 0.386 , c3 = 0.21 , c4 = 0.63;
+ const double c21 = -0.5668800000000000e+01;
+ const double a21 = 0.1544000000000000e+01;
+ const double c31 = -0.2430093356833875e+01 , c32 = -0.2063599157091915e+00;
+ const double a31 = 0.9466785280815826e+00 , a32 = 0.2557011698983284e+00;
+ const double c41 = -0.1073529058151375e+00 , c42 = -0.9594562251023355e+01 , c43 = -0.2047028614809616e+02;
+ const double a41 = 0.3314825187068521e+01 , a42 = 0.2896124015972201e+01 , a43 = 0.9986419139977817e+00;
+ const double c51 = 0.7496443313967647e+01 , c52 = -0.1024680431464352e+02 , c53 = -0.3399990352819905e+02 , c54 = 0.1170890893206160e+02;
+ const double a51 = 0.1221224509226641e+01 , a52 = 0.6019134481288629e+01 , a53 = 0.1253708332932087e+02 , a54 = -0.6878860361058950e+00 ;
+ const double c61 = 0.8083246795921522e+01 , c62 = -0.7981132988064893e+01 , c63 = -0.3152159432874371e+02 , c64 = 0.1631930543123136e+02 , c65 = -0.6058818238834054e+01;
const size_t n = x.size();
matrix_type jac( n , n );
@@ -159,6 +171,11 @@
return sqrt( err / time_type( n ) );
}
+ time_type last_error( void ) const
+ {
+ return m_err_old;
+ }
+
template< class System , class Jacobi >
boost::numeric::odeint::controlled_step_result
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4_stepper.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4_stepper.cpp 2011-01-11 09:22:39 EST (Tue, 11 Jan 2011)
@@ -0,0 +1,140 @@
+/*
+ * rosenbrock4_stepper.cpp
+ *
+ * Created on: Jan 9, 2011
+ * Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+#include <tr1/array>
+
+#include "rosenbrock4.hpp"
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+const static size_t dim = 3;
+typedef double time_type;
+typedef rosenbrock4< time_type > stepper_type;
+typedef stepper_type::state_type state_type;
+typedef stepper_type::matrix_type matrix_type;
+typedef rosenbrock4_controller< time_type > controlled_stepper_type;
+
+//template< class StateType >
+//void system( const StateType &x , StateType &dxdt , time_type t )
+//{
+// dxdt[0] = -0.013 * x[0] - 1000.0 * x[0] * x[2];
+// dxdt[1] = -2500.0 * x[1] * x[2];
+// dxdt[2] = -0.013 * x[0] - 1000.0 * x[0] * x[2] - 2500.0 * x[1] * x[2];
+//}
+//
+//void jacobi( const state_type &x , matrix_type &J , time_type t , state_type &dfdt )
+//{
+// J( 0 , 0 ) = -0.013 - 1000.0 * x[2];
+// J( 0 , 1 ) = 0.0;
+// J( 0 , 2 ) = -1000.0 * x[0];
+// J( 1 , 0 ) = 0.0;
+// J( 1 , 1 ) = -2500.0 * x[2];
+// J( 1 , 2 ) = -2500.0 * x[1];
+// J( 2 , 0 ) = -0.013 - 1000.0 * x[2];
+// J( 2 , 1 ) = -2500.0 * x[2];
+// J( 2 , 2 ) = -1000.0 * x[0] - 2500.0 * x[1];
+//
+// dfdt[0] = 0.0;
+// dfdt[1] = 0.0;
+// dfdt[2] = 0.0;
+//}
+
+
+const time_type sigma = 10.0;
+const time_type R = 28.0;
+const time_type b = 8.0 / 3.0;
+
+template< class StateType >
+void system( const StateType &x , StateType &dxdt , time_type 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];
+}
+
+void jacobi( const state_type &x , matrix_type &J , time_type t , state_type &dfdt )
+{
+ J( 0 , 0 ) = -sigma;
+ J( 0 , 1 ) = sigma;
+ J( 0 , 2 ) = 0.0;
+ J( 1 , 0 ) = R - x[2];
+ J( 1 , 1 ) = -1.0;
+ J( 1 , 2 ) = -x[0];
+ J( 2 , 0 ) = x[1];
+ J( 2 , 1 ) = x[0];
+ J( 2 , 2 ) = -b;
+
+ dfdt[0] = 0.0;
+ dfdt[1] = 0.0;
+ dfdt[2] = 0.0;
+}
+
+
+
+int main( int argc , char **argv )
+{
+ const double dt = 0.01;
+ size_t steps = 1000;
+ double x0 = -12.0 , y0 = -12.0 , z0 = 20.0;
+ if( true )
+ {
+ state_type x( dim ) , xerr( dim );
+ time_type t = 0.0;
+
+ stepper_type stepper;
+ x[0] = x0 ; x[1] = y0 ; x[2] = z0;
+
+ ofstream fout( "dat/ross.dat" );
+ fout.precision( 14 );
+ size_t count = 0;
+ while( count < steps )
+ {
+ fout << t << "\t";
+ fout << x[0] << "\t" << x[1] << "\t" << x[2] << "\t";
+ fout << xerr[0] << "\t" << xerr[1] << "\t" << xerr[2] << "\t";
+ fout <<std::endl;
+
+ stepper.do_step( system< state_type > , jacobi , x , t , dt , xerr );
+ ++count;
+ t += dt;
+ }
+ }
+
+
+
+
+
+ if( true )
+ {
+ typedef std::tr1::array< time_type , 3 > state_type2;
+ typedef explicit_error_rk54_ck< state_type2 > stepper_type2;
+ stepper_type2 rk_stepper;
+ state_type2 x = {{ x0 , y0 , z0 }} , xerr = {{ 0.0 , 0.0 , 0.0 }};
+ time_type t = 0.0;
+
+ ofstream fout( "dat/rk.dat" );
+ fout.precision( 14 );
+ size_t count = 0;
+ while( count < steps )
+ {
+ fout << t << "\t";
+ fout << x[0] << "\t" << x[1] << "\t" << x[2] << "\t";
+ fout << xerr[0] << "\t" << xerr[1] << "\t" << xerr[2] << "\t";
+ fout <<std::endl;
+
+ rk_stepper.do_step( system< state_type2 > , x , t , dt , xerr );
+ ++count;
+ t += dt;
+ }
+ }
+
+ 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