|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r68580 - in sandbox/odeint/branches/karsten: . boost/numeric/odeint/stepper boost/numeric/odeint/util libs/numeric/odeint/regression_test libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2011-01-31 07:18:09
Author: karsten
Date: 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
New Revision: 68580
URL: http://svn.boost.org/trac/boost/changeset/68580
Log:
* finishing rosenbrock4
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_permutation_matrix_resize.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/rosenbrock4.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/rosenbrock4.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/.cproject | 93 ++++++++--------
sandbox/odeint/branches/karsten/TODO | 1
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp | 7 -
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp | 217 +++++++++++++++++++++++++++++----------
sandbox/odeint/branches/karsten/boost/numeric/odeint/util/is_resizeable.hpp | 2
sandbox/odeint/branches/karsten/boost/numeric/odeint/util/size_adjuster.hpp | 71 -------------
sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_resize.hpp | 12 +-
sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile | 4
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/implicit_euler.cpp | 2
10 files changed, 221 insertions(+), 189 deletions(-)
Modified: sandbox/odeint/branches/karsten/.cproject
==============================================================================
--- sandbox/odeint/branches/karsten/.cproject (original)
+++ sandbox/odeint/branches/karsten/.cproject 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -5,52 +5,51 @@
<storageModule moduleId="org.eclipse.cdt.core.settings">
<cconfiguration id="0.1427786045">
<storageModule buildSystemId="org.eclipse.cdt.managedbuilder.core.configurationDataProvider" id="0.1427786045" moduleId="org.eclipse.cdt.core.settings" name="Default">
- <externalSettings/>
- <extensions>
- <extension id="org.eclipse.cdt.core.VCErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
- <extension id="org.eclipse.cdt.core.GCCErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
- <extension id="org.eclipse.cdt.core.GASErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
- <extension id="org.eclipse.cdt.core.GLDErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
- <extension id="org.eclipse.cdt.core.GmakeErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
- <extension id="org.eclipse.cdt.core.CWDLocator" point="org.eclipse.cdt.core.ErrorParser"/>
- </extensions>
- </storageModule>
+<externalSettings/>
+<extensions>
+<extension id="org.eclipse.cdt.core.VCErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+<extension id="org.eclipse.cdt.core.GCCErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+<extension id="org.eclipse.cdt.core.GASErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+<extension id="org.eclipse.cdt.core.GLDErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+<extension id="org.eclipse.cdt.core.MakeErrorParser" point="org.eclipse.cdt.core.ErrorParser"/>
+</extensions>
+</storageModule>
<storageModule moduleId="cdtBuildSystem" version="4.0.0">
- <configuration buildProperties="" description="" id="0.1427786045" name="Default" parent="org.eclipse.cdt.build.core.prefbase.cfg">
- <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 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">
- <outputEntries>
- <entry flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="outputPath" name=""/>
- </outputEntries>
- </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">
- <option id="org.eclipse.cdt.build.core.settings.holder.incpaths.558519945" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
- <listOptionValue builtIn="false" value="/home/karsten/boost/boost_1_45_0"/>
- </option>
- <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"/>
- </tool>
- <tool id="org.eclipse.cdt.build.core.settings.holder.780742962" name="GNU C++" superClass="org.eclipse.cdt.build.core.settings.holder">
- <option id="org.eclipse.cdt.build.core.settings.holder.incpaths.1009256344" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
- <listOptionValue builtIn="false" value="/home/karsten/boost/boost_1_45_0"/>
- </option>
- <inputType id="org.eclipse.cdt.build.core.settings.holder.inType.157834134" languageId="org.eclipse.cdt.core.g++" languageName="GNU C++" sourceContentType="org.eclipse.cdt.core.cxxSource,org.eclipse.cdt.core.cxxHeader" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
- </tool>
- <tool id="org.eclipse.cdt.build.core.settings.holder.421502995" name="GNU C" superClass="org.eclipse.cdt.build.core.settings.holder">
- <option id="org.eclipse.cdt.build.core.settings.holder.incpaths.388283580" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
- <listOptionValue builtIn="false" value="/home/karsten/boost/boost_1_45_0"/>
- </option>
- <inputType id="org.eclipse.cdt.build.core.settings.holder.inType.1344941182" languageId="org.eclipse.cdt.core.gcc" languageName="GNU C" sourceContentType="org.eclipse.cdt.core.cSource,org.eclipse.cdt.core.cHeader" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
- </tool>
- </toolChain>
- </folderInfo>
- <sourceEntries>
- <entry excluding="mkl" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
- </sourceEntries>
- </configuration>
- </storageModule>
+<configuration buildProperties="" description="" id="0.1427786045" name="Default" parent="org.eclipse.cdt.build.core.prefbase.cfg">
+<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 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">
+<outputEntries>
+<entry flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="outputPath" name=""/>
+</outputEntries>
+</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">
+<option id="org.eclipse.cdt.build.core.settings.holder.incpaths.558519945" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
+<listOptionValue builtIn="false" value=""${BOOST_ROOT}""/>
+</option>
+<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"/>
+</tool>
+<tool id="org.eclipse.cdt.build.core.settings.holder.780742962" name="GNU C++" superClass="org.eclipse.cdt.build.core.settings.holder">
+<option id="org.eclipse.cdt.build.core.settings.holder.incpaths.1009256344" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
+<listOptionValue builtIn="false" value=""${BOOST_ROOT}""/>
+</option>
+<inputType id="org.eclipse.cdt.build.core.settings.holder.inType.157834134" languageId="org.eclipse.cdt.core.g++" languageName="GNU C++" sourceContentType="org.eclipse.cdt.core.cxxSource,org.eclipse.cdt.core.cxxHeader" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
+</tool>
+<tool id="org.eclipse.cdt.build.core.settings.holder.421502995" name="GNU C" superClass="org.eclipse.cdt.build.core.settings.holder">
+<option id="org.eclipse.cdt.build.core.settings.holder.incpaths.388283580" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
+<listOptionValue builtIn="false" value=""${BOOST_ROOT}""/>
+</option>
+<inputType id="org.eclipse.cdt.build.core.settings.holder.inType.1344941182" languageId="org.eclipse.cdt.core.gcc" languageName="GNU C" sourceContentType="org.eclipse.cdt.core.cSource,org.eclipse.cdt.core.cHeader" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
+</tool>
+</toolChain>
+</folderInfo>
+<sourceEntries>
+<entry excluding="mkl" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
+</sourceEntries>
+</configuration>
+</storageModule>
<storageModule moduleId="scannerConfiguration">
<autodiscovery enabled="true" problemReportingEnabled="true" selectedProfileId="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile"/>
<profile id="org.eclipse.cdt.make.core.GCCStandardMakePerProjectProfile">
@@ -555,6 +554,6 @@
</cconfiguration>
</storageModule>
<storageModule moduleId="cdtBuildSystem" version="4.0.0">
- <project id="odeint.null.1351650296" name="odeint"/>
- </storageModule>
+<project id="odeint.null.1351650296" name="odeint"/>
+</storageModule>
</cproject>
Modified: sandbox/odeint/branches/karsten/TODO
==============================================================================
--- sandbox/odeint/branches/karsten/TODO (original)
+++ sandbox/odeint/branches/karsten/TODO 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -10,6 +10,7 @@
* discuss
* split check_concepts into check_stepper_concept, check_error_stepper_concept, check_controlled_stepper_concept
* include test/thrust in jam system, use system from
+* implicit euler, include dfdt
* include rosenbrock4 in trunk
DIFFICULT * finishing change of controlled_stepper to units
* check if rosenbrock controller and controlled_stepper can both be used with the explicit steppers
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/implicit_euler.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -17,6 +17,7 @@
#include <boost/numeric/odeint/util/size_adjuster.hpp>
#include <boost/numeric/odeint/util/matrix_vector_adjust_size.hpp>
#include <boost/numeric/odeint/util/ublas_resize.hpp>
+#include <boost/numeric/odeint/util/ublas_permutation_matrix_resize.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
@@ -27,12 +28,6 @@
namespace odeint {
-template< class T >
-struct is_resizeable< boost::numeric::ublas::permutation_matrix< T > >
-{
- struct type : public boost::true_type { };
- const static bool value = type::value;
-};
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -9,6 +9,7 @@
#define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_HPP_
#include <boost/ref.hpp>
+#include <boost/noncopyable.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
@@ -16,55 +17,146 @@
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/util/size_adjuster.hpp>
+#include <boost/numeric/odeint/util/matrix_vector_adjust_size.hpp>
+#include <boost/numeric/odeint/util/ublas_resize.hpp>
+#include <boost/numeric/odeint/util/ublas_permutation_matrix_resize.hpp>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+
namespace boost {
namespace numeric {
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
- * 6. roll out matrix_vector_adjuster
+ * 2. Interfacing for odeint, check if controlled_error_stepper can be used
+ * 3. dense output
*/
+
+
+
template< class Value >
+struct default_rosenbrock_coefficients : boost::noncopyable
+{
+ typedef Value value_type;
+
+ default_rosenbrock_coefficients( void )
+ : gamma ( 0.25 ) ,
+ d1 ( 0.25 ) , d2 ( -0.1043 ) , d3 ( 0.1035 ) , d4 ( 0.3620000000000023e-01 ) ,
+ c2 ( 0.386 ) , c3 ( 0.21 ) , c4 ( 0.63 ) ,
+ c21 ( -0.5668800000000000e+01 ) ,
+ a21 ( 0.1544000000000000e+01 ) ,
+ c31 ( -0.2430093356833875e+01 ) , c32 ( -0.2063599157091915e+00 ) ,
+ a31 ( 0.9466785280815826e+00 ) , a32 ( 0.2557011698983284e+00 ) ,
+ c41 ( -0.1073529058151375e+00 ) , c42 ( -0.9594562251023355e+01 ) , c43 ( -0.2047028614809616e+02 ) ,
+ a41 ( 0.3314825187068521e+01 ) , a42 ( 0.2896124015972201e+01 ) , a43 ( 0.9986419139977817e+00 ) ,
+ c51 ( 0.7496443313967647e+01 ) , c52 ( -0.1024680431464352e+02 ) , c53 ( -0.3399990352819905e+02 ) , c54 ( 0.1170890893206160e+02 ) ,
+ a51 ( 0.1221224509226641e+01 ) , a52 ( 0.6019134481288629e+01 ) , a53 ( 0.1253708332932087e+02 ) , a54 ( -0.6878860361058950e+00 ) ,
+ c61 ( 0.8083246795921522e+01 ) , c62 ( -0.7981132988064893e+01 ) , c63 ( -0.3152159432874371e+02 ) , c64 ( 0.1631930543123136e+02 ) , c65 ( -0.6058818238834054e+01 )
+ {}
+
+ const value_type gamma;
+ const value_type d1 , d2 , d3 , d4;
+ const value_type c2 , c3 , c4;
+ const value_type c21 ;
+ const value_type a21;
+ const value_type c31 , c32;
+ const value_type a31 , a32;
+ const value_type c41 , c42 , c43;
+ const value_type a41 , a42 , a43;
+ const value_type c51 , c52 , c53 , c54;
+ const value_type a51 , a52 , a53 , a54;
+ const value_type c61 , c62 , c63 , c64 , c65;
+};
+
+
+
+template< class Value , class Coefficients = default_rosenbrock_coefficients< Value > , class AdjustSizePolicy = adjust_size_initially_tag >
class rosenbrock4
{
+private:
+
+ void initialize( void )
+ {
+ m_matrix_adjuster.register_state( 0 , m_jac );
+
+ m_pmatrix_adjuster.register_state( 0 , m_pm );
+
+ m_state_adjuster.register_state( 0 , m_dfdt );
+ m_state_adjuster.register_state( 1 , m_dxdt );
+ m_state_adjuster.register_state( 2 , m_g1 );
+ m_state_adjuster.register_state( 3 , m_g2 );
+ m_state_adjuster.register_state( 4 , m_g3 );
+ m_state_adjuster.register_state( 5 , m_g4 );
+ m_state_adjuster.register_state( 6 , m_g5 );
+ m_state_adjuster.register_state( 7 , m_xtmp );
+ m_state_adjuster.register_state( 8 , m_dxdtnew );
+ }
+
+ void copy( const rosenbrock4 &rb )
+ {
+ m_jac = rb.m_jac;
+ m_pm = rb.m_pm;
+ m_dfdt =rb.m_dfdt;
+ m_dxdt = rb.m_dxdt;
+ m_g1 = rb.m_g1;
+ m_g2 = rb.m_g2;
+ m_g3 = rb.m_g3;
+ m_g4 = rb.m_g4;
+ m_g5 = rb.m_g5;
+ m_xtmp = rb.m_xtmp;
+ m_dxdtnew = rb.m_dxdtnew;
+ }
+
public:
typedef Value value_type;
typedef boost::numeric::ublas::vector< value_type > state_type;
typedef state_type deriv_type;
typedef value_type time_type;
- typedef boost::numeric::ublas::matrix< time_type > matrix_type;
+ typedef boost::numeric::ublas::matrix< value_type > matrix_type;
typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
+ typedef AdjustSizePolicy adjust_size_policy;
+ typedef Coefficients rosenbrock_coefficients;
+
rosenbrock4( void )
+ : m_state_adjuster() , m_matrix_adjuster() , m_pmatrix_adjuster() ,
+ m_jac() , m_pm( 1 ) ,
+ m_dfdt() , m_dxdt() ,
+ m_g1() , m_g2() , m_g3() , m_g4() , m_g5() ,
+ m_xtmp() , m_dxdtnew() ,
+ m_coef()
+ {
+ initialize();
+ }
+
+ rosenbrock4( const rosenbrock4 &rb )
+ : m_state_adjuster() , m_matrix_adjuster() , m_pmatrix_adjuster() ,
+ m_jac() , m_pm( 1 ) ,
+ m_dfdt() , m_dxdt() ,
+ m_g1() , m_g2() , m_g3() , m_g4() , m_g5() ,
+ m_xtmp() , m_dxdtnew() ,
+ m_coef()
{
+ initialize();
+ copy( rb );
+ }
+
+ rosenbrock4& operator=( const rosenbrock4 &rb )
+ {
+ copy( rb );
}
template< class System >
void do_step( System system , const state_type &x , time_type t , state_type &xout , time_type dt , state_type &xerr )
{
- const value_type gamma = 0.25;
- const value_type d1 = 0.25 , d2 = -0.1043 , d3 = 0.1035 , d4 = 0.3620000000000023e-01;
- const value_type c2 = 0.386 , c3 = 0.21 , c4 = 0.63;
- const value_type c21 = -0.5668800000000000e+01;
- const value_type a21 = 0.1544000000000000e+01;
- const value_type c31 = -0.2430093356833875e+01 , c32 = -0.2063599157091915e+00;
- const value_type a31 = 0.9466785280815826e+00 , a32 = 0.2557011698983284e+00;
- const value_type c41 = -0.1073529058151375e+00 , c42 = -0.9594562251023355e+01 , c43 = -0.2047028614809616e+02;
- const value_type a41 = 0.3314825187068521e+01 , a42 = 0.2896124015972201e+01 , a43 = 0.9986419139977817e+00;
- const value_type c51 = 0.7496443313967647e+01 , c52 = -0.1024680431464352e+02 , c53 = -0.3399990352819905e+02 , c54 = 0.1170890893206160e+02;
- const value_type a51 = 0.1221224509226641e+01 , a52 = 0.6019134481288629e+01 , a53 = 0.1253708332932087e+02 , a54 = -0.6878860361058950e+00 ;
- const value_type c61 = 0.8083246795921522e+01 , c62 = -0.7981132988064893e+01 , c63 = -0.3152159432874371e+02 , c64 = 0.1631930543123136e+02 , c65 = -0.6058818238834054e+01;
-
+ // get the systen and jacobi function
typedef typename boost::unwrap_reference< System >::type system_type;
typedef typename boost::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
typedef typename boost::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
@@ -72,67 +164,68 @@
deriv_func_type &deriv_func = sys.first;
jacobi_func_type &jacobi_func = sys.second;
+ const size_t n = x.size();
- const size_t n = x.size();
- matrix_type jac( n , n );
- pmatrix_type pm( n );
- state_type dfdt( n ) , dxdt( n );
- deriv_func( x , dxdt , t );
- jacobi_func( x , jac , t , dfdt );
+ // adjust size
+ m_matrix_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
+ m_pmatrix_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
+ for( size_t i=0 ; i<n ; ++i )
+ m_pm( i ) = i;
+ m_state_adjuster.adjust_size_by_policy( x , adjust_size_policy() );
- state_type g1( n ) , g2( n ) , g3( n ) , g4( n ) , g5( n );
- state_type xtmp( n ) , dxdtnew( n );
+ deriv_func( x , m_dxdt , t );
+ jacobi_func( x , m_jac , t , m_dfdt );
- jac *= -1.0;
- jac += 1.0 / gamma / dt * boost::numeric::ublas::identity_matrix< time_type >( n );
- boost::numeric::ublas::lu_factorize( jac , pm );
+ m_jac *= -1.0;
+ m_jac += 1.0 / m_coef.gamma / dt * boost::numeric::ublas::identity_matrix< value_type >( n );
+ boost::numeric::ublas::lu_factorize( m_jac , m_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 );
+ m_g1[i] = m_dxdt[i] + dt * m_coef.d1 * m_dfdt[i];
+ boost::numeric::ublas::lu_substitute( m_jac , m_pm , m_g1 );
for( size_t i=0 ; i<n ; ++i )
- xtmp[i] = x[i] + a21 * g1[i];
- deriv_func( xtmp , dxdtnew , t + c2 * dt );
+ m_xtmp[i] = x[i] + m_coef.a21 * m_g1[i];
+ deriv_func( m_xtmp , m_dxdtnew , t + m_coef.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 );
+ m_g2[i] = m_dxdtnew[i] + dt * m_coef.d2 * m_dfdt[i] + m_coef.c21 * m_g1[i] / dt;
+ boost::numeric::ublas::lu_substitute( m_jac , m_pm , m_g2 );
for( size_t i=0 ; i<n ; ++i )
- xtmp[i] = x[i] + a31 * g1[i] + a32 * g2[i];
- deriv_func( xtmp , dxdtnew , t + c3 * dt );
+ m_xtmp[i] = x[i] + m_coef.a31 * m_g1[i] + m_coef.a32 * m_g2[i];
+ deriv_func( m_xtmp , m_dxdtnew , t + m_coef.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 );
+ m_g3[i] = m_dxdtnew[i] + dt * m_coef.d3 * m_dfdt[i] + ( m_coef.c31 * m_g1[i] + m_coef.c32 * m_g2[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( m_jac , m_pm , m_g3 );
for( size_t i=0 ; i<n ; ++i )
- xtmp[i] = x[i] + a41 * g1[i] + a42 * g2[i] + a43 * g3[i];
- deriv_func( xtmp , dxdtnew , t + c4 * dt );
+ m_xtmp[i] = x[i] + m_coef.a41 * m_g1[i] + m_coef.a42 * m_g2[i] + m_coef.a43 * m_g3[i];
+ deriv_func( m_xtmp , m_dxdtnew , t + m_coef.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 );
+ m_g4[i] = m_dxdtnew[i] + dt * m_coef.d4 * m_dfdt[i] + ( m_coef.c41 * m_g1[i] + m_coef.c42 * m_g2[i] + m_coef.c43 * m_g3[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( m_jac , m_pm , m_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];
- deriv_func( xtmp , dxdtnew , t + dt );
+ m_xtmp[i] = x[i] + m_coef.a51 * m_g1[i] + m_coef.a52 * m_g2[i] + m_coef.a53 * m_g3[i] + m_coef.a54 * m_g4[i];
+ deriv_func( m_xtmp , m_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 );
+ m_g5[i] = m_dxdtnew[i] + ( m_coef.c51 * m_g1[i] + m_coef.c52 * m_g2[i] + m_coef.c53 * m_g3[i] + m_coef.c54 * m_g4[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( m_jac , m_pm , m_g5 );
for( size_t i=0 ; i<n ; ++i )
- xtmp[i] += g5[i];
- deriv_func( xtmp , dxdtnew , t + dt );
+ m_xtmp[i] += m_g5[i];
+ deriv_func( m_xtmp , m_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 );
+ xerr[i] = m_dxdtnew[i] + ( m_coef.c61 * m_g1[i] + m_coef.c62 * m_g2[i] + m_coef.c63 * m_g3[i] + m_coef.c64 * m_g4[i] + m_coef.c65 * m_g5[i] ) / dt;
+ boost::numeric::ublas::lu_substitute( m_jac , m_pm , xerr );
for( size_t i=0 ; i<n ; ++i )
- xout[i] = xtmp[i] + xerr[i];
+ xout[i] = m_xtmp[i] + xerr[i];
}
template< class System >
@@ -142,10 +235,20 @@
}
+
+ template< class StateType >
+ void adjust_size( const StateType &x )
+ {
+ m_state_adjuster.adjust_size( x );
+ m_matrix_adjuster.adjust_size( x );
+ m_pmatrix_adjuster.adjust_size( x );
+ }
+
+
private:
size_adjuster< state_type , 9 > m_state_adjuster;
- size_adjuster< matrix_type , 1 , matrix_vector_adjuster> m_matrix_adjuster;
+ size_adjuster< matrix_type , 1 , matrix_vector_adjust_size > m_matrix_adjuster;
size_adjuster< pmatrix_type , 1 > m_pmatrix_adjuster;
matrix_type m_jac;
@@ -154,7 +257,7 @@
state_type m_g1 , m_g2 , m_g3 , m_g4 , m_g5;
state_type m_xtmp , m_dxdtnew;
-
+ rosenbrock_coefficients m_coef;
};
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -0,0 +1,215 @@
+/*
+ * rosenbrock4_controller.hpp
+ *
+ * Created on: Jan 31, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_
+#define BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_
+
+#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+template< class Stepper >
+class rosenbrock4_controller
+{
+public:
+
+ typedef Stepper stepper_type;
+ typedef typename stepper_type::value_type value_type;
+ typedef typename stepper_type::state_type state_type;
+ typedef typename stepper_type::time_type time_type;
+ typedef typename stepper_type::deriv_type deriv_type;
+
+
+ rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() )
+ : m_rb4() , m_atol( atol ) , m_rtol( rtol ) ,
+ m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
+ m_last_rejected( false )
+ {
+ }
+
+ value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
+ {
+ const size_t n = x.size();
+ value_type err = 0.0 , sk = 0.0;
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ sk = m_atol + m_rtol * std::max( std::abs( xold[i] ) , std::abs( x[i] ) );
+ err += xerr[i] * xerr[i] / sk / sk;
+ }
+ return sqrt( err / value_type( n ) );
+ }
+
+ value_type last_error( void ) const
+ {
+ return m_err_old;
+ }
+
+
+ template< class System >
+ boost::numeric::odeint::controlled_step_result
+ try_step( System sys , state_type &x , value_type &t , value_type &dt )
+ {
+ static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
+
+ const size_t n = x.size();
+ state_type xnew( n ) , xerr( n );
+ m_rb4.do_step( sys , x , t , xnew , dt , xerr );
+ value_type err = error( xnew , x , xerr );
+
+ value_type fac = std::max( fac2 ,std::min( fac1 , std::pow( err , 0.25 ) / safe ) );
+ value_type dt_new = dt / fac;
+ if ( err <= 1.0 )
+ {
+ if( m_first_step )
+ {
+ m_first_step = false;
+ }
+ else
+ {
+ value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
+ fac_pred = std::max( fac2 , std::min( fac1 , fac_pred ) );
+ fac = std::max( fac , fac_pred );
+ dt_new = dt / fac;
+ }
+
+ m_dt_old = dt;
+ m_err_old = std::max( 0.01 , err );
+ if( m_last_rejected )
+ dt_new = ( dt >= 0.0 ? std::min( dt_new , dt ) : std::max( dt_new , dt ) );
+ t += dt;
+ dt = dt_new;
+ m_last_rejected = false;
+ x = xnew;
+ return success_step_size_increased;
+ }
+ else
+ {
+ dt = dt_new;
+ m_last_rejected = true;
+ return step_size_decreased;
+ }
+ }
+
+
+
+private:
+
+ stepper_type m_rb4;
+ value_type m_atol , m_rtol;
+ bool m_first_step;
+ value_type m_err_old , m_dt_old;
+ bool m_last_rejected;
+};
+
+
+
+
+//template< class Value >
+//class rosenbrock4_controller
+//{
+//public:
+//
+// typedef Value value_type;
+// typedef boost::numeric::ublas::vector< value_type > state_type;
+// typedef state_type deriv_type;
+// typedef value_type time_type;
+// typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+// typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
+//
+// rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 )
+// : m_rb4() , m_atol( atol ) , m_rtol( rtol ) ,
+// m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
+// m_last_rejected( false )
+// {
+// }
+//
+// value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
+// {
+// const size_t n = x.size();
+// value_type err = 0.0 , sk = 0.0;
+// for( size_t i=0 ; i<n ; ++i )
+// {
+// sk = m_atol + m_rtol * std::max( std::abs( xold[i] ) , std::abs( x[i] ) );
+// err += xerr[i] * xerr[i] / sk / sk;
+// }
+// return sqrt( err / value_type( n ) );
+// }
+//
+// value_type last_error( void ) const
+// {
+// return m_err_old;
+// }
+//
+//
+// template< class System >
+// boost::numeric::odeint::controlled_step_result
+// try_step( System sys , state_type &x , value_type &t , value_type &dt )
+// {
+// static const value_type safe = 0.9 , fac1 = 5.0 , fac2 = 1.0 / 6.0;
+//
+// const size_t n = x.size();
+// state_type xnew( n ) , xerr( n );
+// m_rb4.do_step( sys , x , t , xnew , dt , xerr );
+// value_type err = error( xnew , x , xerr );
+//
+// value_type fac = std::max( fac2 ,std::min( fac1 , std::pow( err , 0.25 ) / safe ) );
+// value_type dt_new = dt / fac;
+// if ( err <= 1.0 )
+// {
+// if( m_first_step )
+// {
+// m_first_step = false;
+// }
+// else
+// {
+// value_type fac_pred = ( m_dt_old / dt ) * pow( err * err / m_err_old , 0.25 ) / safe;
+// fac_pred = std::max( fac2 , std::min( fac1 , fac_pred ) );
+// fac = std::max( fac , fac_pred );
+// dt_new = dt / fac;
+// }
+//
+// m_dt_old = dt;
+// m_err_old = std::max( 0.01 , err );
+// if( m_last_rejected )
+// dt_new = ( dt >= 0.0 ? std::min( dt_new , dt ) : std::max( dt_new , dt ) );
+// t += dt;
+// dt = dt_new;
+// m_last_rejected = false;
+// x = xnew;
+// return success_step_size_increased;
+// }
+// else
+// {
+// dt = dt_new;
+// m_last_rejected = true;
+// return step_size_decreased;
+// }
+// }
+//
+//
+//
+//private:
+//
+// rosenbrock4< value_type > m_rb4;
+// value_type m_atol , m_rtol;
+// bool m_first_step;
+// value_type m_err_old , m_dt_old;
+// bool m_last_rejected;
+//};
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif /* BOOST_NUMERIC_ODEINT_STEPPER_ROSENBROCK4_CONTROLLER_HPP_ */
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/util/is_resizeable.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/util/is_resizeable.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/util/is_resizeable.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -10,7 +10,7 @@
#include <vector>
-#include <boost/type_traits/integral_constant.hpp> //for true_type and false_type
+#include <boost/type_traits/integral_constant.hpp>
namespace boost {
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/util/size_adjuster.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/util/size_adjuster.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/util/size_adjuster.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -121,77 +121,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*
- * really old stuff
- */
-//template< class State , class AdjustSizeImpl >
-//class size_adjuster
-//{
-//public:
-//
-// size_adjuster( AdjustSizeImpl &adjust_size_impl ) : m_is_initialized( false ) , m_adjust_size_impl( adjust_size_impl ) { }
-//
-// void adjust_size( const State &x )
-// {
-// adjust_size_by_resizeability( x , typename is_resizeable< State >::type() );
-// }
-//
-// void adjust_size_by_policy( const State &x , adjust_size_manually_tag )
-// {
-// }
-//
-// void adjust_size_by_policy( const State &x , adjust_size_initially_tag )
-// {
-// if( !m_is_initialized )
-// {
-// adjust_size_by_resizeability( x , typename is_resizeable< State >::type() );
-// m_is_initialized = true;
-// }
-// }
-//
-// void adjust_size_by_policy( const State &x , adjust_size_always_tag )
-// {
-// adjust_size_by_resizeability( x , typename is_resizeable< State >::type() );
-// }
-//
-//
-//private:
-//
-//
-// void adjust_size_by_resizeability( const State &x , boost::true_type )
-// {
-// m_adjust_size_impl( x );
-// }
-//
-// void adjust_size_by_resizeability( const State &x , boost::false_type )
-// {
-// }
-//
-//
-//private :
-//
-// bool m_is_initialized;
-// AdjustSizeImpl &m_adjust_size_impl;
-//};
-
-
} // odeint
} // numeric
} // boost
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_permutation_matrix_resize.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_permutation_matrix_resize.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -0,0 +1,33 @@
+/*
+ * ublas_permutation_matrix_resize.hpp
+ *
+ * Created on: Jan 31, 2011
+ * Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_UTIL_UBLAS_PERMUTATION_MATRIX_RESIZE_HPP_
+#define BOOST_NUMERIC_ODEINT_UTIL_UBLAS_PERMUTATION_MATRIX_RESIZE_HPP_
+
+#include <boost/type_traits/integral_constant.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+template< class T , class A >
+struct is_resizeable< boost::numeric::ublas::permutation_matrix< T , A > >
+{
+ struct type : public boost::true_type { };
+ const static bool value = type::value;
+};
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif /* BOOST_NUMERIC_ODEINT_UTIL_UBLAS_PERMUTATION_MATRIX_RESIZE_HPP_ */
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_resize.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_resize.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/util/ublas_resize.hpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -9,8 +9,8 @@
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
-#ifndef BOOST_NUMERIC_ODEINT_UBLAS_RESIZE_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_UBLAS_RESIZE_HPP_INCLUDED
+#ifndef BOOST_NUMERIC_ODEINT_UTIL_UBLAS_RESIZE_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_UTIL_UBLAS_RESIZE_HPP_INCLUDED
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
@@ -64,8 +64,8 @@
-} // odeint
-} // numeric
-} // boost
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
-#endif /* BOOST_NUMERIC_ODEINT_UBLAS_RESIZE_HPP_INCLUDED */
+#endif /* BOOST_NUMERIC_ODEINT_UTIL_UBLAS_RESIZE_HPP_INCLUDED */
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/Jamfile 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -23,4 +23,8 @@
exe dense_output_controlled_explicit_fsal
: dense_output_controlled_explicit_fsal.cpp
+ ;
+
+exe rosenbrock4
+ : rosenbrock4.cpp
;
\ No newline at end of file
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/rosenbrock4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/regression_test/rosenbrock4.cpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -0,0 +1,235 @@
+/*
+ * rosenbrock4.cpp
+ *
+ * Created on: Jan 9, 2011
+ * Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+#include <utility>
+#include <tr1/array>
+
+#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
+#include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
+#include <boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp>
+#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+
+struct lorenz
+{
+ template< class StateType >
+ void operator()( const StateType &x , StateType &dxdt , const 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];
+ }
+};
+
+struct lorenz_jacobi
+{
+ template< class State , class Matrix >
+ void operator()( const State &x , Matrix &J , const double &t , State &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;
+ }
+};
+
+
+struct stiff_system
+{
+ template< class State >
+ void operator()( const State &x , State &dxdt , double t )
+ {
+ dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ];
+ dxdt[ 1 ] = x[ 0 ];
+ }
+};
+
+struct stiff_system_jacobi
+{
+ template< class State , class Matrix >
+ void operator()( const State &x , Matrix &J , const double &t , State &dfdt )
+ {
+ J( 0 , 0 ) = -101.0;
+ J( 0 , 1 ) = -100.0;
+ J( 1 , 0 ) = 1.0;
+ J( 1 , 1 ) = 0.0;
+ dfdt[0] = 0.0;
+ dfdt[1] = 0.0;
+ }
+};
+
+
+
+
+
+void test_rosenbrock_stepper_with_lorenz( void )
+{
+ const static size_t dim = 3;
+ typedef rosenbrock4< double > stepper_type;
+ typedef stepper_type::state_type state_type;
+ typedef stepper_type::matrix_type matrix_type;
+
+ const double dt = 0.01;
+ size_t steps = 1000;
+ double x0 = -12.0 , y0 = -12.0 , z0 = 20.0;
+ state_type x( dim ) , xerr( dim );
+ double t = 0.0;
+
+ stepper_type stepper;
+ x[0] = x0 ; x[1] = y0 ; x[2] = z0;
+
+ ofstream fout( "rosenbrock_stepper_lorenz.dat" );
+ fout.precision( 14 );
+ size_t count = 0;
+ while( count < steps )
+ {
+ fout << t << " ";
+ fout << x[0] << " " << x[1] << " " << x[2] << " ";
+ fout << xerr[0] << " " << xerr[1] << " " << xerr[2] << " ";
+ fout <<std::endl;
+
+ stepper.do_step( make_pair( lorenz() , lorenz_jacobi() ) , x , t , dt , xerr );
+ ++count;
+ t += dt;
+ }
+}
+
+void rk54_ck_stepper_with_lorenz( void )
+{
+ typedef std::tr1::array< double , 3 > state_type2;
+ typedef explicit_error_rk54_ck< state_type2 > stepper_type2;
+ stepper_type2 rk_stepper;
+
+ double x0 = -12.0 , y0 = -12.0 , z0 = 20.0;
+ state_type2 x = {{ x0 , y0 , z0 }} , xerr = {{ 0.0 , 0.0 , 0.0 }};
+ double t = 0.0;
+
+ ofstream fout( "rk54_ck_stepper_with_lorenz.dat" );
+ fout.precision( 14 );
+ size_t count = 0;
+ size_t steps = 1000;
+ const double dt = 0.01;
+ 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( lorenz() , x , t , dt , xerr );
+ ++count;
+ t += dt;
+ }
+}
+
+void test_controlled_rosenbrock_with_stiff_system( void )
+{
+ const static size_t dim = 3;
+ typedef rosenbrock4< double > stepper_type;
+ typedef stepper_type::state_type state_type;
+ typedef stepper_type::matrix_type matrix_type;
+ typedef rosenbrock4_controller< stepper_type > controlled_stepper_type;
+
+ state_type x( dim ) , xerr( dim );
+ double t = 0.0 , dt = 0.00001;
+
+ controlled_stepper_type controlled_stepper;
+
+ x[0] = 1.0 ; x[1] = 1.0 ; x[2] = 0.0;
+
+ ofstream fout( "rosenbrock_controller_stiff.dat" );
+ size_t count = 0;
+ while( t < 50.0 )
+ {
+ fout << t << "\t" << dt << "\t" << controlled_stepper.last_error() << "\t";
+ fout << x[0] << "\t" << x[1] << "\t" << x[2] << "\t";
+ fout <<std::endl;
+
+ size_t trials = 0;
+ while( trials < 100 )
+ {
+ if( controlled_stepper.try_step( make_pair( stiff_system() , stiff_system_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;
+}
+
+void rk54_ck_controlled_with_stiff_system( void )
+{
+ typedef std::tr1::array< double , 3 > state_type2;
+ typedef explicit_error_rk54_ck< state_type2 > stepper_type2;
+ typedef controlled_error_stepper< stepper_type2 > controlled_stepper_type2;
+ controlled_stepper_type2 stepper;
+
+ state_type2 x = {{ 1.0 , 1.0 , 0.0 }};
+ double t = 0.0 , dt = 0.00001;
+ ofstream fout( "rk54_ck_controller_stiff.dat" );
+ size_t count = 0;
+ while( t < 50.0 )
+ {
+ 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 )
+ {
+ if( stepper.try_step( stiff_system() , 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;
+}
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+ test_rosenbrock_stepper_with_lorenz();
+ rk54_ck_stepper_with_lorenz();
+ test_controlled_rosenbrock_with_stiff_system();
+ rk54_ck_controlled_with_stiff_system();
+
+ return 0;
+}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/Jamfile 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -27,5 +27,6 @@
[ run stepper_with_units.cpp ]
[ run stepper_copying.cpp ]
[ run stepper_with_ranges.cpp ]
+ [ run rosenbrock4.cpp ]
: <testing.launcher>valgrind
;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/implicit_euler.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/implicit_euler.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/implicit_euler.cpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -17,7 +17,7 @@
#include <boost/test/unit_test.hpp>
#include <boost/numeric/odeint/stepper/implicit_euler.hpp>
-#include <boost/numeric/odeint/algebra/external/ublas_resize.hpp>
+#include <boost/numeric/odeint/util/ublas_resize.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/rosenbrock4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/rosenbrock4.cpp 2011-01-31 07:18:05 EST (Mon, 31 Jan 2011)
@@ -0,0 +1,99 @@
+/* Boost check_implicit_euler.cpp test file
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ This file tests the use of the euler stepper
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE odeint_rosenbrock4
+
+#include <utility>
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef double value_type;
+typedef boost::numeric::ublas::vector< value_type > state_type;
+typedef boost::numeric::ublas::matrix< value_type > matrix_type;
+
+
+void sys( const state_type &x , state_type &dxdt , const value_type &t )
+{
+ dxdt( 0 ) = x( 0 ) + 2 * x( 1 );
+ dxdt( 1 ) = x( 1 );
+}
+
+void jacobi( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt )
+{
+ jacobi( 0 , 0 ) = 1;
+ jacobi( 0 , 1 ) = 2;
+ jacobi( 1 , 0 ) = 0;
+ jacobi( 1 , 1 ) = 1;
+ dfdt( 0 ) = 0.0;
+ dfdt( 1 ) = 0.0;
+}
+
+BOOST_AUTO_TEST_SUITE( rosenbrock4_test )
+
+BOOST_AUTO_TEST_CASE( test_rosenbrock4_stepper )
+{
+ typedef rosenbrock4< value_type > stepper_type;
+ stepper_type stepper;
+
+ typedef stepper_type::state_type state_type;
+ typedef stepper_type::value_type stepper_value_type;
+ typedef stepper_type::deriv_type deriv_type;
+ typedef stepper_type::time_type time_type;
+
+ state_type x( 2 ) , xerr( 2 );
+ x(0) = 0.0; x(1) = 1.0;
+
+
+
+ stepper.do_step( std::make_pair( sys , jacobi ) , x , 0.0 , 0.1 , xerr );
+
+// using std::abs;
+// value_type eps = 1E-12;
+//
+// // compare with analytic solution of above system
+// BOOST_CHECK_SMALL( abs( x(0) - 20.0/81.0 ) , eps );
+// BOOST_CHECK_SMALL( abs( x(1) - 10.0/9.0 ) , eps );
+
+}
+
+BOOST_AUTO_TEST_CASE( test_rosenbrock4_controller )
+{
+
+}
+
+//BOOST_AUTO_TEST_CASE( test_euler )
+//{
+// implicit_euler< value_type > stepper;
+// state_type x( 2 );
+// x(0) = 0.0; x(1) = 1.0;
+//
+// value_type eps = 1E-12;
+//
+// stepper.do_step( std::make_pair( sys , jacobi ) , x , 0.0 , 0.1 );
+//
+// using std::abs;
+//
+// // compare with analytic solution of above system
+// BOOST_CHECK_MESSAGE( abs( x(0) - 20.0/81.0 ) < eps , x(0) - 20.0/81.0 );
+// BOOST_CHECK_MESSAGE( abs( x(1) - 10.0/9.0 ) < eps , x(0) - 10.0/9.0 );
+//
+//}
+
+BOOST_AUTO_TEST_SUITE_END()
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