[Boost-commit] svn:boost r57750 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/stuff/test

Subject: [Boost-commit] svn:boost r57750 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/stuff/test
From: mario.mulansky_at_[hidden]
Date: 2009-11-18 11:48:07

Author: mariomulansky
Date: 2009-11-18 11:48:06 EST (Wed, 18 Nov 2009)
New Revision: 57750
URL: http://svn.boost.org/trac/boost/changeset/57750

Log:
added midpoint stepper - not complete yet...
Added:
   sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp   (contents, props changed)
   sandbox/odeint/libs/numeric/odeint/stuff/test/
   sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/boost/numeric/odeint.hpp |     1 +                                       
   1 files changed, 1 insertions(+), 0 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint.hpp	2009-11-18 11:48:06 EST (Wed, 18 Nov 2009)
@@ -21,6 +21,7 @@
 #include <boost/numeric/odeint/stepper_rk5_ck.hpp>
 #include <boost/numeric/odeint/stepper_rk_generic.hpp>
 #include <boost/numeric/odeint/stepper_half_step.hpp>
+#include <boost/numeric/odeint/stepper_midpoint.hpp>
 
 #include <boost/numeric/odeint/stepsize_controller_standard.hpp>
 

Added: sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/stepper_midpoint.hpp	2009-11-18 11:48:06 EST (Wed, 18 Nov 2009)
@@ -0,0 +1,142 @@
+/* Boost odeint/stepper_midpoint.hpp header file
+ 
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ 
+ This file includes the explicit midpoint solver for
+ ordinary differential equations.
+
+ It solves any ODE dx/dt = f(x,t) via
+ x(t+dt) = x(t) + dt*f(x,t)
+
+ 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)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_MIDPOINT_HPP
+#define BOOST_NUMERIC_ODEINT_STEPPER_MIDPOINT_HPP
+
+#include <boost/concept_check.hpp>
+
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/resizer.hpp>
+
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+    template<
+        class Container ,
+        class Time = double ,
+        class Resizer = resizer< Container >
+        >
+    class stepper_midpoint
+    {
+        // provide basic typedefs
+    public:
+
+        typedef Container container_type;
+        typedef Resizer resizer_type;
+        typedef Time time_type;
+        typedef const unsigned short order_type;
+        typedef typename container_type::value_type value_type;
+        typedef typename container_type::iterator iterator;
+
+
+
+
+        // check the concept of the ContainerType
+    private:
+
+        BOOST_CLASS_REQUIRE( container_type ,
+			     boost::numeric::odeint, Container );
+
+        
+        // private memebers
+    private:
+        resizer_type m_resizer;
+
+        container_type m_x0;
+        container_type m_x1;
+        container_type m_x2;
+        container_type m_dxdt;
+
+    public:
+        
+        order_type order() const { return 2; }
+
+        template< class DynamicalSystem >
+        void next_step( 
+                DynamicalSystem &system ,
+                container_type &x ,
+                container_type &dxdt ,
+                time_type t ,
+                time_type dt ,
+                unsigned short n = 2 )
+        {
+            if( n < 2 ) return;
+
+            const time_type h = dt/static_cast<time_type>( n );
+            const time_type h2 = static_cast<time_type>( 2.0 )*h;
+            const time_type t_1 = static_cast<time_type>( 1.0 );
+            const time_type t_05 = static_cast<time_type>( 0.5 );
+            time_type th = t + h;
+
+            m_resizer.adjust_size(x, m_x0);
+            m_resizer.adjust_size(x, m_x1);
+            m_resizer.adjust_size(x, m_x2);
+
+            using namespace detail::it_algebra;
+
+            scale_sum( m_x1.begin(), m_x1.end(),
+                       t_1, x.begin(),
+                       h, dxdt.begin() );
+            system( m_x1, m_x2, th );
+
+            m_x1 = x;
+
+            unsigned short i = 2;
+            while( i != n )
+            {   // general step
+                m_x0 = m_x1;
+                scale_sum( m_x1.begin(), m_x1.end(),
+                           t_1, m_x0.begin(),
+                           h2, m_x2.begin() );
+                th += h;
+                system( m_x1, m_x2, th);
+                i++;
+            }
+            
+            // last step
+            scale_sum( x.begin(), x.end(),
+                       t_05, m_x0.begin(),
+                       t_05, m_x1.begin(),
+                       t_05*h, m_x2.begin() );
+        }
+
+
+        template< class DynamicalSystem >
+        void next_step( 
+                DynamicalSystem &system ,
+                container_type &x ,
+                time_type t ,
+                time_type dt ,
+                unsigned short n = 2 )
+        {
+            m_resizer.adjust_size(x, m_dxdt);
+            system( x, m_dxdt, t );
+            next_step( system , x, m_dxdt, t, dt, n );
+        }
+            
+
+    };
+
+}
+}
+}
+
+#endif

Added: sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/stuff/test/stepper_midpoint.cpp	2009-11-18 11:48:06 EST (Wed, 18 Nov 2009)
@@ -0,0 +1,56 @@
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+#include <cmath>
+
+#include <boost/numeric/odeint.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+typedef std::tr1::array< double , 2 > state_type;
+
+
+const double gam = .15;
+
+void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t)
+{
+    dxdt[0] = x[1];
+    dxdt[1] = -x[0] - gam*x[1];
+}
+
+
+int main( int argc , char **argv )
+{
+    const double dt = 0.01;
+    const size_t olen = 10000;
+    state_type x1 = {{ 1.0, 0.0 }};
+    state_type x2 = {{ 1.0, 0.0 }};
+    state_type x3 = {{ 1.0, 0.0 }};
+    state_type x4 = {{ 1.0, 0.0 }};
+
+    stepper_midpoint< state_type > stepper_mp;
+    stepper_rk4< state_type > stepper_rk4;
+
+    double t = 0.0;
+    for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
+    {
+        stepper_mp.next_step( harmonic_oscillator , x1 , t , dt , 2 );
+        stepper_mp.next_step( harmonic_oscillator , x2 , t , dt , 4 );
+        stepper_mp.next_step( harmonic_oscillator , x3 , t , dt , 8 );
+        stepper_rk4.next_step( harmonic_oscillator , x4 , t , dt );
+        cout<< t << tab << x1[0]*x1[0] + x1[1]*x1[1];
+        cout<< tab << x2[0]*x2[0] + x2[1]*x2[1];
+        cout<< tab << x3[0]*x3[0] + x3[1]*x3[1];
+        cout<< tab << x4[0]*x4[0] + x4[1]*x4[1] << endl;
+    }
+}
+
+/*
+  Compile with
+  g++ -Wall -O3 -I$BOOST_ROOT -I../../../../../ lorenz_stepper_cmp.cpp
+*/