[Boost-commit] svn:boost r57966 - in sandbox/odeint: boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/examples

Subject: [Boost-commit] svn:boost r57966 - in sandbox/odeint: boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-11-27 03:05:00

Author: karsten
Date: 2009-11-27 03:04:59 EST (Fri, 27 Nov 2009)
New Revision: 57966
URL: http://svn.boost.org/trac/boost/changeset/57966

Log:
rk78 implementation

Added:
   sandbox/odeint/libs/numeric/odeint/examples/dnls_stepper_compare.cpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp |   186 +++++++++++++++++++++++++++++++++++++++ 
   sandbox/odeint/boost/numeric/odeint/stepper_rk4.hpp             |     2                                         
   sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp   |   186 ++++++++++++++++++++++++++++++++++++++- 
   sandbox/odeint/libs/numeric/odeint/examples/Jamfile             |     1                                         
   4 files changed, 361 insertions(+), 14 deletions(-)

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-27 03:04:59 EST (Fri, 27 Nov 2009)
@@ -102,15 +102,16 @@
     // computes y = alpha1*x1 + alpha2*x2
     template <
         class OutputIterator ,
-        class InputIterator ,
+        class InputIterator1 ,
+        class InputIterator2 ,
         class T
         >
     inline void scale_sum( OutputIterator y_begin ,
                            OutputIterator y_end ,
                            T alpha1 ,
-                           InputIterator x1_begin ,
+                           InputIterator1 x1_begin ,
                            T alpha2 ,
-                           InputIterator x2_begin )
+                           InputIterator2 x2_begin )
     {
         while( y_begin != y_end )
             (*y_begin++) = alpha1 * (*x1_begin++) + alpha2 * (*x2_begin++);
@@ -193,7 +194,8 @@
     }
 
 
-    // computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5 + alpha6*x6
+    // computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
+    // + alpha6*x6
     template <
         class OutputIterator ,
         class InputIterator ,
@@ -224,6 +226,182 @@
                 alpha6 * (*x6_begin++);
     }
 
+
+    // computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
+    // + alpha6*x6 + alpha7*x7
+    template <
+        class OutputIterator ,
+        class InputIterator ,
+        class T
+        >
+    inline void scale_sum( OutputIterator y_begin ,
+                           OutputIterator y_end ,
+                           T alpha1 ,
+                           InputIterator x1_begin ,
+                           T alpha2 ,
+                           InputIterator x2_begin ,
+                           T alpha3 ,
+                           InputIterator x3_begin ,
+                           T alpha4 ,
+                           InputIterator x4_begin ,
+                           T alpha5 ,
+                           InputIterator x5_begin ,
+                           T alpha6 ,
+                           InputIterator x6_begin ,
+                           T alpha7 ,
+                           InputIterator x7_begin )
+
+    {
+        while( y_begin != y_end )
+            (*y_begin++) = 
+                alpha1 * (*x1_begin++) + 
+                alpha2 * (*x2_begin++) + 
+                alpha3 * (*x3_begin++) +
+                alpha4 * (*x4_begin++) +
+                alpha5 * (*x5_begin++) +
+                alpha6 * (*x6_begin++) +
+                alpha7 * (*x7_begin++) ;
+    }
+
+
+
+    // computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
+    // + alpha6*x6 + alpha7*x7 + alpha8*x8
+    template <
+        class OutputIterator ,
+        class InputIterator ,
+        class T
+        >
+    inline void scale_sum( OutputIterator y_begin ,
+                           OutputIterator y_end ,
+                           T alpha1 ,
+                           InputIterator x1_begin ,
+                           T alpha2 ,
+                           InputIterator x2_begin ,
+                           T alpha3 ,
+                           InputIterator x3_begin ,
+                           T alpha4 ,
+                           InputIterator x4_begin ,
+                           T alpha5 ,
+                           InputIterator x5_begin ,
+                           T alpha6 ,
+                           InputIterator x6_begin ,
+                           T alpha7 ,
+                           InputIterator x7_begin ,
+                           T alpha8 ,
+                           InputIterator x8_begin )
+
+
+    {
+        while( y_begin != y_end )
+            (*y_begin++) = 
+                alpha1 * (*x1_begin++) + 
+                alpha2 * (*x2_begin++) + 
+                alpha3 * (*x3_begin++) +
+                alpha4 * (*x4_begin++) +
+                alpha5 * (*x5_begin++) +
+                alpha6 * (*x6_begin++) +
+                alpha7 * (*x7_begin++) +
+                alpha8 * (*x8_begin++);
+    }
+
+    // computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
+    // + alpha6*x6 + alpha7*x7 + alpha8*x8 + alpha9*x9
+    template <
+        class OutputIterator ,
+        class InputIterator ,
+        class T
+        >
+    inline void scale_sum( OutputIterator y_begin ,
+                           OutputIterator y_end ,
+                           T alpha1 ,
+                           InputIterator x1_begin ,
+                           T alpha2 ,
+                           InputIterator x2_begin ,
+                           T alpha3 ,
+                           InputIterator x3_begin ,
+                           T alpha4 ,
+                           InputIterator x4_begin ,
+                           T alpha5 ,
+                           InputIterator x5_begin ,
+                           T alpha6 ,
+                           InputIterator x6_begin ,
+                           T alpha7 ,
+                           InputIterator x7_begin ,
+                           T alpha8 ,
+                           InputIterator x8_begin ,
+                           T alpha9 ,
+                           InputIterator x9_begin )
+
+
+    {
+        while( y_begin != y_end )
+            (*y_begin++) = 
+                alpha1 * (*x1_begin++) + 
+                alpha2 * (*x2_begin++) + 
+                alpha3 * (*x3_begin++) +
+                alpha4 * (*x4_begin++) +
+                alpha5 * (*x5_begin++) +
+                alpha6 * (*x6_begin++) +
+                alpha7 * (*x7_begin++) +
+                alpha8 * (*x8_begin++) +
+                alpha9 * (*x9_begin++);
+    }
+
+
+
+
+    // computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
+    // + alpha6*x6 + alpha7*x7 + alpha8*x8 + alpha9*x9 + alpha10*x10
+    template <
+        class OutputIterator ,
+        class InputIterator ,
+        class T
+        >
+    inline void scale_sum( OutputIterator y_begin ,
+                           OutputIterator y_end ,
+                           T alpha1 ,
+                           InputIterator x1_begin ,
+                           T alpha2 ,
+                           InputIterator x2_begin ,
+                           T alpha3 ,
+                           InputIterator x3_begin ,
+                           T alpha4 ,
+                           InputIterator x4_begin ,
+                           T alpha5 ,
+                           InputIterator x5_begin ,
+                           T alpha6 ,
+                           InputIterator x6_begin ,
+                           T alpha7 ,
+                           InputIterator x7_begin ,
+                           T alpha8 ,
+                           InputIterator x8_begin ,
+                           T alpha9 ,
+                           InputIterator x9_begin ,
+                           T alpha10 ,
+                           InputIterator x10_begin 
+        )
+    {
+        while( y_begin != y_end )
+            (*y_begin++) = 
+                alpha1 * (*x1_begin++) + 
+                alpha2 * (*x2_begin++) + 
+                alpha3 * (*x3_begin++) +
+                alpha4 * (*x4_begin++) +
+                alpha5 * (*x5_begin++) +
+                alpha6 * (*x6_begin++) +
+                alpha7 * (*x7_begin++) +
+                alpha8 * (*x8_begin++) +
+                alpha9 * (*x9_begin++) +
+                alpha10 * (*x10_begin++);
+    }
+
+
+
+
+
+
+
     // generic version for n values
     template <
         class OutputIterator ,

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-27 03:04:59 EST (Fri, 27 Nov 2009)
@@ -121,7 +121,7 @@
                        dt, m_dxm.begin() );
 
             // dt * m_dxh = k4
