Subject: [Boost-commit] svn:boost r67869 - in sandbox/odeint/branches/karsten: . boost/numeric/odeint/stepper libs/numeric/odeint/ideas/rosenbrock4
From: karsten.ahnert_at_[hidden]
Date: 2011-01-09 12:33:19
Author: karsten
Date: 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
New Revision: 67869
URL: http://svn.boost.org/trac/boost/changeset/67869
Log:
* small template parameter change in controlled_error_stepper
* introducing ideas/rosenbrock4
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/.cproject | 2 +-
sandbox/odeint/branches/karsten/.project | 8 ++++----
sandbox/odeint/branches/karsten/Jamroot | 20 ++++++++++++++------
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/controlled_error_stepper.hpp | 21 ++++++++++++---------
4 files changed, 31 insertions(+), 20 deletions(-)
Modified: sandbox/odeint/branches/karsten/.cproject
==============================================================================
--- sandbox/odeint/branches/karsten/.cproject (original)
+++ sandbox/odeint/branches/karsten/.cproject 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -20,7 +20,7 @@
<folderInfo id="0.1427786045." name="/" resourcePath="">
<toolChain id="org.eclipse.cdt.build.core.prefbase.toolchain.1701227309" name="No ToolChain" resourceTypeBasedDiscovery="false" superClass="org.eclipse.cdt.build.core.prefbase.toolchain">
<targetPlatform id="org.eclipse.cdt.build.core.prefbase.toolchain.1701227309.1923743366" name=""/>
- <builder id="org.eclipse.cdt.build.core.settings.default.builder.1866451942" keepEnvironmentInBuildfile="false" managedBuildOn="false" name="Gnu Make Builder" superClass="org.eclipse.cdt.build.core.settings.default.builder"/>
+ <builder cleanBuildTarget="--clean" command="bjam" id="org.eclipse.cdt.build.core.settings.default.builder.1866451942" incrementalBuildTarget="" keepEnvironmentInBuildfile="false" managedBuildOn="false" name="Gnu Make Builder" superClass="org.eclipse.cdt.build.core.settings.default.builder"/>
<tool id="org.eclipse.cdt.build.core.settings.holder.libs.295828857" name="holder for library settings" superClass="org.eclipse.cdt.build.core.settings.holder.libs"/>
<tool id="org.eclipse.cdt.build.core.settings.holder.1728088817" name="Assembly" superClass="org.eclipse.cdt.build.core.settings.holder">
<inputType id="org.eclipse.cdt.build.core.settings.holder.inType.1379843062" languageId="org.eclipse.cdt.core.assembly" languageName="Assembly" sourceContentType="org.eclipse.cdt.core.asmSource" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
Modified: sandbox/odeint/branches/karsten/.project
==============================================================================
--- sandbox/odeint/branches/karsten/.project (original)
+++ sandbox/odeint/branches/karsten/.project 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -29,11 +29,11 @@
</dictionary>
<dictionary>
<key>org.eclipse.cdt.make.core.buildCommand</key>
- <value>make</value>
+ <value>bjam</value>
</dictionary>
<dictionary>
<key>org.eclipse.cdt.make.core.cleanBuildTarget</key>
- <value>clean</value>
+ <value>--clean</value>
</dictionary>
<dictionary>
<key>org.eclipse.cdt.make.core.contents</key>
@@ -53,7 +53,7 @@
</dictionary>
<dictionary>
<key>org.eclipse.cdt.make.core.fullBuildTarget</key>
- <value>all</value>
+ <value></value>
</dictionary>
<dictionary>
<key>org.eclipse.cdt.make.core.stopOnError</key>
@@ -61,7 +61,7 @@
</dictionary>
<dictionary>
<key>org.eclipse.cdt.make.core.useDefaultBuildCmd</key>
- <value>true</value>
+ <value>false</value>
</dictionary>
</arguments>
</buildCommand>
Modified: sandbox/odeint/branches/karsten/Jamroot
==============================================================================
--- sandbox/odeint/branches/karsten/Jamroot (original)
+++ sandbox/odeint/branches/karsten/Jamroot 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -5,16 +5,24 @@
import os ;
import modules ;
-import path ;
+import path ;
+
+path-constant BOOST_ROOT : [ os.environ BOOST_ROOT ] ;
+
+project
+ : requirements
+ <include>$(BOOST_ROOT) ;
-project odeint
- : requirements ;
-# <include>$BOOST_ROOT ;
-build-project libs/numeric/odeint/examples ;
build-project libs/numeric/odeint/test ;
+build-project libs/numeric/odeint/examples ;
+
+
+
+# ideas
build-project libs/numeric/odeint/ideas/butcher ;
build-project libs/numeric/odeint/ideas/generic_stepper ;
+build-project libs/numeric/odeint/ideas/rosenbrock4 ;
# additional tests with external libraries :
build-project libs/numeric/odeint/test/gmp ;
@@ -26,7 +34,7 @@
-path-constant BOOST_ROOT : [ os.environ BOOST_ROOT ] ;
+
###### The following is copied from another sandbox project #####
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-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -39,9 +39,10 @@
template<
class ErrorStepper ,
class ErrorChecker = error_checker_standard< typename ErrorStepper::state_type ,
- typename ErrorStepper::time_type ,
- typename ErrorStepper::algebra_type ,
- typename ErrorStepper::operations_type > ,
+ typename ErrorStepper::time_type ,
+ typename ErrorStepper::algebra_type ,
+ typename ErrorStepper::operations_type > ,
+ class AdjustSizePolicy = typename ErrorStepper::adjust_size_policy ,
class ErrorStepperCategory = typename ErrorStepper::stepper_category
>
class controlled_error_stepper { };
@@ -54,9 +55,10 @@
*/
template<
class ErrorStepper ,
- class ErrorChecker
+ class ErrorChecker ,
+ class AdjustSizePolicy
>
-class controlled_error_stepper< ErrorStepper , ErrorChecker , explicit_error_stepper_tag > : boost::noncopyable
+class controlled_error_stepper< ErrorStepper , ErrorChecker , AdjustSizePolicy , explicit_error_stepper_tag > : boost::noncopyable
{
public:
@@ -64,7 +66,7 @@
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::order_type order_type;
- typedef typename stepper_type::adjust_size_policy adjust_size_policy;
+ typedef AdjustSizePolicy adjust_size_policy;
typedef ErrorChecker error_checker_type;
@@ -220,9 +222,10 @@
*/
template<
class ErrorStepper ,
- class ErrorChecker
+ class ErrorChecker ,
+ class AdjustSizePolicy
>
-class controlled_error_stepper< ErrorStepper , ErrorChecker , explicit_error_stepper_fsal_tag > : boost::noncopyable
+class controlled_error_stepper< ErrorStepper , ErrorChecker , AdjustSizePolicy , explicit_error_stepper_fsal_tag > : boost::noncopyable
{
public:
@@ -230,7 +233,7 @@
typedef typename stepper_type::state_type state_type;
typedef typename stepper_type::time_type time_type;
typedef typename stepper_type::order_type order_type;
- typedef typename stepper_type::adjust_size_policy adjust_size_policy;
+ typedef AdjustSizePolicy adjust_size_policy;
typedef ErrorChecker error_checker_type;
controlled_error_stepper(
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -0,0 +1,21 @@
+# (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)
+
+import os ;
+import modules ;
+import path ;
+
+path-constant HOME : [ os.environ HOME ] ;
+path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
+
+project
+ : requirements
+ <define>BOOST_ALL_NO_LIB=1
+ <include>../../../../..
+ <include>$(BOOST_ROOT)
+ ;
+
+exe rosenbrock4
+ : rosenbrock4.cpp
+ ;
\ No newline at end of file
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.cpp 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -0,0 +1,123 @@
+/*
+ * rosenbrock4.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;
+}
+
+
+int main( int argc , char **argv )
+{
+ if( true )
+ {
+ state_type x( dim ) , xerr( dim );
+ time_type t = 0.0 , dt = 0.00001;
+
+ stepper_type stepper;
+ stepper.do_step( system< state_type > , jacobi , x , t , dt , xerr );
+ controlled_stepper_type controlled_stepper;
+
+ x[0] = 1.0 ; x[1] = 1.0 ; x[2] = 0.0;
+
+ ofstream fout( "dat/ross.dat" );
+ size_t count = 0;
+ while( t < 50.0 )
+ {
+ fout << t << "\t" << dt << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << std::endl;
+ size_t trials = 0;
+ while( trials < 100 )
+ {
+ if( controlled_stepper.try_step( system< state_type > , jacobi , x , t , dt ) != step_size_decreased )
+ break;
+ ++trials;
+ }
+ if( trials == 100 )
+ {
+ cerr << "Error : stepper did not converge! " << endl;
+ break;
+ }
+ ++count;
+ }
+ clog << "Rosenbrock : " << count << endl;
+ }
+
+
+
+
+
+ if( true )
+ {
+ typedef std::tr1::array< time_type , 3 > state_type2;
+ typedef explicit_error_rk54_ck< state_type2 > stepper_type2;
+ typedef controlled_error_stepper< stepper_type2 > controlled_stepper_type2;
+ stepper_type2 rk_stepper;
+ controlled_stepper_type2 stepper( rk_stepper );
+
+ state_type2 x = {{ 1.0 , 1.0 , 0.0 }};
+ time_type t = 0.0 , dt = 0.00001;
+ ofstream fout( "dat/rk.dat" );
+ size_t count = 0;
+ while( t < 50.0 )
+ {
+ fout << t << "\t" << dt << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << std::endl;
+ size_t trials = 0;
+ while( trials < 100 )
+ {
+ if( stepper.try_step( system< state_type2 > , x , t , dt ) != step_size_decreased )
+ break;
+ ++trials;
+ }
+ if( trials == 100 )
+ {
+ cerr << "Error : stepper did not converge! " << endl;
+ break;
+ }
+ ++count;
+ }
+ clog << "RK 54 : " << count << endl;
+ }
+
+ return 0;
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/rosenbrock4.hpp 2011-01-09 12:33:16 EST (Sun, 09 Jan 2011)
@@ -0,0 +1,216 @@
+/*
+ * rosenbrock4.hpp
+ *
+ * Created on: Jan 9, 2011
+ * Author: karsten
+ */
+
+#ifndef ROSENBROCK4_HPP_
+#define ROSENBROCK4_HPP_
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+
+#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+
+template< class Value >
+class rosenbrock4
+{
+public:
+
+ typedef Value time_type;
+ typedef boost::numeric::ublas::vector< time_type > state_type;
+ typedef boost::numeric::ublas::matrix< time_type > matrix_type;
+ typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
+
+ rosenbrock4( void )
+ {
+ }
+
+
+ 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 size_t n = x.size();
+ matrix_type jac( n , n );
+ pmatrix_type pm( n );
+ state_type dfdt( n ) , dxdt( n );
+ system( x , dxdt , t );
+ jacobi( x , jac , t , dfdt );
+
+ state_type g1( n ) , g2( n ) , g3( n ) , g4( n ) , g5( n );
+ state_type xtmp( n ) , dxdtnew( n );
+
+
+ jac *= -1.0;
+ jac += 1.0 / gamma / dt * boost::numeric::ublas::identity_matrix< time_type >( n );
+ boost::numeric::ublas::lu_factorize( jac , pm );
+
+ for( size_t i=0 ; i<n ; ++i )
+ g1[i] = dxdt[i] + dt * d1 * dfdt[i];
+ boost::numeric::ublas::lu_substitute( jac , pm , g1 );
+
+
+ for( size_t i=0 ; i<n ; ++i )
+ xtmp[i] = x[i] + a21 * g1[i];
+ system( xtmp , dxdtnew , t + c2 * dt );
+ for( size_t i=0 ; i<n ; ++i )
+ g2[i] = dxdtnew[i] + dt * d2 * dfdt[i] + c21 * g1[i] / dt;
+ boost::numeric::ublas::lu_substitute( jac , pm , g2 );
+
+
+ for( size_t i=0 ; i<n ; ++i )
+ xtmp[i] = x[i] + a31 * g1[i] + a32 * g2[i];
+ system( xtmp , dxdtnew , t + c3 * dt );
+ for( size_t i=0 ; i<n ; ++i )
+ g3[i] = dxdtnew[i] + dt * d3 * dfdt[i] + ( c31 * g1[i] + c32 * g2[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( jac , pm , g3 );
+
+
+ for( size_t i=0 ; i<n ; ++i )
+ xtmp[i] = x[i] + a41 * g1[i] + a42 * g2[i] + a43 * g3[i];
+ system( xtmp , dxdtnew , t + c4 * dt );
+ for( size_t i=0 ; i<n ; ++i )
+ g4[i] = dxdtnew[i] + dt * d4 * dfdt[i] + ( c41 * g1[i] + c42 * g2[i] + c43 * g3[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( jac , pm , g4 );
+
+
+ for( size_t i=0 ; i<n ; ++i )
+ xtmp[i] = x[i] + a51 * g1[i] + a52 * g2[i] + a53 * g3[i] + a54 * g4[i];
+ system( xtmp , dxdtnew , t + dt );
+ for( size_t i=0 ; i<n ; ++i )
+ g5[i] = dxdtnew[i] + ( c51 * g1[i] + c52 * g2[i] + c53 * g3[i] + c54 * g4[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( jac , pm , g5 );
+
+ for( size_t i=0 ; i<n ; ++i )
+ xtmp[i] += g5[i];
+ system( xtmp , dxdtnew , t + dt );
+ for( size_t i=0 ; i<n ; ++i )
+ xerr[i] = dxdtnew[i] + ( c61 * g1[i] + c62 * g2[i] + c63 * g3[i] + c64 * g4[i] + c65 * g5[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( jac , pm , xerr );
+
+ for( size_t i=0 ; i<n ; ++i )
+ xout[i] = xtmp[i] + xerr[i];
+ }
+
+ template< class System , class Jacobi >
+ void do_step( System &system , Jacobi &jacobi , state_type &x , time_type t , time_type dt , state_type &xerr )
+ {
+ state_type out( x.size() );
+ do_step( system , jacobi , x , t , out , dt , xerr );
+ x = out;
+ }
+
+
+private:
+
+};
+
+
+
+template< class Value >
+class rosenbrock4_controller
+{
+public:
+
+ typedef Value time_type;
+ typedef boost::numeric::ublas::vector< time_type > state_type;
+ typedef boost::numeric::ublas::matrix< time_type > matrix_type;
+ typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
+
+ rosenbrock4_controller( time_type atol = 1.0e-6 , time_type rtol = 1.0e-6 )
+ : m_rb4() , m_atol( atol ) , m_rtol( rtol )
+ {
+ }
+
+ time_type error( const state_type &x , const state_type &xold , const state_type &xerr )
+ {
+ const size_t n = x.size();
+ time_type err = 0.0 , sk = 0.0;
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ sk = m_atol + m_rtol * std::max( std::abs( x[i] ) , std::abs( xold[i] ) );
+ err += xerr[i] * xerr[i] / sk / sk;
+ }
+ return sqrt( err / time_type( n ) );
+ }
+
+
+ template< class System , class Jacobi >
+ boost::numeric::odeint::controlled_step_result
+ try_step( System &sys , Jacobi &jacobi , state_type &x , time_type &t , time_type &dt )
+ {
+ // hier gehts weiter
+// const size_t n = x.size();
+ return success_step_size_unchanged;
+ }
+
+// template <class D> stepperross.h
+// StepperRoss<D>::Controller::Controller() : reject(false), first_step(true) {}
+// Step size controller for fourth-order Rosenbrock method.
+// template <class D>
+// bool StepperRoss<D>::Controller::success(Doub err, Doub &h) {
+// Returns true if err 1, false otherwise. If step was successful, sets hnext to the estimated
+// optimal stepsize for the next step. If the step failed, reduces h appropriately for another try.
+// static const Doub safe=0.9,fac1=5.0,fac2=1.0/6.0;
+// Doub fac=MAX(fac2,MIN(fac1,pow(err,0.25)/safe));
+// Doub hnew=h/fac; Ensure 1=fac1 hnew=h 1=fac2.
+// if (err <= 1.0) { Step succeeded.
+// if (!first_step) { Predictive control.
+// Doub facpred=(hold/h)*pow(err*err/errold,0.25)/safe;
+// facpred=MAX(fac2,MIN(fac1,facpred));
+// fac=MAX(fac,facpred);
+// hnew=h/fac;
+// }
+// first_step=false;
+// hold=h;
+// errold=MAX(0.01,err);
+// if (reject) Don’t let step increase if last one was rejected.
+// hnew=(h >= 0.0 ? MIN(hnew,h) : MAX(hnew,h));
+// hnext=hnew;
+// reject=false;
+// return true;
+// } else { Truncation error too large, reduce stepsize.
+// h=hnew;
+// reject=true;
+// return false;
+// }
+// }
+
+
+private:
+
+ rosenbrock4< time_type > m_rb4;
+ time_type m_atol , m_rtol;
+
+};
+
+
+}
+}
+}
+
+
+#endif /* ROSENBROCK4_HPP_ */