[Boost-commit] svn:boost r68530 - sandbox/odeint/branches/karsten/libs/numeric/odeint/doc

Subject: [Boost-commit] svn:boost r68530 - sandbox/odeint/branches/karsten/libs/numeric/odeint/doc
From: karsten.ahnert_at_[hidden]
Date: 2011-01-28 07:18:17

Author: karsten
Date: 2011-01-28 07:18:15 EST (Fri, 28 Jan 2011)
New Revision: 68530
URL: http://svn.boost.org/trac/boost/changeset/68530

Log:
* tutorial qbks 
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_harmonic_oscillator.qbk   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk   (contents, props changed)

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_harmonic_oscillator.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_harmonic_oscillator.qbk	2011-01-28 07:18:15 EST (Fri, 28 Jan 2011)
@@ -0,0 +1,126 @@
+[section Harmonic oscillator]
+
+[section Define the ODE]
+First of all, you have to specify the datatype that represents a state of your
+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 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 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 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.
+
+To integrate a differential equation numerically, one has to define the rhs of
+the equation ['x' = f(x)]. In odeint you supply this function in terms of
+an object that implements the ()-operator with a certain parameter structure.
+Hence, the straight forward way would be to just define a function, e.g:
+
+[rhs_function]
+
+The parameters of the function must follow the example above where [*x] is the
+current state, [*dxdt] is the derivative ['x'] and should be filled by the 
+function with ['f(x)] and [*t] is the current time.
+
+A more sophisticated approach is to implement the system as a class where the
+rhs function is defined as the ()-operator of the class with the same parameter
+structure as above:
+
+[rhs_class]
+
+odeint can deal with instances of such classes instead of pure functions which
+allows for cleaner code.
+
+[endsect]
+
+[section Stepper Types]
+
+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 
+accurate up to term ['~dt^q] that means the error in ['x] made by such a step 
+is ['~dt^(q+1)]. odeint provides several steppers of different orders from which
+you can choose. There are three types of steppers: [*Stepper], [*ErrorStepper]
+and [*ControlledStepper]. 
+
+[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]]
+]
+
+[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
+enough to have reasonable small errors. However, you should apply some sort of
+validity check of your results (such as observing conserved quantities) becasue
+you have no other control of the error. The following example defines a basic 
+stepper based on the classical Runge-Kutta scheme of 4th order.
+
+[define_const_stepper]
+
+The declaration of the stepper requires the state type as template parameter.
+The integration can now be done by using the `integrate_const( Stepper, System,
+ state, start_time, end_time, step_size )` function from odeint:
+
+[integrate_const]
+
+This call integrates the system defined by `harmonic_oscillator` using the rk4
+method from t=0 to 10 with a stepsize dt=0.01 and the initial condition given
+in `x`. The result, ['x(t=10)] is stored in `x` (in-place).
+
+[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 
+error estimation of the applied step.
+Odeint provides a number of such *ErrorSteppers* and we will show their usage on
+the example of stepper_rk5_ck -- a 5th order Runge-Kutta method with 4th order
+error estimation and coefficients introduced by Cash-Karp.
+
+[define_adapt_stepper]
+
+Given the error stepper, one still needs an instance that checks the error and 
+adjusts the step size accordingly.
+In odeint, this is done by *ControlledSteppers*.
+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 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| ) ]
+by decreasesing the step size. Note, that the stepsize is also increased if the error
+gets too small compared to the rhs of the above relation.
+Now we have everything needed to integrate the harmonic oscillator using an adaptive 
+step size method.
+Similar to the case with constant step size above, there exists a `integrate_adaptive` 
+function with a similar parameter structure, but it requires the controlled stepper 
+create by `make_controlled_stepper_standard`.
+
+[integrate_adapt]
+
+As above, this integrates the system defined by `harmonic_oscillator` using an adaptive
+step size method based on the rk5_ck scheme from t=0 to 10 with an initial step size of 
+dt=0.01 (will be adjusted) and the initial condition given in x. The result, x(t=10), will 
+also be stored in x (in-place).
+
+[endsect]
+
+[endsect]

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk	2011-01-28 07:18:15 EST (Fri, 28 Jan 2011)
@@ -0,0 +1,85 @@
+[section Solar system]
+
+[section Gravitation and energy conservation]
+
+The next example in this tutorial is a simulation of the solar system. In the
+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
+
+F_ij = -gamma m_i m_j (q_i-q_j)/|q_i-q_j|^3
+
+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_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
+
+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 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=mi pi dpi=f(qi).
+
+[endsect]
+
+
+[section Define the system function]
+
+To implement this system we define a point type which will represent the space
+as well as the velocity. Therefore, we use the operators from
+<boost/operator.hpp>:
+
+[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:
+
+
+[momentum_function]
+
+[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, 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]