-            system( m_xt , m_dxh , value_type( t + dt ) );  
+            system( m_xt , m_dxh , t + dt );  
             //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
             scale_sum( x.begin(), x.end(),
                        val1, x.begin(),

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk78_fehlberg.hpp	2009-11-27 03:04:59 EST (Fri, 27 Nov 2009)
@@ -44,7 +44,7 @@
         container_type m_dxdt;
         container_type m_xt;
         container_type m_k02 , m_k03 , m_k04 , m_k05 , m_k06 , m_k07 ,
-            m_k08 , m_k09 , m_k10 , m_k11 , m_k12;
+            m_k08 , m_k09 , m_k10 , m_k11 , m_k12 , m_k13;
 
         resizer_type m_resizer;
 
@@ -64,7 +64,7 @@
         template< class DynamicalSystem >
         void do_step( DynamicalSystem &system ,
                         container_type &x ,
-                        const container_type &dxdt ,
+                        container_type &dxdt ,
                         time_type t ,
                         time_type dt )
         {
@@ -87,7 +87,6 @@
             const time_type c10 = static_cast<time_type>( c09 );
             const time_type c12 = static_cast<time_type>( 41.0 / 840.0 );
             const time_type c13 = static_cast<time_type>( c12 );
-            const time_type c11 = static_cast<time_type>( 0. );
 
             // the coefficients for each step
             const time_type b02_01 = static_cast<time_type>( 2.0 / 27.0 );
@@ -188,38 +187,207 @@
             m_resizer.adjust_size( x , m_k10 );
             m_resizer.adjust_size( x , m_k11 );
             m_resizer.adjust_size( x , m_k12 );
+            m_resizer.adjust_size( x , m_k13 );
 
 
             // k1, the first system call has allready been evaluated
+
+            // k2 step
             scale_sum( m_xt.begin() , m_xt.end() , 
                        val1 , x.begin() , 
                        dt * b02_01 , dxdt.begin() );
-
-            // k2 step
             system( m_xt , m_k02 , m_t02 );
+
+            // k3 step
             scale_sum( m_xt.begin() , m_xt.end() ,
                        val1 , x.begin() ,
                        dt * b03_01 , dxdt.begin() ,
-                       dt * b03_02 , m_k02 );
+                       dt * b03_02 , m_k02.begin() );
+            system( m_xt , m_k03 , m_t03 );
+
+
+            // k4 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b04_01 , dxdt.begin() ,
+                       dt * b04_03 , m_k03.begin() );
+            system( m_xt , m_k04 , m_t04 );
+
+
+            // k5 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b05_01 , dxdt.begin() ,
+                       dt * b05_03 , m_k03.begin() ,
+                       dt * b05_04 , m_k04.begin() );
+            system( m_xt , m_k05 , m_t05 );
+
+
+            // k6 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b06_01 , dxdt.begin() ,
+                       dt * b06_04 , m_k04.begin() ,
+                       dt * b06_05 , m_k05.begin() );
+            system( m_xt , m_k06 , m_t06 );
+
+
+            // k7 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b07_01 , dxdt.begin() ,
+                       dt * b07_04 , m_k04.begin() ,
+                       dt * b07_05 , m_k05.begin() ,
+                       dt * b07_06 , m_k06.begin() );
+            system( m_xt , m_k07 , m_t07 );
 
 
+            // k8 step 
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b08_01 , dxdt.begin() ,
+                       dt * b08_05 , m_k05.begin() ,
+                       dt * b08_06 , m_k06.begin() ,
+                       dt * b08_07 , m_k07.begin() );
+            system( m_xt , m_k08 , m_t08 );
+
+
+            // k9 step 
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b09_01 , dxdt.begin() ,
+                       dt * b09_04 , m_k04.begin() ,
+                       dt * b09_05 , m_k05.begin() ,
+                       dt * b09_06 , m_k06.begin() ,
+                       dt * b09_07 , m_k07.begin() ,
+                       dt * b09_08 , m_k08.begin() );
+            system( m_xt , m_k09 , m_t09 );
 
             
+            // k10 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b10_01 , dxdt.begin() ,
+                       dt * b10_04 , m_k04.begin() ,
+                       dt * b10_05 , m_k05.begin() ,
+                       dt * b10_06 , m_k06.begin() ,
+                       dt * b10_07 , m_k07.begin() ,
+                       dt * b10_08 , m_k08.begin() ,
+                       dt * b10_09 , m_k09.begin() );
+            system( m_xt , m_k10 , m_t10 );
 
 
+            // k11 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b11_01 , dxdt.begin() ,
+                       dt * b11_04 , m_k04.begin()  ,
+                       dt * b11_05 , m_k05.begin() ,
+                       dt * b11_06 , m_k06.begin() ,
+                       dt * b11_07 , m_k07.begin() ,
+                       dt * b11_08 , m_k08.begin() ,
+                       dt * b11_09 , m_k09.begin() ,
+                       dt * b11_10 , m_k10.begin() );
+            system( m_xt , m_k11 , m_t11 );
+
+
+            // k12 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b12_01 , dxdt.begin() ,
+                       dt * b12_06 , m_k06.begin() ,
+                       dt * b12_07 , m_k07.begin() ,
+                       dt * b12_08 , m_k08.begin() ,
+                       dt * b12_09 , m_k09.begin() ,
+                       dt * b12_10 , m_k10.begin() );
+            system( m_xt , m_k12 , m_t12 );
+
+
+            // k13 step
+            scale_sum( m_xt.begin() , m_xt.end() ,
+                       val1 , x.begin() ,
+                       dt * b13_01 , dxdt.begin() ,
+                       dt * b13_04 , m_k04.begin() ,
+                       dt * b13_05 , m_k05.begin() ,
+                       dt * b13_06 , m_k06.begin() ,
+                       dt * b13_07 , m_k07.begin() ,
+                       dt * b13_08 , m_k08.begin() ,
+                       dt * b13_09 , m_k09.begin() ,
+                       dt * b13_10 , m_k10.begin() ,
+                       dt * b13_12 , m_k12.begin() );
+            system( m_xt , m_k13 , m_t13 );
+
+            scale_sum( x.begin() , x.end() ,
+                       val1 , x.begin() ,
+                       dt * c06 , m_k06.begin() ,
+                       dt * c07 , m_k07.begin() ,
+                       dt * c08 , m_k08.begin() ,
+                       dt * c09 , m_k09.begin() ,
+                       dt * c10 , m_k10.begin() ,
+                       dt * c12 , m_k12.begin() ,
+                       dt * c13 , m_k13.begin() );
         }
 
 
         template< class DynamicalSystem >
         void do_step( DynamicalSystem &system ,
-                        container_type &x ,
-                        time_type t ,
-                        time_type dt )
+                      container_type &x ,
+                      time_type t ,
+                      time_type dt )
         {
             m_resizer.adjust_size( x , m_dxdt );
             system( x , m_dxdt , t );
             do_step( system , x , m_dxdt , t , dt );
         }
+
+
+/*
+
+        template< class DynamicalSystem >
+        void do_step( DynamicalSystem &system ,
+                      container_type &x ,
+                      container_type &dxdt ,
+                      time_type t ,
+                      time_type dt ,
+                      container_type &xerr )
+        {
+            const time_type cc01 = static_cast<time_type>( 41.0 / 840.0 );
+            const time_type cc06 = static_cast<time_type>( 34.0 / 105.0 );
+            const time_type cc07 = static_cast<time_type>( 9.0 / 35.0 );
+            const time_type cc08 = static_cast<time_type>( cc08 );
+            const time_type cc09 = static_cast<time_type>( 9.0 / 280.0 );
+            const time_type cc10 = static_cast<time_type>( cc09 );
+            const time_type cc11 = static_cast<time_type>( cc01 );
+
+            xerr = x;
+            do_step( system , xerr , dxdt , t , dt );
+
+            // now, k1-k13 are calculated and stored in m_k01 - m_k13
+            scale_sum( x.begin() , x.end() ,
+                       static_cast<time_type>(1.0) , x.begin(),
+                       dt * cc01 , dxdt ,
+                       dt * cc06 , m_k06 ,
+                       dt * cc07 , m_k07 ,
+                       dt * cc08 , m_k08 ,
+                       dt * cc09 , m_k09 ,
+                       dt * cc10 , m_k10 ,
+                       dt * cc11 , m_k11 );
+
+            increment( xerr.begin() , xerr.end() , x.begin() , static_cast<time_type>(-1.0) );
+        }
+
+        template< class DynamicalSystem >
+        void do_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 );
+            do_step( system , x , m_dxdt , t , dt , xerr );
+        }
+*/
     };
 
 } // namespace odeint

