|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r68477 - in sandbox/odeint/branches/karsten: . libs/numeric/odeint/doc libs/numeric/odeint/doc/html libs/numeric/odeint/examples libs/numeric/odeint/ideas/butcher libs/numeric/odeint/ideas/generic_stepper libs/numeric/odeint/ideas/rosenbrock4 libs/numeric/odeint/ideas/units
From: karsten.ahnert_at_[hidden]
Date: 2011-01-27 12:11:07
Author: karsten
Date: 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
New Revision: 68477
URL: http://svn.boost.org/trac/boost/changeset/68477
Log:
* clean up and preparing the documentation
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/concepts.qbk (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/extend.qbk (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference.qbk (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference_integrate_functions.qbk (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference_steppers.qbk (contents, props changed)
Removed:
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/container_traits.qbk
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/integrate_functions.qbk
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/steppers.qbk
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Makefile
Text files modified:
sandbox/odeint/branches/karsten/Jamroot | 9 -
sandbox/odeint/branches/karsten/TODO | 6 +
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/Jamfile | 24 ++++
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk | 9 -
sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial.qbk | 167 ++++++++++++++++++++++++++++++++-------
sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/Jamfile | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp | 34 ++++----
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp | 28 +++---
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/units/Jamfile | 1
12 files changed, 201 insertions(+), 81 deletions(-)
Modified: sandbox/odeint/branches/karsten/Jamroot
==============================================================================
--- sandbox/odeint/branches/karsten/Jamroot (original)
+++ sandbox/odeint/branches/karsten/Jamroot 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -14,6 +14,7 @@
<include>$(BOOST_ROOT) ;
+# tests, regression tests and examples
build-project libs/numeric/odeint/test ;
build-project libs/numeric/odeint/examples ;
build-project libs/numeric/odeint/regression_test ;
@@ -25,17 +26,15 @@
build-project libs/numeric/odeint/test_external/gsl ;
-
-
# ideas
-# build-project libs/numeric/odeint/ideas/butcher ;
-# build-project libs/numeric/odeint/ideas/generic_stepper ;
+build-project libs/numeric/odeint/ideas/butcher ;
+build-project libs/numeric/odeint/ideas/generic_stepper ;
build-project libs/numeric/odeint/ideas/rosenbrock4 ;
build-project libs/numeric/odeint/ideas/units ;
# docs:
-# build-project libs/numeric/odeint/doc ;
+build-project libs/numeric/odeint/doc ;
Modified: sandbox/odeint/branches/karsten/TODO
==============================================================================
--- sandbox/odeint/branches/karsten/TODO (original)
+++ sandbox/odeint/branches/karsten/TODO 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -5,6 +5,7 @@
* test vector_space_algebra, maybe with some proto lib
* test copying
* test gsl
+ * test explicit stepper with ranges
OK * change standard_operations::rel_error in order to word with units and test it
OK * include implicit euler
OK * call via std::pair< deriv , jacobi >
@@ -30,10 +31,15 @@
* split resizing and copy/destruct/construct in different files
OK * change resizing concept, in order to word within the implicit steppers
OK * in all tests and regression test do not include odeint.hpp, only include the headers which are really needed
+* start new doc or cleanup the old project
* Integrate functions
* skript for setting the include defines according to the position in file system an writing a general copyright comment at the beginning
+* Documentation
+
+
+
* Adaptoren
* GMP
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/Jamfile 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -1,14 +1,30 @@
+# (C) Copyright 2010 : Karsten Ahnert, 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 quickbook ;
+project
+ : requirements
+ ;
+
+
xml odeint :
odeint.qbk
;
-boostbook standalone :
+boostbook standalone
+ :
odeint
+ :
+ <xsl:param>boost.root=../../../..
+ <xsl:param>generate.section.toc.level=1
+ <xsl:param>chunk.section.depth=1
+ <xsl:param>chunk.first.sections=1
;
-install html :
- /boost//doc/html/boostbook.css
- ;
+# install html :
+# /boost//doc/src/boostbook.css
+# ;
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/concepts.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/concepts.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -0,0 +1,255 @@
+[section Concepts]
+
+The odeint library defines three concepts for stepping objects.
+
+[/
+
+[section Dynamical system]
+
+The dynamical system represents the right hand side of the differential equation:
+
+[* x'(t) = f(x(t) , t)]
+
+In this odeint library the dynamical system is realized by a callable object,
+e.g. a function pointer, a functor or an instance of an object with an
+appropriate `()`-operator.
+The parameter structure has to be the following:
+
+`(const container_type &x, container_type &dxdt, const time_type t)`
+
+`x` is the current state of the system and `t` is the current time.
+The result *x' = f(x,t)* is then returned in `dxdt`.
+See the tutorial for examples on how to define the dynamical system.
+
+[endsect]
+
+]
+
+
+
+[section Basic stepper]
+
+Basic steppers execute one timestep of a specific order with a given stepsize.
+They usually allocate internal memory to store intermediate function call
+results.
+If state types with variable size are used (e.g. `vector`), it has to be assured
+that the stepper gets informed about any change of the state size by calling its
+`adjust_size` method.
+
+[*Associated Types]
+[table
+ [[] [] [Description]]
+ [[Time] [Stepper::time_type] [Type of the time variable, e.g. `double`]]
+ [[Container] [Stepper::container_type] [Type of the system state, e.g. `vector<double>`]]
+ [[Value] [Stepper::value_type] [Value type of the state, e.g. `double`]]
+ [[Order Type] [Stepper::order_type] [Type of the order parameter, usually `unsigned short`]]
+]
+
+[*Methods]
+
+*`Stepper()` Constructor.
+
+*`Stepper( container_type &x )` Constructor that allocates internal memory to
+store intermediate results of the same size as `x`.
+
+*`void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt )`
+
+Executes one timestep with the given parameters:
+[table
+ [[Parameter] [Type] [Description]]
+ [[system] [DynamicalSystem] [Function (callable object) that computes the rhs of the ode]]
+ [[x] [container_type] [The current state of the system *x(t)*]]
+ [[t] [time_type] [The current time *t*]]
+ [[dt] [time_type] [Length of the timestep to be executed]]
+]
+The result of this method is the (approximate) state of the system *x(t+dt)* and
+is stored in the variable `x` (in-place).
+Note, that the time `t` is not automatically increased by this method.
+
+*`void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ const container_type &dxdt ,
+ time_type t ,
+ time_type dt )`
+
+The same as above but with the additional parameter `dxdt` that represents the
+derivative [*x'(t) = f(x,t)] at the time *t*.
+
+*`void adjust_size( const container_type &x )`
+Adjusts the internal memory to store intermediate results of the same size as
+`x`.
+This function /must/ be called whenever the system size changes during the
+integration.
+
+*`order_type order_step()`
+Returns the order of the algorithm. If *n* is the order of a method, then the
+result of one iteration with the timestep *dt* is accurate up to *dt^n*. That
+means the error made by the time discretization is of order *dt^(n+1)*.
+
+[*Stepper that model this concept]
+
+*`stepper_euler`
+*`stepper_rk4`
+*`stepper_rk78_fehlberg`
+
+[endsect]
+
+[section Error stepper]
+
+Error steppers execute one timestep of a specific order with a given stepsize.
+Additionally, an error estimation of the obtained result is computed that can
+be used to control the error introduced by the time discretization.
+Like the basic steppers, error steppers usually allocate internal memory to store
+intermediate function call results. If state types with variable size are used
+(e.g. `vector`), it has to be assured that the stepper gets informed about any
+change of the state size by calling its `adjust_size` method.
+
+[*Associated Types]
+
+Same as for /basic steppers/ above.
+
+[*Methods]
+
+*`Error_Stepper()` Constructor.
+
+*`Error_Stepper( container_type &x )` Constructor that allocates internal memory to
+store intermediate results of the same size as `x`.
+
+*`void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt ,
+ container_type &xerr)`
+
+Executes one timestep with the given parameters:
+[table
+ [[Parameter] [Type] [Description]]
+ [[system] [DynamicalSystem] [Function (callable object) that computes the rhs of the ode]]
+ [[x] [container_type] [The current state of the system *x(t)*]]
+ [[t] [time_type] [The current time *t*]]
+ [[dt] [time_type] [Length of the timestep to be executed]]
+ [[xerr] [container_type] [Used by the method to return the error estimation of this computation]]
+]
+The result of this method is the (approximate) state of the system *x(t+dt)*,
+which is returned in the variable `x` (in-place), and the corresponding error
+estimation returned in `xerr`.
+Note, that the time `t` is not automatically increased by this method.
+
+*`void do_step( DynamicalSystem &system ,
+ container_type &x ,
+ const container_type &dxdt ,
+ time_type t ,
+ time_type dt ,
+ container_type &xerr)`
+
+The same as above but with the additional parameter `dxdt` that represents the
+derivative [*x'(t) = f(x,t)] at the time *t*.
+
+*`void adjust_size( const container_type &x )`
+Adjusts the internal memory to store intermediate results of the same size as
+`x`.
+This function /must/ be called whenever the system size changes during the
+integration.
+
+*`order_type order_error_step()`
+Returns the order of the result *x(t+dt)* of the algorithm.
+
+*`order_type order_error()`
+Returns the order of the error estimation of the algorithm.
+
+[*Stepper that model this concept]
+
+*`stepper_rk5_ck`
+*`stepper_rk78_fehlberg`
+*`stepper_half_step`
+
+[endsect]
+
+[section Controlled stepper]
+
+Controlled steppers try to execute a timestep with a given error threshold.
+If the estimated error of the obtained solution is too big, the result is
+rejected and a new stepsize is proposed.
+If the error is small enough the timestep is accepted and possibly an increased
+stepsize is proposed.
+
+[*Associated Types]
+
+Same as for /basic steppers/ above.
+
+[*Methods]
+
+*`Controlled_Stepper( time_type abs_err, time_type rel_err,
+ time_type factor_x, time_type factor_dxdt )`
+
+Constructor that initializes the controlled stepper with several parameters of
+the error control. The controlled stepper assures that the error done by each
+individual timestep yields:
+
+[* xerr < 1.1 ( eps_abs + eps_rel * (factor_x |x| + factor_dxdt h |x'|) ) ]
+
+The factor 1.1 is for safety to avoid unnecessary many stepsize adjustings.
+The above inequality should be understand to hold for /all/ components of the
+possibly more dimensional vectors *x*, *x'* and *xerr*.
+If the estimated error is too large, a reduced stepsize will be suggested.
+If the estimated error is less than half of the desired error, an increased
+stepsize will be suggested.
+
+*`Controlled_Stepper( container_type &x, time_type abs_err, time_type rel_err,
+ time_type factor_x, time_type factor_dxdt )`
+Same as above, but with additional allocation of the internal memory to store
+intermediate results of the same size as `x`.
+
+*`controlled_step_result try_step(
+ DynamicalSystem &system,
+ container_type &x,
+ time_type &t,
+ time_type &dt )`
+
+Tries one timestep with the given parameters
+[table
+ [[Parameter] [Type] [Description]]
+ [[system] [DynamicalSystem] [Function (callable object) that computes the rhs of the ode]]
+ [[x] [container_type] [The current state of the system *x(t)*]]
+ [[t] [time_type] [The current time *t*]]
+ [[dt] [time_type] [Length of the timestep to be executed]]
+]
+This method has three possible outcomes represented by the returned value `result`:
+If `result = success` the step has been applied and x contains the new state
+*x(t)* where the time has also been increased *t += dt*.
+If `result = step_size_increased` the step has also been accomplished, but the
+estimated error was so small that a new stepsize is proposed in the variable
+`dt`.
+If `result = step_size_decreased` the step has been rejected due to a too big
+error. `x` and `t` remain unchanged and `dt` now containes the suggested reduced
+stepsize that should give an error below the desired level.
+
+*`controlled_step_result try_step(
+ DynamicalSystem &system,
+ container_type &x,
+ const container_type &dxdt,
+ time_type &t,
+ time_type &dt )`
+Same as above but with the additional parameter `dxdt` that that represents the
+derivative *x'(t) = f(x,t)* at the time *t*.
+
+*`void adjust_size( const container_type &x )`
+Adjusts the internal memory to store intermediate results of the same size as
+`x`.
+This function /must/ be called whenever the system size changes during the
+integration.
+
+*`order_type order_error_step()`
+Returns the order of the result *x(t+dt)* of the algorithm.
+
+[*Stepper that model this concept]
+
+*`controlled_stepper_standard`
+*`controlled_stepper_bs`
+
+[endsect]
+
+[endsect]
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/container_traits.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/container_traits.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
+++ (empty file)
@@ -1,9 +0,0 @@
-[section Container traits]
-
-definitions
-
-[section Define your own traits]
-
-[endsect]
-
-[endsect]
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/extend.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/extend.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -0,0 +1,19 @@
+[section Extend odeint]
+
+[section Write own steppers]
+
+[endsect]
+
+[section Adapt your own containers]
+
+gsl_vector, gsl_matrix, ublas::matrix, blitz::matrix, thrust
+
+[endsect]
+
+[section Adapt your own operations]
+
+gsl_complex, complex, thrust
+
+[endsect]
+
+[endsect]
\ No newline at end of file
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/integrate_functions.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/integrate_functions.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
+++ (empty file)
@@ -1,11 +0,0 @@
-[section Integration functions]
-
-[section Constant step-size functions]
-
-[endsect]
-
-[section Adaptive step-size functions]
-
-[endsect]
-
-[endsect]
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -96,12 +96,11 @@
[endsect]
-[import ../examples/doc_integrate.cpp]
-
[include tutorial.qbk]
-[include steppers.qbk]
+[include extend.qbk]
+
+[include concepts.qbk]
-[include integrate_functions.qbk]
+[include reference.qbk]
-[include container_traits.qbk]
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -0,0 +1,32 @@
+[section Reference]
+
+
+[include reference_steppers.qbk]
+
+[include reference_integrate_functions.qbk]
+
+[section Algebras]
+
+[endsect]
+
+[section Operations]
+
+[endsect]
+
+[section Resizing]
+
+[endsect]
+
+
+
+
+
+
+
+
+
+
+
+
+
+[endsect]
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference_integrate_functions.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference_integrate_functions.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -0,0 +1,11 @@
+[section Integration functions]
+
+[section Constant step-size functions]
+
+[endsect]
+
+[section Adaptive step-size functions]
+
+[endsect]
+
+[endsect]
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference_steppers.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/reference_steppers.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -0,0 +1,118 @@
+[section Stepper classes]
+
+[table Stepper Algorithms
+ [[Method] [Class] [Order] [Error Estimation]]
+ [[Euler] [stepper_euler] [1] [No]]
+ [[Runge-Kutta 4] [stepper_rk4] [4] [No]]
+ [[Runge-Kutta Cash-Karp] [stepper_rk5_ck] [5] [Yes (Order 4)]]
+ [[Runge-Kutta Fehlberg] [stepper_rk78_fehlberg] [7] [Yes (Order 8)]]
+ [[Midpoint] [stepper_midpoint] [variable] [No]]
+ [[Bulirsch-Stoer] [controlled_stepper_bs] [variable] [Controlled]]
+]
+
+[/
+
+[section stepper_euler]
+
+[*Concept:] [link boost_sandbox_numeric_odeint.stepper.concepts.basic_stepper Basic stepper]
+
+[*Accuracy:] The produced solution is accurate up to order 1.
+
+[*Complexity:] This method uses 1 function evaluation.
+
+[*Template Parameters:]
+[template tpl_parameter_table_[]
+[table
+ [[Parameter] [Default Value] [Description]]
+ [[class Container] [no default value] [Type of container that stores the state of the system]]
+ [[class Time] [double] [Type of the time variable in the ode]]
+ [[class Traits] [container_traits< Container >] [container_traits used by the stepper]]
+]]
+[tpl_parameter_table_]
+
+The Euler algorithm used in this class performs one integration
+step according to the formula:
+[: x(t+dt) = x(t) + dt*f(x,t)]
+The Euler stepper is the most simple algorithm to integrate a
+differential equation. It's only purpose in this library is for educational
+reasons - use more sophisticated steppers for real life problems.
+
+[endsect]
+
+[section stepper_rk4]
+
+[*Concept:] [link boost_sandbox_numeric_odeint.stepper.concepts.basic_stepper Basic stepper]
+
+[*Accuracy:] The produced solution is accurate up to order 4.
+
+[*Complexity:] This method uses 4 function evaluations.
+
+[*Template Parameters:]
+
+[tpl_parameter_table_]
+
+The Runge-Kutta 4 algorithm is /the/ classical method to integrate odes.
+The [@http://en.wikipedia.org/wiki/Runge%e2%80%93Kutta_methods#Explicit_Runge.E2.80.93Kutta_methods
+Butcher Tableau] for this method is as follows:
+
+[pre
+Butcher Tableau for Runge-Kutta 4
+ 0 |
+ 1/2 | 1/2
+ 1/2 | 0 1/2
+ 1 | 0 0 1
+ ------------------------
+ | 1/6 1/3 1/3 1/6
+]
+This method produces fast, though not too accurate solutions. Use it for quick
+preliminary results and then switch to error controlled methods like Cash-Karp
+for more reliable calculations.
+
+[endsect]
+
+[section stepper_rk5_ck]
+
+[*Concept:] [link boost_sandbox_numeric_odeint.stepper.concepts.error_stepper Error stepper]
+
+[*Accuracy:] The produced solution is accurate up to order 5 and the estimated
+error is also of order 5.
+
+[*Complexity:] This method uses 6 function evaluation.
+
+[*Template Parameters:]
+
+[tpl_parameter_table_]
+
+The Runge-Kutta method of order 5 with Cash-Karp coefficients is a good trade-off
+between numerical effort and obtained accuracy with error estimation.
+It is the all-round method and suitable for most problems.
+See [@http://en.wikipedia.org/wiki/Cash%e2%80%93Karp_method wikipedia] for details and
+the Butcher tableau.
+
+[endsect]
+
+[section stepper_rk78_fehlberg]
+
+[*Concept:] [link boost_sandbox_numeric_odeint.stepper.concepts.error_stepper Error Stepper]
+
+[*Accuracy:] The produced solution is accurate up to order 7 and the estimated
+error is of order 8.
+
+[*Complexity:] This method uses 13 function evaluations.
+
+[*Template Parameters:]
+
+[tpl_parameter_table_]
+
+The Runge-Kutta 78 Fehlberg method is a high-order algorithm that produces very
+good accuracy but for the cost of high numerical effort.
+Whenever extra-ordinarily high precision is required, you can fall back to this
+method.
+
+[endsect]
+
+
+]
+
+[endsect]
+
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/steppers.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/steppers.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
+++ (empty file)
@@ -1,45 +0,0 @@
-[section Stepper]
-
-
-[section Concepts]
-
-; huhu
-
-[section Basic stepper]
-
-definitions
-
-which class models this concept
-
-[endsect]
-
-[section Error stepper]
-
-definitions
-
-which class models this concept
-
-
-[endsect]
-
-[section Controlled stepper]
-
-definitions
-
-which class models this concept
-
-
-[endsect]
-
-[endsect]
-
-
-[section Stepper classes]
-
-overview table
-
-all stepper classes
-
-[endsect]
-
-[endsect]
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial.qbk (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial.qbk 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -7,13 +7,12 @@
system [*x]. Mathematically, this usually is an n-dimensional vector with
real numbers or complex numbers as scalar objects. For odeint the most natural
way is to use `vector< double >` or `vector< complex< double > >` to represent
-the system state. However, odeint can deal with other container types than
-vector as well, e.g. `tr1/array< double , N >` as long as it is able to obtain
+the system state. However, odeint can deal with other container types as well, e.g. `tr1/array< double , N >` as long as it is able to obtain
a ForwardIterator going through all of the container's elements. The scalar type
-must have several operators ( +, -, +=, -= ) and the `abs()`-fcuntion defined .
+must have several operations ( + , - , += , -= ) and the `abs()`-function defined.
Furthermore, one can choose the datatype of the time (that is, the parameter to
which respect the differentiation is done). The standard type for this is
-`double`, but one might be able to choose, for example, `complex< double >` as
+`double`, but it should be possible to use, for example, `complex< double >` as
well (untested). It must be possible to multiply the time type and the scalar
type of the vector. For the most cases `vector< double >` as state type and the
standard `double` for the time type should be sufficient.
@@ -41,7 +40,8 @@
[endsect]
[section Stepper Types]
-Numerical integration works iteratevly, that means you start at a state ['x(t)]
+
+Numerical integration works iteratively, that means you start at a state ['x(t)]
and performs a timestep of length ['dt] to obtain the approximate state
['x(t+dt)]. There exist many different methods to perform such a timestep each
of which has a certain order ['q]. If the order of a method is ['q] than it is
@@ -63,6 +63,7 @@
[endsect]
[section Integration with Constant Step Size]
+
The basic stepper just performs one timestep and doesn't give you any
information about the error that was made (except that you know it is of order
q+1). Such steppers are used with constant step size that should be chosen small
@@ -86,6 +87,7 @@
[endsect]
[section Integration with Adaptive Step Size]
+
To improve the numerical results and additionally minimize the computational
effort, the application of a step size control is advisable.
Step size control is realized via stepper algorithms that additionally provide an
@@ -102,7 +104,7 @@
The usual way to create a controlled stepper is via the
`make_controlled_stepper_standard( ErrorStepper , eps_abs , eps_rel , a_x , a_dxdt )`
function that takes an error stepper as parameter and four values defining the maximal
-absolute and relative error allowed for on integration step.
+absolute and relative error allowed for one integration step.
The standard controlled stepper created by this method ensures that the error ['err]
of the solution fulfills
['err < eps_abs + eps_rel * ( a_x * |x| + a_dxdt * dt * |dxdt| ) ]
@@ -147,35 +149,39 @@
solar system each planet, and of course also the sun will be represented by
mass points. The interaction force between each object is the gravitational
force which can be written as
-Fij = - gamma m M / (qi-qj)3 * (qi-qj)
-where gamma is the gravitational constant, mi and mj are the masses and qi and
-qj are the locations of the two objects.
+F_ij = -gamma m_i m_j (q_i-q_j)/|q_i-q_j|^3
-dqi = pi
-dpi = 1/m sum ji Fij
+where gamma is the gravitational constant, m_i and m_j are the masses and q_i
+and q_j are the locations of the two objects. The equations of motion are then
+
+dq_i/dt = p_i
+dp_i/dt = 1/m_i sum_ji F_ij
where pi is the momenta of object i. The equations of motion can also be
derived from the Hamiltonian
-H = sum over i pi^2 / 2 + sum j V( qi , qj )
+H = sum_i p_i^2/(2m_i) + sum_j V( qi , qj )
+
+with the interaction potential V(q_i,q_j). The Hamiltonian equations give the
+equations of motion
-via dqi = dH / dpi, dpi = - dH / dq1. V(qi,qj) is the interaction
-potential.
+dq_i/dt = dH/dp_i
+
+dp_i = -dH/dq_i.
In time independent Hamiltonian system the energy is conserved and special
integration methods have to be applied in order to ensure energy
-conservation. The odeint library provides two stepper classes for Hamiltonian
-systems, which are separable and can be written in the form H = sum pi^2/2 +
-Hq.
+conservation. The odeint library provides classes for Hamiltonian
+systems, which are separable and can be written in the form H = sum p_i^2/2m_i +
+Hq(q), where Hq(q) only depends on the coordinates.
hamiltonian_stepper_euler
hamiltonian_stepper_rk
Alltough this functional form might look a bit arbitrary it covers nearly all
classical mechanical systems with inertia and without dissipation, or where
-the equations of motion can be written in the form dqi=pi dpi=f(qi).
-
+the equations of motion can be written in the form dqi=mi pi dpi=f(qi).
[endsect]
@@ -186,36 +192,135 @@
as well as the velocity. Therefore, we use the operators from
<boost/operator.hpp>:
-show the code
+[import ../examples/point_type.hpp]
+[point_type]
+
+
+The next step is to define a container type storing the values of q and p and
+to define systems (derivative) functions. As container type we use
+std::tr1::array and the state type is then simply
+
+[import ../examples/solar_system.cpp]
+[state_type_definition]
+and represents all space coordinates q or all momenta coordinates p. As system
+function we have to provide f(p) = dq and f(q) = -dp, which acts only on p or q:
-The next step is to define the state type and the system (derivative)
-function. As state type we use std::tr1::array and a state type represents all
-space coordinates q or all momenta coordinates p. As system function we have
-to provide f(q)
-show the code
+[momentum_function]
-Note, that we have allready define the masses of all planets in the solar
-system.
+[coordinate_function]
In general a three body-system is chaotic, hence we can not expect that
arbitray initial conditions of the system will lead to a dynamic which is
comparable with the solar system. That is we have to define proper initial
-conditions.
-
-show the code
+conditions, which are taken from the book of Hairer, Wannier, Lubich.
Now, we use the rk stepper to integrate the solar system. To visualize the
motion we save the trajectory of each planet in a circular buffer. The output
can be piped directly into gnuplot and a very nice visualization of the motion
appears.
+[integration_solar_system]
+
+[endsect]
+
+[endsect]
+
+
+
+
+
+[section Stiff systems]
+
+blah blah
+
+[endsect]
+
+[section Lyapunov exponents]
+
+blah blah
+
+[endsect]
+
+[section Using boost::units]
+
+blah blah
+
+[endsect]
+
+[section Using Cuda and Thrust]
+
+blah blah
+
+[endsect]
+
+[section Using matrices as state types]
+
+Expanding resizing
+
[endsect]
-usage of the steppers
+
+[section Special topics]
+
+[section Pass by value or by reference]
+
+blah blah
[endsect]
+[section Using boost::range]
+
+blah blah
+
+[endsect]
+
+[endsect]
+
+
+
+
+
+
+
+
+
+
+
+
+[section References]
+
+[*General informations about numerical integration of ordinary differential equations:]
+
+[1] Press William H et al., Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd ed. (Cambridge University Press, 2007).
+
+[2] Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed. (Springer, Berlin, 2009).
+
+[3] Ernst Hairer and Gerhard Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed. (Springer, Berlin, 2010).
+
+
+[*Symplectic integration of numerical integration:]
+
+[4] Ernst Hairer, Gerhard Wanner, and Christian Lubich, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed. (Springer-Verlag Gmbh, 2006).
+
+[5] Leimkuhler Benedict and Reich Sebastian, Simulating Hamiltonian Dynamics (Cambridge University Press, 2005).
+
+
+[*Special symplectic methods:]
+
+[6] Haruo Yoshida, “Construction of higher order symplectic integrators,†Physics Letters A 150, no. 5 (November 12, 1990): 262-268.
+
+[7] Robert I. McLachlan, “On the numerical integration of ordinary differential equations by symmetric composition methods,†SIAM J. Sci. Comput. 16, no. 1 (1995): 151-168.
+
+[endsect]
+
+
+
+
+
+
+
+
[endsect]
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 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -7,7 +7,6 @@
project
: requirements
<include>../../../..
- <include>$(BOOST_ROOT)
<define>BOOST_ALL_NO_LIB=1
: build-dir .
;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -13,7 +13,6 @@
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
- <include>$(BOOST_ROOT)
<include>$(CHRONO_ROOT)
;
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Makefile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Makefile 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
+++ (empty file)
@@ -1,22 +0,0 @@
-CC = g++
-
-RATIO_ROOT = $(HOME)/boost/chrono
-
-CXXFLAGS = -O3 -I$(BOOST_ROOT) -I$(ODEINT_ROOT) -I$(RATIO_ROOT)
-
-HEADERS = algebra.hpp convert_value.hpp diagnostic.hpp explicit_runge_kutta.hpp predefined_steppers.hpp
-
-all : butcher_performance butcher_lorenz butcher_vdp
-
-butcher_performance : butcher_performance.o
-butcher_performance.o : butcher_performance.cpp $(HEADERS)
-
-butcher_vdp : butcher_vdp.o
-butcher_vdp.o : butcher_vdp.cpp $(HEADERS)
-
-butcher_lorenz : butcher_lorenz.o
-butcher_lorenz.o : butcher_lorenz.cpp $(HEADERS)
-
-
-clean :
- -rm *.o *~ butcher_performance butcher_vdp butcher_lorenz
\ No newline at end of file
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -11,8 +11,8 @@
#include <boost/mpl/int.hpp>
#include <boost/mpl/at.hpp>
-#include <boost/numeric/odeint/algebra/standard_algebra.hpp>
-#include <boost/numeric/odeint/algebra/standard_operations.hpp>
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include "convert_value.hpp"
@@ -29,7 +29,7 @@
struct algebra< state_type , coef_type , mpl::int_< 0 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -44,8 +44,8 @@
// *first1++ = *first2++ + a1 * *first3++ ;
- std_algebra::for_each3( x_tmp , x , k_vector[0] ,
- std_op::scale_sum2( 1.0 , a1 ) );
+ std_algebra::for_each3()( x_tmp , x , k_vector[0] ,
+ std_op::scale_sum2<>( 1.0 , a1 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
@@ -59,7 +59,7 @@
struct algebra< state_type , coef_type , mpl::int_< 1 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -74,8 +74,8 @@
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++;
- std_algebra::for_each4( x_tmp , x , k_vector[0] , k_vector[1] ,
- std_op::scale_sum3( 1.0 , a1 , a2 ) );
+ std_algebra::for_each4()( x_tmp , x , k_vector[0] , k_vector[1] ,
+ std_op::scale_sum3<>( 1.0 , a1 , a2 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
@@ -90,7 +90,7 @@
struct algebra< state_type , coef_type , mpl::int_< 2 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -106,8 +106,8 @@
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++;
- std_algebra::for_each5( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
- std_op::scale_sum4( 1.0 , a1 , a2 , a3 ) );
+ std_algebra::for_each5()( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
+ std_op::scale_sum4<>( 1.0 , a1 , a2 , a3 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
@@ -123,7 +123,7 @@
struct algebra< state_type , coef_type , mpl::int_< 3 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -141,8 +141,8 @@
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++ + a4 * *first6++;
- std_algebra::for_each6( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
- std_op::scale_sum5( 1.0 , a1 , a2 , a3 , a4 ) );
+ std_algebra::for_each6()( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
+ std_op::scale_sum5<>( 1.0 , a1 , a2 , a3 , a4 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
@@ -159,7 +159,7 @@
struct algebra< state_type , coef_type , mpl::int_< 4 > >
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -178,8 +178,8 @@
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++ + a4 * *first6++ + a5 * *first7++;
- std_algebra::for_each7( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] , k_vector[4] ,
- std_op::scale_sum6( 1.0 , a1 , a2 , a3 , a4 , a5 ) );
+ std_algebra::for_each7()( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] , k_vector[4] ,
+ std_op::scale_sum6<>( 1.0 , a1 , a2 , a3 , a4 , a5 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -15,7 +15,6 @@
<include>../../../../..
<include>../butcher
<include>$(CHRONO_ROOT)
- <include>$(BOOST_ROOT)
;
#exe test
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -12,15 +12,15 @@
#include <boost/array.hpp>
-#include <boost/numeric/odeint/algebra/standard_algebra.hpp>
-#include <boost/numeric/odeint/algebra/standard_operations.hpp>
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
template< size_t n >
struct fusion_algebra
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
template< class state_type >
inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , n > &a ,
@@ -34,13 +34,13 @@
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
template< class state_type >
inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 1 > &a ,
const state_type *k_vector , const double dt )
{
- std_algebra::for_each3( x_tmp , x , k_vector[0] , std_op::scale_sum2( 1.0 , a[0]*dt ) );
+ std_algebra::for_each3()( x_tmp , x , k_vector[0] , std_op::scale_sum2<>( 1.0 , a[0]*dt ) );
}
};
@@ -51,14 +51,14 @@
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
template< class state_type >
inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 2 > &a ,
const state_type *k_vector , const double dt )
{
- std_algebra::for_each4( x_tmp , x , k_vector[0] , k_vector[1] ,
- std_op::scale_sum3( 1.0 , a[0]*dt , a[1]*dt ) );
+ std_algebra::for_each4()( x_tmp , x , k_vector[0] , k_vector[1] ,
+ std_op::scale_sum3<>( 1.0 , a[0]*dt , a[1]*dt ) );
}
};
@@ -69,14 +69,14 @@
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
template< class state_type >
inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 3 > &a ,
const state_type *k_vector , const double dt )
{
- std_algebra::for_each5( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
- std_op::scale_sum4( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt) );
+ std_algebra::for_each5()( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
+ std_op::scale_sum4<>( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt) );
}
};
@@ -88,14 +88,14 @@
{
typedef boost::numeric::odeint::range_algebra std_algebra;
- typedef boost::numeric::odeint::default_operations< double > std_op;
+ typedef boost::numeric::odeint::default_operations std_op;
template< class state_type >
inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 4 > &a ,
const state_type *k_vector , const double dt )
{
- std_algebra::for_each6( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
- std_op::scale_sum5( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt , a[3]*dt ) );
+ std_algebra::for_each6()( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
+ std_op::scale_sum5<>( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt , a[3]*dt ) );
}
};
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/rosenbrock4/Jamfile 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -13,7 +13,6 @@
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
- <include>$(BOOST_ROOT)
;
exe rosenbrock4
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/units/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/units/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/units/Jamfile 2011-01-27 12:10:58 EST (Thu, 27 Jan 2011)
@@ -13,7 +13,6 @@
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
- <include>$(BOOST_ROOT)
;
exe test_units
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk