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

Subject: [Boost-commit] svn:boost r57632 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples libs/numeric/odeint/stuff/gsl_compare
From: mario.mulansky_at_[hidden]
Date: 2009-11-13 07:31:30

Author: mariomulansky
Date: 2009-11-13 07:31:29 EST (Fri, 13 Nov 2009)
New Revision: 57632
URL: http://svn.boost.org/trac/boost/changeset/57632

Log:
added RK4 Cash Karp method
updated gsl compare example
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/observer.hpp                            |     2                                         
   sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp                         |   134 ++++++++++++++++++++++++++++++++++++++++
   sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp    |     5 +                                       
   sandbox/odeint/libs/numeric/odeint/stuff/gsl_compare/lorenz_stepper_cmp.cpp |    38 +++++++++-                              
   4 files changed, 172 insertions(+), 7 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint/observer.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/observer.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/observer.hpp	2009-11-13 07:31:29 EST (Fri, 13 Nov 2009)
@@ -21,7 +21,7 @@
 
 
     template< class Time , class Container , class System >
-    void do_nothing_observer( Time , Container& , System& )
+    inline void do_nothing_observer( Time , Container& , System& )
     {
     }
     

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-13 07:31:29 EST (Fri, 13 Nov 2009)
@@ -27,6 +27,22 @@
 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 ,
@@ -66,6 +82,7 @@
         container_type m_dxt;
         container_type m_dxm;
         container_type m_xt;
+        container_type m_x4, m_x5, m_x6;
         resizer_type m_resizer;
 
         
@@ -128,6 +145,123 @@
         }
 
 
+        /* RK4 step with error estimation to 5th order according to Cash Karp */
+        template< class DynamicalSystem >
+        void next_step( DynamicalSystem &system ,
+                        container_type &x ,
+                        const 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) );
+
+            system( m_xt , m_dxt , t + dt*time_type(rk4_a2) ); // m_dxt = nr_ak2
+            iterator x_i = x.begin();
+            iterator m_xt_i = m_xt.begin();
+            typename container_type::const_iterator dxdt_i = dxdt.begin();
+            iterator m_dxt_i = m_dxt.begin();
+            while( x_i != x.end() ) {
+                *m_xt_i++ = (*x_i++) + dt*( time_type(rk4_b31)*(*dxdt_i++) + 
+                                            time_type(rk4_b32)*(*m_dxt_i++) );
+            }
+
+            system( m_xt , m_dxm , t + dt*time_type(rk4_a3) ); // m_dxm = nr_ak3
+            x_i = x.begin();
+            m_xt_i = m_xt.begin();
+            dxdt_i = dxdt.begin();
+            m_dxt_i = m_dxt.begin();
+            iterator m_dxm_i = m_dxm.begin();
+            while( x_i != x.end() ) {
+                *m_xt_i++ = (*x_i++) + dt*( time_type(rk4_b41)*(*dxdt_i++) + 
+                                            time_type(rk4_b42)*(*m_dxt_i++) + 
+                                            time_type(rk4_b43)*(*m_dxm_i++) );
+            }
+
+            system( m_xt, m_x4 , t + dt*time_type(rk4_a4) ); // m_x4 = nr_ak4
+            x_i = x.begin();
+            m_xt_i = m_xt.begin();
+            dxdt_i = dxdt.begin();
+            m_dxt_i = m_dxt.begin();
+            m_dxm_i = m_dxm.begin();
+            iterator m_x4_i = m_x4.begin();
+            while( x_i != x.end() ) {
+                *m_xt_i++ = (*x_i++) + dt*( time_type(rk4_b51)*(*dxdt_i++) + 
+                                            time_type(rk4_b52)*(*m_dxt_i++) +
+                                            time_type(rk4_b53)*(*m_dxm_i++) + 
+                                            time_type(rk4_b54)*(*m_x4_i++) ) ;
+            }
+
+            system( m_xt , m_x5 , t + dt*time_type(rk4_a5) ); // m_x5 = nr_ak5
+            x_i = x.begin();
+            m_xt_i = m_xt.begin();
+            dxdt_i = dxdt.begin();
+            m_dxt_i = m_dxt.begin();
+            m_dxm_i = m_dxm.begin();
+            m_x4_i = m_x4.begin();
+            iterator m_x5_i = m_x5.begin();
+            while( x_i != x.end() ) {
+                *m_xt_i++ = (*x_i++) + dt*( time_type(rk4_b61)*(*dxdt_i++) +
+                                          time_type(rk4_b62)*(*m_dxt_i++) +
+                                          time_type(rk4_b63)*(*m_dxm_i++) +
+                                          time_type(rk4_b64)*(*m_x4_i++) +
+                                          time_type(rk4_b65)*(*m_x5_i++) );
+            }
+
+            system( m_xt , m_x6 , t + dt*time_type(rk4_a6) ); // m_x6 = nr_ak6
+            x_i = x.begin();
+            dxdt_i = dxdt.begin();
+            m_dxm_i = m_dxm.begin();
+            m_x4_i = m_x4.begin();
+            iterator m_x6_i = m_x6.begin();
+            while( x_i != x.end() ) {
+                (*x_i++) += dt*( time_type(rk4_c1)*(*dxdt_i++) +
+                                 time_type(rk4_c3)*(*m_dxm_i++) +
+                                 time_type(rk4_c4)*(*m_x4_i++) +
+                                 time_type(rk4_c6)*(*m_x6_i++) );
+            }
+
+            // error estimate
+            iterator xerr_i = xerr.begin();
+            dxdt_i = dxdt.begin();
+            m_dxm_i = m_dxm.begin();
+            m_x4_i = m_x4.begin();
+            m_x5_i = m_x5.begin();
+            m_x6_i = m_x6.begin();
+            while( xerr_i != xerr.end() ) {
+                *xerr_i++ = dt*( time_type(rk4_dc1)*(*dxdt_i++) +
+                                  time_type(rk4_dc3)*(*m_dxm_i++) +
+                                  time_type(rk4_dc4)*(*m_x4_i++) +
+                                  time_type(rk4_dc5)*(*m_x5_i++) +
+                                  time_type(rk4_dc6)*(*m_x6_i++) );
+            }
+        }
+
+        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 );
+        }
+
+
+
     };
 
 } // namespace odeint

Modified: sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp	2009-11-13 07:31:29 EST (Fri, 13 Nov 2009)
@@ -6,7 +6,10 @@
  Shows the usage of odeint by integrating the equations of a 
  pendulum with horizontally vibrating pivot:
 
- 
+ d^2x/dt^2 + sin(x) + alpha*x = a*omega^2*sin(omega*t)*cos(x)
+
+ for large enough omega >sim 20 two new fixpoints (of the 
+ slow dynamics) arise, that can be seen in the simulations as well
 
  Distributed under the Boost Software License, Version 1.0.
 (See accompanying file LICENSE_1_0.txt or

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-13 07:31:29 EST (Fri, 13 Nov 2009)
@@ -105,7 +105,7 @@
 
 
 
-    // odeint method
+    // odeint method 1 - rk4 with half stepsize error estimation
     state_type x1 = x_start , x1_err;
     stepper_half_step< stepper_rk4< state_type > > stepper;
 
@@ -114,12 +114,20 @@
     for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
         stepper.next_step( lorenz , x1 , t , dt , x1_err );
     end = clock();
-    clog << "odeint array : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
-
+    clog << "odeint half step array : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
 
+    // odeint method 2 - rk45 with Cash Karp error estimation
+    x1 = x_start;
+    stepper_rk4< state_type > stepper_cash_karp;
+    start = clock();
+    t = 0.0;
+    for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
+        stepper_cash_karp.next_step( lorenz , x1 , t , dt , x1_err );
+    end = clock();
+    clog << "odeint Cash-Karp array : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
 
 
-    // gsl method
+    // gsl method rk4
     gsl_odeiv_step *s = gsl_odeiv_step_alloc( gsl_odeiv_step_rk4 , 3);
     gsl_odeiv_system sys = { lorenz_gsl , 0 , 3 , 0 };
     double x2[3] = { x_start[0] , x_start[1] , x_start[2] } , x2_err[3];
@@ -130,6 +138,19 @@
         gsl_odeiv_step_apply ( s , t , dt , x2 , x2_err , 0 , 0 , &sys );
     end = clock();
     clog << "gsl rk4 : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+
+
+    // gsl method rk4 cash-karp
+    gsl_odeiv_step *s_rkck = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck , 3);
+    //    sys = { lorenz_gsl , 0 , 3 , 0 };
+    copy( x_start.begin() , x_start.end() , x2 );
+
+    start= clock();
+    t = 0.0;
+    for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
+        gsl_odeiv_step_apply ( s_rkck , t , dt , x2 , x2_err , 0 , 0 , &sys );
+    end = clock();
+    clog << "gsl rk4 Cash-Karp : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
      
 
 
@@ -138,6 +159,9 @@
     x1 = x_start;
     copy( x_start.begin() , x_start.end() , x2 );
     double x3[3] = { x_start[0] , x_start[1] , x_start[2] };
+
+    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.flags( ios::scientific );
 
@@ -148,7 +172,11 @@
         gsl_odeiv_step_apply ( s , t , dt , x2 , x2_err , 0 , 0 , &sys );
         rk4_lorenz( x3 , 0.5*dt );
         rk4_lorenz( x3 , 0.5*dt );
-        cout << t << tab << x1[0] << tab << x2[0] << tab << x3[0] << endl;
+        cout << t << tab << x1[0] << tab <<  x2[0] << tab << x3[0] << tab;
+        // compare cash-karp methods
+        stepper_cash_karp.next_step( lorenz , x_odeint_rkck , t , dt , x1_err );
+        gsl_odeiv_step_apply ( s_rkck , t , dt , x_gsl_rkck , x2_err , 0 , 0 , &sys );
+        cout << x_odeint_rkck[0] << tab << x_gsl_rkck[0] << endl;
     }
 
     gsl_odeiv_step_free (s);