Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r57656 - in sandbox/odeint: boost/numeric boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/stuff/gsl_compare
From: mario.mulansky_at_[hidden]
Date: 2009-11-14 13:00:11


Author: mariomulansky
Date: 2009-11-14 13:00:10 EST (Sat, 14 Nov 2009)
New Revision: 57656
URL: http://svn.boost.org/trac/boost/changeset/57656

Log:
moved cash karp solver to separate class (it's an 5th order solver)
Added:
   sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp (contents, props changed)
Text files modified:
   sandbox/odeint/boost/numeric/odeint.hpp | 1
   sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp | 2
   sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp | 104 ---------------------------------------
   sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp | 4
   4 files changed, 5 insertions(+), 106 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2009-11-14 13:00:10 EST (Sat, 14 Nov 2009)
@@ -17,6 +17,7 @@
 
 #include <boost/numeric/odeint/stepper_euler.hpp>
 #include <boost/numeric/odeint/stepper_rk4.hpp>
+#include <boost/numeric/odeint/stepper_rk5_ck.hpp>
 #include <boost/numeric/odeint/stepper_half_step.hpp>
 
 #include <boost/numeric/odeint/stepsize_controller_standard.hpp>

Modified: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp 2009-11-14 13:00:10 EST (Sat, 14 Nov 2009)
@@ -99,7 +99,7 @@
                 ( (*first2++) + (*first3++) + alpha2*(*first4++) );
     }
 
- // computes y = x1 + alpha2*x2
+ // computes y = alpha1*x1 + alpha2*x2
     template <
         class OutputIterator ,
         class InputIterator ,

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp 2009-11-14 13:00:10 EST (Sat, 14 Nov 2009)
@@ -15,35 +15,17 @@
 */
 
 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_RK4_HPP
-#define BOOST_NUMERIC_ODEINT_STEPPER_RK_4_HPP
+#define BOOST_NUMERIC_ODEINT_STEPPER_RK4_HPP
 
 #include <boost/concept_check.hpp>
 
 #include <boost/numeric/odeint/concepts/state_concept.hpp>
 #include <boost/numeric/odeint/resizer.hpp>
-#include <iostream>
 
 namespace boost {
 namespace numeric {
 namespace odeint {
 
-namespace detail {
-namespace constants {
- // constants for rk4 cash karp
- static const double rk4_a2 = 0.2, rk4_a3=0.3, rk4_a4 = 0.6, rk4_a5 = 1.0,
- rk4_a6 = 0.875, rk4_b21 = 0.2, rk4_b31 = 3.0/40.0, rk4_b32 = 9.0/40.0,
- rk4_b41 = 0.3, rk4_b42 = -0.9, rk4_b43 = 1.2, rk4_b51 = -11.0/54.0,
- rk4_b52 = 2.5, rk4_b53 = -70.0/27.0, rk4_b54 = 35.0/27.0,
- rk4_b61 = 1631.0/55296.0, rk4_b62 = 175.0/512.0, rk4_b63 = 575.0/13824.0,
- rk4_b64 = 44275.0/110592.0, rk4_b65 = 253.0/4096.0,
- rk4_c1 = 37.0/378.0, rk4_c3 = 250.0/621.0, rk4_c4 = 125.0/594.0,
- rk4_c6 = 512.0/1771.0, rk4_dc5 = -277.0/14336.0,
- rk4_dc1 = rk4_c1 - 2825.0/27648, rk4_dc3 = rk4_c3 - 18575.0/48384.0,
- rk4_dc4 = rk4_c4 - 13525.0/55296.0, rk4_dc6 = rk4_c6 - 0.25;
-}
-}
-
-
     template<
         class Container ,
         class Time = double ,
@@ -146,90 +128,6 @@
 
 
         /* RK4 step with error estimation to 5th order according to Cash Karp */
- template< class DynamicalSystem >
- void next_step( DynamicalSystem &system ,
- container_type &x ,
- container_type &dxdt ,
- time_type t ,
- time_type dt ,
- container_type &xerr )
- {
- using namespace detail::it_algebra;
- using namespace detail::constants;
- m_resizer.adjust_size( x , m_dxt );
- m_resizer.adjust_size( x , m_dxm );
- m_resizer.adjust_size( x , m_xt );
- m_resizer.adjust_size( x , m_x4 );
- m_resizer.adjust_size( x , m_x5 );
- m_resizer.adjust_size( x , m_x6 );
-
- //m_xt = x + dt*b21*dxdt
- assign_sum( m_xt.begin() , m_xt.end() , x.begin() , dxdt.begin() ,
- dt*time_type(rk4_b21) );
-
- time_type t_1 = time_type(1.0);
-
- system( m_xt , m_dxt , t + dt*time_type(rk4_a2) ); // m_dxt = nr_ak2
- // m_xt = x + dt*rk4_b31*dxt + dt*rk4_b32*m_dxt
- scale_sum( m_xt.begin(), m_xt.end(),
- t_1, x.begin(),
- dt*time_type(rk4_b31), dxdt.begin(),
- dt*time_type(rk4_b32), m_dxt.begin() );
-
- system( m_xt , m_dxm , t + dt*time_type(rk4_a3) ); // m_dxm = nr_ak3
- scale_sum( m_xt.begin(), m_xt.end(),
- t_1, x.begin(),
- dt*time_type(rk4_b41), dxdt.begin(),
- dt*time_type(rk4_b42), m_dxt.begin(),
- dt*time_type(rk4_b43), m_dxm.begin() );
-
- system( m_xt, m_x4 , t + dt*time_type(rk4_a4) ); // m_x4 = nr_ak4
- scale_sum( m_xt.begin(), m_xt.end(),
- t_1, x.begin(),
- dt*time_type(rk4_b51), dxdt.begin(),
- dt*time_type(rk4_b52), m_dxt.begin(),
- dt*time_type(rk4_b53), m_dxm.begin(),
- dt*time_type(rk4_b54), m_x4.begin() );
-
- system( m_xt , m_x5 , t + dt*time_type(rk4_a5) ); // m_x5 = nr_ak5
- scale_sum( m_xt.begin(), m_xt.end(),
- t_1, x.begin(),
- dt*time_type(rk4_b61), dxdt.begin(),
- dt*time_type(rk4_b62), m_dxt.begin(),
- dt*time_type(rk4_b63), m_dxm.begin(),
- dt*time_type(rk4_b64), m_x4.begin(),
- dt*time_type(rk4_b65), m_x5.begin() );
-
- system( m_xt , m_x6 , t + dt*time_type(rk4_a6) ); // m_x6 = nr_ak6
- scale_sum( x.begin(), x.end(),
- t_1, x.begin(),
- dt*time_type(rk4_c1), dxdt.begin(),
- dt*time_type(rk4_c3), m_dxm.begin(),
- dt*time_type(rk4_c4), m_x4.begin(),
- dt*time_type(rk4_c6), m_x6.begin() );
-
- // error estimate
- scale_sum(xerr.begin(), xerr.end(),
- dt*time_type(rk4_dc1), dxdt.begin(),
- dt*time_type(rk4_dc3), m_dxm.begin(),
- dt*time_type(rk4_dc4), m_x4.begin(),
- dt*time_type(rk4_dc5), m_x5.begin(),
- dt*time_type(rk4_dc6), m_x6.begin() );
- }
-
- template< class DynamicalSystem >
- void next_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt ,
- container_type &xerr )
- {
- m_resizer.adjust_size( x , m_dxdt );
- system( x , m_dxdt , t );
- next_step( system , x , m_dxdt , t , dt , xerr );
- }
-
-
 
     };
 

Added: sandbox/odeint/boost/numeric/odeint/stepper_rk5_ck.hpp
==============================================================================
Binary file. No diff available.

Modified: sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp 2009-11-14 13:00:10 EST (Sat, 14 Nov 2009)
@@ -118,7 +118,7 @@
 
     // odeint method 2 - rk45 with Cash Karp error estimation
     x1 = x_start;
- stepper_rk4< state_type > stepper_cash_karp;
+ stepper_rk5_ck< state_type > stepper_cash_karp;
     start = clock();
     t = 0.0;
     for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
@@ -162,7 +162,7 @@
 
     state_type x_odeint_rkck = x_start;
     double x_gsl_rkck[3] = { x_start[0] , x_start[1] , x_start[2] };
- cout.precision( 14 );
+ cout.precision( 16 );
     cout.flags( ios::scientific );
 
     t = 0.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