Subject: [Boost-commit] svn:boost r61797 - in sandbox/odeint/branches/karsten: . boost/numeric boost/numeric/odeint boost/numeric/odeint/concepts boost/numeric/odeint/container_traits boost/numeric/odeint/detail boost/numeric/odeint/integrate_functions boost/numeric/odeint/steppers boost/numeric/odeint/steppers/detail libs/numeric/odeint/examples libs/numeric/odeint/test
From: karsten.ahnert_at_[hidden]
Date: 2010-05-05 15:52:10
Author: karsten
Date: 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
New Revision: 61797
URL: http://svn.boost.org/trac/boost/changeset/61797
Log:
structurized my branch
Added:
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits/
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits/container_traits.hpp
- copied unchanged from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits/container_traits_tr1_array.hpp
- copied unchanged from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_tr1_array.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits/container_traits_ublas_matrix.hpp
- copied unchanged from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_ublas_matrix.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_adaptive_stepsize.hpp
- copied, changed from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_adaptive_stepsize.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_constant_stepsize.hpp
- copied, changed from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_constant_stepsize.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/observer.hpp
- copied unchanged from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/observer.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers.hpp (contents, props changed)
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/controlled_stepper_standard.hpp
- copied, changed from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_standard.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/detail/
- copied from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/detail/
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/error_checker_standard.hpp
- copied unchanged from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/error_checker_standard.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/hamiltonian_stepper_euler.hpp
- copied, changed from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_euler.hpp
- copied, changed from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_half_step.hpp
- copied, changed from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_half_step.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/point_type.hpp
- copied unchanged from r61775, /sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp
Removed:
sandbox/odeint/branches/karsten/Jamroot
sandbox/odeint/branches/karsten/boost/numeric/odeint/concepts/
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_blitz_array.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_mtl4_dense2d.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_tr1_array.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_ublas_matrix.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_bs.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_standard.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/detail/
sandbox/odeint/branches/karsten/boost/numeric/odeint/error_checker_standard.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/gram_schmitt.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_rk.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_adaptive_stepsize.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_constant_stepsize.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/observer.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_base.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_half_step.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_midpoint.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk4.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk4_classical.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk5_ck.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk78_fehlberg.hpp
sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk_generic.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/dnls.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/harmonic_oscillator.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_controlled.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_integrator.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_stepper.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp
Text files modified:
sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp | 26 ---
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_adaptive_stepsize.hpp | 7
sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_constant_stepsize.hpp | 2
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/controlled_stepper_standard.hpp | 6
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/hamiltonian_stepper_euler.hpp | 6
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_euler.hpp | 9 -
sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_half_step.hpp | 3
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile | 14 -
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/doc_integrate.cpp | 11
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/solar_system.cpp | 272 ++++++++++++++++-----------------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp | 23 +--
11 files changed, 143 insertions(+), 236 deletions(-)
Deleted: sandbox/odeint/branches/karsten/Jamroot
==============================================================================
--- sandbox/odeint/branches/karsten/Jamroot 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,30 +0,0 @@
-# Copyright 2009 Karsten Ahnert and Mario Mulansky.
-# 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)
-
-import os ;
-import modules ;
-import path ;
-
-project odeint
- : requirements
- <include>$BOOST_ROOT ;
-
-build-project libs/numeric/odeint/examples ;
-build-project libs/numeric/odeint/test ;
-# build-project libs/numeric/odeint/doc ;
-
-
-
-path-constant BOOST_ROOT : [ os.environ BOOST_ROOT ] ;
-
-
-###### The following is copied from another sandbox project #####
-###### to get the quickbook and boostbook working ... #####
-
-local boost-root = [ modules.peek : BOOST_ROOT ] ;
-local explore-header-include = $(top)/../.. ;
-use-project /boost/regex : $(boost-root)/libs/regex/build ;
-
-##################################################################
Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -18,29 +18,7 @@
#include <boost/config.hpp>
#include <boost/numeric/odeint/container_traits.hpp>
-#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
-
-#include <boost/numeric/odeint/stepper_euler.hpp>
-#include <boost/numeric/odeint/stepper_half_step.hpp>
-#include <boost/numeric/odeint/stepper_rk4.hpp>
-#include <boost/numeric/odeint/stepper_rk4_classical.hpp>
-#include <boost/numeric/odeint/stepper_rk5_ck.hpp>
-//#include <boost/numeric/odeint/stepper_rk_generic.hpp>
-#include <boost/numeric/odeint/stepper_midpoint.hpp>
-#include <boost/numeric/odeint/stepper_rk78_fehlberg.hpp>
-
-#include <boost/numeric/odeint/hamiltonian_stepper_euler.hpp>
-#include <boost/numeric/odeint/hamiltonian_stepper_rk.hpp>
-
-
-#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
-#include <boost/numeric/odeint/controlled_stepper_bs.hpp>
-
-#include <boost/numeric/odeint/error_checker_standard.hpp>
-
-#include <boost/numeric/odeint/integrator_adaptive_stepsize.hpp>
-#include <boost/numeric/odeint/integrator_constant_stepsize.hpp>
-
-#include <boost/numeric/odeint/gram_schmitt.hpp>
+#include <boost/numeric/odeint/steppers.hpp>
+#include <boost/numeric/odeint/integrate_functions.hpp>
#endif // BOOST_NUMERIC_ODEINT_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,81 +0,0 @@
-/* Boost odeint/container_traits.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- This file includes container_traits functionality for containers
-
- container_traits
-
- 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_CONTAINER_TRAITS_HPP
-#define BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_HPP
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template< class Container >
- struct container_traits
- {
-
- typedef Container container_type;
-
- typedef typename container_type::value_type value_type;
- typedef typename container_type::iterator iterator;
- typedef typename container_type::const_iterator const_iterator;
-
-
- static void resize( const container_type &x , container_type &dxdt )
- {
- dxdt.resize( x.size() );
- }
-
- static bool same_size(
- const container_type &x1 ,
- const container_type &x2
- )
- {
- return (x1.size() == x2.size());
- }
-
- static void adjust_size(
- const container_type &x1 ,
- container_type &x2
- )
- {
- if( !same_size( x1 , x2 ) ) resize( x1 , x2 );
- }
-
- static iterator begin( container_type &x )
- {
- return x.begin();
- }
-
- static const_iterator begin( const container_type &x )
- {
- return x.begin();
- }
-
- static iterator end( container_type &x )
- {
- return x.end();
- }
-
- static const_iterator end( const container_type &x )
- {
- return x.end();
- }
- };
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_blitz_array.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_blitz_array.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,85 +0,0 @@
-/*
- boost header: numeric/odeint/blitz_container_traits.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- Container traits for blitz::Array.
-
- 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_CONTAINER_TRAITS_BLITZ_ARRAY_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_BLITZ_ARRAY_HPP_INCLUDED
-
-#include <cstdlib>
-#include "container_traits.hpp"
-#include <blitz/array.h>
-
-#include<iostream>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template< typename T , int n >
- struct container_traits< blitz::Array< T , n > >
- {
-
- typedef blitz::Array< T , n > container_type;
- typedef T value_type;
- typedef typename container_type::iterator iterator;
- typedef typename container_type::const_iterator const_iterator;
-
-
-
- static void resize( const container_type &x , container_type &dxdt )
- {
- //dxdt.resize( x.shape() );
- dxdt.resizeAndPreserve( x.shape() );
- }
-
- static bool same_size( const container_type &x1 , const container_type &x2 )
- {
- for( int d=0; d<x1.dimensions(); d++ )
- if( x1.extent(d) != x2.extent(d) )
- return false;
- return true;
- }
-
- static void adjust_size( const container_type &x1 , container_type &x2 )
- {
- //if( !same_size( x1 , x2 ) ) resize( x1 , x2 );
- x2.resizeAndPreserve( x1.shape() );
- }
-
- static iterator begin( container_type &x )
- {
- return x.begin();
- }
-
- static const_iterator begin( const container_type &x )
- {
- return x.begin();
- }
-
- static iterator end( container_type &x )
- {
- return x.end();
- }
-
- static const_iterator end( const container_type &x )
- {
- return x.end();
- }
- };
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif //BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_BLITZ_ARRAY_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_mtl4_dense2d.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_mtl4_dense2d.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,170 +0,0 @@
-/*
- boost header: numeric/odeint/mtl4_dense2d_container_traits.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- Container traits for mtl4::dense2D matrices
-
- 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_CONTAINER_TRAITS_MTL4_DENSE2D_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_MTL4_DENSE2D_HPP_INCLUDED
-
-#include <boost/numeric/odeint/container_traits.hpp>
-#include <boost/numeric/mtl/mtl.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template <typename Value, typename Parameters >
- struct container_traits< mtl::dense2D< Value , Parameters > >
- {
-
- typedef mtl::matrix::dense2D< Value , Parameters > container_type;
-
- typedef typename container_type::value_type value_type;
-
- typedef typename mtl::traits::range_generator<
- mtl::tag::all, container_type >::type cursor;
- typedef typename mtl::traits::value< container_type >::type value_map;
- typedef typename mtl::traits::const_value<
- container_type >::type const_value_map;
- typedef mtl::utilities::iterator_adaptor<
- value_map, cursor, value_type > iterator;
- typedef mtl::utilities::iterator_adaptor<
- const_value_map, cursor, value_type > const_iterator;
-
- static void resize( const container_type &x , container_type &dxdt )
- {
- dxdt = container_type( x.num_rows(), x.num_cols() );
- }
-
- static bool same_size(
- const container_type &x1 ,
- const container_type &x2
- )
- {
- return ( (x1.num_rows() == x2.num_rows()) && (x1.num_cols() == x2.num_cols()));
- }
-
- static void adjust_size(
- const container_type &x1 ,
- container_type &x2
- )
- {
- if( !same_size( x1 , x2 ) ) resize( x1 , x2 );
- }
-
- static iterator begin( container_type &x )
- {
- cursor c = mtl::begin< mtl::tag::all >(x);
- value_map v(x);
- return iterator(v, c);
- }
-
- static const_iterator begin( const container_type &x )
- {
- cursor c = mtl::begin< mtl::tag::all >(x);
- const_value_map v(x);
- return const_iterator(v, c);
- }
-
- static iterator end( container_type &x )
- {
- cursor c = mtl::end< mtl::tag::all >(x);
- value_map v(x);
- return iterator(v, c);
- }
-
- static const_iterator end( const container_type &x )
- {
- cursor c = mtl::end< mtl::tag::all >(x);
- const_value_map v(x);
- return const_iterator(v, c);
- }
-
- };
-
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-
-/* Template Specialization to provide += operator for iterator return type */
-namespace mtl { namespace utilities { namespace detail {
-
-template <typename Cursor, typename ValueType>
-struct iterator_proxy<typename mtl::traits::value< mtl::matrix::dense2D< ValueType > >::type, Cursor, ValueType>
-{
- typedef iterator_proxy self;
- typedef typename mtl::traits::value< mtl::matrix::dense2D< ValueType > >::type property_map;
-
- iterator_proxy(property_map map, Cursor cursor)
- : map(map), cursor(cursor) {}
-
- operator ValueType() const
- {
- return static_cast<ValueType>(map(*cursor));
- }
-
- self& operator=(ValueType const& value)
- {
- map(*cursor, value);
- return *this;
- }
-
- self& operator+=(ValueType const& value)
- {
- map( *cursor, value + static_cast<ValueType>(map(*cursor)) );
- return *this;
- }
-
- property_map map;
- Cursor cursor;
-};
-
-}}} // namespace mtl::utilities::detail
-
-
-// define iterator traits for dense2D iterators
-namespace std {
-
-template < typename Cursor , typename Value , typename Parameters>
-struct iterator_traits<
- mtl::utilities::iterator_adaptor< mtl::detail::direct_value< mtl::dense2D< Value , Parameters > >, Cursor, Value >
->
-{
- typedef Value value_type;
- typedef int difference_type; // ?
- typedef Value* pointer;
- typedef Value& reference;
- typedef std::random_access_iterator_tag iterator_category;
-};
-
-
-template < typename Cursor , typename Value , typename Parameters>
-struct iterator_traits<
- mtl::utilities::iterator_adaptor< mtl::detail::direct_const_value< mtl::dense2D< Value , Parameters > >, Cursor, Value >
->
-{
- typedef Value value_type;
- typedef int difference_type; // ?
- typedef const Value* pointer;
- typedef const Value& reference;
- typedef std::random_access_iterator_tag iterator_category;
-};
-
-
-}
-
-
-#endif //BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_MTL4_DENSE2D_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_tr1_array.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_tr1_array.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,83 +0,0 @@
-/*
- boost header: numeric/odeint/container_traits_tr1_array.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_CONTAINER_TRAITS_TR1_ARRAY_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_TR1_ARRAY_HPP_INCLUDED
-
-#include <tr1/array>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- // Template Specialization for fixed size array - no resizing can happen
- template< class T , size_t N >
- struct container_traits< std::tr1::array< T , N > >
- {
- public:
-
- typedef std::tr1::array< T , N > container_type;
- typedef typename container_type::value_type value_type;
- typedef typename container_type::iterator iterator;
- typedef typename container_type::const_iterator const_iterator;
-
-
- static void resize( const container_type &x , container_type &dxdt )
- {
- throw; // should never be called
- }
-
- static const bool same_size(
- const container_type &x1 ,
- const container_type &x2
- )
- {
- return true; // if this was false, the code wouldn't compile
- }
-
- static void adjust_size(
- const container_type &x1 ,
- container_type &x2
- )
- {
- if( !same_size( x1 , x2 ) ) throw;
- }
-
- static iterator begin( container_type &x )
- {
- return x.begin();
- }
-
- static const_iterator begin( const container_type &x )
- {
- return x.begin();
- }
-
- static iterator end( container_type &x )
- {
- return x.end();
- }
-
- static const_iterator end( const container_type &x )
- {
- return x.end();
- }
- };
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif //BOOST_NUMERIC_ODEINT_CONTAINER_TRAITS_TR1_ARRAY_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_ublas_matrix.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/container_traits_ublas_matrix.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,82 +0,0 @@
-/*
- boost header: numeric/odeint/ublas_matrix_container_traits.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_UBLAS_MATRIX_CONTAINER_TRAITS_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_UBLAS_MATRIX_CONTAINER_TRAITS_HPP_INCLUDED
-
-#include "container_traits.hpp"
-#include <boost/numeric/ublas/matrix.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- template<class T, class L, class A>
- struct container_traits< boost::numeric::ublas::matrix< T , L , A > >
- {
-
- typedef boost::numeric::ublas::matrix< T , L , A > container_type;
- typedef typename container_type::value_type value_type;
- typedef typename container_type::array_type::iterator iterator;
- typedef typename container_type::array_type::const_iterator const_iterator;
-
-
-
- static void resize( const container_type &x , container_type &dxdt )
- {
- dxdt.resize( x.size1() , x.size2() );
- }
-
- static bool same_size(
- const container_type &x1 ,
- const container_type &x2
- )
- {
- return ( ( x1.size1() == x2.size1() ) && ( x1.size2() == x2.size2() ) );
- }
-
- static void adjust_size(
- const container_type &x1 ,
- container_type &x2
- )
- {
- if( !same_size( x1 , x2 ) ) resize( x1 , x2 );
- }
-
- static iterator begin( container_type &x )
- {
- return x.data().begin();
- }
-
- static const_iterator begin( const container_type &x )
- {
- return x.data().begin();
- }
-
- static iterator end( container_type &x )
- {
- return x.data().end();
- }
-
- static const_iterator end( const container_type &x )
- {
- return x.data().end();
- }
- };
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif //BOOST_NUMERIC_ODEINT_UBLAS_MATRIX_CONTAINER_TRAITS_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_bs.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_bs.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,369 +0,0 @@
-/* Boost odeint/controlled_stepper_bs.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- This file includes the explicit Burlisch Stoer
- 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_CONTROLLED_STEPPER_BS_HPP
-#define BOOST_NUMERIC_ODEINT_CONTROLLED_STEPPER_BS_HPP
-
-#include <limits>
-#include <vector>
-
-#include <boost/numeric/odeint/stepper_midpoint.hpp>
-#include <boost/numeric/odeint/error_checker_standard.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class controlled_stepper_bs
- {
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
-
-
- // private memebers
- private:
-
- stepper_midpoint< container_type, time_type, traits_type > m_stepper_mp;
- error_checker_standard< container_type, time_type , traits_type > m_error_checker;
-
- const unsigned short m_k_max;
-
- const time_type m_safety1;
- const time_type m_safety2;
- const time_type m_max_dt_factor;
- const time_type m_min_step_scale;
- const time_type m_max_step_scale;
-
- bool m_continuous_calls;
- bool m_decreased_step_during_last_call;
-
- time_type m_dt_last;
- time_type m_t_last;
- time_type m_current_eps;
-
- unsigned short m_current_k_max;
- unsigned short m_current_k_opt;
-
- container_type m_x0;
- container_type m_xerr;
- container_type m_x_mp;
- container_type m_x_scale;
- container_type m_dxdt;
-
- typedef std::vector< time_type > value_vector;
- typedef std::vector< std::vector< time_type > > value_matrix;
- typedef std::vector< unsigned short > us_vector;
-
- value_vector m_error; // errors of repeated midpoint steps and extrapolations
- value_vector m_a; // stores the work (number of f calls) required for the orders
- value_matrix m_alpha; // stores convergence factor for stepsize adjustment
- us_vector m_interval_sequence;
-
- value_vector m_times;
- std::vector< container_type > m_d;
- container_type m_c;
-
- // public functions
- public:
-
- // constructor
- controlled_stepper_bs(
- time_type abs_err, time_type rel_err,
- time_type factor_x, time_type factor_dxdt )
- : m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
- m_k_max(8),
- m_safety1(0.25), m_safety2(0.7),
- m_max_dt_factor( 0.1 ), m_min_step_scale(5E-5), m_max_step_scale(0.7),
- m_continuous_calls(false), m_decreased_step_during_last_call( false ),
- m_dt_last( 1.0E30),
- m_current_eps( -1.0 )
- {
- m_error.resize(m_k_max);
- m_a.resize(m_k_max+1);
- m_alpha.resize(m_k_max); // k_max * k_max matrix
- typename value_matrix::iterator it = m_alpha.begin();
- while( it != m_alpha.end() )
- (*it++).resize(m_k_max);
- m_interval_sequence.resize(m_k_max+1);
- for( unsigned short i = 1; i <= m_k_max+1; i++ )
- m_interval_sequence[i-1] = (2*i);
-
- m_times.resize(m_k_max);
- m_d.resize(m_k_max);
- }
-
-
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size(x, m_xerr);
- traits_type::adjust_size(x, m_x_mp);
- traits_type::adjust_size(x, m_x_scale);
- traits_type::adjust_size(x, m_dxdt);
- m_stepper_mp.adjust_size( x );
- }
-
-
- template< class DynamicalSystem >
- controlled_step_result try_step(
- DynamicalSystem &system ,
- container_type &x ,
- container_type &dxdt ,
- time_type &t ,
- time_type &dt )
- {
-
- // get error scale
- m_error_checker.fill_scale(x, dxdt, dt, m_x_scale);
-
- //std::clog << " x scale: " << m_x_scale[0] << '\t' << m_x_scale[1] << std::endl;
-
- m_x0 = x; // save starting state
- time_type max_eps = m_error_checker.get_epsilon();
- if( m_current_eps != max_eps )
- { // error changed -> recalculate tableau
- initialize_step_adjust_tableau( max_eps );
- m_current_eps = max_eps;
- m_dt_last = m_t_last = std::numeric_limits< value_type >::max(); // unrealistic
- }
- // if t and dt are exactly the parameters from last step we are called continuously
- bool continuous_call = ((t == m_t_last) || (dt == m_dt_last ));
- if( !continuous_call )
- m_current_k_opt = m_current_k_max;
-
- bool converged = false;
- value_type step_scale = 1.0;
- unsigned short k_conv = 0;
-
- for( unsigned short k=0; k<=m_current_k_max; k++ )
- { // loop through interval numbers
- unsigned short step_number = m_interval_sequence[k];
- //out-of-place midpoint step
- m_stepper_mp.set_step_number(step_number);
- m_stepper_mp.midpoint_step(system, m_x0, dxdt, t, dt, m_x_mp);
- //std::clog << "x_mp: " << k << '\t' << m_x_mp[0] << '\t' << m_x_mp[1] << std::endl;
- time_type t_est = (dt/step_number)*(dt/step_number);
- extrapolate(k, t_est, m_x_mp, x, m_xerr);
- //std::clog << "Error: " << k << '\t' << m_xerr[0] << '\t' << m_xerr[1] << std::endl;
- if( k != 0 )
- {
- time_type max_err = m_error_checker.get_max_error_ratio(m_xerr, m_x_scale);
- m_error[k-1] = std::pow( max_err/m_safety1, 1.0/(2*k+1) );
- if( (k >= m_current_k_opt-1) || !continuous_call )
- { //we're in the order window where convergence is expected
- if( max_err < 1.0 )
- {
- k_conv = k;
- converged = true;
- //std::clog << "Converged at: " << k << '\t' << max_err << std::endl;
- break;
- } else {
- converged = false;
- if( (k == m_current_k_max) || (k == m_current_k_opt+1) )
- {
- step_scale = m_safety2/m_error[k-1];
- break;
- } else if( (k == m_current_k_opt) &&
- (m_alpha[m_current_k_opt-1][m_current_k_opt] < m_error[k-1] ) )
- {
- step_scale = static_cast<value_type>(1.0)/m_error[k-1];
- break;
- } else if( (m_current_k_opt == m_current_k_max) &&
- (m_alpha[k-1][m_current_k_max-1] < m_error[k-1]) )
- {
- step_scale = m_alpha[k-1][m_current_k_max-1]*m_safety2/m_error[k-1];
- //std::clog << " Case 3 " << m_alpha[k-1][m_current_k_max-1] << '\t' << m_error[k-1] << '\t' << step_scale << std::endl;
- break;
- } else if( m_alpha[k-1][m_current_k_opt] < m_error[k-1] )
- {
- step_scale = m_alpha[k-1][m_current_k_opt-1]/m_error[k-1];
- break;
- }
- }
- }
- }
- }
- if( !converged ) { // dt was too large - no convergence up to order k_max
-
- // step_scale is always > 1
- step_scale = std::max(step_scale, m_min_step_scale); // at least m_min ...
- step_scale = std::min(step_scale, m_max_step_scale); // at most m_max ...
- dt *= step_scale;
- m_dt_last = dt;
- m_t_last = t;
- m_decreased_step_during_last_call = true;
- x = m_x0; // copy the state back
- return step_size_decreased;
-
- } else { //converged
-
- t += dt; // step accomplished
-
- time_type work_min = std::numeric_limits< time_type >::max();
- for( unsigned short k=0; k < k_conv ; k++ )
- { // compute optimal convergence order and corresponding stepsize
- time_type factor = std::max(m_error[k], m_max_dt_factor);
- time_type work = factor * m_a[k+1];
- if( work < work_min )
- {
- step_scale = factor;
- work_min = work;
- m_current_k_opt = k+1;
- }
- }
- m_dt_last = dt/step_scale;
-
- //std::clog << "Internal: " << dt << '\t' << k_conv-1 << '\t' << m_current_k_opt << '\t' << m_current_k_max << '\t' << m_decreased_step_during_last_call << std::endl;
-
- if( (m_current_k_opt >= k_conv) &&
- (m_current_k_opt != m_current_k_max) &&
- !m_decreased_step_during_last_call )
- { // check possible order increase, if step was not decreased before
- time_type factor = std::max(
- step_scale/m_alpha[m_current_k_opt-1][m_current_k_opt],
- m_max_dt_factor );
- if( m_a[m_current_k_opt+1]*factor <= work_min )
- {
- m_dt_last = dt/factor;
- m_current_k_opt++;
- }
- }
- dt = m_dt_last;
- m_t_last = t;
- m_decreased_step_during_last_call = false;
- //std::clog << "Internal: " << dt << '\t' << m_current_k_opt << std::endl;
- return step_size_increased;
- }
- }
-
- template< class System >
- controlled_step_result try_step(
- System &system,
- container_type &x,
- time_type &t,
- time_type &dt )
- {
- system(x, m_dxdt, t);
- return try_step(system, x, m_dxdt, t, dt );
- }
-
-
- private: // private functions
-
- void initialize_step_adjust_tableau( time_type eps )
- {
- m_a[0] = m_interval_sequence[0]+1;
- for( unsigned short k=0; k<m_k_max; k++ )
- m_a[k+1] = m_a[k] + m_interval_sequence[k+1];
- for( unsigned short i=1; i<m_k_max; i++ )
- {
- for( unsigned short k=0; k<i; k++ )
- {
- m_alpha[k][i] = std::pow(
- m_safety1*eps,
- (m_a[k+1]-m_a[i+1])/
- ((m_a[i+1]-m_a[0]+1.0)*(2*k+3))
- );
- }
- }
- m_current_k_opt = m_k_max-1;
- for( unsigned short k=1; k<m_k_max-1; k++ )
- {
- if( m_a[k+1] > m_a[k]*m_alpha[k-1][k] )
- {
- m_current_k_opt = k;
- break;
- }
- }
- m_current_k_max = m_current_k_opt;
- }
-
- void extrapolate(
- unsigned short k_est,
- time_type t_est,
- container_type &x_est,
- container_type &x,
- container_type &x_err )
- {
- //traits_type::adjust_size(x, m_c);
- //std::vector< container_type > m_d_iter = m_d.begin();
- //while( m_d_iter != m_d.end() )
- // traits_type::adjust_size(x, (*m_d_iter++));
-
- m_times[k_est] = t_est;
- x_err = x = x_est;
-
- const iterator x_end = traits_type::end(x);
-
- if( k_est == 0 )
- {
- m_d[0] = x_est;
- }
- else
- {
- m_c = x_est;
- for( unsigned short k=0; k<k_est; k++ )
- {
- value_type delta = static_cast<value_type>(1.0) /
- (m_times[k_est-k-1]-static_cast<value_type>(t_est));
- value_type val1 = static_cast<value_type>(t_est)*delta;
- value_type val2 = m_times[k_est-k-1]*delta;
-
- //std::clog << " values: " << delta << '\t' << val1 << '\t' << val2 << std::endl;
-
- iterator x_iter = traits_type::begin(x);
- iterator x_err_iter = traits_type::begin(x_err);
- iterator d_k_iter = m_d[k].begin();
- iterator c_iter = m_c.begin();
- while( x_iter != x_end )
- {
- //std::clog << " extrapolate: " << '\t' << *x_iter << '\t' << *x_err_iter << '\t' << *d_k_iter << std::endl;
- value_type q = *d_k_iter;
- *d_k_iter++ = *x_err_iter;
- delta = *c_iter - q;
- *x_err_iter = val1 * delta;
- *c_iter++ = val2 * delta;
- *x_iter++ += *x_err_iter++;
- }
- //std::clog << " extrapolate: " << '\t' << x[1] << '\t' << x_err[1] << '\t' << m_d[k][1] << std::endl;
- m_d[k_est] = x_err;
- }
- }
- }
-
- };
-}
-}
-}
-
-#endif
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_standard.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_standard.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,212 +0,0 @@
-/* Boost odeint/controlled_stepper_standard.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- This file includes the standard controlled stepper that should be
- used with runge kutta routines.
-
- 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_CONTROLLED_STEPPER_STANDARD_HPP
-#define BOOST_NUMERIC_ODEINT_CONTROLLED_STEPPER_STANDARD_HPP
-
-#include <cmath> // for pow( ) and abs()
-#include <algorithm>
-#include <complex>
-
-#include <boost/numeric/odeint/error_checker_standard.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-
-
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- typedef enum
- {
- success ,
- step_size_decreased ,
- step_size_increased
- } controlled_step_result;
-
-
- /*
- The initial state is given in x.
-
- The stepsize is adjust such that the following maximal relative error is
- small enough for each step:
- R = max( x_err_n / [eps_abs + eps_rel*( a_x * |x_n| + a_dxdt * |dxdt_n| )] )
- where the max refers to the componentwise maximum the expression.
-
- if R > 1.1 the stepsize is decreased:
- dt = dt*S*R^(-1/(q-1))
-
- if R < 0.5 the stepsize is increased:
- dt = dt*S*R^(-1/q))
-
- q is the order of the error term provided by the stepper.order_error() function
- (e.g. 2 for simple euler and 5 for rk5_ck) and S is a safety factor set to
- S = 0.9. See e.g. Numerical Recipes for more details on this strategy.
-
- To avoid extensive chages in dt, the decrease factor is limited to 0.2 and
- the increase factor to 5.0.
- */
-
- template< class ErrorStepper >
- class controlled_stepper_standard
- {
-
- public:
-
- // forward types from ErrorStepper
- typedef ErrorStepper stepper_type;
- typedef typename stepper_type::order_type order_type;
- typedef typename stepper_type::container_type container_type;
- typedef typename stepper_type::time_type time_type;
- typedef typename stepper_type::traits_type traits_type;
- typedef typename stepper_type::value_type value_type;
-// typedef typename stepper_type::iterator iterator;
-// typedef typename stepper_type::const_iterator const_iterator;
-
- typedef error_checker_standard<
- container_type, time_type , traits_type
- > error_checker_type;
-
-
- // private members
- private:
-
- stepper_type m_stepper;
- error_checker_type m_error_checker;
-
- time_type m_eps_abs;
- time_type m_eps_rel;
- time_type m_a_x;
- time_type m_a_dxdt;
-
- container_type m_dxdt;
- container_type m_x_tmp;
- container_type m_x_err;
- container_type m_x_scale;
-
-
- // public functions
- public:
-
- order_type order_error_step() { return m_stepper.order_error_step(); }
-
- controlled_stepper_standard(
- time_type abs_err, time_type rel_err,
- time_type factor_x, time_type factor_dxdt )
- : m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
- m_eps_abs(abs_err),
- m_eps_rel(rel_err),
- m_a_x(factor_x),
- m_a_dxdt(factor_dxdt)
- {
- }
-
- controlled_stepper_standard(
- const container_type &x ,
- time_type abs_err, time_type rel_err,
- time_type factor_x, time_type factor_dxdt )
- : m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
- m_eps_abs(abs_err),
- m_eps_rel(rel_err),
- m_a_x(factor_x),
- m_a_dxdt(factor_dxdt)
- {
- adjust_size( x );
- }
-
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size( x , m_x_err );
- traits_type::adjust_size( x , m_x_scale );
- traits_type::adjust_size( x , m_dxdt );
- m_stepper.adjust_size( x );
- }
-
-
-
-
- /* Tries a controlled step with the given stepsize dt. If dt is too large,
- x remains unchanged, an appropriate stepsize is assigned to dt and
- step_size_decreased is returned. If dt is too small, the small step is
- accomplished, a larger stepsize is assigned to dt and step_size_increased
- is returned. If dt is appropriate, the step is accomplished and success
- is returned.
- */
- template< class DynamicalSystem >
- controlled_step_result try_step(
- DynamicalSystem &system,
- container_type &x,
- const container_type &dxdt,
- time_type &t,
- time_type &dt )
- {
- m_error_checker.fill_scale( x , dxdt , dt , m_x_scale );
-
- m_x_tmp = x;
- m_stepper.do_step( system , x , dxdt, t , dt , m_x_err );
-
- time_type max_rel_err = m_error_checker.get_max_error_ratio(m_x_err, m_x_scale);
-
- if( max_rel_err > 1.1 )
- {
- // error too large - decrease dt
- // limit scaling factor to 0.2
- dt *= std::max( 0.9 * pow(max_rel_err , -1.0/(m_stepper.order_error()-1.0)),
- 0.2 );
-
- // reset state
- x = m_x_tmp;
- return step_size_decreased;
- }
- else
- {
- if( max_rel_err < 0.5 )
- {
- //error too small - increase dt and keep the evolution
- t += dt;
- // limit scaling factor to 5.0
- dt *= std::min( 0.9*pow(max_rel_err , -1.0/m_stepper.order_error_step()), 5.0 );
- return step_size_increased;
- }
- else
- {
- t += dt;
- return success;
- }
- }
- }
-
- template< class DynamicalSystem >
- controlled_step_result try_step(
- DynamicalSystem &system,
- container_type &x,
- time_type &t,
- time_type &dt )
- {
- system( x , m_dxdt , t );
- return try_step( system , x , m_dxdt , t , dt );
- }
-
- };
-
-
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_STEPSIZE_CONTROLLER_STANDARD_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/error_checker_standard.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/error_checker_standard.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,95 +0,0 @@
-/* Boost odeint/error_checker_standard.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- This file includes the standard error checker to be used with
- controlled steppers. It's purpose is to provide a method
- that calculates the error ration of a given error-estimation
- with respect to some user defined tolerance.
-
- 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_ERROR_CHECKER_STANDARD_HPP
-#define BOOST_NUMERIC_ODEINT_ERROR_CHECKER_STANDARD_HPP
-
-#include <cmath>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template< class Container,
- class Time,
- class Traits = container_traits< Container > >
- class error_checker_standard
- {
-
-
- public:
-
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
-
-
- private:
-
- time_type m_eps_abs;
- time_type m_eps_rel;
- time_type m_a_x;
- time_type m_a_dxdt;
-
- public:
-
- // constructor
- error_checker_standard(
- time_type abs_err, time_type rel_err,
- time_type factor_x, time_type factor_dxdt )
- : m_eps_abs( abs_err ), m_eps_rel( rel_err ),
- m_a_x( factor_x ), m_a_dxdt( factor_dxdt )
- {
- }
-
-
- void fill_scale(
- const container_type &x,
- const container_type &dxdt,
- time_type dt,
- container_type &scale ) const
- {
- detail::it_algebra::weighted_scale( traits_type::begin(scale),
- traits_type::end(scale),
- traits_type::begin(x),
- traits_type::begin(dxdt),
- m_eps_abs,
- m_eps_rel,
- m_a_x,
- m_a_x*dt );
- }
-
- time_type get_max_error_ratio(
- const container_type &x_err,
- const container_type &scale) const
- {
- return detail::it_algebra::max_ratio( traits_type::begin(x_err),
- traits_type::end(x_err),
- traits_type::begin(scale),
- static_cast<time_type>(0.0) );
- }
-
- const time_type get_epsilon() { return std::max(m_eps_rel, m_eps_abs); }
- };
-
-}
-}
-}
-
-#endif
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/gram_schmitt.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/gram_schmitt.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,84 +0,0 @@
-/*
- boost header: numeric/odeint/gram_schmitt.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_GRAM_SCHMITT_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
-
-#include <iterator>
-#include <algorithm>
-#include <numeric>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template< class Iterator , class T >
- void normalize( Iterator first , Iterator last , T norm )
- {
- while( first != last ) *first++ /= norm;
- }
-
- template< class Iterator , class T >
- void substract_vector( Iterator first1 , Iterator last1 ,
- Iterator first2 , T val )
- {
- while( first1 != last1 )
- *first1++ -= val * *first2++;
- }
-
- template< class StateType , class LyapType >
- void gram_schmitt( StateType &x , LyapType &lyap ,
- size_t n , size_t num_of_lyap )
- {
- if( !num_of_lyap ) return;
- if( ptrdiff_t( ( num_of_lyap + 1 ) * n ) != std::distance( x.begin() , x.end() ) )
- throw std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." );
-
- typedef typename StateType::value_type value_type;
- typedef typename StateType::iterator iterator;
-
- value_type norm[num_of_lyap] , tmp[num_of_lyap];
- iterator first = x.begin() + n;
- iterator beg1 = first , end1 = first + n ;
-
- std::fill( norm , norm+num_of_lyap , 0.0 );
-
- // normalize first vector
- norm[0] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
- normalize( beg1 , end1 , norm[0] );
-
- beg1 += n ;
- end1 += n;
-
- for( size_t j=1 ; j<num_of_lyap ; ++j , beg1+=n , end1+=n )
- {
- for( size_t k=0 ; k<j ; ++k )
- tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 );
-
- for( size_t k=0 ; k<j ; ++k )
- substract_vector( beg1 , end1 , first + k*n , tmp[k] );
-
- // nromalize j-th vector
- norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
- normalize( beg1 , end1 , norm[j] );
- }
-
- for( size_t j=0 ; j<num_of_lyap ; j++ )
- lyap[j] += log( norm[j] );
- }
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-#endif //BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,88 +0,0 @@
-/*
- boost header: numeric/odeint/hamiltonian_stepper_euler.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_HAMILTONIAN_STEPPER_EULER_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_EULER_HPP_INCLUDED
-
-#include <stdexcept>
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class hamiltonian_stepper_euler
- {
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
- private:
-
- container_type m_dpdt;
-
-
-
- public:
-
- template< class CoordinateFunction >
- void do_step( CoordinateFunction &qfunc ,
- container_type &q ,
- container_type &p ,
- time_type dt )
- {
- if( !traits_type::same_size( q , p ) )
- {
- std::string msg( "hamiltonian_stepper_euler::do_step(): " );
- msg += "q and p have different sizes";
- throw std::invalid_argument( msg );
- }
-
- traits_type::adjust_size( p , m_dpdt );
-
- detail::it_algebra::increment( traits_type::begin(q) ,
- traits_type::end(q) ,
- traits_type::begin(p) ,
- dt );
- qfunc( q , m_dpdt );
- detail::it_algebra::increment( traits_type::begin(p) ,
- traits_type::end(p) ,
- traits_type::begin(m_dpdt) ,
- dt );
- }
-
- };
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-
-
-
-#endif //BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_EULER_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_rk.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,126 +0,0 @@
-/*
- boost header: numeric/odeint/hamiltonian_stepper_euler.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_HAMILTONIAN_STEPPER_RK_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_RK_HPP_INCLUDED
-
-#include <stdexcept>
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class hamiltonian_stepper_rk
- {
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
- private:
-
- container_type m_dpdt;
-
-/*
- rk_a[0]=0.40518861839525227722;
- rk_a[1]=-0.28714404081652408900;
- rk_a[2]=0.5-(rk_a[0]+rk_a[1]);
- rk_a[3]=rk_a[2];
- rk_a[4]=rk_a[1];
- rk_a[5]=rk_a[0];
-
- rk_b[0]=-3.0/73.0;
- rk_b[1]=17.0/59.0;
- rk_b[2]=1.0-2.0*(rk_b[0]+rk_b[1]);
- rk_b[3]=rk_b[1];
- rk_b[4]=rk_b[0];
- rk_b[5]=0.0;
-*/
-
- public:
-
- template< class CoordinateFunction >
- void do_step( CoordinateFunction &qfunc ,
- container_type &q ,
- container_type &p ,
- time_type dt )
- {
- const size_t order = 6;
-
- const time_type rk_a[order] = {
- static_cast<time_type>( 0.40518861839525227722 ) ,
- static_cast<time_type>( -0.28714404081652408900 ) ,
- static_cast<time_type>( 0.3819554224212718118 ) ,
- static_cast<time_type>( 0.3819554224212718118 ) ,
- static_cast<time_type>( -0.28714404081652408900 ) ,
- static_cast<time_type>( 0.40518861839525227722 )
- };
- const time_type rk_b[order] = {
- static_cast<time_type>( -3.0/73.0 ) ,
- static_cast<time_type>( 17.0/59.0 ) ,
- static_cast<time_type>( 0.50592059438123984212 ) ,
- static_cast<time_type>( 17.0/59.0 ) ,
- static_cast<time_type>( -3.0/73.0 ) ,
- static_cast<time_type>( 0.0 )
- };
-
-
-
- if( !traits_type::same_size( q , p ) )
- {
- std::string msg( "hamiltonian_stepper_rk::do_step(): " );
- msg += "q and p have different sizes";
- throw std::invalid_argument( msg );
- }
-
- traits_type::adjust_size( p , m_dpdt );
-
- for( size_t l=0 ; l<order ; ++l )
- {
- detail::it_algebra::increment( traits_type::begin(q) ,
- traits_type::end(q) ,
- traits_type::begin(p) ,
- rk_a[l]*dt );
- qfunc( q , m_dpdt );
- detail::it_algebra::increment( traits_type::begin(p) ,
- traits_type::end(p) ,
- traits_type::begin(m_dpdt) ,
- rk_b[l]*dt );
- }
- }
-
- };
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-
-
-
-#endif //BOOST_NUMERIC_ODEINT_HAMILTONIAN_STEPPER_RK_HPP_INCLUDED
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -0,0 +1,19 @@
+/*
+ boost header: boost/numeric/odeint/integrate_functions.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ 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_BOOST_NUMERIC_ODEINT_INTEGRATE_FUNCTIONS_HPP_INCLUDED
+#define BOOST_BOOST_NUMERIC_ODEINT_INTEGRATE_FUNCTIONS_HPP_INCLUDED
+
+#include <boost/numeric/odeint/integrate_functions/integrator_constant_stepsize.hpp>
+#include <boost/numeric/odeint/integrate_functions/integrator_adaptive_stepsize.hpp>
+
+#endif //BOOST_BOOST_NUMERIC_ODEINT_INTEGRATE_FUNCTIONS_HPP_INCLUDED
Copied: sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_adaptive_stepsize.hpp (from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_adaptive_stepsize.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_adaptive_stepsize.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_adaptive_stepsize.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -13,8 +13,8 @@
#ifndef BOOST_NUMERIC_ODEINT_INTEGRATOR_ADAPTIVE_STEPSIZE_HPP
#define BOOST_NUMERIC_ODEINT_INTEGRATOR_ADAPTIVE_STEPSIZE_HPP
-#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
-#include <boost/numeric/odeint/observer.hpp>
+#include <boost/numeric/odeint/steppers/controlled_stepper_standard.hpp>
+#include <boost/numeric/odeint/integrate_functions/observer.hpp>
#include <vector>
#include <limits>
@@ -154,7 +154,8 @@
)
{
// we use cash karp stepper as base stepper
- typedef stepper_rk5_ck< ContainerType , T > stepper_type;
+// typedef stepper_rk5_ck< ContainerType , T > stepper_type;
+ typedef stepper_half_step< stepper_euler< ContainerType , T > > stepper_type;
// we use the standard controller for this adaptive integrator
controlled_stepper_standard< stepper_type >
Copied: sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_constant_stepsize.hpp (from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_constant_stepsize.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_constant_stepsize.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/integrate_functions/integrator_constant_stepsize.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -13,7 +13,7 @@
#ifndef BOOST_NUMERIC_ODEINT_INTEGRATOR_CONSTANT_STEPSIZE_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_INTEGRATOR_CONSTANT_STEPSIZE_HPP_INCLUDED
-#include <boost/numeric/odeint/observer.hpp>
+#include <boost/numeric/odeint/integrate_functions/observer.hpp>
namespace boost {
namespace numeric {
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_adaptive_stepsize.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_adaptive_stepsize.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,174 +0,0 @@
-/* Boost odeint/integrator_adaptive.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- This file includes integration methods with adaptive stepsize.
-
- 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_INTEGRATOR_ADAPTIVE_STEPSIZE_HPP
-#define BOOST_NUMERIC_ODEINT_INTEGRATOR_ADAPTIVE_STEPSIZE_HPP
-
-#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
-#include <boost/numeric/odeint/observer.hpp>
-#include <vector>
-#include <limits>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- template<
- class ControlledStepper,
- class DynamicalSystem,
- class Observer
- >
- size_t integrate_adaptive(
- ControlledStepper &stepper,
- DynamicalSystem &system,
- typename ControlledStepper::container_type &state,
- typename ControlledStepper::time_type start_time,
- typename ControlledStepper::time_type end_time,
- typename ControlledStepper::time_type dt,
- Observer &observer )
- {
- typedef typename ControlledStepper::time_type time_type;
-
- stepper.adjust_size( state );
-
- controlled_step_result result;
- size_t iterations = 0;
- time_type t = start_time , dt_ = dt;
-
- observer(t, state, system);
-
- while( t < end_time )
- {
- if( (end_time - t) < dt ) dt = end_time - t;
- // do a controlled step
- result = stepper.try_step( system, state, t, dt_ );
-
- if( result != step_size_decreased )
- { // we actually did a step forward (dt was small enough)
- observer(t, state, system);
- iterations++;
- }
-
- if( !( t+dt_ > t) )
- throw; // we've reached machine precision with dt - no advancing in t
- }
- return iterations;
- }
-
- template<
- class ControlledStepper,
- class DynamicalSystem
- >
- size_t integrate_adaptive(
- ControlledStepper &stepper,
- DynamicalSystem &system,
- typename ControlledStepper::container_type &state,
- typename ControlledStepper::time_type start_time,
- typename ControlledStepper::time_type end_time,
- typename ControlledStepper::time_type dt )
- {
- return integrate_adaptive(
- stepper, system ,
- state, start_time , end_time, dt,
- do_nothing_observer<
- typename ControlledStepper::time_type ,
- typename ControlledStepper::container_type ,
- DynamicalSystem >
- );
- }
-
-
-
-
- /* Integration of an ode with adaptive stepsize methods.
- Integrates an ode give by system using the integration scheme stepper and the
- step-size controller controller.
- The initial state is given in x.
- t is an vector including the times at which the state will be written into
- the vector x_vec.
- x_vec must provide enough space to hold times.size() states.
- dt is the initial step size (will be adjusted according to the controller).
- This function returns the total number of steps required to integrate the
- whole intervale times.begin() - times.end().
- Note that the values in times don't influence the stepsize, but only the
- time points at which the state is stored into x_vec.
- */
- template<
- class ControlledStepper,
- class DynamicalSystem,
- class TimeInsertIterator,
- class StateInsertIterator
- >
- size_t integrate(
- ControlledStepper &stepper,
- DynamicalSystem &system,
- typename ControlledStepper::container_type &state,
- typename ControlledStepper::time_type start,
- typename ControlledStepper::time_type end,
- typename ControlledStepper::time_type dt,
- TimeInsertIterator time_inserter,
- StateInsertIterator state_inserter)
- {
- state_copy_observer<TimeInsertIterator, StateInsertIterator>
- observer(time_inserter, state_inserter);
- return integrate_adaptive(stepper, system, state,
- start , end, dt , observer);
- }
-
-
- /*
- Integration of an ode with adaptive stepsize methods.
- Integrates an ode give by system using the integration scheme stepper and the
- a standard step-size controller that ensures the error being below the values
- given below.
- */
- template<
- class DynamicalSystem,
- class ContainerType,
- class TimeInsertIterator,
- class StateInsertIterator,
- class T
- >
- size_t integrate(
- DynamicalSystem &system,
- ContainerType &state,
- T start_time ,
- T end_time ,
- TimeInsertIterator time_inserter,
- StateInsertIterator state_inserter,
- T dt = 1E-4,
- T eps_abs = 1E-6,
- T eps_rel = 1E-7,
- T a_x = 1.0 ,
- T a_dxdt = 1.0
- )
- {
- // we use cash karp stepper as base stepper
- typedef stepper_rk5_ck< ContainerType , T > stepper_type;
-
- // we use the standard controller for this adaptive integrator
- controlled_stepper_standard< stepper_type >
- controlled_stepper( eps_abs, eps_rel, a_x, a_dxdt );
- // initialized with values from above
-
- // call the normal integrator
- return integrate(controlled_stepper, system, state,
- start_time, end_time, dt, time_inserter, state_inserter);
- }
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-#endif
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_constant_stepsize.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/integrator_constant_stepsize.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,138 +0,0 @@
-/*
- boost header: numeric/odeint/integrator_constant_step.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_INTEGRATOR_CONSTANT_STEPSIZE_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_INTEGRATOR_CONSTANT_STEPSIZE_HPP_INCLUDED
-
-#include <boost/numeric/odeint/observer.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- template<
- class Stepper ,
- class DynamicalSystem ,
- class Observer
- >
- size_t integrate_const(
- Stepper &stepper ,
- DynamicalSystem &system ,
- typename Stepper::container_type &state ,
- typename Stepper::time_type start_time ,
- typename Stepper::time_type end_time ,
- typename Stepper::time_type dt ,
- Observer &observer
- )
- {
- stepper.adjust_size( state );
- size_t iteration = 0;
- while( start_time < end_time )
- {
- observer( start_time , state , system );
- stepper.do_step( system , state , start_time , dt );
- start_time += dt;
- ++iteration;
- }
- observer( start_time , state , system );
-
- return iteration;
- }
-
-
-
- template<
- class Stepper ,
- class DynamicalSystem
- >
- size_t integrate_const(
- Stepper &stepper ,
- DynamicalSystem &system ,
- typename Stepper::container_type &state ,
- typename Stepper::time_type start_time ,
- typename Stepper::time_type end_time ,
- typename Stepper::time_type dt
- )
- {
- return integrate_const(
- stepper , system , state, start_time , end_time , dt ,
- do_nothing_observer<
- typename Stepper::time_type ,
- typename Stepper::container_type ,
- DynamicalSystem >
- );
- }
-
-
-
- template<
- class Stepper ,
- class DynamicalSystem ,
- class Observer
- >
- typename Stepper::time_type integrate_const_steps(
- Stepper &stepper ,
- DynamicalSystem &system ,
- typename Stepper::container_type &state ,
- typename Stepper::time_type start_time ,
- typename Stepper::time_type dt ,
- size_t num_of_steps ,
- Observer &observer
- )
- {
- stepper.adjust_size( state );
- size_t iteration = 0;
- while( iteration < num_of_steps )
- {
- observer( start_time , state , system );
- stepper.do_step( system , state , start_time , dt );
- start_time += dt;
- ++iteration;
- }
- observer( start_time , state , system );
-
- return start_time;
- }
-
-
- template<
- class Stepper ,
- class DynamicalSystem
- >
- typename Stepper::time_type integrate_const_steps(
- Stepper &stepper ,
- DynamicalSystem &system ,
- typename Stepper::container_type &state ,
- typename Stepper::time_type start_time ,
- typename Stepper::time_type dt ,
- size_t num_of_steps
- )
- {
- return integrate_const_steps(
- stepper , system , state , start_time , dt , num_of_steps ,
- do_nothing_observer<
- typename Stepper::time_type ,
- typename Stepper::container_type ,
- DynamicalSystem >
- );
- }
-
-
-
-
-
-} // odeint
-} // numeric
-} // boost
-
-#endif //BOOST_NUMERIC_ODEINT_INTEGRATOR_CONSTANT_STEPSIZE_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/observer.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/observer.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,63 +0,0 @@
-/*
- boost header: numeric/odeint/observer.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_OBSERVER_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_OBSERVER_HPP_INCLUDED
-
-#include<vector>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- template< class Time, class Container , class System >
- inline void do_nothing_observer( Time , Container& , System& )
- {
- }
-
-
-
- template< class TimeInsertIterator, class StateInsertIterator >
- class state_copy_observer
- {
-
- private:
-
- TimeInsertIterator m_time_inserter;
- StateInsertIterator m_state_inserter;
-
- public:
-
- state_copy_observer( TimeInsertIterator time_inserter ,
- StateInsertIterator state_inserter )
- : m_time_inserter(time_inserter),
- m_state_inserter(state_inserter)
- { }
-
- void reset( void ) { }
-
- template< class Time, class Container, class System >
- void operator () (Time t, Container &state, System &system )
- {
- *m_time_inserter++ = t; // insert time
- *m_state_inserter++ = state; // insert state
- }
- };
-
-
-} // odeint
-} // numeric
-} // boost
-
-
-#endif //BOOST_NUMERIC_ODEINT_OBSERVER_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_base.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_base.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,156 +0,0 @@
-/*
- boost header: numeric/odeint/stepper_base.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_BASE_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_STEPPER_BASE_HPP_INCLUDED
-
-#include <boost/noncopyable.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- typedef unsigned char order_type;
-
- template<
- class Stepper ,
- class Container ,
- order_type Order ,
- class Time ,
- class Traits
- >
- class explicit_stepper_base : private boost::noncopyable
- {
- public:
-
- // some typedef
-
- typedef Stepper stepper_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
- public:
-
- // provide some functions
-
- explicit_stepper_base( stepper_type &stepper )
- : m_stepper( stepper ) { }
-
- explicit_stepper_base( stepper_type &stepper , container_type &x )
- : m_stepper( stepper )
- {
- traits_type::adjust_size( x , m_dxdt );
- }
-
-// order_type order() const { return m_order; }
- const static order_type order = Order;
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- system( x , m_dxdt , t );
- stepper.do_step( system , x , m_dxdt , t , dt );
- }
-
- void adjust_size( container_type &x )
- {
- traits_type::adjust_size( x , m_dxdt );
- }
-
- private:
-
- container_type m_dxdt;
- stepper_type &m_stepper;
- };
-
-
-
-
- template<
- class ErrorStepper ,
- class Container ,
- order_type StepperOrder ,
- order_type ErrorOrder
- class Time ,
- class Traits
- >
- class explicit_error_stepper_base : private boost::noncopyable
- {
- public:
-
- // some typedef
-
- typedef ErrorStepper stepper_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
- public:
-
- // provide some functions
-
- explicit_stepper_base( stepper_type &stepper )
- : m_stepper( stepper ) { }
-
- explicit_stepper_base( stepper_type &stepper , container_type &x )
- : m_stepper( stepper )
- {
- traits_type::adjust_size( x , m_dxdt );
- }
-
-// order_type order() const { return m_order; }
- const static order_type stepper_order = StepperOrder;
- const static order_type error_order = ErrorOrder;
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt ,
- container_type &err )
- {
- system( x , m_dxdt , t );
- stepper.do_step_with_deriv( system , x , m_dxdt , t , dt , err );
- }
-
- void adjust_size( container_type &x )
- {
- traits_type::adjust_size( x , m_dxdt );
- }
-
- private:
-
- container_type m_dxdt;
- stepper_type &m_stepper;
- };
-
-}
-}
-}
-
-#endif //BOOST_NUMERIC_ODEINT_STEPPER_BASE_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,130 +0,0 @@
-/* Boost odeint/stepper_euler.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- This file includes the explicit euler 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_EULER_HPP
-#define BOOST_NUMERIC_ODEINT_STEPPER_EULER_HPP
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-// #include <boost/numeric/odeint/stepper_base.hpp>
-
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_euler
- {
- //
- // provide basic typedefs
- //
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
-
-
- //
- // private members
- //
- private:
-
- container_type m_dxdt;
-
-
-
- //
- // public interface
- //
- public:
-
- // return the order of the stepper
- order_type order_step() const { return 1; }
-
-
- // standard constructor, m_dxdt is not adjusted
- stepper_euler( void )
- {
- }
-
-
-
- // contructor, which adjusts m_dxdt
- stepper_euler( const container_type &x )
- {
- adjust_size( x );
- }
-
-
-
- // adjust the size of m_dxdt
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size( x , m_dxdt );
- }
-
-
-
- // performs one step with the knowledge of dxdt(t)
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- // x = x + dt*dxdt
- detail::it_algebra::increment( traits_type::begin(x) ,
- traits_type::end(x) ,
- traits_type::begin(dxdt) ,
- dt );
- }
-
-
-
- // performs one step
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- system( x , m_dxdt , t );
- do_step( system , x , m_dxdt , t , dt );
- }
-
- };
-
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_STEPPER_EULER_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_half_step.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,179 +0,0 @@
-/* Boost odeint/stepper_half_stepr.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- This file includes a stepper which calculates the
- error during one step from performing two steps with
- the halt stepsize. It works with arbitray steppers
-
- 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_HALF_STEP_HPP
-#define BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
-
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-
-#include <iostream>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- template< class Stepper >
- class stepper_half_step
- {
- //
- // provide basic typedefs
- //
- public:
-
- typedef Stepper stepper_type;
- typedef typename stepper_type::container_type container_type;
- typedef typename stepper_type::traits_type traits_type;
- typedef typename stepper_type::time_type time_type;
- typedef typename stepper_type::order_type order_type;
- typedef typename stepper_type::value_type value_type;
-// typedef typename stepper_type::iterator iterator;
-// typedef typename stepper_type::const_iterator const_iterator;
-
-
-
- //
- // private members
- //
- private:
-
- container_type m_dxdt;
- container_type m_xtemp;
- stepper_type m_stepper;
-
-
- //
- // public interface
- //
- public:
-
-
- // the order of the step if a normal step is performed
- order_type order_step( void ) const
- {
- return m_stepper.order_step();
- }
-
-
- // the order of the step if an error step is performed
- order_type order_error_step( void ) const
- {
- return m_stepper.order_step();
- }
-
-
- // the order of the error term if the error step is performed
- order_type order_error( void ) const
- {
- return m_stepper.order_step() + 1;
- }
-
-
- // standard constructor
- stepper_half_step( void )
- {
- }
-
-
- // contructor, which adjust the size of internal containers
- stepper_half_step( const container_type &x )
- {
- adjust_size( x );
- }
-
-
- // adjust the size of m_dxdt , m_xtemp und m_stepper
- void adjust_size( const container_type &x )
- {
- m_stepper.adjust_size( x );
- traits_type::adjust_size( x , m_dxdt );
- traits_type::adjust_size( x , m_xtemp );
- }
-
-
-
- // performs a normal step, without error calculation
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- m_stepper.do_step( system , x , dxdt , t , dt );
- }
-
-
- // performs a normal step, without error calculation
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- m_stepper.do_step( system , x , t , dt );
- }
-
-
-
- // performs a error step with error calculation
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt ,
- container_type &xerr )
- {
- time_type dt2 = static_cast<time_type>(0.5) * dt;
-
- m_xtemp = x;
-
- do_step( system , m_xtemp , dxdt , t , dt );
- do_step( system , x , dxdt , t , dt2 );
- do_step( system , x , t+dt2 , dt2 );
-
- detail::it_algebra::scale_sum( traits_type::begin(xerr) ,
- traits_type::end(xerr) ,
- static_cast< value_type >(1.0),
- traits_type::begin(m_xtemp) ,
- static_cast< value_type >(-1.0),
- traits_type::begin(x) );
- }
-
-
-
-
- // performs a error step with error calculation
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt ,
- container_type &xerr )
- {
- system( x , m_dxdt , t );
- do_step( system , x , m_dxdt , t , dt , xerr );
- }
- };
-
-
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_midpoint.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_midpoint.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,193 +0,0 @@
-/* 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/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-#include <iostream>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_midpoint
- {
- //
- // provide basic typedefs
- //
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
-
- // private memebers
- private:
-
- unsigned short m_step_number;
-
- container_type m_x0;
- container_type m_x1;
- container_type m_dxdt;
-
-
-
-
- public:
-
- order_type order_step() const { return 2; }
-
- // standard constructor, the size of the internal container is not set
- stepper_midpoint( unsigned short step_number = 2 )
- {
- set_step_number( step_number );
- }
-
- // constructor, which adjusts the size of the internal containers
- stepper_midpoint( const container_type &x , unsigned short step_number = 2 )
- {
- adjust_size( x );
- set_step_number( step_number );
- }
-
- // adjusts the size of the internal containers
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size( x , m_x0 );
- traits_type::adjust_size( x , m_x1 );
- traits_type::adjust_size( x , m_dxdt );
- }
-
- void set_step_number( unsigned short step_number )
- {
- if( step_number > 1 )
- m_step_number = step_number;
- }
-
- unsigned short get_step_number() const
- {
- return m_step_number;
- }
-
-
-
- // performs a midpoint step
- template< class DynamicalSystem >
- void midpoint_step(
- DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt ,
- container_type &x_out )
- {
- using namespace detail::it_algebra;
-
- const time_type t_1 = static_cast<time_type>( 1.0 );
- const time_type t_05 = static_cast<time_type>( 0.5 );
-
- const time_type h = dt / static_cast<time_type>( m_step_number );
- const time_type h2 = static_cast<time_type>( 2.0 ) * h;
-
-
- time_type th = t + h;
-
- // m_x1 = x + h*dxdt
- scale_sum( traits_type::begin(m_x1),
- traits_type::end(m_x1),
- t_1, traits_type::begin(x),
- h, traits_type::begin(dxdt) );
- system( m_x1, m_dxdt, th );
-
- m_x0 = x;
-
- unsigned short i = 1;
- while( i != m_step_number )
- {
- // general step
- //tmp = m_x1; m_x1 = m_x0 + h2*m_dxdt; m_x0 = tmp
- scale_sum_swap( traits_type::begin(m_x1),
- traits_type::end(m_x1),
- traits_type::begin(m_x0),
- h2, traits_type::begin(m_dxdt) );
- th += h;
- system( m_x1, m_dxdt, th);
- i++;
- }
-
- // last step
- // x = 0.5*( m_x0 + m_x1 + h*m_dxdt )
- scale_sum( traits_type::begin(x_out), traits_type::end(x_out),
- t_05, traits_type::begin(m_x0),
- t_05, traits_type::begin(m_x1),
- t_05*h, traits_type::begin(m_dxdt) );
- }
-
-
-
-
-
- template< class DynamicalSystem >
- void do_step(
- DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- midpoint_step( system, x, dxdt, t, dt, x );
- }
-
-
-
-
-
- template< class DynamicalSystem >
- void do_step(
- DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- system( x, m_dxdt, t );
- do_step( system , x, m_dxdt, t, dt );
- }
-
-
- };
-
-}
-}
-}
-
-#endif
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk4.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk4.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,164 +0,0 @@
-/* Boost odeint/stepper_rk4.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- This file includes the explicit 4th order runge kutta
- solver for ordinary differential equations.
-
- It solves any ODE dx/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_RK4_HPP
-#define BOOST_NUMERIC_ODEINT_STEPPER_RK4_HPP
-
-#include <boost/numeric/odeint/container_traits.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_rk4
- {
-
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
-
-
- // private members
- private:
-
- container_type m_dxdt;
- container_type m_dxt;
- container_type m_dxm;
- container_type m_dxh;
- container_type m_xt;
-
-
- // private member functions
- private:
-
-
- // public interface
- public:
-
-
- order_type order_step() const { return 4; }
-
- // standard constructor, internal containers are not initialized
- stepper_rk4( void )
- {
- }
-
- // constructor, which adjusts the internal containers
- stepper_rk4( const container_type &x )
- {
- adjust_size( x );
- }
-
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size( x , m_dxdt );
- traits_type::adjust_size( x , m_dxt );
- traits_type::adjust_size( x , m_dxm );
- traits_type::adjust_size( x , m_xt );
- traits_type::adjust_size( x , m_dxh );
- }
-
-
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- using namespace detail::it_algebra;
-
- const time_type val1 = static_cast<time_type>( 1.0 );
-
- time_type dh = static_cast<time_type>( 0.5 ) * dt;
- time_type th = t + dh;
-
- // dt * dxdt = k1
- // m_xt = x + dh*dxdt
- scale_sum( traits_type::begin(m_xt),
- traits_type::end(m_xt),
- val1, traits_type::begin(x),
- dh, traits_type::begin(dxdt) );
-
-
- // dt * m_dxt = k2
- system( m_xt , m_dxt , th );
- // m_xt = x + dh*m_dxt
- scale_sum( traits_type::begin(m_xt) ,
- traits_type::end(m_xt) ,
- val1, traits_type::begin(x) ,
- dh, traits_type::begin(m_dxt) );
-
-
- // dt * m_dxm = k3
- system( m_xt , m_dxm , th );
- //m_xt = x + dt*m_dxm
- scale_sum( traits_type::begin(m_xt), traits_type::end(m_xt),
- val1, traits_type::begin(x) ,
- dt, traits_type::begin(m_dxm) );
-
-
- // dt * m_dxh = k4
- system( m_xt , m_dxh , t + dt );
- //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
- time_type dt6 = dt / static_cast<time_type>( 6.0 );
- time_type dt3 = dt / static_cast<time_type>( 3.0 );
- scale_sum_inplace( traits_type::begin(x) , traits_type::end(x),
- dt6 , traits_type::begin(dxdt),
- dt3 , traits_type::begin(m_dxt),
- dt3 , traits_type::begin(m_dxm),
- dt6 , traits_type::begin(m_dxh) );
- }
-
-
-
-
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- system( x , m_dxdt , t );
- do_step( system , x , m_dxdt , t , dt );
- }
-
- };
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_STEPPER_RK4_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk4_classical.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,154 +0,0 @@
-/* Boost odeint/stepper_rk4.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- This file includes the explicit runge kutta solver for
- ordinary differential equations.
-
- It solves any ODE dx/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_RK4_CLASSICAL_HPP
-#define BOOST_NUMERIC_ODEINT_STEPPER_RK4_CLASSICAL_HPP
-
-#include <boost/numeric/odeint/container_traits.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_rk4_classical
- {
-
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
-
-
- // private members
- private:
-
- container_type m_dxdt;
- container_type m_dxt;
- container_type m_dxm;
- container_type m_xt;
-
-
-
-
-
- // public interface
- public:
-
- order_type order_step() const { return 4; }
-
- // standard constructor, internal containers are not initialized
- stepper_rk4_classical( void )
- {
- }
-
- // constructor, which adjusts the internal containers
- stepper_rk4_classical( const container_type &x )
- {
- adjust_size( x );
- }
-
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size( x , m_dxdt );
- traits_type::adjust_size( x , m_dxt );
- traits_type::adjust_size( x , m_dxm );
- traits_type::adjust_size( x , m_xt );
- }
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- using namespace detail::it_algebra;
-
- const time_type val2 = time_type( 2.0 );
-
-
- time_type dh = time_type( 0.5 ) * dt;
- time_type th = t + dh;
-
- //m_xt = x + dh*dxdt
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- static_cast< time_type >(1.0) ,
- traits_type::begin(x) ,
- dh,
- traits_type::begin(dxdt) );
-
-
- //m_xt = x + dh*m_dxdt
- system( m_xt , m_dxt , th );
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- static_cast< time_type >(1.0),
- traits_type::begin(x) ,
- dh ,
- traits_type::begin(m_dxt) );
-
-
- //m_xt = x + dt*m_dxm ; m_dxm += m_dxt
- system( m_xt , m_dxm , th );
- assign_sum_increment( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- traits_type::begin(x) , traits_type::begin(m_dxm) ,
- traits_type::begin(m_dxt) , dt );
-
-
- //x = dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
- system( m_xt , m_dxt , t + dt );
- increment_sum_sum( traits_type::begin(x) , traits_type::end(x) ,
- traits_type::begin(dxdt) ,
- traits_type::begin(m_dxt) ,
- traits_type::begin(m_dxm) ,
- dt / time_type( 6.0 ) , val2 );
- }
-
-
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- system( x , m_dxdt , t );
- do_step( system , x , m_dxdt , t , dt );
- }
-
-
- };
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif // BOOST_NUMERIC_ODEINT_STEPPER_RK4_CLASSICAL_HPP
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk5_ck.hpp
==============================================================================
Binary file. No diff available.
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk78_fehlberg.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk78_fehlberg.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,415 +0,0 @@
-/*
- boost header: numeric/odeint/stepper_rk78_fehlberg.hpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- 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_RK78_FEHLBERG_HPP_INCLUDED
-#define BOOST_NUMERIC_ODEINT_STEPPER_RK78_FEHLBERG_HPP_INCLUDED
-
-#include <boost/numeric/odeint/container_traits.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_rk78_fehlberg
- {
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
- // private members
- private:
-
- 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_k13;
-
-
- // the times at which system is called
- time_type m_t02 , m_t03 , m_t04 , m_t05 , m_t06 , m_t07 , m_t08 ,
- m_t09 , m_t10 , m_t11 , m_t12 , m_t13;
-
-
-
-
- // public interface
- public:
-
- order_type order_step( void ) const { return 8; }
-
- order_type order_error_step( void ) const { return 7; }
-
- order_type order_error( void ) const { return 8; }
-
- // standard constructor, leaves the internal containers uninitialized
- stepper_rk78_fehlberg( void )
- {
- }
-
-
- // constructor, which adjusts the internal containers
- stepper_rk78_fehlberg( const container_type &x )
- {
- adjust_size( x );
- }
-
-
- void adjust_size( const container_type &x )
- {
- traits_type::adjust_size( x , m_dxdt );
- traits_type::adjust_size( x , m_xt );
- traits_type::adjust_size( x , m_k02 );
- traits_type::adjust_size( x , m_k03 );
- traits_type::adjust_size( x , m_k04 );
- traits_type::adjust_size( x , m_k05 );
- traits_type::adjust_size( x , m_k06 );
- traits_type::adjust_size( x , m_k07 );
- traits_type::adjust_size( x , m_k08 );
- traits_type::adjust_size( x , m_k09 );
- traits_type::adjust_size( x , m_k10 );
- traits_type::adjust_size( x , m_k11 );
- traits_type::adjust_size( x , m_k12 );
- traits_type::adjust_size( x , m_k13 );
- }
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- // the time constant
- const time_type a02 = static_cast<time_type>( 2.0 / 27.0 );
- const time_type a03 = static_cast<time_type>( 1.0 / 9.0 );
- const time_type a04 = static_cast<time_type>( 1.0 / 6.0 );
- const time_type a05 = static_cast<time_type>( 5.0 / 12.0 );
- const time_type a06 = static_cast<time_type>( 0.50 );
- const time_type a07 = static_cast<time_type>( 5.0 / 6.0 );
- const time_type a08 = static_cast<time_type>( 1.0 / 6.0 );
- const time_type a09 = static_cast<time_type>( 2.0 / 3.0 );
- const time_type a10 = static_cast<time_type>( 1.0 / 3.0 );
-
- // the weights for each step
- const time_type c06 = static_cast<time_type>( 34.0 / 105.0 );
- const time_type c07 = static_cast<time_type>( 9.0 / 35.0 );
- const time_type c08 = static_cast<time_type>( c07 );
- const time_type c09 = static_cast<time_type>( 9.0 / 280.0 );
- 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 );
-
- // the coefficients for each step
- const time_type b02_01 = static_cast<time_type>( 2.0 / 27.0 );
-
- const time_type b03_01 = static_cast<time_type>( 1.0 / 36.0 );
- const time_type b03_02 = static_cast<time_type>( 1.0 / 12.0 );
-
- const time_type b04_01 = static_cast<time_type>( 1.0 / 24.0 );
- const time_type b04_03 = static_cast<time_type>( 1.0 / 8.0 );
-
- const time_type b05_01 = static_cast<time_type>( 5.0 / 12.0 );
- const time_type b05_03 = static_cast<time_type>( -25.0 / 16.0 );
- const time_type b05_04 = static_cast<time_type>( -b05_03 );
-
- const time_type b06_01 = static_cast<time_type>( 0.050 );
- const time_type b06_04 = static_cast<time_type>( 0.250 );
- const time_type b06_05 = static_cast<time_type>( 0.20 );
-
- const time_type b07_01 = static_cast<time_type>( -25.0 / 108.0 );
- const time_type b07_04 = static_cast<time_type>( 125.0 / 108.0 );
- const time_type b07_05 = static_cast<time_type>( -65.0 / 27.0 );
- const time_type b07_06 = static_cast<time_type>( 125.0 / 54.0 );
-
- const time_type b08_01 = static_cast<time_type>( 31.0 / 300.0 );
- const time_type b08_05 = static_cast<time_type>( 61.0 / 225.0 );
- const time_type b08_06 = static_cast<time_type>( -2.0 / 9.0 );
- const time_type b08_07 = static_cast<time_type>( 13.0 / 900.0 );
-
- const time_type b09_01 = static_cast<time_type>( 2.0 );
- const time_type b09_04 = static_cast<time_type>( -53.0 / 6.0 );
- const time_type b09_05 = static_cast<time_type>( 704.0 / 45.0 );
- const time_type b09_06 = static_cast<time_type>( -107.0 / 9.0 );
- const time_type b09_07 = static_cast<time_type>( 67.0 / 90.0 );
- const time_type b09_08 = static_cast<time_type>( 3.0 );
-
- const time_type b10_01 = static_cast<time_type>( -91.0 / 108.0 );
- const time_type b10_04 = static_cast<time_type>( 23.0 / 108.0 );
- const time_type b10_05 = static_cast<time_type>( -976.0 / 135.0 );
- const time_type b10_06 = static_cast<time_type>( 311.0 / 54.0 );
- const time_type b10_07 = static_cast<time_type>( -19.0 / 60.0 );
- const time_type b10_08 = static_cast<time_type>( 17.0 / 6.0 );
- const time_type b10_09 = static_cast<time_type>( -1.0 / 12.0 );
-
- const time_type b11_01 = static_cast<time_type>( 2383.0 / 4100.0 );
- const time_type b11_04 = static_cast<time_type>( -341.0 / 164.0 );
- const time_type b11_05 = static_cast<time_type>( 4496.0 / 1025.0 );
- const time_type b11_06 = static_cast<time_type>( -301.0 / 82.0 );
- const time_type b11_07 = static_cast<time_type>( 2133.0 / 4100.0 );
- const time_type b11_08 = static_cast<time_type>( 45.0 / 82.0 );
- const time_type b11_09 = static_cast<time_type>( 45.0 / 164.0 );
- const time_type b11_10 = static_cast<time_type>( 18.0 / 41.0 );
-
- const time_type b12_01 = static_cast<time_type>( 3.0 / 205.0 );
- const time_type b12_06 = static_cast<time_type>( -6.0 / 41.0 );
- const time_type b12_07 = static_cast<time_type>( -3.0 / 205.0 );
- const time_type b12_08 = static_cast<time_type>( -3.0 / 41.0 );
- const time_type b12_09 = static_cast<time_type>( 3.0 / 41.0 );
- const time_type b12_10 = static_cast<time_type>( 6.0 / 41.0 );
-
- const time_type b13_01 = static_cast<time_type>( -1777.0 / 4100.0 );
- const time_type b13_04 = static_cast<time_type>( b11_04 );
- const time_type b13_05 = static_cast<time_type>( b11_05 );
- const time_type b13_06 = static_cast<time_type>( -289.0 / 82.0 );
- const time_type b13_07 = static_cast<time_type>( 2193.0 / 4100.0 );
- const time_type b13_08 = static_cast<time_type>( 51.0 / 82.0 );
- const time_type b13_09 = static_cast<time_type>( 33.0 / 164.0 );
- const time_type b13_10 = static_cast<time_type>( 12.0 / 41.0 );
- const time_type b13_12 = static_cast<time_type>( 1.0 );
-
- const time_type val1 = static_cast<time_type>( 1.0 );
-
- using namespace detail::it_algebra;
-
- // compute the times at which system is evaluated
- m_t02 = t + a02 * dt;
- m_t03 = t + a03 * dt;
- m_t04 = t + a04 * dt;
- m_t05 = t + a05 * dt;
- m_t06 = t + a06 * dt;
- m_t07 = t + a07 * dt;
- m_t08 = t + a08 * dt;
- m_t09 = t + a09 * dt;
- m_t10 = t + a10 * dt;
- m_t11 = t + dt;
- m_t12 = t;
- m_t13 = t + dt;
-
- // k1, the first system call has allready been evaluated
-
- // k2 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b02_01 , traits_type::begin(dxdt) );
- system( m_xt , m_k02 , m_t02 );
-
- // k3 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b03_01 , traits_type::begin(dxdt) ,
- dt * b03_02 , traits_type::begin(m_k02) );
- system( m_xt , m_k03 , m_t03 );
-
-
- // k4 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b04_01 , traits_type::begin(dxdt) ,
- dt * b04_03 , traits_type::begin(m_k03) );
- system( m_xt , m_k04 , m_t04 );
-
-
- // k5 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b05_01 , traits_type::begin(dxdt) ,
- dt * b05_03 , traits_type::begin(m_k03) ,
- dt * b05_04 , traits_type::begin(m_k04) );
- system( m_xt , m_k05 , m_t05 );
-
-
- // k6 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b06_01 , traits_type::begin(dxdt) ,
- dt * b06_04 , traits_type::begin(m_k04) ,
- dt * b06_05 , traits_type::begin(m_k05) );
- system( m_xt , m_k06 , m_t06 );
-
-
- // k7 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b07_01 , traits_type::begin(dxdt) ,
- dt * b07_04 , traits_type::begin(m_k04) ,
- dt * b07_05 , traits_type::begin(m_k05) ,
- dt * b07_06 , traits_type::begin(m_k06) );
- system( m_xt , m_k07 , m_t07 );
-
-
- // k8 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b08_01 , traits_type::begin(dxdt) ,
- dt * b08_05 , traits_type::begin(m_k05) ,
- dt * b08_06 , traits_type::begin(m_k06) ,
- dt * b08_07 , traits_type::begin(m_k07) );
- system( m_xt , m_k08 , m_t08 );
-
-
- // k9 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b09_01 , traits_type::begin(dxdt) ,
- dt * b09_04 , traits_type::begin(m_k04) ,
- dt * b09_05 , traits_type::begin(m_k05) ,
- dt * b09_06 , traits_type::begin(m_k06) ,
- dt * b09_07 , traits_type::begin(m_k07) ,
- dt * b09_08 , traits_type::begin(m_k08) );
- system( m_xt , m_k09 , m_t09 );
-
-
- // k10 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b10_01 , traits_type::begin(dxdt) ,
- dt * b10_04 , traits_type::begin(m_k04) ,
- dt * b10_05 , traits_type::begin(m_k05) ,
- dt * b10_06 , traits_type::begin(m_k06) ,
- dt * b10_07 , traits_type::begin(m_k07) ,
- dt * b10_08 , traits_type::begin(m_k08) ,
- dt * b10_09 , traits_type::begin(m_k09) );
- system( m_xt , m_k10 , m_t10 );
-
-
- // k11 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b11_01 , traits_type::begin(dxdt) ,
- dt * b11_04 , traits_type::begin(m_k04) ,
- dt * b11_05 , traits_type::begin(m_k05) ,
- dt * b11_06 , traits_type::begin(m_k06) ,
- dt * b11_07 , traits_type::begin(m_k07) ,
- dt * b11_08 , traits_type::begin(m_k08) ,
- dt * b11_09 , traits_type::begin(m_k09) ,
- dt * b11_10 , traits_type::begin(m_k10) );
- system( m_xt , m_k11 , m_t11 );
-
-
- // k12 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b12_01 , traits_type::begin(dxdt) ,
- dt * b12_06 , traits_type::begin(m_k06) ,
- dt * b12_07 , traits_type::begin(m_k07) ,
- dt * b12_08 , traits_type::begin(m_k08) ,
- dt * b12_09 , traits_type::begin(m_k09) ,
- dt * b12_10 , traits_type::begin(m_k10) );
- system( m_xt , m_k12 , m_t12 );
-
-
- // k13 step
- scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
- val1 , traits_type::begin(x) ,
- dt * b13_01 , traits_type::begin(dxdt) ,
- dt * b13_04 , traits_type::begin(m_k04) ,
- dt * b13_05 , traits_type::begin(m_k05) ,
- dt * b13_06 , traits_type::begin(m_k06) ,
- dt * b13_07 , traits_type::begin(m_k07) ,
- dt * b13_08 , traits_type::begin(m_k08) ,
- dt * b13_09 , traits_type::begin(m_k09) ,
- dt * b13_10 , traits_type::begin(m_k10) ,
- dt * b13_12 , traits_type::begin(m_k12) );
- system( m_xt , m_k13 , m_t13 );
-
- scale_sum( traits_type::begin(x) , traits_type::end(x) ,
- val1 , traits_type::begin(x) ,
- dt * c06 , traits_type::begin(m_k06) ,
- dt * c07 , traits_type::begin(m_k07) ,
- dt * c08 , traits_type::begin(m_k08) ,
- dt * c09 , traits_type::begin(m_k09) ,
- dt * c10 , traits_type::begin(m_k10) ,
- dt * c12 , traits_type::begin(m_k12) ,
- dt * c13 , traits_type::begin(m_k13) );
- }
-
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- 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>( cc07 );
- 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 );
-
- using namespace detail::it_algebra;
-
- xerr = x;
- do_step( system , xerr , dxdt , t , dt );
-
- // now, k1-k13 are calculated and stored in m_k01 - m_k13
- scale_sum( traits_type::begin(x) , traits_type::end(x) ,
- static_cast<time_type>(1.0) , traits_type::begin(x) ,
- dt * cc01 , traits_type::begin(dxdt) ,
- dt * cc06 , traits_type::begin(m_k06) ,
- dt * cc07 , traits_type::begin(m_k07) ,
- dt * cc08 , traits_type::begin(m_k08) ,
- dt * cc09 , traits_type::begin(m_k09) ,
- dt * cc10 , traits_type::begin(m_k10) ,
- dt * cc11 , traits_type::begin(m_k11) );
-
- increment( traits_type::begin(xerr) , traits_type::end(xerr) ,
- traits_type::begin(x) , 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 )
- {
- system( x , m_dxdt , t );
- do_step( system , x , m_dxdt , t , dt , xerr );
- }
- };
-
-} // namespace odeint
-} // namespace numeric
-} // namespace boost
-
-
-#endif //BOOST_NUMERIC_ODEINT_STEPPER_RK78_FEHLBERG_HPP_INCLUDED
Deleted: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk_generic.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_rk_generic.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,438 +0,0 @@
-/* Boost odeint/stepper_rk_generic.hpp header file
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
- Copyright 2009 Andre Bergner
-
- This file includes generic Runge Kutta solver
- for ordinary differential equations.
-
- It solves any ODE dx/dt = f(x,t) using a general Runge Kutta scheme.
-
- 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_RK_GENERIC_HPP
-#define BOOST_NUMERIC_ODEINT_STEPPER_RK_GENERIC_HPP
-
-#include <vector>
-#include <exception>
-#include <cmath> // for pow( ) and abs()
-#include <limits>
-
-#include <boost/numeric/odeint/container_traits.hpp>
-
-
-namespace boost {
-namespace numeric {
-namespace odeint {
-
-
- class butcher_tableau_consistency_exception : public std::exception {
-
- virtual const char* what() const throw()
- {
- return "Consistency Check of Butcher Tableau failed!";
- }
- };
-
- class butcher_tableau_order_condition_exception : public std::exception {
-
- virtual const char* what() const throw()
- {
- return "Order Condition Check of Butcher Tableau failed!";
- }
- };
-
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_rk_generic
- {
-
-
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
-
-
-
-
- // private variables
- private:
-
- typedef std::vector< container_type > container_vector;
- typedef std::vector< iterator > container_iterator_vector;
- typedef std::vector< time_type > parameter_vector;
- typedef std::vector< parameter_vector > parameter_matrix;
-
- container_vector m_xvec;
- container_iterator_vector m_xiter_vec;
- container_type m_xtmp;
- const parameter_vector m_a;
- const parameter_matrix m_b;
- const parameter_vector m_c;
-
- order_type m_q;
-
-
-
- // private member functions
- private:
-
- void reset_iter(typename container_iterator_vector::iterator xiter_iter)
- {
- typename container_vector::iterator x_iter = m_xvec.begin();
- while( x_iter != m_xvec.end() ) {
- (*xiter_iter++) = traits_type::begin(*x_iter++);
- }
- }
-
- void check_consitency()
- {
- typename parameter_vector::const_iterator a_iter = m_a.begin();
- typename parameter_vector::const_iterator c_iter = m_c.begin();
- typename parameter_matrix::const_iterator b_iter = m_b.begin();
- typename parameter_vector::const_iterator b_iter_iter;
-
- // check 1: a_i = sum_j b_ij
- while( a_iter != m_a.end() ) {
- time_type tmp = 0.0;
- b_iter_iter = (*b_iter).begin();
- while( b_iter_iter != (*b_iter).end() ) {
- tmp += *b_iter_iter++;
- }
- b_iter++;
- if( *a_iter++ != tmp )
- throw butcher_tableau_consistency_exception();
- }
-
- // check 2: sum_i c_i * (a_i)^(k-1) = 1/k k = 1..q
- for( unsigned short k = 1; k <= m_q; k++ ) {
- time_type tmp = 0.0;
- a_iter = m_a.begin();
- c_iter = m_c.begin();
- if( k == 1 ) // special treatment for 0^0 = 1
- tmp += *c_iter++;
- else
- c_iter++;
- while( a_iter != m_a.end() ) {
- //std::clog<<pow( *a_iter , k-1 )<< ", ";
- tmp += (*c_iter++) * pow( *a_iter++ , k-1 );
- }
- //std::clog << tmp << " = " << time_type(1.0)/time_type(k) << "?" << std::endl;
- if( std::abs(time_type(tmp - time_type(1.0)/time_type(k))) >
- std::numeric_limits<time_type>::epsilon() ) {
- //std::clog << std::abs(time_type(tmp - time_type(1.0)/time_type(k))) << std::endl;
- throw butcher_tableau_order_condition_exception();
- }
- }
- }
-
-
- public:
-
-
-
- /* Constructor
-
- a,b,c are vectors providing the butcher tableau for the Runge Kutta scheme
-
- 0 |
- a_1 | b_21
- a_2 | b_31 b_32
- ... | ... ...
- a_q-1 | b_q1 b_q2 ... b_qq-1
- -------------------------------
- | c_1 c_2 ... c_q-1 c_q
-
- The size of c is determining the order of the scheme q.
- a is of length q-1
- b is of length q-1, b[0] of length 1, b[1] of length 2 and so on
- c is of length q (which defines q).
-
- There are 2 conditions that these parameters have to fullfill:
- Consitency:
-
- a_i = sum_j b_ij for i = 1 ... q-1
-
- Condition on the order of the method (all error terms dt^k , k<q+1 cancel out):
-
- sum_i c_i * (a_(i+1))^(k-1) = 1/k k = 1 ... q
-
- Note, that a_0 = 1 (implicitely) and 0^0 = 1
- so this means sum_i c_i = 1 at k=1
- */
- stepper_rk_generic( parameter_vector &a,
- parameter_matrix &b,
- parameter_vector &c)
- : m_a(a), m_b(b), m_c(c), m_q(c.size())
- {
- m_xvec.resize(m_q);
- m_xiter_vec.resize(m_q);
-
- check_consitency();
- }
-
- order_type order() const { return m_q; }
-
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- using namespace detail::it_algebra;
- typename container_vector::iterator x_iter = m_xvec.begin();
- typename container_iterator_vector::iterator xiter_iter = m_xiter_vec.begin();
-
- (*x_iter) = dxdt;
- (*xiter_iter++) = traits_type::begin(*x_iter++);
-
- while( x_iter != m_xvec.end() )
- {
- traits_type::adjust_size(x, (*x_iter));
- (*xiter_iter++) = traits_type::begin(*x_iter++);
- }
- traits_type::adjust_size(x, m_xtmp);
-
- x_iter = m_xvec.begin()+1;
-
- typename parameter_vector::const_iterator a_iter = m_a.begin();
- typename parameter_matrix::const_iterator b_iter = m_b.begin();
- while( x_iter != m_xvec.end() )
- {
- reset_iter(m_xiter_vec.begin());
- scale_sum_generic( traits_type::begin(m_xtmp), traits_type::end(m_xtmp),
- (*b_iter).begin(), (*b_iter).end(), dt,
- traits_type::begin(x), m_xiter_vec.begin() );
- system( m_xtmp, *x_iter++ , t + dt*(*a_iter++) );
- b_iter++;
- }
-
- reset_iter(m_xiter_vec.begin());
- typename parameter_vector::const_iterator c_iter = m_c.begin();
- scale_sum_generic( traits_type::begin(x), traits_type::end(x),
- m_c.begin(), m_c.end(), dt,
- traits_type::begin(x), m_xiter_vec.begin() );
- }
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- traits_type::adjust_size(x, m_xtmp);
- system(x, m_xtmp, t);
- do_step( system, x, m_xtmp, t, dt);
- }
-
- };
-
-
-
- /* ############################################################
- ############################################################
-
- C-Array Version of a,b,c handling
-
- ############################################################
- ############################################################
- */
-
-
-
-
- template<
- class Container ,
- class Time = double ,
- class Traits = container_traits< Container >
- >
- class stepper_rk_generic_ptr
- {
-
-
- // provide basic typedefs
- public:
-
- typedef unsigned short order_type;
- typedef Time time_type;
- typedef Traits traits_type;
- typedef typename traits_type::container_type container_type;
- typedef typename traits_type::value_type value_type;
- typedef typename traits_type::iterator iterator;
- typedef typename traits_type::const_iterator const_iterator;
-
- // private variables
- private:
-
- typedef std::vector< container_type > container_vector;
- typedef std::vector< iterator > container_iterator_vector;
- typedef std::vector< time_type > parameter_vector;
- typedef std::vector< parameter_vector > parameter_matrix;
-
- container_vector m_xvec;
- container_iterator_vector m_xiter_vec;
- container_type m_xtmp;
- const time_type* m_a;
- const time_type* m_b;
- const time_type* m_c;
-
- order_type m_q;
-
-
- private:
- void reset_iter(typename container_iterator_vector::iterator xiter_iter)
- {
- typename container_vector::iterator x_iter = m_xvec.begin();
- while( x_iter != m_xvec.end() ) {
- (*xiter_iter++) = traits_type::begin(*x_iter++);
- }
- }
-
- void check_consitency()
- {
- const time_type* a_iter = &m_a[0];
- const time_type* b_iter = &m_b[0];
- const time_type* c_iter = &m_c[0];
-
- unsigned short b_len = 1;
- // check 1: a_i = sum_j b_ij
- while( a_iter != &m_a[+m_q-1] ) {
- time_type tmp = 0.0;
- const time_type* b_end = b_iter + b_len;
- while( b_iter != b_end ) {
- tmp += *b_iter++;
- }
- b_len++;
- if( *a_iter++ != tmp )
- throw butcher_tableau_consistency_exception();
- }
-
- // check 2: sum_i c_i * (a_i)^(k-1) = 1/k k = 1..q
- for( unsigned short k = 1; k <= m_q; k++ ) {
- time_type tmp = 0.0;
- a_iter = &m_a[0];
- c_iter = &m_c[0];
- if( k == 1 ) // special treatment for 0^0 = 1
- tmp += *c_iter++;
- else
- c_iter++;
- while( a_iter != &m_a[m_q-1] ) {
- tmp += (*c_iter++) * pow( *a_iter++ , k-1 );
- }
- if( std::abs(time_type(tmp - time_type(1.0)/time_type(k))) >
- std::numeric_limits<time_type>::epsilon() ) {
- //std::clog << std::abs(time_type(tmp - time_type(1.0)/time_type(k))) << std::endl;
- throw butcher_tableau_order_condition_exception();
- }
- }
- }
-
-
- public:
-
- /* Constructor
-
- Same as for the stepper_rk_generic class, but with a, b and c provided
- as time_type*. The order q of the integration scheme is given separately
- as parameter. a is a pointer to an array of time_type[] of length q-1,
- b is of length 1 + 2 + ... q-1 = (q-1)(q-2)/2 and ordered as follows:
- b[0] = b_21, b[1] = b_31, b[2] = b_32, b[3] = b_41, b[4] = b42 ...
- c has length q.
-
- */
- stepper_rk_generic_ptr( const time_type* a,
- const time_type* b,
- const time_type* c,
- const unsigned short q)
- : m_a(a), m_b(b), m_c(c), m_q(q)
- {
- m_xvec.resize(m_q);
- m_xiter_vec.resize(m_q);
-
- check_consitency();
- }
-
- order_type order() const { return m_q; }
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- container_type &dxdt ,
- time_type t ,
- time_type dt )
- {
- using namespace detail::it_algebra;
- typename container_vector::iterator x_iter = m_xvec.begin();
- typename container_iterator_vector::iterator xiter_iter = m_xiter_vec.begin();
-
- (*x_iter) = dxdt;
- (*xiter_iter++) = traits_type::begin(*x_iter++);
-
- while( x_iter != m_xvec.end() )
- {
- traits_type::adjust_size(x, (*x_iter));
- (*xiter_iter++) = traits_type::begin(*x_iter++);
- }
- traits_type::adjust_size(x, m_xtmp);
-
- x_iter = m_xvec.begin()+1;
-
- const time_type* a_iter = &m_a[0];
- const time_type* b_iter = &m_b[0];
- unsigned short b_len= 1;
- while( x_iter != m_xvec.end() )
- {
- reset_iter(m_xiter_vec.begin());
- const time_type* b_end = b_iter + b_len;
- scale_sum_generic( traits_type::begin(m_xtmp), traits_type::end(m_xtmp),
- b_iter, b_end, dt,
- traits_type::begin(x), m_xiter_vec.begin() );
- system( m_xtmp, *x_iter++ , t + dt*(*a_iter++) );
- b_iter = b_end;
- b_len++;
- }
-
- reset_iter(m_xiter_vec.begin());
- scale_sum_generic( traits_type::begin(x), traits_type::end(x),
- &m_c[0], &m_c[m_q], dt,
- traits_type::begin(x), m_xiter_vec.begin() );
- }
-
-
-
- template< class DynamicalSystem >
- void do_step( DynamicalSystem &system ,
- container_type &x ,
- time_type t ,
- time_type dt )
- {
- traits_type::adjust_size(x, m_xtmp);
- system(x, m_xtmp, t);
- do_step( system, x, m_xtmp, t, dt);
- }
-
- };
-
-}
-}
-}
-
-#endif
Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -0,0 +1,33 @@
+/*
+ boost header: boost/numeric/odeint/steppers.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ 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_BOOST_NUMERIC_ODEINT_STEPPERS_HPP_INCLUDED
+#define BOOST_BOOST_NUMERIC_ODEINT_STEPPERS_HPP_INCLUDED
+
+// steppers
+#include <boost/numeric/odeint/steppers/stepper_euler.hpp>
+
+
+// error steppers
+#include <boost/numeric/odeint/steppers/stepper_half_step.hpp>
+
+
+// hamiltonian steppers
+#include <boost/numeric/odeint/steppers/hamiltonian_stepper_euler.hpp>
+
+
+// controlled steppers
+#include <boost/numeric/odeint/steppers/error_checker_standard.hpp>
+#include <boost/numeric/odeint/steppers/controlled_stepper_standard.hpp>
+
+
+#endif //BOOST_BOOST_NUMERIC_ODEINT_STEPPERS_HPP_INCLUDED
Copied: sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/controlled_stepper_standard.hpp (from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_standard.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/boost/numeric/odeint/controlled_stepper_standard.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/controlled_stepper_standard.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -18,10 +18,10 @@
#include <algorithm>
#include <complex>
-#include <boost/numeric/odeint/error_checker_standard.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
+#include <boost/numeric/odeint/steppers/error_checker_standard.hpp>
+#include <boost/numeric/odeint/container_traits/container_traits.hpp>
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/steppers/detail/iterator_algebra.hpp>
Copied: sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/hamiltonian_stepper_euler.hpp (from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/boost/numeric/odeint/hamiltonian_stepper_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/hamiltonian_stepper_euler.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -15,8 +15,8 @@
#include <stdexcept>
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
+#include <boost/numeric/odeint/steppers/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/container_traits/container_traits.hpp>
namespace boost {
namespace numeric {
@@ -49,7 +49,7 @@
public:
template< class CoordinateFunction >
- void do_step( CoordinateFunction &qfunc ,
+ void do_step( CoordinateFunction qfunc ,
container_type &q ,
container_type &p ,
time_type dt )
Copied: sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_euler.hpp (from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_euler.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_euler.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -18,11 +18,8 @@
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_EULER_HPP
#define BOOST_NUMERIC_ODEINT_STEPPER_EULER_HPP
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
-#include <boost/numeric/odeint/container_traits.hpp>
-
-// #include <boost/numeric/odeint/stepper_base.hpp>
-
+#include <boost/numeric/odeint/steppers/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/container_traits/container_traits.hpp>
namespace boost {
namespace numeric {
@@ -45,8 +42,6 @@
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
-// typedef typename traits_type::iterator iterator;
-// typedef typename traits_type::const_iterator const_iterator;
//
Copied: sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_half_step.hpp (from r61775, /sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_half_step.hpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper_half_step.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/steppers/stepper_half_step.hpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -16,7 +16,8 @@
#ifndef BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
#define BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
-#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/container_traits/container_traits.hpp>
+#include <boost/numeric/odeint/steppers/detail/iterator_algebra.hpp>
#include <iostream>
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -8,22 +8,10 @@
: requirements
<include>../../../..
<include>$BOOST_ROOT
- <include>$BLITZ_ROOT
- <include>$MTL4_ROOT
<define>BOOST_ALL_NO_LIB=1
: build-dir .
;
-exe harmonic_oscillator : harmonic_oscillator.cpp ;
-
-
-exe lorenz_cmp_rk4_rk_generic : lorenz_cmp_rk4_rk_generic.cpp ;
-exe lorenz_controlled : lorenz_controlled.cpp ;
-exe lorenz_integrate_constant_step : lorenz_integrate_constant_step.cpp ;
-exe lorenz_integrator : lorenz_integrator.cpp ;
-exe lorenz_stepper : lorenz_stepper.cpp ;
-exe pendulum_vibrating_pivot : pendulum_vibrating_pivot.cpp ;
-exe dnls : dnls.cpp ;
-
exe doc_harm_osc : doc_harm_osc.cpp ;
exe doc_integrate : doc_integrate.cpp ;
+exe solar_system : solar_system.cpp point_type.hpp ;
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/dnls.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/dnls.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,115 +0,0 @@
-/* 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;
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/doc_integrate.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/doc_integrate.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/doc_integrate.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -33,20 +33,21 @@
harm_osc harmonic_oscillator(0.15);
//[ define_const_stepper
- stepper_rk4< state_type > rk4;
+ stepper_euler< state_type > euler;
//]
//[ integrate_const
- integrate_const( rk4, harmonic_oscillator, x , 0.0, 10.0 , 0.01);
+ integrate_const( euler, harmonic_oscillator, x , 0.0, 10.0 , 0.01);
//]
//[ define_adapt_stepper
- stepper_rk5_ck< state_type > rk5;
+ stepper_half_step< stepper_euler< state_type > > half_stepper;
//]
//[ define_conntrolled_stepper
- controlled_stepper_standard< stepper_rk5_ck< state_type > >
+ controlled_stepper_standard<
+ stepper_half_step< stepper_euler< state_type > > >
controlled_rk5( 1E-6 , 1E-7 , 1.0 , 1.0 );
//]
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/harmonic_oscillator.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/harmonic_oscillator.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,294 +0,0 @@
-/* Boost numeric/odeint/examples/harm_osc.cpp
-
- Copyright 2009 A. Constantine Ashford
-
- Demonstrates the use of odeint to estimate the (known) solution to a simple,
- undriven harmonic oscillator:
-
- y'' = -y - gamma*y'
-
- This example uses the Euler solver as well as he fourth-order Runge Kutta.
- The output of each is printed in the output. The accuracy can be compared, or
- the program can be profiled for a speed comparison.
-
- The analytical solution to this ODE is (with initial displacement y and no
- initial velocity y'):
-
- Y = exp( -.5*gamma*t ) * cos( sqrt(1-gamma^2)*t )
-
- This example calculates the analytical solution using the standard math
- library and outputs the values at the same time points for comparison and
- plotting.
-
- On a unix system, a plot of the output could be created easily using gnuplot
- and the following commands:
-
- $ plot "harm_osc.out" using 1:2 title 'Euler' with dots,\
- $ "harm_osc.out" using 1:5 title 'Runge Kutta 4' with dots,\
- $ "harm_osc.out" using 1:8 title 'Exact' with dots
-
- The step size, number of steps to integrate, output file name, and gamma can
- be specified on the command line IFF you have the boost library files installed
- *which are at the same version as this example is compiled against*. You can
- enable the command line options by compiling with the preprocessor variable
- BOOST_PROGRAM_OPTIONS_INSTALLED defined.
-
-
- Compile with:
- g++ -Wall -O3 -march=k8 -I$BOOST_DIR -I../../../../ harmonic_oscillator.cpp
-
- 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 <fstream>
-#include <iomanip>
-#include <vector>
-#include <cmath>
-
-#include <boost/numeric/odeint.hpp>
-
-#ifdef BOOST_PROGRAM_OPTIONS_INSTALLED
-#include <boost/program_options.hpp>
-namespace opts=boost::program_options;
-#endif
-
-namespace odeint=boost::numeric::odeint;
-
-// Namespace for variables specific to this example
-namespace harm_osc
-{
- double gamma = .15;
- double dt = 0.01;
- size_t olen = 10000;
- std::ofstream fout;
-
- // The type of container used to hold the state vector
- typedef std::vector<double> state_type;
-
- /*
- void harmonic_oscillator(state_type &x, state_type &dxdt, double t)
-
- Calculates the derivatives of the state vector for the harmonic oscillator
- system. The system equation is described in the opening comments.
-
- The state vector is:
- x[0] = Y (the output)
- x[1] = Y' (dY/dt)
-
- Arguments:
- state_type &x: reference to the current state vector [ Y Y' ]
- state_type &dxdt: reference to the vector of calculated derivitives
- [ Y' Y'']
- double t: the current integration time
- */
- void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t)
- {
- dxdt[0] = x[1];
- dxdt[1] = -x[0] - gamma*x[1];
- }
-
- /*
- Optionally included function to parse the command line. See detailed
- comments and implementation at the bottom of this file
- */
- int parse_command_line(int ac, char ** av);
-}
-
-// Main
-int main(int argc, char **argv)
-{
- using namespace std;
-
- // Save the original destination of cout in case we redirect to a file
- // as a result of --filename command line option
- streambuf * cout_buf = cout.rdbuf();
-
- // Parse the command line options iff boost::program_options is installed
- #ifdef BOOST_PROGRAM_OPTIONS_INSTALLED
- int parse_ret = harm_osc::parse_command_line(argc, argv);
- if(parse_ret)
- return parse_ret;
- #endif
-
- // Declare state vectors for the Euler system, and the Runge Kutta system
- harm_osc::state_type x_euler(2);
- harm_osc::state_type x_rk4(2);
-
- // Set identical initial conditions for the Euler solver and the Runge Kutta
- x_euler[0] = 1.0;
- x_euler[1] = 0.0;
- x_rk4[0] = 1.0;
- x_rk4[1] = 0.0;
-
- // Initialize the solvers
- odeint::stepper_euler<harm_osc::state_type> euler;
- odeint::stepper_rk4<harm_osc::state_type> rk4;
-
- euler.adjust_size( x_euler );
- rk4.adjust_size( x_rk4 );
-
-
- // Write the output header. The '#' symbol creates a comment in gnuplot
- cout << "# Output from the simple harmonic oscillator example." << endl;
- cout << "#" << endl;
- cout << "# Columns:" << endl;
- cout << "#" << endl;
- cout << "# 1: time (s)" << endl;
- cout << "# 2: Euler output" << endl;
- cout << "# 3: Euler output dY/dt" << endl;
- cout << "# 4: Euler output error" << endl;
- cout << "# 5: 4th Order Runge Kutta output" << endl;
- cout << "# 6: 4th Order Runge Kutta dY/dt" << endl;
- cout << "# 7: 4th Order Runge Kutta output error" << endl;
- cout << "# 8: Analytical (calculated) solution" <<endl;
- cout << "#" << endl;
- cout << "# A. Constantine Ashford" << endl;
- cout << "#" << endl;
-
-// Initialize the time variable, and the analytical solution variable
- double t = 0.0;
- double reference;
- double e_error = 0, r_error = 0;
- double e_error_max = 0, r_error_max = 0;
- double e_error_rms = 0, r_error_rms = 0;
- int data_column_width = 15;
-
- // Integrate and write the approximate and analytical results to a table
- for(size_t oi=0; oi < harm_osc::olen; ++oi, t += harm_osc::dt)
- {
- // Calculate the analytical solution for this time
- reference = exp(-.5*harm_osc::gamma*t) *
- cos(sqrt(1-pow(harm_osc::gamma,2))*t);
-
- // Calculate absolute errors for the data table
- e_error = abs((reference - x_euler[0])/reference);
- r_error = abs((reference - x_rk4[0])/reference);
-
- // Write all the results to a row in the table
- cout << setw(5) << t
- << setw(data_column_width) << x_euler[0]
- << setw(data_column_width) << x_euler[1]
- << setw(data_column_width) << e_error
- << setw(data_column_width) << x_rk4[0]
- << setw(data_column_width) << x_rk4[1]
- << setw(data_column_width) << r_error
- << setw(data_column_width) << reference
- << endl;
-
- // Aggregate statistics
- if( (e_error = abs(e_error)) > e_error_max)
- e_error_max = e_error;
-
- if((r_error = abs(r_error)) > r_error_max)
- r_error_max = r_error;
-
- e_error_rms += pow(e_error, 2);
- r_error_rms += pow(r_error, 2);
-
- // Advance the solvers
- euler.do_step(harm_osc::harmonic_oscillator,
- x_euler, t, harm_osc::dt);
-
- rk4.do_step(harm_osc::harmonic_oscillator,
- x_rk4, t, harm_osc::dt);
-
- }
-
- e_error_rms = sqrt(e_error_rms/harm_osc::olen);
- r_error_rms = sqrt(r_error_rms/harm_osc::olen);
-
- // Restore cout to direct to standard out in case we redirected to a file
- cout.rdbuf(cout_buf);
- cout << "#Integration complete." << endl;
- cout << "#Errors for the Euler solver: MAX: " << e_error_max << " RMS: "
- << e_error_rms << endl;
- cout << "#Errors for the Runge Kutta solver: MAX: " << r_error_max << " RMS: "
- << r_error_rms << endl;
-
-
- // Success
- return 0;
-}
-
-#ifdef BOOST_PROGRAM_OPTIONS_INSTALLED
-/*
-parse_command_line(int ac, char ** av)
-
-Uses boost::program_options to parse the command line.
-Allows the example to use runtime options such as:
---step=.001 or
---step .001
-
-Arguments:
- int ac: the number of command line arguments. Usually just pass in argc
- from main.
-
- char ** av: an array of c-style strings holding the command line arguments.
- Usually just pass in argv from main.
-
-Return value: 1 if the user chose the --help option, and does *not* wish to
-run an integration. Otherwise 0.
-
-*/
-int harm_osc::parse_command_line(int ac, char ** av)
-{
- using namespace std;
-
- // Declare supported options
- opts::options_description opts_desc("Allowed options");
- opts_desc.add_options()
- ("help", "Display this message.")
- ("file", opts::value<string>(), "Output file for data.(Default: stdout)")
- ("step", opts::value<double>(), "Set step size (in seconds)")
- ("duration", opts::value<size_t>(), "Set integration duration (in steps)")
- ("gamma", opts::value<double>(), "Set gamma equation parameter");
-
- // Let Boost Program Options parse the command line
- opts::variables_map opts_map;
- opts::store(opts::parse_command_line(ac, av, opts_desc), opts_map);
- opts::notify(opts_map);
-
- // Print help message then exit
- if(opts_map.count("help"))
- {
- cout << opts_desc << endl;
- return 1;
- }
-
- // Alter step size in seconds
- if(opts_map.count("step"))
- {
- harm_osc::dt = opts_map["step"].as<double>();
- cout << "Step size set to " << harm_osc::dt << " seconds" << endl;
- }
-
- // Change integration duration (in steps)
- if(opts_map.count("duration"))
- {
- harm_osc::olen = opts_map["duration"].as<size_t>();
- cout << "Integration duration set to " << harm_osc::olen << " steps."
- << endl;
- }
-
- // Change damping constant
- if(opts_map.count("gamma"))
- {
- harm_osc::gamma = opts_map["gamma"].as<double>();
- cout << "Parameter GAMMA set to " << harm_osc::gamma << "." << endl;
- }
-
- // Redirect output from standard out to file
- if(opts_map.count("file"))
- {
- cout << "Writing output to file " << opts_map["file"].as<string>()
- << endl;
- fout.open(opts_map["file"].as<string>().c_str());
- cout.rdbuf(fout.rdbuf());
- }
-
- return 0;
-}
-#endif
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,151 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- Lorenz Eqquations:
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- with sigma = 10, r=28, b = 8/3
-
- 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 <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/stepper_rk_generic.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::numeric::odeint;
-
-
-typedef std::tr1::array< double , 3 > state_type;
-
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-int main( int argc , char **argv )
-{
- const double dt = 0.01;
- const size_t olen = 100000000;
-
- state_type x1 = {{1.0, 0.0, 0.0}};
- state_type x2 = {{1.0, 0.0, 0.0}};
- state_type x3 = {{1.0, 0.0, 0.0}};
- state_type x4 = {{1.0, 0.0, 0.0}};
-
- stepper_rk4< state_type > stepper_rk4;
- stepper_rk4_classical< state_type > stepper_rk4_classical;
-
- vector< double > a(3);
-
- a[0] = 0.5;
- a[1] = 0.5;
- a[2] = 1;
-
- vector< vector<double> > b(3);
- b[0].resize(1);
- b[1].resize(2);
- b[2].resize(3);
-
- b[0][0] = 0.5;
- b[1][0] = 0.0; b[1][1] = 0.5;
- b[2][0] = 0.0; b[2][1] = 0.0; b[2][2] = 1.0;
-
- vector< double > c(4);
- c[0] = 1.0/6.0;
- c[1] = 1.0/3.0;
- c[2] = 1.0/3.0;
- c[3] = 1.0/6.0;
-
- const double a2[3] = {0.5, 0.5, 1.0 };
- const double b2[6] = {0.5, 0.0, 0.5, 0.0, 0.0, 1.0};
- const double c2[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
-
- stepper_rk_generic< state_type > stepper_generic4(a, b, c);
-
- stepper_rk_generic_ptr< state_type > stepper_generic4_ptr(a2, b2, c2, 4);
-
- clock_t start , end;
- double t;
-
- cout.precision(16);
-
- cout << "Hand written Runge Kutta 4"<<endl;
-
- t = 0.0;
- start= clock();
- for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
- stepper_rk4.do_step( lorenz , x1 , t , dt );
- if( oi < 5 )
- cout << "x after step "<<oi<<": "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
-
- }
- end = clock();
- cout << "x after "<<olen<<" steps: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
- cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) <<"s"<< endl;
-
- cout << endl << "###########################" << endl << endl;
-
- cout << "Runge Kutta 4 via generic Runge Kutta implementation"<<endl;
- t = 0.0;
- start= clock();
- for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
- stepper_generic4.do_step( lorenz , x2 , t , dt );
- if( oi < 5 )
- cout << "x after step "<<oi<<": "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
- }
- end = clock();
- cout << "x after "<<olen<<" steps: "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
- cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
-
- cout << endl << "###########################" << endl << endl;
- cout << "Runge Kutta 4 via C-Array styled generic Runge Kutta implementation"<<endl;
-
- t = 0.0;
- start= clock();
- for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
- stepper_generic4_ptr.do_step( lorenz , x3 , t , dt );
- if( oi < 5 )
- cout << "x after step "<<oi<<": "<<x3[0]<<tab<<x3[1]<<tab<<x3[2]<<endl;
- }
- end = clock();
- cout << "x after "<<olen<<" steps: "<<x3[0]<<tab<<x3[1]<<tab<<x3[2]<<endl;
- cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
-
-
- cout << endl << "###########################" << endl << endl;
-
- cout << "Runge Kutta 4 with classical NR-style Runge Kutta implementation"<<endl;
- t = 0.0;
- start= clock();
- for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
- stepper_rk4_classical.do_step( lorenz , x4 , t , dt );
- if( oi < 5 )
- cout << "x after step "<<oi<<": "<<x4[0]<<tab<<x4[1]<<tab<<x4[2]<<endl;
- }
- end = clock();
- cout << "x after "<<olen<<" steps: "<<x4[0]<<tab<<x4[1]<<tab<<x4[2]<<endl;
- cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
-
- return 0;
-}
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_controlled.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_controlled.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,117 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_controlled.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- Shows the usage of odeint by integrating the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- mit sigma = 10, r=28, b = 8/3
-
- 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 <fstream>
-#include <vector>
-#include <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/container_traits_tr1_array.hpp>
- // #include <boost/numeric/odeint/controlled_stepper_bs.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace std::tr1;
-using namespace boost::numeric::odeint;
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-typedef array<double, 3> state_type;
-
-size_t function_calls = 0;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
- function_calls++;
-}
-
-
-class output_observer
-{
-
- ofstream m_file;
-
-public:
- output_observer( const char* file_name )
- : m_file( file_name )
- {
- m_file.precision(10);
- //m_file.setf(ios::fixed,ios::floatfield);
- }
-
- ~output_observer()
- {
- m_file.close();
- }
-
- void operator()( double t, state_type &x, ... )
- {
- m_file << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
- }
-};
-
-int main( int argc , char **argv )
-{
- const double end_time = 25.0;
-
- const double eps_abs = 1E-10;
- const double eps_rel = 1E-10;;
-
- state_type x;
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 20.0;
-
- stepper_rk5_ck< state_type > rk5;
- controlled_stepper_standard< stepper_rk5_ck< state_type > > controlled_rk5( eps_abs , eps_rel, 1.0, 1.0 );
- output_observer rk5_obs("lorenz_rk5.dat");
- size_t steps = integrate_adaptive( controlled_rk5, lorenz, x, 0.0, end_time, 1E-2, rk5_obs );
-
- clog << "RK5: " << steps << " steps. (" << function_calls << " function calls)" << endl;
-
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 20.0;
-
- function_calls = 0;
-
- controlled_stepper_bs< state_type > controlled_bs(eps_abs, eps_rel, 1.0, 1.0);
-
- output_observer bs_obs("lorenz_bs.dat");
- steps = integrate_adaptive( controlled_bs, lorenz, x, 0.0, end_time, 1E-2, bs_obs );
-
- clog << "BS: " << steps << " steps. (" << function_calls << " function calls)" << endl;
-
-
-
- return 0;
-}
-
-/*
- Compile with
- g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_controlled.cpp
-*/
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_integrate_constant_step.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,98 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_integrate_constant_step.cpp
-
- Copyright 2009 Karsten Ahnert
-
- Shows the usage of odeint, and integrates the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- mit sigma = 10, r=28, b = 8/3
-
- 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 <list>
-#include <tr1/array>
-
-#include <boost/lambda/lambda.hpp>
-#include <boost/lambda/bind.hpp>
-#include <boost/lambda/if.hpp>
-#include <boost/lambda/loops.hpp>
-#include <boost/lambda/switch.hpp>
-#include <boost/lambda/construct.hpp>
-#include <boost/lambda/casts.hpp>
-#include <boost/lambda/exceptions.hpp>
-#include <boost/lambda/numeric.hpp>
-#include <boost/lambda/algorithm.hpp>
-
-#include <boost/numeric/odeint.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::lambda;
-using namespace boost::numeric::odeint;
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-const double dt = 0.01;
-const size_t olen = 10000;
-
-//typedef std::tr1::array< double , 3 > state_type;
-typedef vector< double> state_type;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-int main( int argc , char **argv )
-{
- state_type x(3);
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 0.0;
-
- stepper_rk4< state_type > rk4;
- stepper_euler< state_type > euler;
- integrate_const( rk4 , lorenz , x , 0.0 , 100.0 , 0.01 ,
- cout << _1 << tab << _2[0] << "\n" );
-
- integrate_const_steps( rk4 , lorenz , x, 0.0 , 0.01 , 100 ,
- cout << _1 << tab << _2[0] << "\n" );
-
- integrate_const( rk4 , lorenz , x , 0.0 , 100.0, 0.01 );
- integrate_const_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 );
-
-/* integrate( euler , lorenz , 0.0 , 0.01 , x , 100.0 ,
- cout << _1 << tab << _2[0] << "\n" );*/
-
-/* vector<double> traj;
- back_insert_iterator< vector<double> > iter(traj);
- integrate( euler , lorenz , 0.0 , 0.01 , x , 1.0 , var(*iter++) = _2[1] );
- copy( traj.begin() , traj.end() ,
- ostream_iterator<double>( cout , "\n" ) ); */
-
-
-
-
- return 0;
-}
-
-
-
-/*
- Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz_integrate_constant_step.cpp
-*/
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_integrator.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_integrator.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,115 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_integrator.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- Shows the usage of odeint integrator by integrating the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- mit sigma = 10, r=28, b = 8/3
-
- 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 <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace std::tr1;
-using namespace boost::numeric::odeint;
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-const double eps_abs = 1E-6;
-const double eps_rel = 1E-7;
-
-const size_t time_points = 101;
-
-typedef array<double, 3> state_type;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-int main( int argc , char **argv )
-{
- state_type x1;
- x1[0] = 1.0;
- x1[1] = 0.0;
- x1[2] = 20.0;
- state_type x2(x1);
- state_type x3(x1);
-
-
- vector<state_type> x1_t_vec;
- vector<double> t1_vec;
- vector<state_type> x2_t_vec;
- vector<double> t2_vec;
- vector<state_type> x3_t_vec;
- vector<double> t3_vec;
-
- stepper_half_step< stepper_euler< state_type > > euler;
- controlled_stepper_standard< stepper_half_step< stepper_euler< state_type > > >
- euler_controlled( eps_abs, eps_rel, 1.0, 1.0);
- size_t steps = integrate( euler_controlled, lorenz, x1,
- 0.0, 10.0, 1E-4,
- back_inserter(t1_vec),
- back_inserter(x1_t_vec));
-
- clog << "Euler Half Step: #steps " << steps << endl;
-
- stepper_half_step< stepper_rk4< state_type > > rk4;
- controlled_stepper_standard< stepper_half_step< stepper_rk4< state_type > > >
- rk4_controlled( eps_abs, eps_rel, 1.0, 1.0);
- steps = integrate( rk4_controlled, lorenz, x2, 0.0, 10.0, 1E-4,
- back_inserter(t2_vec),
- back_inserter(x2_t_vec));
-
- clog << "RK4 Half Step: #steps " << steps << endl;
-
-
- stepper_rk5_ck< state_type > rk5;
- controlled_stepper_standard< stepper_rk5_ck< state_type > >
- rk5_controlled( eps_abs, eps_rel, 1.0, 1.0);
- steps = integrate( rk5_controlled, lorenz, x3, 0.0, 10.0, 1E-4,
- back_inserter(t3_vec),
- back_inserter(x3_t_vec));
-
- clog << "RK5 Cash-Karp: #steps: " << steps << endl;
-
- cout.precision(16);
- cout.setf(ios::fixed,ios::floatfield);
-
- cout << t1_vec.size() << tab << t2_vec.size() << tab << t3_vec.size() << endl;
-
-
- //cout << "current state: " ;
- cout << (x1_t_vec.back())[0] << tab << (x1_t_vec.back())[1] << tab << (x1_t_vec.back())[2] << tab;
- cout << x2_t_vec.back()[0] << tab << x2_t_vec.back()[1] << tab << x2_t_vec.back()[2] << tab;
- cout << x3_t_vec.back()[0] << tab << x3_t_vec.back()[1] << tab << x3_t_vec.back()[2] << endl;
-
- return 0;
-}
-
-/*
- Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz_integrator.cpp
-*/
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_stepper.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/lorenz_stepper.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,103 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_stepper.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- Shows the usage of odeint, and integrates the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- with sigma = 10, r=28, b = 8/3
-
- Furthermore, the usage of std::tr1::array and std::vector in odeint is
- shown and the performance of both containers is compared.
-
- 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 <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::numeric::odeint;
-
-
-typedef std::vector< double > state_type1;
-typedef std::tr1::array< double , 3 > state_type2;
-
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-void lorenz1( state_type1 &x , state_type1 &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-void lorenz2( state_type2 &x , state_type2 &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-
-
-
-int main( int argc , char **argv )
-{
- const double dt = 0.01;
- const size_t olen = 100000000;
-
- state_type1 x1(3);
- x1[0] = 1.0;
- x1[1] = 0.0;
- x1[2] = 0.0;
- state_type2 x2 = {{ 1.0 , 0.0 , 0.0 }};
-
- stepper_rk4< state_type1 > stepper1( x1 );
- stepper_rk4< state_type2 > stepper2( x2 );
-
- clock_t start , end;
- double t;
-
- start= clock();
- t = 0.0;
- for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- stepper1.do_step( lorenz1 , x1 , t , dt );
- end = clock();
- cout << "vector : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
- cout << "x: "<<x1[0]<<tab<<x1[1]<<tab<<x1[2]<<endl;
-
-
- start= clock();
- t = 0.0;
- for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- stepper2.do_step( lorenz2 , x2 , t , dt );
- end = clock();
- cout << "array : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
- cout << "x: "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
-
- return 0;
-}
-
-
-
-/*
- Compile with
- g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_stepper.cpp
-*/
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
+++ (empty file)
@@ -1,82 +0,0 @@
-/* Boost numeric/odeint/examples/coupled_vdp.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- 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
- copy at http://www.boost.org/LICENSE_1_0.txt)
-*/
-
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <list>
-#include <tr1/array>
-
-#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 alpha = 0.1;
-const double omega = 20;
-const double a = 0.1;
-
-/*
- Defines the right hand side f(x,t) of the dynamical equations dx/dt = f(x,t)
- x consists of x=(x, dx/dt) and f has explicit time dependence
-*/
-void my_system( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = x[1];
- dxdt[1] = -sin(x[0]) - alpha*x[1] + a*omega*omega*sin(omega*t)*cos(x[0]);
-}
-
-int main( int argc , char **argv )
-{
- state_type x = {{ 1.0, 0.0 }};
-
- vector<double> times;
- vector<state_type> x_t_vec;
-
- stepper_half_step< stepper_rk4< state_type > > stepper;
-
- controlled_stepper_standard
- < stepper_half_step< stepper_rk4< state_type > >
- > controlled_stepper( 1E-6 , 1E-7 , 1.0 , 1.0 );
-
- size_t steps = integrate( controlled_stepper, my_system, x,
- 0.0, 100.0,1E-4,
- back_inserter(times),
- back_inserter(x_t_vec));
-
- clog << "Steps: " << steps << endl;
-
- cout.precision(5);
- cout.setf(ios::fixed,ios::floatfield);
-
-
- for( size_t i=0; i<times.size(); i++ )
- {
- //cout << "current state: " ;
- cout << times[i] << tab;
- cout << x_t_vec[i][0] << tab << x_t_vec[i][1] << endl;
- }
-
- return 0;
-
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/solar_system.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/solar_system.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/solar_system.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -1,114 +1,42 @@
+/* Boost libs/numeric/odeint/examples/solar_system.cpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+
+ solar system example for Hamiltonian stepper
+
+ 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 <bits/stdc++.h>
#include <bits/stdtr1c++.h>
-#include <boost/operators.hpp>
#include <boost/circular_buffer.hpp>
+#include <boost/ref.hpp>
#include <boost/numeric/odeint.hpp>
+#include "point_type.hpp"
+
#define tab "\t"
using namespace std;
using namespace boost::numeric::odeint;
-template< class T >
-class point :
- boost::additive1< point< T > ,
- boost::additive2< point< T > , T ,
- boost::multiplicative2< point< T > , T
- > > >
-{
-public:
-
- typedef T value_type;
- typedef point< value_type > point_type;
-
- value_type x , y , z;
-
- point( void ) : x(0.0) , y(0.0) , z(0.0) { }
-
- point( value_type val ) : x(val) , y(val) , z(val) { }
-
- point( value_type _x , value_type _y , value_type _z )
- : x(_x) , y(_y) , z(_z) { }
-
-
- point_type& operator+=( const point_type& p )
- {
- x += p.x ; y += p.y ; z += p.z;
- return *this;
- }
-
- point_type& operator-=( const point_type& p )
- {
- x -= p.x ; y -= p.y ; z -= p.z;
- return *this;
- }
-
- point_type& operator+=( const value_type& val )
- {
- x += val ; y += val ; z += val;
- return *this;
- }
+// we simulate 5 planets and the sun
+const size_t n = 6;
- point_type& operator-=( const value_type& val )
- {
- x -= val ; y -= val ; z -= val;
- return *this;
- }
-
- point_type& operator*=( const value_type &val )
- {
- x *= val ; y *= val ; z *= val;
- return *this;
- }
-
- point_type& operator/=( const value_type &val )
- {
- x /= val ; y /= val ; z /= val;
- return *this;
- }
-};
-
-template< class T >
-T inner_product( const point< T > &p1 , const point< T > &p2 )
-{
- return p1.x*p2.x + p1.y*p2.y + p1.z*p2.z;
-}
-
-template< class T >
-T norm( const point< T > &p )
-{
- return inner_product( p , p );
-}
-
-template< class T >
-T abs( const point< T > &p )
-{
- return sqrt( norm( p ) );
-}
-
-template< class T >
-ostream& operator<<( ostream &out , const point< T > &p )
-{
- out << p.x << tab << p.y << tab << p.z;
- return out;
-}
-
-
-
-const size_t n = 3;
-typedef point< double > point_type;
+typedef point< double , 3 > point_type;
typedef std::tr1::array< point_type , n > state_type;
typedef std::tr1::array< double , n > mass_type;
+typedef hamiltonian_stepper_euler< state_type > stepper_type;
-typedef hamiltonian_stepper_rk< state_type > stepper_type;
-typedef boost::circular_buffer< point_type > buffer_type;
+const double gravitational_constant = 2.95912208286e-4;
-const mass_type masses = {{ 1.0e9 , 1.0e9 , 1.0e12}};
-const double gravitational_constant = 6.657e-20;
-
ostream& operator<<( ostream &out , const state_type &x )
{
typedef state_type::value_type value_type;
@@ -117,99 +45,121 @@
return out;
}
-point_type get_mean( const state_type &x )
-{
- point_type mean( 0.0 );
- if( x.empty() ) return mean;
- for( size_t i=0 ; i<x.size() ; ++i ) mean += x[i];
- mean /= double( x.size() );
- return mean;
-}
-point_type get_center_of_mass( const state_type &x ,
- const mass_type &m )
+point_type center_of_mass( const state_type &x , const mass_type &m )
{
- point_type mean( 0.0 );
- if( x.empty() ) return mean;
double overall_mass = 0.0;
+ point_type mean( 0.0 );
for( size_t i=0 ; i<x.size() ; ++i )
{
overall_mass += m[i];
mean += m[i] * x[i];
}
- mean /= overall_mass;
+ if( x.size() != 0 ) mean /= overall_mass;
return mean;
-
}
-void center_system( state_type &x , point_type mean )
-{
- for( size_t i=0 ; i<x.size() ; ++i ) x[i] -= mean;
-}
-
-void solar_system( state_type &q , state_type &dpdt )
+double energy( const state_type &q , const state_type &p ,
+ const mass_type &masses )
{
- point_type diff , tmp;
- fill( dpdt.begin() , dpdt.end() , 0.0 );
+ const size_t n = q.size();
+ double en = 0.0;
for( size_t i=0 ; i<n ; ++i )
{
- for( size_t j=i+1 ; j<n ; ++j )
- {
- diff = q[j] - q[i];
- tmp = gravitational_constant * diff / pow( abs( diff ) , 3.0 );
- dpdt[i] += tmp * masses[j];
- dpdt[j] -= tmp * masses[i];
- }
+ en += 0.5 * norm( p[i] ) / masses[i];
+ for( size_t j=0 ; j<i ; ++j )
+ {
+ double diff = abs( q[i] - q[j] );
+ en -= gravitational_constant * masses[j] * masses[i] / diff;
+ }
+ }
+ return en;
+}
+
+struct solar_system
+{
+ mass_type &m_masses;
+
+ solar_system( mass_type &masses ) : m_masses( masses ) { }
+
+ void operator()( const state_type &q , state_type &dpdt )
+ {
+ const size_t n = q.size();
+ fill( dpdt.begin() , dpdt.end() , 0.0 );
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ for( size_t j=i+1 ; j<n ; ++j )
+ {
+ point_type diff = q[j] - q[i];
+ double d = abs( diff );
+ diff = gravitational_constant * diff / d / d / d;
+ dpdt[i] += diff * m_masses[j];
+ dpdt[j] -= diff * m_masses[i];
+ }
+ }
}
-}
+};
int main( int argc , char **argv )
{
- state_type q , p;
- stepper_type stepper;
-
- fill( q.begin() , q.end() , 0.0 );
- fill( p.begin() , p.end() , 0.0 );
- q[0] = point_type( 0.0 , 1.0 , 0.0 );
- p[0] = point_type( 0.00001 , 0.0 , 0.0 );
- q[2] = point_type( 1.0 , 0.0 , 0.0 );
-
- center_system( q , get_center_of_mass( q , masses ) );
- center_system( p , get_center_of_mass( p , masses ) );
-
- const double dt = 1.0;
+ mass_type masses = {{
+ 1.00000597682 , // sun
+ 0.000954786104043 , // jupiter
+ 0.000285583733151 , // saturn
+ 0.0000437273164546 , // uranus
+ 0.0000517759138449 , // neptune
+ 1.0 / ( 1.3e8 ) // pluto
+ }};
+
+ state_type q = {{
+ point_type( 0.0 , 0.0 , 0.0 ) , // sun
+ point_type( -3.5023653 , -3.8169847 , -1.5507963 ) , // jupiter
+ point_type( 9.0755314 , -3.0458353 , -1.6483708 ) , // saturn
+ point_type( 8.3101420 , -16.2901086 , -7.2521278 ) , // uranus
+ point_type( 11.4707666 , -25.7294829 , -10.8169456 ) , // neptune
+ point_type( -15.5387357 , -25.2225594 , -3.1902382 ) // pluto
+ }};
+
+ state_type p = {{
+ point_type( 0.0 , 0.0 , 0.0 ) , // sun
+ point_type( 0.00565429 , -0.00412490 , -0.00190589 ) , // jupiter
+ point_type( 0.00168318 , 0.00483525 , 0.00192462 ) , // saturn
+ point_type( 0.00354178 , 0.00137102 , 0.00055029 ) , // uranus
+ point_type( 0.00288930 , 0.00114527 , 0.00039677 ) , // neptune
+ point_type( 0.00276725 , -0.00170702 , -0.00136504 ) // pluto
+ }};
+
+
+ point_type qmean = center_of_mass( q , masses );
+ point_type pmean = center_of_mass( p , masses );
+ for( size_t i=0 ; i<n ; ++i ) { q[i] -= qmean ; p[i] -= pmean; }
- buffer_type buffers[n];
- for( size_t i=0 ; i<n ; ++i ) buffers[i].set_capacity( 100 );
+ stepper_type stepper;
- cout << "unset key\n";
- const size_t ilen = 1000;
+ const double dt = 100.0;
double t = 0.0;
- while( true )
+ while( t < 10000000.0 )
{
- clog << get_mean( p ) << tab << get_mean( q ) << endl;
- for( size_t i=0 ; i<n ; ++i ) buffers[i].push_back( q[i] );
-
- cout << "p [-20:20][-20:20] '-' u 1:2 w l\n";
- for( size_t i=0 ; i<n ; ++i )
- {
- copy( buffers[i].begin() , buffers[i].end() ,
- ostream_iterator<point_type>( cout , "\n" ) );
- cout << "\n";
- }
- cout << "e" << endl;
-
-// getchar();
-
-// cout << "p [-2:2][-2:2] '-' u 1:2 pt 7 ps 5 \n" << q << "e" << endl;
-
- for( size_t ii=0 ; ii<ilen ; ++ii,t+=dt )
- {
- stepper.do_step( solar_system , q , p , dt );
- }
+ clog << t << tab << energy( q , p , masses ) << tab;
+ clog << center_of_mass( q , masses ) << tab;
+ clog << center_of_mass( p , masses ) << endl;
+
+ cout << t;
+ for( size_t i=0 ; i<n ; ++i ) cout << tab << q[i];
+ cout << endl;
+
+ for( size_t i=0 ; i<1 ; ++i,t+=dt )
+ stepper.do_step( solar_system( masses ) , q , p , dt );
+ t += dt;
}
return 0;
}
+
+
+/*
+Plot with gnuplot:
+p "solar_system.dat" u 2:4 w l,"solar_system.dat" u 5:7 w l,"solar_system.dat" u 8:10 w l,"solar_system.dat" u 11:13 w l,"solar_system.dat" u 14:16 w l,"solar_system.dat" u 17:19 w l
+*/
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp 2010-05-05 15:52:04 EDT (Wed, 05 May 2010)
@@ -15,15 +15,9 @@
#include <boost/test/unit_test.hpp>
-#include <boost/numeric/odeint/stepper_euler.hpp>
-#include <boost/numeric/odeint/stepper_half_step.hpp>
-#include <boost/numeric/odeint/stepper_midpoint.hpp>
-#include <boost/numeric/odeint/stepper_rk4_classical.hpp>
-#include <boost/numeric/odeint/stepper_rk4.hpp>
-#include <boost/numeric/odeint/stepper_rk5_ck.hpp>
-#include <boost/numeric/odeint/stepper_rk78_fehlberg.hpp>
-
-#include <boost/numeric/odeint/controlled_stepper_standard.hpp>
+#include <boost/numeric/odeint/steppers/stepper_euler.hpp>
+#include <boost/numeric/odeint/steppers/stepper_half_step.hpp>
+#include <boost/numeric/odeint/steppers/controlled_stepper_standard.hpp>
using namespace boost::unit_test;
using namespace boost::numeric::odeint;
@@ -149,6 +143,7 @@
check_error_stepper_concept( stepper , 1 , 2 );
}
+/*
void test_midpoint_concept()
{
stepper_midpoint< std::vector< double > > stepper;
@@ -182,16 +177,16 @@
check_stepper_concept( stepper , 8 );
check_error_stepper_concept( stepper , 7 , 8 );
}
+*/
void test_controlled_stepper_standard_concept()
{
- typedef stepper_rk5_ck< std::vector< double > > stepper_type;
+ typedef stepper_euler< std::vector< double > > stepper_type;
typedef controlled_stepper_standard< stepper_type > controlled_stepper_type;
controlled_stepper_type stepper( 1.0 , 1.0 , 1.0 , 1.0 );
check_controlled_stepper_concept( stepper );
-}
-
+ }
@@ -201,11 +196,11 @@
test->add( BOOST_TEST_CASE( &test_euler_concept ) );
test->add( BOOST_TEST_CASE( &test_half_step_euler_concept ) );
- test->add( BOOST_TEST_CASE( &test_midpoint_concept ) );
+/* test->add( BOOST_TEST_CASE( &test_midpoint_concept ) );
test->add( BOOST_TEST_CASE( &test_rk4_classical_concept ) );
test->add( BOOST_TEST_CASE( &test_rk4_concept ) );
test->add( BOOST_TEST_CASE( &test_rk5_ck_concept ) );
- test->add( BOOST_TEST_CASE( &test_rk78_fehlberg_concept ) );
+ test->add( BOOST_TEST_CASE( &test_rk78_fehlberg_concept ) );*/
test->add( BOOST_TEST_CASE( &test_controlled_stepper_standard_concept ) );
return test;