[Boost-commit] svn:boost r57656 - in sandbox/odeint: boost/numeric boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/stuff/gsl_compare

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;