Modified: sandbox/odeint/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/Jamfile	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/Jamfile	2009-11-27 03:04:59 EST (Fri, 27 Nov 2009)
@@ -19,3 +19,4 @@
 exe lorenz_integrator : lorenz_integrator.cpp ;
 exe lorenz_stepper : lorenz_stepper.cpp ;
 exe pendulum_vibrating_pivot : pendulum_vibrating_pivot.cpp ;
+exe dnls_stepper_compare : dnls_stepper_compare.cpp ;

Added: sandbox/odeint/libs/numeric/odeint/examples/dnls_stepper_compare.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/examples/dnls_stepper_compare.cpp	2009-11-27 03:04:59 EST (Fri, 27 Nov 2009)
@@ -0,0 +1,115 @@
+/* Boost numeric/odeint/examples/lorenz_integrator.cpp
+ 
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ compares the steppers rk4 and rk78
+ the system is the dnls, which is complex and Hamiltonian
+
+ 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)
+*/
+
+#include <iostream>
+#include <vector>
+#include <iterator>
+#include <list>
+#include <algorithm>
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace std::tr1;
+using namespace boost::numeric::odeint;
+
+
+const size_t n = 64;
+const double beta = 1.0;
+
+typedef array< complex<double> , n > state_type;
+
+const complex<double> II( 0.0 , -1.0 );
+
+void dnls( state_type &x , state_type &dxdt , double t )
+{
+    dxdt[0] = II * ( beta * norm( x[0] ) * x[0] + x[1] + x[n-1] );
+    for( size_t i=1 ; i<n-1 ; ++i )
+        dxdt[i] = II * ( beta * norm( x[i] ) * x[i] + x[i+1] + x[i-1] );
+    dxdt[n-1] = II * ( beta * norm( x[n-1] ) * x[n-1] + x[0] + x[n-2] );
+}
+
+double norm( const state_type &x )
+{
+    double nn = 0.0;
+    state_type::const_iterator iter = x.begin() ;
+    while( iter != x.end() ) nn += norm( *iter++ );
+    return nn;
+}
+
+double energy( const state_type &x )
+{
+    double e = 0.0 , nn;
+    for( size_t i=0 ; i<n-1 ; ++i )
+    {
+        nn = norm( x[i] );
+        e += 0.5 * beta * nn * nn + 2.0 * ( x[i]*conj(x[i+1]) ).real();
+    }
+    nn = norm( x[n-1] );
+    e += 0.5 * beta * nn * nn + 2.0 * ( x[n-1]*conj(x[0]) ).real();
+    return e;
+}
+
+ostream& operator<<( ostream &out , const state_type &x )
+{
+    state_type::const_iterator iter = x.begin() ;
+    while( iter != x.end() )
+    {
+        const complex<double> &y = *iter++;
+        out << y.real() << tab << y.imag() << tab << norm( y ) << "\n";
+    }
+    return out;
+}
+
+int main( int argc , char **argv )
+{
+    state_type x;
+
+    generate( x.begin() , x.end() , drand48 );
+
+    state_type x1( x ) , x2( x );
+    stepper_rk4< state_type > rk4;
+    stepper_rk78_fehlberg< state_type > rk78;
+
+    double norm0 = norm( x1 ) , energy0 = energy( x1 );
+
+    const size_t olen = 10000 , ilen = 100;
+    const double dt = 0.01;
+    double t = 0.0;
+    cout.precision(14);
+    cout.flags( ios::scientific );
+    for( size_t oi=0 ; oi<olen ; ++oi )
+    {
+        double norm1 = norm( x1 ) , norm2 = norm( x2 );
+        double energy1 = energy( x1 ) , energy2 = energy( x2 );
+
+        cout << t << tab;
+        cout << norm1 << tab << norm2 << tab;
+        cout << energy1 << tab << energy2 << tab;
+        cout << norm1 - norm0 << tab << norm2 - norm0 << tab;
+        cout << energy1 - energy0 << tab << energy2 - energy0 << endl;
+
+        for( size_t ii=0 ; ii<ilen ; ++ii,t+=dt )
+        {
+            rk4.do_step( dnls , x1 , t , dt );
+            rk78.do_step( dnls , x2 , t , dt );
+        }
+    }
+
+
+
+    return 0;
+}