mirror of
https://github.com/boostorg/odeint.git
synced 2025-05-09 23:24:01 +00:00
spell checking the docs and some minor changes
This commit is contained in:
Karsten Ahnert
parent
61b8f1a307
commit
e7fd6f3a27
@ -31,7 +31,7 @@
|
||||
</h5>
|
||||
<p>
|
||||
A controlled stepper following this Controlled Stepper concept provides the
|
||||
possibilty to perform one step of the solution <span class="emphasis"><em>x(t)</em></span>
|
||||
possibility to perform one step of the solution <span class="emphasis"><em>x(t)</em></span>
|
||||
of an ODE with step-size <span class="emphasis"><em>dt</em></span> to obtain <span class="emphasis"><em>x(t+dt)</em></span>
|
||||
with a given step-size <span class="emphasis"><em>dt</em></span>. Depending on an error estimate
|
||||
of the solution the step might be rejected and a smaller step-size is suggested.
|
||||
@ -135,7 +135,7 @@
|
||||
<td>
|
||||
<p>
|
||||
Tries one step of step size <code class="computeroutput"><span class="identifier">dt</span></code>.
|
||||
If the step was succesfull, <code class="computeroutput"><span class="identifier">success</span></code>
|
||||
If the step was successful, <code class="computeroutput"><span class="identifier">success</span></code>
|
||||
is returned, the resulting state is written to <code class="computeroutput"><span class="identifier">x</span></code>,
|
||||
the new time is stored in <code class="computeroutput"><span class="identifier">t</span></code>
|
||||
and <code class="computeroutput"><span class="identifier">dt</span></code> now contains
|
||||
|
||||
@ -21,7 +21,7 @@
|
||||
Output Stepper</a>
|
||||
</h3></div></div></div>
|
||||
<p>
|
||||
This concept specifies the interface a dense ouput stepper has to fulfill
|
||||
This concept specifies the interface a dense output stepper has to fulfill
|
||||
to be used within <a class="link" href="../odeint_in_detail/integrate_functions.html" title="Integrate functions">integrate
|
||||
functions</a>.
|
||||
</p>
|
||||
@ -30,7 +30,7 @@
|
||||
<span><a name="boost_sandbox_numeric_odeint.concepts.dense_output_stepper.description"></a></span><a class="link" href="dense_output_stepper.html#boost_sandbox_numeric_odeint.concepts.dense_output_stepper.description">Description</a>
|
||||
</h5>
|
||||
<p>
|
||||
A dense ouput stepper following this Dense Output Stepper concept provides
|
||||
A dense output stepper following this Dense Output Stepper concept provides
|
||||
the possibility to perform a single step of the solution <span class="emphasis"><em>x(t)</em></span>
|
||||
of an ODE to obtain <span class="emphasis"><em>x(t+dt)</em></span>. The step-size <code class="computeroutput"><span class="identifier">dt</span></code> might be adjusted automatically due
|
||||
to error control. Dense output steppers also can interpolate the solution
|
||||
|
||||
@ -119,8 +119,8 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Calls the observer which can do some analyzation of the state
|
||||
<span class="emphasis"><em>x</em></span> at time <span class="emphasis"><em>t</em></span>.
|
||||
Calls the observer which can do some analysis or output of the
|
||||
state <span class="emphasis"><em>x</em></span> at time <span class="emphasis"><em>t</em></span>.
|
||||
</p>
|
||||
</td>
|
||||
</tr></tbody>
|
||||
|
||||
@ -35,7 +35,7 @@
|
||||
</tr>
|
||||
<tr><td align="left" valign="top"><p>
|
||||
The following does not apply to implicit steppers like implicit_euler or
|
||||
rosenbrock 4 as there the <code class="computeroutput"><span class="identifier">state_type</span></code>
|
||||
Rosenbrock 4 as there the <code class="computeroutput"><span class="identifier">state_type</span></code>
|
||||
can not be changed from <code class="computeroutput"><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">vector</span></code>
|
||||
and no algebra/operations are used.
|
||||
</p></td></tr>
|
||||
@ -672,7 +672,7 @@
|
||||
<td>
|
||||
<p>
|
||||
Anything that defines operators + within itself and * with scalar
|
||||
(Mathematically spoken, anyhting that is a vector space).
|
||||
(Mathematically spoken, anything that is a vector space).
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
|
||||
@ -20,7 +20,7 @@
|
||||
<a name="boost_sandbox_numeric_odeint.concepts.stepper"></a><a class="link" href="stepper.html" title="Stepper">Stepper</a>
|
||||
</h3></div></div></div>
|
||||
<p>
|
||||
This concepts specifies the interface a simple stepper has to fullfill to
|
||||
This concepts specifies the interface a simple stepper has to fulfill to
|
||||
be used within the <a class="link" href="../odeint_in_detail/integrate_functions.html" title="Integrate functions">integrate
|
||||
functions</a>.
|
||||
</p>
|
||||
|
||||
@ -27,7 +27,7 @@
|
||||
The System concept models the algorithmic implementation of the rhs. of the
|
||||
ODE <span class="emphasis"><em>x' = f(x,t)</em></span>. The only requirement for this concept
|
||||
is that it should be callable with a specific parameter syntax (see below).
|
||||
A System is typically implemented as a function or a functor. Systems fullfilling
|
||||
A System is typically implemented as a function or a functor. Systems fulfilling
|
||||
this concept are required by all Runge-Kutta steppers as well as the Bulirsch-Stoer
|
||||
steppers. However, symplectic and implicit steppers work with other system
|
||||
concepts, see <a class="link" href="symplectic_system.html" title="Symplectic System">Symplectic
|
||||
|
||||
@ -33,7 +33,7 @@
|
||||
differential equations. Mathematically, these problems are formulated as
|
||||
follows: <span class="emphasis"><em>x'(t) = f(x,t)</em></span>, <span class="emphasis"><em>x(0) = x0</em></span>.
|
||||
<span class="emphasis"><em>x</em></span> and <span class="emphasis"><em>f</em></span> can be vectors and the
|
||||
solution is some function <span class="emphasis"><em>x(t)</em></span> fullfilling both equations
|
||||
solution is some function <span class="emphasis"><em>x(t)</em></span> fulfilling both equations
|
||||
above. In the following we will refer to <span class="emphasis"><em>x'(t)</em></span> also
|
||||
<code class="computeroutput"><span class="identifier">dxdt</span></code> which is also our notation
|
||||
for the derivative in the source code.
|
||||
@ -388,7 +388,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Standard method with error control and dense ouput, to be used
|
||||
Standard method with error control and dense output, to be used
|
||||
in controlled_error_stepper and in dense_output_controlled_explicit_fsal.
|
||||
</p>
|
||||
</td>
|
||||
@ -676,9 +676,9 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Dense ouput for <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
Dense output for <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
and <a class="link" href="../concepts/error_stepper.html" title="Error Stepper">Error
|
||||
Stepper</a> from above if they provide dense ouput functionality
|
||||
Stepper</a> from above if they provide dense output functionality
|
||||
(like <code class="computeroutput"><span class="identifier">euler</span></code> and
|
||||
<code class="computeroutput"><span class="identifier">runge_kutta_dopri5</span></code>).
|
||||
Order depends on the given stepper.
|
||||
@ -778,8 +778,8 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Stepper with step size and order control as well as dense ouput.
|
||||
Very good if high precision and dense ouput is required.
|
||||
Stepper with step size and order control as well as dense output.
|
||||
Very good if high precision and dense output is required.
|
||||
</p>
|
||||
</td>
|
||||
</tr>
|
||||
@ -936,12 +936,12 @@
|
||||
<tr>
|
||||
<td>
|
||||
<p>
|
||||
Dense Ouput Rosenbrock 4
|
||||
Dense Output Rosenbrock 4
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
<code class="computeroutput"><span class="identifier">rosenbrock4_dense_ouput</span></code>
|
||||
<code class="computeroutput"><span class="identifier">rosenbrock4_dense_output</span></code>
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
@ -1090,7 +1090,7 @@
|
||||
differential equations. Be sure, you will find them in every discipline.
|
||||
They also occur if partial differential equations (PDEs) are discretized
|
||||
in one coordinate. Then, a system of coupled ordinary differential occurs,
|
||||
sometimes also refered as lattices ODEs.
|
||||
sometimes also referred as lattices ODEs.
|
||||
</p>
|
||||
</div>
|
||||
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
|
||||
|
||||
@ -70,7 +70,7 @@
|
||||
For the integration itself we'll use the <code class="computeroutput">integrate</code>
|
||||
function, which is a convenient way to get quick results. It is based on
|
||||
the error-controlled <code class="computeroutput">runge_kutta_rk5_ck</code>
|
||||
stepper (5th order) and uses adaptive stepsize.
|
||||
stepper (5th order) and uses adaptive step-size.
|
||||
</p>
|
||||
<p>
|
||||
</p>
|
||||
@ -84,7 +84,7 @@
|
||||
above, the initial state <code class="computeroutput"><span class="identifier">x</span></code>,
|
||||
the start-and end-time of the integration as well as the initial time step.
|
||||
Note, that <code class="computeroutput">integrate</code>
|
||||
uses an adaptive stepsize during the integration steps so the time points
|
||||
uses an adaptive step-size during the integration steps so the time points
|
||||
will not be equally spaced. The integration returns the number of steps that
|
||||
were applied.
|
||||
</p>
|
||||
|
||||
@ -108,7 +108,7 @@
|
||||
<p>
|
||||
Of course, all controllers and dense-output steppers in odeint can be used
|
||||
with these mechanisms. In the table below you will find, which steppers is
|
||||
contructed from <code class="computeroutput"><span class="identifier">make_controlled</span></code>
|
||||
constructed from <code class="computeroutput"><span class="identifier">make_controlled</span></code>
|
||||
or <code class="computeroutput"><span class="identifier">make_dense_output</span></code> if applied
|
||||
on a stepper from odeint:
|
||||
</p>
|
||||
|
||||
@ -41,7 +41,7 @@
|
||||
observer calls</a>
|
||||
</h5>
|
||||
<p>
|
||||
If observer calls at equdistant time intervals <span class="emphasis"><em>dt</em></span> are
|
||||
If observer calls at equidistant time intervals <span class="emphasis"><em>dt</em></span> are
|
||||
needed, the <code class="computeroutput"><span class="identifier">integrate_const</span></code>
|
||||
function should be used.
|
||||
</p>
|
||||
@ -271,10 +271,10 @@
|
||||
<th align="left">Warning</th>
|
||||
</tr>
|
||||
<tr><td align="left" valign="top"><p>
|
||||
This function works only with steppers fullfilling the <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
This function works only with steppers fulfilling the <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
or <a class="link" href="../concepts/error_stepper.html" title="Error Stepper">Error
|
||||
Stepper</a> concept! Providing a <a class="link" href="../concepts/controlled_stepper.html" title="Controlled Stepper">Controlled
|
||||
Stepper</a> or __dense_out_stepper results in a compilationerror.
|
||||
Stepper</a> or __dense_out_stepper results in a compilation error.
|
||||
</p></td></tr>
|
||||
</table></div>
|
||||
<p>
|
||||
@ -301,7 +301,7 @@
|
||||
starting at <span class="emphasis"><em>x<sub>0</sub></em></span> and <span class="emphasis"><em>t<sub>0</sub></em></span>. If provided,
|
||||
<code class="computeroutput"><span class="identifier">observer</span></code> is called after
|
||||
every step and at the beginning with <code class="computeroutput"><span class="identifier">t0</span></code>.
|
||||
This function works only with steppers fullfilling the <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
This function works only with steppers fulfilling the <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
or <a class="link" href="../concepts/error_stepper.html" title="Error Stepper">Error
|
||||
Stepper</a> concept! It performs exactly <code class="computeroutput"><span class="identifier">n</span></code>
|
||||
steps with fixed step size <code class="computeroutput"><span class="identifier">dt</span></code>
|
||||
@ -319,7 +319,7 @@
|
||||
<p>
|
||||
Additionally to the sophisticated integrate function above odeint also provides
|
||||
a simple <code class="computeroutput"><span class="identifier">integrate</span></code> routine
|
||||
which uses a dense ouput stepper based on <code class="computeroutput"><span class="identifier">runge_kutta_dopri5</span></code>
|
||||
which uses a dense output stepper based on <code class="computeroutput"><span class="identifier">runge_kutta_dopri5</span></code>
|
||||
with standard error bounds <span class="emphasis"><em>10<sup>-6</sup></em></span> for the steps.
|
||||
</p>
|
||||
<p>
|
||||
|
||||
@ -33,7 +33,7 @@
|
||||
the user full control over the state type and the mathematical operations
|
||||
for this state type. Technically, this is done by introducing three concepts:
|
||||
StateType, Algebra, Operations. Most of the steppers in odeint expect three
|
||||
classes types fullfilling these concepts as template parameters. Note that
|
||||
classes types fulfilling these concepts as template parameters. Note that
|
||||
these concepts are not fully independent of each other but rather a valid
|
||||
combination must be provided in order to make the steppers work. In the following
|
||||
we will give some examples on reasonable state_type-algebra-operations combinations.
|
||||
@ -755,7 +755,7 @@
|
||||
<p>
|
||||
A similar class exists for the <code class="computeroutput"><span class="keyword">const</span></code>
|
||||
version of the iterator. Then we have a function returning the end iterator
|
||||
(similarily for <code class="computeroutput"><span class="keyword">const</span></code> again):
|
||||
(similarly for <code class="computeroutput"><span class="keyword">const</span></code> again):
|
||||
</p>
|
||||
<pre class="programlisting"><span class="identifier">gsl_vector_iterator</span> <span class="identifier">end_iterator</span><span class="special">(</span> <span class="identifier">gsl_vector</span> <span class="special">*</span><span class="identifier">x</span> <span class="special">)</span>
|
||||
<span class="special">{</span>
|
||||
@ -786,7 +786,7 @@
|
||||
Again with similar definitions for the <code class="computeroutput"><span class="keyword">const</span></code>
|
||||
versions. This eventually makes odeint work with gsl vectors as state
|
||||
types. The full code for these bindings is found in <a href="https://github.com/headmyshoulder/odeint-v2/tree/master/boost/numeric/odeint/external/gsl/gsl_wrapper.hpp" target="_top">gsl_wrapper.hpp</a>.
|
||||
It might look rather complicated but keep in mind that gsl is a pre compiled
|
||||
It might look rather complicated but keep in mind that gsl is a pre-compiled
|
||||
C library so it is quite an achievement that it's possible to get it
|
||||
working at all in the first place.
|
||||
</p>
|
||||
@ -934,7 +934,7 @@
|
||||
</tbody>
|
||||
</table></div>
|
||||
<p>
|
||||
Definign these operators makes your state type work with any basic Runge-Kutta
|
||||
Defining these operators makes your state type work with any basic Runge-Kutta
|
||||
stepper. However, if you want to use step-size control, some more functionality
|
||||
is required. Specifically, operations like <span class="emphasis"><em>max<sub>i</sub>( |err<sub>i</sub>| / (alpha
|
||||
* |s<sub>i</sub>|) )</em></span> have to be performed. <span class="emphasis"><em>err</em></span> and
|
||||
@ -992,7 +992,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Calculates the elementwise division 'x/y'
|
||||
Calculates the element-wise division 'x/y'
|
||||
</p>
|
||||
</td>
|
||||
</tr>
|
||||
@ -1074,7 +1074,7 @@
|
||||
<p>
|
||||
As an example for the employment of the <code class="computeroutput"><span class="identifier">vector_space_algebra</span></code>
|
||||
we will adopt <code class="computeroutput"><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">vector</span></code> from <a href="http://www.boost.org/doc/libs/release/libs/numeric/ublas/index.html" target="_top">Boost.UBlas</a>
|
||||
to work as a state type in odeint. This is particularily easy because
|
||||
to work as a state type in odeint. This is particularly easy because
|
||||
<code class="computeroutput"><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">vector</span></code> supports vector-vector addition
|
||||
and scalar-vector multiplication described above as well as <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">size</span></code>. It also has a resize member function
|
||||
so all that has to be done in this case is to declare resizability:
|
||||
@ -1230,7 +1230,7 @@
|
||||
use controlled steppers. For simple steppers definition of the simple
|
||||
<code class="computeroutput"><span class="special">+=</span></code> and <code class="computeroutput"><span class="special">*=</span></code>
|
||||
operators are sufficient. Having defined such a point type, we can easily
|
||||
perform the integration on a lorenz system by using the <code class="computeroutput"><span class="identifier">vector_space_algebra</span></code> again:
|
||||
perform the integration on a Lorenz system by using the <code class="computeroutput"><span class="identifier">vector_space_algebra</span></code> again:
|
||||
</p>
|
||||
<p>
|
||||
</p>
|
||||
|
||||
@ -40,7 +40,7 @@
|
||||
own steppers</a></span></dt>
|
||||
</dl></div>
|
||||
<p>
|
||||
Solving ordinary differential equation numerically is ususally done iteratively,
|
||||
Solving ordinary differential equation numerically is usually done iteratively,
|
||||
that is a given state of an ordinary differential equation is iterated forward
|
||||
<span class="emphasis"><em>x(t) -> x(t+dt) -> x(t+2dt)</em></span>. The steppers in odeint
|
||||
perform one single step. The most general stepper type is described by the
|
||||
@ -50,7 +50,7 @@
|
||||
we briefly present the mathematical and numerical details of the steppers.
|
||||
The <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
has two versions of the <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
method, one with an in-place transform of the currrent state and one with
|
||||
method, one with an in-place transform of the current state and one with
|
||||
an out-of-place transform:
|
||||
</p>
|
||||
<p>
|
||||
@ -90,7 +90,7 @@
|
||||
</p>
|
||||
<p>
|
||||
Other types of system function represent Hamiltonian systems or system which
|
||||
also compute the Jacobian needed in implicit steppers. For informations which
|
||||
also compute the Jacobian needed in implicit steppers. For information which
|
||||
stepper uses which system function see the stepper table below. It might
|
||||
be possible, that odeint will introduce new system types in near future.
|
||||
Since the system function is strongly related to the stepper type, such an
|
||||
@ -171,7 +171,7 @@
|
||||
This method also fills the derivative at time <span class="emphasis"><em>t+dt</em></span>
|
||||
into <code class="computeroutput"><span class="identifier">dxdtout</span></code>. Of course,
|
||||
the performance gain of such FSAL steppers only appears when combining
|
||||
with intergrate error estimation, like in the Runge-Kutta-Dopri5 stepper.
|
||||
with integrate error estimation, like in the Runge-Kutta-Dopri5 stepper.
|
||||
The FSAL-trick is sometimes also referred as the Fehlberg-Trick. An example
|
||||
how the FSAL steppers can be used is
|
||||
</p>
|
||||
@ -221,11 +221,11 @@
|
||||
A special class of symplectic systems are separable systems which can be
|
||||
written in the form <span class="emphasis"><em>dqdt/dt = f1(p)</em></span>, <span class="emphasis"><em>dpdt/dt
|
||||
= f2(q)</em></span>, where <span class="emphasis"><em>(q,p)</em></span> are the state of system.
|
||||
The space of <span class="emphasis"><em>(q,p)</em></span> is sometimes refered as the phase
|
||||
The space of <span class="emphasis"><em>(q,p)</em></span> is sometimes referred as the phase
|
||||
space and <span class="emphasis"><em>q</em></span> and <span class="emphasis"><em>p</em></span> are said the
|
||||
be the phase space variables. Symplectic systems in this special form occur
|
||||
widely in nature. For example the complete classical mechanics as written
|
||||
down by Newton, Lagrange and Hamilton can be forumulated in this framework.
|
||||
down by Newton, Lagrange and Hamilton can be formulated in this framework.
|
||||
Of course, the separability of the system depends on the specific choice
|
||||
of coordinates.
|
||||
</p>
|
||||
@ -233,7 +233,7 @@
|
||||
Integrable symplectic systems can be solved by odeint by means of the symplectic_euler
|
||||
stepper and a symplectic Runge-Kutta-Nystrom method of sixth-order. These
|
||||
stepper assume that the system is autonomous, hence the time will not explicitly
|
||||
occur. Furhter they fullfil in principle the default Stepper concept, but
|
||||
occur. Further they fulfill in principle the default Stepper concept, but
|
||||
they expect the system to be a pair of callable objects. The first entry
|
||||
of this pair calculates <span class="emphasis"><em>f1(p)</em></span> while the second calculates
|
||||
<span class="emphasis"><em>f2(q)</em></span>. The syntax is <code class="computeroutput"><span class="identifier">sys</span><span class="special">.</span><span class="identifier">first</span><span class="special">(</span><span class="identifier">p</span><span class="special">,</span><span class="identifier">dqdt</span><span class="special">)</span></code>
|
||||
@ -314,7 +314,7 @@
|
||||
form. This kind of system can be used in the symplectic solvers, by simply
|
||||
passing <span class="emphasis"><em>f(p)</em></span> to the <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
method, again <span class="emphasis"><em>f(p)</em></span> will be represented by a simple
|
||||
C-function or a functor. Here, the above example of the harmonic oscaillator
|
||||
C-function or a functor. Here, the above example of the harmonic oscillator
|
||||
can be written as
|
||||
</p>
|
||||
<p>
|
||||
@ -331,7 +331,7 @@
|
||||
is exactly the same function as in the above examples.
|
||||
</p>
|
||||
<p>
|
||||
Note, that the state of the ODE must not be contructed explicitly via
|
||||
Note, that the state of the ODE must not be constructed explicitly via
|
||||
<code class="computeroutput"><span class="identifier">pair</span><span class="special"><</span>
|
||||
<span class="identifier">vector_type</span> <span class="special">,</span>
|
||||
<span class="identifier">vector_type</span> <span class="special">></span>
|
||||
@ -369,7 +369,7 @@
|
||||
system possesses two or more time scales of very different order. Solvers
|
||||
for stiff systems are usually implicit, meaning that they solve equations
|
||||
like <span class="emphasis"><em>x(t+dt) = x(t) + dt * f(x(t+1))</em></span>. This particular
|
||||
scheme is the implicit euler method. Implicit methods usually solve the
|
||||
scheme is the implicit Euler method. Implicit methods usually solve the
|
||||
system of equations by a root finding algorithm like the Newton method
|
||||
and therefore need to know the Jacobian of the system <span class="emphasis"><em>J<sub>​ij</sub> = df<sub>​i</sub> /
|
||||
dx<sub>​j</sub></em></span>.
|
||||
@ -392,15 +392,15 @@
|
||||
<p>
|
||||
Another large class of solvers are multi-step method. They save a small
|
||||
part of the history of the solution and compute the next step with the
|
||||
help of this history. Since multistep methods know a part of their history
|
||||
help of this history. Since multi-step methods know a part of their history
|
||||
they do not need to compute the system function very often, usually it
|
||||
is only computed once. This makes multistep methods preferable if a call
|
||||
is only computed once. This makes multi-step methods preferable if a call
|
||||
of the system function is expensive. Examples are ODEs defined on networks,
|
||||
where the computation of the interaction is usually where expensive (and
|
||||
might be of order O(N^2)).
|
||||
</p>
|
||||
<p>
|
||||
Multistep methods differ from the normal steppers. They safe a part of
|
||||
Multi-step methods differ from the normal steppers. They safe a part of
|
||||
their history and this part has to be explicitly calculated and initialized.
|
||||
In the following example an Adams-Bashforth-stepper with a history of 5
|
||||
steps is instantiated and initialized;
|
||||
@ -429,7 +429,7 @@
|
||||
<p>
|
||||
</p>
|
||||
<p>
|
||||
Many multistep methods are also explicit steppers, hence the parameter
|
||||
Many multi-step methods are also explicit steppers, hence the parameter
|
||||
of <code class="computeroutput"><span class="identifier">do_step</span></code> method do not
|
||||
differ from the explicit steppers.
|
||||
</p>
|
||||
@ -439,7 +439,7 @@
|
||||
<th align="left">Caution</th>
|
||||
</tr>
|
||||
<tr><td align="left" valign="top"><p>
|
||||
The multistep methods have some internal variables which depend on the
|
||||
The multi-step methods have some internal variables which depend on the
|
||||
explicit solution. Hence you can not exchange the system of the state
|
||||
between two consecutive calls of <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
since then the internal variable do not correspond with the ODE and the
|
||||
@ -456,7 +456,7 @@
|
||||
</h4></div></div></div>
|
||||
<p>
|
||||
Many of the above introduced steppers possess the possibility to use adaptive
|
||||
stepsize control. Adaptive step size integration works in principle as
|
||||
step-size control. Adaptive step size integration works in principle as
|
||||
follows:
|
||||
</p>
|
||||
<div class="orderedlist"><ol class="orderedlist" type="1">
|
||||
@ -469,8 +469,8 @@
|
||||
</li>
|
||||
<li class="listitem">
|
||||
This error is compared against some predefined error tolerances. Are
|
||||
the tolerance violated the step is reject and the stepsize is decreases.
|
||||
Otherwise the step is accepted and possibly the stepsize is increased.
|
||||
the tolerance violated the step is reject and the step-size is decreases.
|
||||
Otherwise the step is accepted and possibly the step-size is increased.
|
||||
</li>
|
||||
</ol></div>
|
||||
<p>
|
||||
@ -478,8 +478,8 @@
|
||||
Stepper</a> concept. They are usually constructed from the underlying
|
||||
error steppers. An example is the controller for the explicit Runge-Kutta
|
||||
steppers. The Runge-Kutta steppers enter the controller as a template argument.
|
||||
Additionally one can pass the Runge-Kutta stepper to the contructor, but
|
||||
this step is not neccessary; the stepper is default-constructed if possible.
|
||||
Additionally one can pass the Runge-Kutta stepper to the constructor, but
|
||||
this step is not necessary; the stepper is default-constructed if possible.
|
||||
</p>
|
||||
<p>
|
||||
Different step size controlling mechanism exist. They all have in common
|
||||
@ -491,7 +491,7 @@
|
||||
in the integration functions.
|
||||
</p>
|
||||
<p>
|
||||
A classical way to decide wether a step is rejected or accepted is
|
||||
A classical way to decide whether a step is rejected or accepted is
|
||||
</p>
|
||||
<p>
|
||||
<span class="emphasis"><em>val = || | err<sub>​i</sub> | / ( ε<sub>​abs</sub> + ε<sub>​rel</sub> * ( a<sub>​x</sub> | x<sub>​i</sub> | + a<sub>​dxdt</sub> | | dxdt<sub>​i</sub> | )||
|
||||
@ -663,8 +663,9 @@
|
||||
<p>
|
||||
To ease to generation of the controlled stepper generation functions exist
|
||||
which take the absolute and relative error tolerances and a predefined
|
||||
error stepper and contruct from this knowledge an appropirate controlled
|
||||
stepper. The generation functions are explained in detail in XYZ.
|
||||
error stepper and construct from this knowledge an appropriate controlled
|
||||
stepper. The generation functions are explained in detail in <a class="link" href="generation_functions.html" title="Generation functions">Generation
|
||||
functions</a>.
|
||||
</p>
|
||||
</div>
|
||||
<div class="section">
|
||||
@ -695,11 +696,11 @@
|
||||
<p>
|
||||
Dense output stepper have their own concept. The main difference is that
|
||||
they control the state by them-self. If you call <code class="computeroutput"><span class="identifier">do_step</span></code>,
|
||||
only the ODE is passed as argument. Furhtermore <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
only the ODE is passed as argument. Furthermore <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
return the last time interval, hence you interpolate the solution between
|
||||
these two times points. Another difference is that they must be initialized
|
||||
with <code class="computeroutput"><span class="identifier">initialize</span></code>, otherwise
|
||||
the internal state of the stepper is default contructed which might produce
|
||||
the internal state of the stepper is default constructed which might produce
|
||||
funny errors or bugs.
|
||||
</p>
|
||||
<p>
|
||||
@ -716,8 +717,8 @@
|
||||
<p>
|
||||
</p>
|
||||
<p>
|
||||
Of course, this statement is also lengthly; it demonstrates how <code class="computeroutput"><span class="identifier">make_dense_output</span></code> can be used with the
|
||||
<code class="computeroutput"><span class="identifier">result_of</span></code> protocoll. The
|
||||
Of course, this statement is also lengthy; it demonstrates how <code class="computeroutput"><span class="identifier">make_dense_output</span></code> can be used with the
|
||||
<code class="computeroutput"><span class="identifier">result_of</span></code> protocol. The
|
||||
parameters to <code class="computeroutput"><span class="identifier">make_dense_output</span></code>
|
||||
are the absolute error tolerance, the relative error tolerance and the
|
||||
stepper. This explicitly assumes that the underlying stepper is a controlled
|
||||
@ -739,8 +740,8 @@
|
||||
steppers</a>
|
||||
</h4></div></div></div>
|
||||
<p>
|
||||
This section contains some general informations about the usage of the
|
||||
steppers in odeint.
|
||||
This section contains some general information about the usage of the steppers
|
||||
in odeint.
|
||||
</p>
|
||||
<p>
|
||||
<span class="bold"><strong>Steppers are copied by value</strong></span>
|
||||
@ -763,13 +764,13 @@
|
||||
</p></td></tr>
|
||||
</table></div>
|
||||
<p>
|
||||
Some steppers require to store some informations about the state of the
|
||||
ODE between two steps. Examples are the multistep method which store a
|
||||
Some steppers require to store some information about the state of the
|
||||
ODE between two steps. Examples are the multi-step method which store a
|
||||
part of the solution during the evolution of the ODE, or the FSAL steppers
|
||||
which store the last derivative at time <span class="emphasis"><em>t+dt</em></span>, to be
|
||||
used in the next step. In both cases the steppers expect that consecutive
|
||||
calls of <code class="computeroutput"><span class="identifier">do_step</span></code> are from
|
||||
the same solution and the same ODE. In this case it is absolutely neccessary
|
||||
the same solution and the same ODE. In this case it is absolutely necessary
|
||||
that you call <code class="computeroutput"><span class="identifier">do_step</span></code> with
|
||||
the same system function and the same state, see also the examples for
|
||||
the FSAL steppers above.
|
||||
@ -789,10 +790,10 @@
|
||||
</p>
|
||||
<p>
|
||||
All these stepper have in common, that they initially fill their internal
|
||||
state by themself. Hence you are not required to call initialize. See how
|
||||
this works for the Adams-Bashforth-Moulton stepper: in the example we instantiate
|
||||
a fourth order Adams-Bashforth-Moulton stepper, meaning that it will store
|
||||
4 internal derivatives of the solution at times <code class="computeroutput"><span class="special">(</span><span class="identifier">t</span><span class="special">-</span><span class="identifier">dt</span><span class="special">,</span><span class="identifier">t</span><span class="special">-</span><span class="number">2d</span><span class="identifier">t</span><span class="special">,</span><span class="identifier">t</span><span class="special">-</span><span class="number">3d</span><span class="identifier">t</span><span class="special">,</span><span class="identifier">t</span><span class="special">-</span><span class="number">4d</span><span class="identifier">t</span><span class="special">)</span></code>.
|
||||
state by themselves. Hence you are not required to call initialize. See
|
||||
how this works for the Adams-Bashforth-Moulton stepper: in the example
|
||||
we instantiate a fourth order Adams-Bashforth-Moulton stepper, meaning
|
||||
that it will store 4 internal derivatives of the solution at times <code class="computeroutput"><span class="special">(</span><span class="identifier">t</span><span class="special">-</span><span class="identifier">dt</span><span class="special">,</span><span class="identifier">t</span><span class="special">-</span><span class="number">2d</span><span class="identifier">t</span><span class="special">,</span><span class="identifier">t</span><span class="special">-</span><span class="number">3d</span><span class="identifier">t</span><span class="special">,</span><span class="identifier">t</span><span class="special">-</span><span class="number">4d</span><span class="identifier">t</span><span class="special">)</span></code>.
|
||||
</p>
|
||||
<p>
|
||||
</p>
|
||||
@ -826,7 +827,7 @@
|
||||
all stepper have a method <code class="computeroutput"><span class="identifier">adjust_size</span></code>
|
||||
which takes a parameter representing a state type and which manually adjusts
|
||||
the size of the internal variables. Before performing the actual resizing
|
||||
odeint always checks if the sizes of the state and the interal variable
|
||||
odeint always checks if the sizes of the state and the internal variable
|
||||
differ and only resizes if they are different.
|
||||
</p>
|
||||
<p>
|
||||
@ -836,7 +837,7 @@
|
||||
system and your state you have to call <code class="computeroutput"><span class="identifier">adjust_size</span></code>
|
||||
by hand. The second resizer is the <code class="computeroutput"><span class="identifier">always_resizer</span></code>
|
||||
which tries to resize the internal variables at every call of <code class="computeroutput"><span class="identifier">do_step</span></code>. Typical use cases for this kind
|
||||
of resizer are self expanding lattics like shown in the tutorial or partial
|
||||
of resizer are self expanding lattices like shown in the tutorial or partial
|
||||
differential equations with an adaptive grid. The third class of resizer
|
||||
is the <code class="computeroutput"><span class="identifier">never_resizer</span></code> which
|
||||
means that the internal variables are never adjusted automatically and
|
||||
@ -844,15 +845,13 @@
|
||||
</p>
|
||||
<p>
|
||||
There is a second mechanism which influences the resizing and which controls
|
||||
if a state type is at least resizeable - a metafunction <code class="computeroutput"><span class="identifier">is_resizeable</span></code>.
|
||||
This metafunction returns a static boolean value if any type is resizable.
|
||||
For example it will return <code class="computeroutput"><span class="keyword">true</span></code>
|
||||
for <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span>
|
||||
<span class="identifier">T</span> <span class="special">></span></code>
|
||||
but <code class="computeroutput"><span class="keyword">false</span></code> for <code class="computeroutput"><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="identifier">T</span> <span class="special">></span></code>.
|
||||
if a state type is at least resizeable - a meta-function <code class="computeroutput"><span class="identifier">is_resizeable</span></code>. This meta-function returns
|
||||
a static boolean value if any type is resizable. For example it will return
|
||||
<code class="computeroutput"><span class="keyword">true</span></code> for <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span> <span class="identifier">T</span> <span class="special">></span></code> but <code class="computeroutput"><span class="keyword">false</span></code>
|
||||
for <code class="computeroutput"><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="identifier">T</span> <span class="special">></span></code>.
|
||||
By default and for unknown types <code class="computeroutput"><span class="identifier">is_resizeable</span></code>
|
||||
returns <code class="computeroutput"><span class="keyword">false</span></code>, so if you have
|
||||
your own type you need to specialize this metafunction. For more details
|
||||
your own type you need to specialize this meta-function. For more details
|
||||
on the resizing mechanism see the section <a class="link" href="state_types__algebras_and_operations.html" title="State types, algebras and operations">Adapt
|
||||
you own state types</a>.
|
||||
</p>
|
||||
@ -875,7 +874,7 @@
|
||||
<li class="listitem">
|
||||
<code class="computeroutput"><span class="identifier">runge_kutta4</span></code> is a good
|
||||
stepper for constant step sizes. It is widely used and very well known.
|
||||
If you need to create artifical time series this stepper should be
|
||||
If you need to create artificial time series this stepper should be
|
||||
the first choice.
|
||||
</li>
|
||||
<li class="listitem">
|
||||
@ -1193,7 +1192,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Standard method with error control and dense ouput, to be used
|
||||
Standard method with error control and dense output, to be used
|
||||
in controlled_error_stepper and in dense_output_controlled_explicit_fsal.
|
||||
</p>
|
||||
</td>
|
||||
@ -1481,9 +1480,9 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Dense ouput for <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
Dense output for <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
and <a class="link" href="../concepts/error_stepper.html" title="Error Stepper">Error
|
||||
Stepper</a> from above if they provide dense ouput functionality
|
||||
Stepper</a> from above if they provide dense output functionality
|
||||
(like <code class="computeroutput"><span class="identifier">euler</span></code> and
|
||||
<code class="computeroutput"><span class="identifier">runge_kutta_dopri5</span></code>).
|
||||
Order depends on the given stepper.
|
||||
@ -1583,8 +1582,8 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Stepper with step size and order control as well as dense ouput.
|
||||
Very good if high precision and dense ouput is required.
|
||||
Stepper with step size and order control as well as dense output.
|
||||
Very good if high precision and dense output is required.
|
||||
</p>
|
||||
</td>
|
||||
</tr>
|
||||
@ -1741,12 +1740,12 @@
|
||||
<tr>
|
||||
<td>
|
||||
<p>
|
||||
Dense Ouput Rosenbrock 4
|
||||
Dense Output Rosenbrock 4
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
<code class="computeroutput"><span class="identifier">rosenbrock4_dense_ouput</span></code>
|
||||
<code class="computeroutput"><span class="identifier">rosenbrock4_dense_output</span></code>
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
@ -1899,12 +1898,12 @@
|
||||
</h4></div></div></div>
|
||||
<p>
|
||||
Of course, one can write own steppers which are fully compatible with odeint.
|
||||
They only have to fullfil one or several of the stepper <a class="link" href="../concepts.html" title="Concepts">Concepts</a>
|
||||
They only have to fulfill one or several of the stepper <a class="link" href="../concepts.html" title="Concepts">Concepts</a>
|
||||
of odeint.
|
||||
</p>
|
||||
<p>
|
||||
We will illustrate how to write your own stepper with the example of the
|
||||
stochastic euler method. This method is suited to solve stochastic differential
|
||||
stochastic Euler method. This method is suited to solve stochastic differential
|
||||
equations (SDEs). A SDE has the form
|
||||
</p>
|
||||
<p>
|
||||
@ -1918,8 +1917,8 @@
|
||||
the SDE is said to have additive noise. It is not possible to solve SDE
|
||||
with the classical solvers for ODEs since the noisy part of the SDE has
|
||||
to be scaled differently then the deterministic part with respect to the
|
||||
time step. BUt there exist many solvers for SDEs. A classical and easy
|
||||
method is the stochastic euler solver. It works by iterating
|
||||
time step. But there exist many solvers for SDEs. A classical and easy
|
||||
method is the stochastic Euler solver. It works by iterating
|
||||
</p>
|
||||
<p>
|
||||
<span class="emphasis"><em>x(t+Δ t) = x(t) + Δ t f(x(t)) + Δ t<sup>1/2</sup> g(x) ξ(t)</em></span>
|
||||
@ -1953,7 +1952,7 @@
|
||||
</p>
|
||||
<p>
|
||||
The types are needed in order to fulfill the stepper concept. As internal
|
||||
state and deriv type we use simple arrays in the stochastic euler, they
|
||||
state and deriv type we use simple arrays in the stochastic Euler, they
|
||||
are needed for the temporaries. The stepper has the order one which is
|
||||
returned from the <code class="computeroutput"><span class="identifier">order</span><span class="special">()</span></code> function.
|
||||
</p>
|
||||
@ -1988,7 +1987,7 @@
|
||||
<p>
|
||||
</p>
|
||||
<p>
|
||||
This is all. It is quite simple and the stochastic euler stepper implement
|
||||
This is all. It is quite simple and the stochastic Euler stepper implement
|
||||
here is quite general. Of course it can be enhanced, for example
|
||||
</p>
|
||||
<div class="itemizedlist"><ul class="itemizedlist" type="disc">
|
||||
@ -2018,7 +2017,7 @@
|
||||
<p>
|
||||
where ξ is Gaussian white noise with standard deviation <span class="emphasis"><em>σ</em></span>.
|
||||
Implementing the Ornstein-Uhlenbeck process is quite simple. We need two
|
||||
functions or funtors - one for the deterministic and one for the stochastic
|
||||
functions or functors - one for the deterministic and one for the stochastic
|
||||
part:
|
||||
</p>
|
||||
<p>
|
||||
|
||||
@ -22,7 +22,7 @@
|
||||
</h3></div></div></div>
|
||||
<p>
|
||||
Most steppers in odeint also accept the state give as a range. A range is
|
||||
sequence of values modelled by a range concept. See <a href="http://www.boost.org/doc/libs/release/libs/range/index.html" target="_top">Boost.Range</a>
|
||||
sequence of values modeled by a range concept. See <a href="http://www.boost.org/doc/libs/release/libs/range/index.html" target="_top">Boost.Range</a>
|
||||
for an overview over existing concepts and examples of ranges. This means
|
||||
that the <code class="computeroutput"><span class="identifier">state_type</span></code> of the
|
||||
stepper must not necessarily used to call the <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
@ -40,7 +40,7 @@
|
||||
<p>
|
||||
Another use case is a system consisting of coupled units where you want to
|
||||
initialize each unit separately with the ODE of the uncoupled unit. An example
|
||||
is a chain of coupled van-der-Pol-oscillators which are initialized uniformily
|
||||
is a chain of coupled van-der-Pol-oscillators which are initialized uniformly
|
||||
from the uncoupled van-der-Pol-oscillator. Then you can use <a href="http://www.boost.org/doc/libs/release/libs/range/index.html" target="_top">Boost.Range</a>
|
||||
to solve only one individual oscillator in the chain.
|
||||
</p>
|
||||
@ -207,8 +207,8 @@
|
||||
</table></div>
|
||||
</div>
|
||||
<br class="table-break"><div class="table">
|
||||
<a name="boost_sandbox_numeric_odeint.odeint_in_detail.using_boost__range.algebra_supporting_boost_range"></a><p class="title"><b>Table 1.11. Algebra supporting Boost.Range</b></p>
|
||||
<div class="table-contents"><table class="table" summary="Algebra supporting Boost.Range">
|
||||
<a name="boost_sandbox_numeric_odeint.odeint_in_detail.using_boost__range.algebras_supporting_boost_range"></a><p class="title"><b>Table 1.11. Algebras supporting Boost.Range</b></p>
|
||||
<div class="table-contents"><table class="table" summary="Algebras supporting Boost.Range">
|
||||
<colgroup><col></colgroup>
|
||||
<thead><tr><th>
|
||||
<p>
|
||||
|
||||
@ -33,7 +33,7 @@
|
||||
<p>
|
||||
</p>
|
||||
<p>
|
||||
This behaviour is suitable for most systems, especially if your system only
|
||||
This behavior is suitable for most systems, especially if your system only
|
||||
does not contain any data or only a few parameters. However, in some cases
|
||||
you might pass some large amount of data with you system function and passing
|
||||
them by value is not desired since the data will be copied.
|
||||
@ -41,8 +41,8 @@
|
||||
<p>
|
||||
In such cases you can easily use <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span></code> (and
|
||||
its relative <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">cref</span></code>) which passes its argument by reference
|
||||
(or const reference). odeint will unpack the arguments and no copying at
|
||||
all of you system will take place:
|
||||
(or constant reference). odeint will unpack the arguments and no copying
|
||||
at all of you system will take place:
|
||||
</p>
|
||||
<p>
|
||||
</p>
|
||||
|
||||
@ -155,7 +155,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
The 2D phase oscillator example shows how a two-dimenstional lattice
|
||||
The 2D phase oscillator example shows how a two-dimensional lattice
|
||||
can used with odeint and how matrix types can be used as state
|
||||
types in odeint.
|
||||
</p>
|
||||
@ -210,7 +210,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
The Lorenz paramaters examples show how ensembles of ordinary differential
|
||||
The Lorenz parameters examples show how ensembles of ordinary differential
|
||||
equations can be solved by means of Thrust to study the dependence
|
||||
of an ODE on some parameters.
|
||||
</p>
|
||||
@ -224,7 +224,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
The MTL-Gauss-packet example shows how the mtl can be easily used
|
||||
The MTL-Gauss-packet example shows how the MTL can be easily used
|
||||
with odeint.
|
||||
</p>
|
||||
</td>
|
||||
|
||||
@ -37,16 +37,16 @@
|
||||
system. In principle <span class="emphasis"><em>n</em></span> of these perturbations exit,
|
||||
they form a hypercube and evolve in the time. The Lyapunov exponents are
|
||||
then defined as logarithmic growth rates of the perturbations. If one Lyapunov
|
||||
exponent is larger then zero nearby trajectories divergy exponentially hence
|
||||
exponent is larger then zero nearby trajectories diverge exponentially hence
|
||||
they are chaotic. If the largest Lyapunov exponent is zero one is usually
|
||||
faced with periodic motion. In the case of a largest Lyapunov exponent smaller
|
||||
then zero the converges the a fixed point. More informations about Lyapunov
|
||||
then zero the converges the a fixed point. More information's about Lyapunov
|
||||
exponents and nonlinear dynamical systems can be found in many textbooks,
|
||||
see for example .
|
||||
</p>
|
||||
<p>
|
||||
To calculate the Lyapunov exponents numerically one usually solves the equations
|
||||
of motion for <span class="emphasis"><em>n</em></span> perturbations und orthonormalized them
|
||||
of motion for <span class="emphasis"><em>n</em></span> perturbations and orthonormalizes them
|
||||
every <span class="emphasis"><em>k</em></span> steps. The Lyapunov exponent is the average
|
||||
the logarithm of the stretching factor of each perturbation.
|
||||
</p>
|
||||
@ -78,7 +78,7 @@
|
||||
<p>
|
||||
We need also need to integrate the set of the perturbations. This is done
|
||||
in parallel to the original system, hence within one system function. Of
|
||||
course, we want to use the above defintion of the Lorenz system, hence the
|
||||
course, we want to use the above definition of the Lorenz system, hence the
|
||||
definition of the system function including the Lorenz system itself and
|
||||
the perturbation could look like:
|
||||
</p>
|
||||
@ -161,7 +161,7 @@
|
||||
The problem is now, that <code class="computeroutput"><span class="identifier">x</span></code>
|
||||
is the full state containing also the perturbations and <code class="computeroutput"><span class="identifier">integrate_n_steps</span></code>
|
||||
does not know that it should only use 3 element. In detail, odeint and its
|
||||
steppers determine the length of the system under consideration be determing
|
||||
steppers determine the length of the system under consideration be determining
|
||||
the length of the state. In the classical solvers from the problem was solved
|
||||
by pointer to the state and an appropriate length, something similar to
|
||||
</p>
|
||||
@ -219,7 +219,7 @@
|
||||
<li class="listitem">
|
||||
We initialize the perturbations. They are stored linearly behind the
|
||||
state of the Lorenz system. The perturbation are initialized such that
|
||||
they fullfill <span class="emphasis"><em>< p <sub>​i</sub> , p<sub>​j</sub> > = δ <sub>​ij</sub> </em></span> where <span class="emphasis"><em><,></em></span>
|
||||
they fulfill <span class="emphasis"><em>< p <sub>​i</sub> , p<sub>​j</sub> > = δ <sub>​ij</sub> </em></span> where <span class="emphasis"><em><,></em></span>
|
||||
is the classical scalar product and <span class="emphasis"><em>δ <sub>​ij</sub></em></span> is the Kronecker
|
||||
symbol.
|
||||
</li>
|
||||
|
||||
@ -36,12 +36,12 @@
|
||||
the ODE</a>
|
||||
</h4></div></div></div>
|
||||
<p>
|
||||
First of all, you have to specify the datatype that represents a state
|
||||
First of all, you have to specify the data type that represents a state
|
||||
of your system <span class="emphasis"><em>x</em></span>. 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 <code class="computeroutput"><span class="identifier">vector</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span></code> or <code class="computeroutput"><span class="identifier">vector</span><span class="special"><</span> <span class="identifier">complex</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="special">></span></code>
|
||||
to represent the system state. However, odeint can deal with other container
|
||||
types as well, e.g. <code class="computeroutput"><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">N</span> <span class="special">></span></code> as long as it is fullfils some requirements
|
||||
types as well, e.g. <code class="computeroutput"><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">N</span> <span class="special">></span></code> as long as it is fulfills some requirements
|
||||
defined below.
|
||||
</p>
|
||||
<p>
|
||||
@ -108,10 +108,10 @@
|
||||
</h4></div></div></div>
|
||||
<p>
|
||||
Numerical integration works iteratively, that means you start at a state
|
||||
<span class="emphasis"><em>x(t)</em></span> and perform a timestep of length <span class="emphasis"><em>dt</em></span>
|
||||
<span class="emphasis"><em>x(t)</em></span> and perform a time-step of length <span class="emphasis"><em>dt</em></span>
|
||||
to obtain the approximate state <span class="emphasis"><em>x(t+dt)</em></span>. There exist
|
||||
many different methods to perform such a timestep each of which has a certain
|
||||
order <span class="emphasis"><em>q</em></span>. If the order of a method is <span class="emphasis"><em>q</em></span>
|
||||
many different methods to perform such a time-step each of which has a
|
||||
certain order <span class="emphasis"><em>q</em></span>. If the order of a method is <span class="emphasis"><em>q</em></span>
|
||||
than it is accurate up to term <span class="emphasis"><em>~dt<sup>q</sup></em></span> that means the
|
||||
error in <span class="emphasis"><em>x</em></span> made by such a step is <span class="emphasis"><em>~dt<sup>q+1</sup></em></span>.
|
||||
odeint provides several steppers of different orders from which you can
|
||||
@ -414,7 +414,7 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Standard method with error control and dense ouput, to be used
|
||||
Standard method with error control and dense output, to be used
|
||||
in controlled_error_stepper and in dense_output_controlled_explicit_fsal.
|
||||
</p>
|
||||
</td>
|
||||
@ -702,9 +702,9 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Dense ouput for <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
Dense output for <a class="link" href="../concepts/stepper.html" title="Stepper">Stepper</a>
|
||||
and <a class="link" href="../concepts/error_stepper.html" title="Error Stepper">Error
|
||||
Stepper</a> from above if they provide dense ouput functionality
|
||||
Stepper</a> from above if they provide dense output functionality
|
||||
(like <code class="computeroutput"><span class="identifier">euler</span></code> and
|
||||
<code class="computeroutput"><span class="identifier">runge_kutta_dopri5</span></code>).
|
||||
Order depends on the given stepper.
|
||||
@ -804,8 +804,8 @@
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
Stepper with step size and order control as well as dense ouput.
|
||||
Very good if high precision and dense ouput is required.
|
||||
Stepper with step size and order control as well as dense output.
|
||||
Very good if high precision and dense output is required.
|
||||
</p>
|
||||
</td>
|
||||
</tr>
|
||||
@ -962,12 +962,12 @@
|
||||
<tr>
|
||||
<td>
|
||||
<p>
|
||||
Dense Ouput Rosenbrock 4
|
||||
Dense Output Rosenbrock 4
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
<p>
|
||||
<code class="computeroutput"><span class="identifier">rosenbrock4_dense_ouput</span></code>
|
||||
<code class="computeroutput"><span class="identifier">rosenbrock4_dense_output</span></code>
|
||||
</p>
|
||||
</td>
|
||||
<td>
|
||||
@ -1124,15 +1124,16 @@
|
||||
with Constant Step Size</a>
|
||||
</h4></div></div></div>
|
||||
<p>
|
||||
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 <span class="emphasis"><em>q+1</em></span>).
|
||||
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.
|
||||
The declaration of the stepper requires the state type as template parameter.
|
||||
The integration can now be done by using the <code class="computeroutput"><span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">Stepper</span><span class="special">,</span> <span class="identifier">System</span><span class="special">,</span> <span class="identifier">state</span><span class="special">,</span> <span class="identifier">start_time</span><span class="special">,</span> <span class="identifier">end_time</span><span class="special">,</span> <span class="identifier">step_size</span>
|
||||
The basic stepper just performs one time-step and doesn't give you any
|
||||
information about the error that was made (except that you know it is of
|
||||
order <span class="emphasis"><em>q+1</em></span>). 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) because 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. The declaration of the stepper requires
|
||||
the state type as template parameter. The integration can now be done by
|
||||
using the <code class="computeroutput"><span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">Stepper</span><span class="special">,</span> <span class="identifier">System</span><span class="special">,</span> <span class="identifier">state</span><span class="special">,</span> <span class="identifier">start_time</span><span class="special">,</span> <span class="identifier">end_time</span><span class="special">,</span> <span class="identifier">step_size</span>
|
||||
<span class="special">)</span></code> function from odeint:
|
||||
</p>
|
||||
<p>
|
||||
@ -1145,7 +1146,7 @@
|
||||
<p>
|
||||
This call integrates the system defined by <code class="computeroutput"><span class="identifier">harmonic_oscillator</span></code>
|
||||
using the RK4 method from <span class="emphasis"><em>t=0</em></span> to <span class="emphasis"><em>10</em></span>
|
||||
with a stepsize <span class="emphasis"><em>dt=0.01</em></span> and the initial condition
|
||||
with a step-size <span class="emphasis"><em>dt=0.01</em></span> and the initial condition
|
||||
given in <code class="computeroutput"><span class="identifier">x</span></code>. The result,
|
||||
<span class="emphasis"><em>x(t=10)</em></span> is stored in <code class="computeroutput"><span class="identifier">x</span></code>
|
||||
(in-place). Each stepper defines a <code class="computeroutput"><span class="identifier">do_step</span></code>
|
||||
@ -1207,7 +1208,7 @@
|
||||
This is a nice technique; however one drawback is that one always needs
|
||||
to define both steppers. Of course, one could also write the instantiation
|
||||
of the controlled stepper into the call of the integrate function but a
|
||||
complete knowledge of the underlying stepper types is still neccessary.
|
||||
complete knowledge of the underlying stepper types is still necessary.
|
||||
Another point is, that the error tolerances for the step size control are
|
||||
not easily included into the controlled stepper. Both issues can be solved
|
||||
by using <code class="computeroutput"><span class="identifier">make_controlled</span></code>:
|
||||
@ -1257,7 +1258,7 @@
|
||||
</p>
|
||||
<p>
|
||||
In the tables below, one can find all steppers which are working with
|
||||
<code class="computeroutput"><span class="identifier">make_controlled</span></code> and <code class="computeroutput"><span class="identifier">make_dense_output</span></code> which is the analogon
|
||||
<code class="computeroutput"><span class="identifier">make_controlled</span></code> and <code class="computeroutput"><span class="identifier">make_dense_output</span></code> which is the analog
|
||||
for the dense output steppers.
|
||||
</p>
|
||||
<div class="table">
|
||||
|
||||
@ -20,7 +20,7 @@
|
||||
<a name="boost_sandbox_numeric_odeint.tutorial.references"></a><a class="link" href="references.html" title="References">References</a>
|
||||
</h3></div></div></div>
|
||||
<p>
|
||||
<span class="bold"><strong>General informations about numerical integration of
|
||||
<span class="bold"><strong>General information about numerical integration of
|
||||
ordinary differential equations:</strong></span>
|
||||
</p>
|
||||
<p>
|
||||
|
||||
@ -81,7 +81,7 @@
|
||||
to ensure these conservation laws. The odeint library provides classes
|
||||
for separable Hamiltonian systems, which can be written in the form <span class="emphasis"><em>H
|
||||
= Σ p<sub>​i</sub><sup>2</sup> / (2m<sub>​i</sub>) + H<sub>​q</sub>(q)</em></span>, where <span class="emphasis"><em>H<sub>​q</sub>(q)</em></span> only
|
||||
depends on the coordinates. Alltough this functional form might look a
|
||||
depends on the coordinates. Although 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 <span class="emphasis"><em>dq<sub>​i</sub> / dt = p<sub>​i</sub></em></span> / m<sub>​i</sub> , <span class="emphasis"><em>dp<sub>​i</sub> / dt =
|
||||
@ -151,7 +151,7 @@
|
||||
<code class="computeroutput"><span class="identifier">pair</span><span class="special"><</span>
|
||||
<span class="identifier">container_type</span> <span class="special">,</span>
|
||||
<span class="identifier">container_type</span> <span class="special">></span></code>
|
||||
since it needs the informations about the coordinates and the momenta.
|
||||
since it needs the information about the coordinates and the momenta.
|
||||
</p>
|
||||
<p>
|
||||
As system function we have to provide <span class="emphasis"><em>f(p)</em></span> and <span class="emphasis"><em>f(q)</em></span>:
|
||||
@ -206,7 +206,7 @@
|
||||
</p>
|
||||
<p>
|
||||
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
|
||||
arbitrary 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.
|
||||
</p>
|
||||
|
||||
@ -84,7 +84,7 @@
|
||||
<p>
|
||||
</p>
|
||||
<p>
|
||||
Of courst, one can also use classical functions to implement the state
|
||||
Of course, one can also use classical functions to implement the state
|
||||
function. In this cast the Stuart-Landau oscillator looks like
|
||||
</p>
|
||||
<p>
|
||||
@ -148,7 +148,7 @@
|
||||
The FPU is solved again by a symplectic solver, but in our case we can
|
||||
speed up the computation because the <span class="emphasis"><em>q</em></span> components
|
||||
trivially reduce to <span class="emphasis"><em>dq<sub>​i</sub> / dt = p<sub>​i</sub></em></span>. odeint is capable
|
||||
of doing this peformance improvement. All you have to do is to call the
|
||||
of doing this performance improvement. All you have to do is to call the
|
||||
symplectic solver with an state function for the <span class="emphasis"><em>p</em></span>
|
||||
components. Here, is how this function looks like
|
||||
</p>
|
||||
@ -281,10 +281,10 @@
|
||||
<span class="emphasis"><em>dφ<sub>​k</sub> / dt = ω<sub>​k</sub> + ε / N Σ<sub>​j</sub> sin( φ<sub>​j</sub> - φ<sub>​k</sub> )</em></span>
|
||||
</p>
|
||||
<p>
|
||||
The natural frequencies <span class="emphasis"><em>ω<sub>​i</sub></em></span> of each oscialltor follow
|
||||
The natural frequencies <span class="emphasis"><em>ω<sub>​i</sub></em></span> of each oscillator follow
|
||||
some distribution and <span class="emphasis"><em>ε</em></span> is the coupling strength. We
|
||||
choose here a Lorentzian distribution for <span class="emphasis"><em>ω<sub>​i</sub></em></span>. Interestingly
|
||||
a phase transition can be observed if the coupling strenght exceeds a critical
|
||||
a phase transition can be observed if the coupling strength exceeds a critical
|
||||
value. Above this value synchronization sets in and some of the oscillators
|
||||
oscillate with the same frequency despite their different natural frequencies.
|
||||
The transition is also called Kuramoto transition. Its behavior can be
|
||||
@ -486,7 +486,7 @@
|
||||
Note, that the <code class="computeroutput"><span class="identifier">state_type</span></code>
|
||||
and the <code class="computeroutput"><span class="identifier">deriv_type</span></code> are
|
||||
now a compile-time fusion sequences. They have to be explicitly defined.
|
||||
Next, we define the ordinary differential equation which is completety
|
||||
Next, we define the ordinary differential equation which is completely
|
||||
equivalent to the example in the first section of this tutorial
|
||||
</p>
|
||||
<p>
|
||||
@ -588,7 +588,7 @@
|
||||
<p>
|
||||
By default, odeint can be used with <code class="computeroutput"><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">matrix</span><span class="special"><</span> <span class="identifier">T</span> <span class="special">></span></code> as state type for matrices. A simple
|
||||
example is a two-dimensional lattice of coupled phase oscillators. We like
|
||||
phas oscillators, they are extremly easy and might serve for different
|
||||
phase oscillators, they are extremely easy and might serve for different
|
||||
demonstration purposes. Other matrix types like <code class="computeroutput"><span class="identifier">mtl</span><span class="special">::</span><span class="identifier">dense_matrix</span></code>
|
||||
or blitz arrays and matrices can used as well but need some kind of activation
|
||||
in order to work with odeint. This activation is described in following
|
||||
@ -678,7 +678,7 @@
|
||||
<code class="computeroutput"><span class="keyword">double</span></code>, <code class="computeroutput"><span class="identifier">complex</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span></code> you can also use arbitrary precision
|
||||
types, like the types from <a href="http://gmplib.org/" target="_top">gmp</a>
|
||||
and <a href="http://www.mpfr.org/" target="_top">mpfr</a>. But you have to be
|
||||
carful about instantiating any numbers.
|
||||
careful about instantiating any numbers.
|
||||
</p>
|
||||
<p>
|
||||
For gmp types you have to set the default precision before any number is
|
||||
@ -766,8 +766,8 @@
|
||||
and adjust it's internal storage accordingly.
|
||||
</p>
|
||||
<p>
|
||||
We exemplarily show this for a Hamiltonian system of nonlinear, disordered
|
||||
oscillators with nonlinear nearest neighbour coupling.
|
||||
We exemplary show this for a Hamiltonian system of nonlinear, disordered
|
||||
oscillators with nonlinear nearest neighbor coupling.
|
||||
</p>
|
||||
<p>
|
||||
The system function is implemented in terms of a class that also provides
|
||||
|
||||
@ -24,7 +24,7 @@
|
||||
An important class of ordinary differential equations are so called stiff
|
||||
system which are characterized by two or more time scales of different order.
|
||||
Examples of such systems are found in chemical systems where reaction rates
|
||||
of individual subreaction might differ over large ranges.
|
||||
of individual sub-reaction might differ over large ranges.
|
||||
</p>
|
||||
<p>
|
||||
CHECK:
|
||||
|
||||
@ -82,7 +82,7 @@
|
||||
perfectly suited for such kinds of problems where one needs a large number
|
||||
of particles (oscillators). We start by defining the state type which is
|
||||
a <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span></code>. The content of this vector
|
||||
lives on the GPU. If you are not familar with this we recommend reading
|
||||
lives on the GPU. If you are not familiar with this we recommend reading
|
||||
the <span class="emphasis"><em>Getting started</em></span> section on the <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a>
|
||||
website.
|
||||
</p>
|
||||
@ -317,7 +317,7 @@
|
||||
In principle we can use all the techniques from the previous phase oscillator
|
||||
ensemble example, but we have to take special care about the coupling of
|
||||
the oscillators. To efficiently implement the coupling you can use a very
|
||||
elegant way employing Thrust's permutation iterator. A permuation iterator
|
||||
elegant way employing Thrust's permutation iterator. A permutation iterator
|
||||
behaves like a normal iterator on a vector but it does not iterate along
|
||||
the usual order of the elements. It rather iterates along some permutation
|
||||
of the elements defined by some index map. To realize the nearest neighbor
|
||||
@ -402,9 +402,9 @@
|
||||
Note, how easy you can obtain the value for the left and right neighboring
|
||||
oscillator in the system functor using the permutation iterators. But,
|
||||
the call of the <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span></code>
|
||||
function looks relativ complicated. Every term of the r.h.s. of the ODE
|
||||
is resembled by one iterator packed in exactly the same way as it is unpacked
|
||||
in the functor above.
|
||||
function looks relatively complicated. Every term of the r.h.s. of the
|
||||
ODE is resembled by one iterator packed in exactly the same way as it is
|
||||
unpacked in the functor above.
|
||||
</p>
|
||||
<p>
|
||||
Now we put everything together. We create random initial conditions and
|
||||
@ -456,8 +456,8 @@
|
||||
and Cuda are parameter studies of relatively small systems. Consider for
|
||||
example the three-dimensional Lorenz system from the chaotic systems example
|
||||
in the previous section which has three parameters. If you want to study
|
||||
the behaviour of this system for different parameters you usually have
|
||||
to integrate the system for many parameter values. Using thrust and odeint
|
||||
the behavior of this system for different parameters you usually have to
|
||||
integrate the system for many parameter values. Using thrust and odeint
|
||||
you can do this integration in parallel, hence you integrate a whole ensemble
|
||||
of Lorenz systems where each individual realization has a different parameter
|
||||
value.
|
||||
@ -465,7 +465,7 @@
|
||||
<p>
|
||||
In the following we will show how you can use <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a>
|
||||
to integrate the above mentioned ensemble of Lorenz systems. We will vary
|
||||
only the parameter <span class="emphasis"><em>β</em></span> but it is straightfoward to vary
|
||||
only the parameter <span class="emphasis"><em>β</em></span> but it is straightforward to vary
|
||||
other parameters or even two or all three parameters. Furthermore, we will
|
||||
use the largest Lyapunov exponent to quantify the behavior of the system
|
||||
(chaoticity).
|
||||
@ -551,7 +551,7 @@
|
||||
The length of the state is <span class="emphasis"><em>3N</em></span> where <span class="emphasis"><em>N</em></span>
|
||||
is the number of systems. The system is encoded into this vector such that
|
||||
all <span class="emphasis"><em>x</em></span> components come first, then every <span class="emphasis"><em>y</em></span>
|
||||
compoents and finally every <span class="emphasis"><em>z</em></span> components. Implementing
|
||||
components and finally every <span class="emphasis"><em>z</em></span> components. Implementing
|
||||
the device function is then a simple task, you only have to decompose the
|
||||
tuple originating from the zip iterators.
|
||||
</p>
|
||||
@ -559,7 +559,7 @@
|
||||
Besides the system without perturbations we furthermore need to calculate
|
||||
the system including linearized equations governing the time evolution
|
||||
of small perturbations. Using the method from above this is straightforward,
|
||||
with a small difficulty that thrust's tuples have a maximal arity of 10.
|
||||
with a small difficulty that Thrust's tuples have a maximal arity of 10.
|
||||
But this is only a small problem since we can create a zip iterator packed
|
||||
with zip iterators. So the top level zip iterator contains one zip iterator
|
||||
for the state, one normal iterator for the parameter, and one zip iterator
|
||||
@ -604,10 +604,10 @@
|
||||
We initialize them such that <span class="emphasis"><em>x=y=z=10</em></span>, <span class="emphasis"><em>dx=1</em></span>,
|
||||
and <span class="emphasis"><em>dy=dz=0</em></span>. We define a stepper type, a controlled
|
||||
Runge-Kutta Dormand-Prince 5 stepper. We start with some integration to
|
||||
overcome the transient behavior. For this, we do not involve the pertubation
|
||||
overcome the transient behavior. For this, we do not involve the perturbation
|
||||
and run the algorithm only on the state <span class="emphasis"><em>x,y,z</em></span> without
|
||||
any observer. Note, how <a href="http://www.boost.org/doc/libs/release/libs/range/index.html" target="_top">Boost.Range</a>
|
||||
is used for partial integration of the state vector without pertubations
|
||||
is used for partial integration of the state vector without perturbations
|
||||
(the first half of the whole state). After the transient, the full system
|
||||
with perturbations is integrated and the Lyapunov exponents are calculated
|
||||
and written to <code class="computeroutput"><span class="identifier">stdout</span></code>.
|
||||
|
||||
109
doc/index.html
109
doc/index.html
@ -1,20 +1,20 @@
|
||||
<html>
|
||||
<head>
|
||||
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
|
||||
<title>Chapter 1. boost.sandbox.numeric.odeint</title>
|
||||
<title>Chapter 1. Boost.Numeric.Odeint</title>
|
||||
<link rel="stylesheet" href="boostbook.css" type="text/css">
|
||||
<meta name="generator" content="DocBook XSL Stylesheets V1.75.2">
|
||||
<link rel="home" href="index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
|
||||
<link rel="next" href="boost_sandbox_numeric_odeint/getting_started.html" title="Getting started">
|
||||
<link rel="home" href="index.html" title="Chapter 1. Boost.Numeric.Odeint">
|
||||
<link rel="next" href="boost_numeric_odeint/getting_started.html" title="Getting started">
|
||||
</head>
|
||||
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
|
||||
<table cellpadding="2" width="100%"><tr><td valign="top"></td></tr></table>
|
||||
<hr>
|
||||
<div class="spirit-nav"><a accesskey="n" href="boost_sandbox_numeric_odeint/getting_started.html"><img src="images/next.png" alt="Next"></a></div>
|
||||
<div class="spirit-nav"><a accesskey="n" href="boost_numeric_odeint/getting_started.html"><img src="images/next.png" alt="Next"></a></div>
|
||||
<div class="chapter">
|
||||
<div class="titlepage"><div>
|
||||
<div><h2 class="title">
|
||||
<a name="odeint"></a>Chapter 1. boost.sandbox.numeric.odeint</h2></div>
|
||||
<a name="odeint"></a>Chapter 1. Boost.Numeric.Odeint</h2></div>
|
||||
<div><div class="author"><h3 class="author">
|
||||
<span class="firstname">Karsten</span> <span class="surname">Ahnert</span>
|
||||
</h3></div></div>
|
||||
@ -33,102 +33,71 @@
|
||||
<div class="toc">
|
||||
<p><b>Table of Contents</b></p>
|
||||
<dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/getting_started.html">Getting started</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/getting_started.html">Getting started</a></span></dt>
|
||||
<dd><dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/getting_started/overview.html">Overview</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/getting_started/usage__compilation__headers.html">Usage,
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/getting_started/overview.html">Overview</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/getting_started/usage__compilation__headers.html">Usage,
|
||||
Compilation, Headers</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/getting_started/short_example.html">Short
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/getting_started/short_example.html">Short
|
||||
Example</a></span></dt>
|
||||
</dl></dd>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial.html">Tutorial</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial.html">Tutorial</a></span></dt>
|
||||
<dd><dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/harmonic_oscillator.html">Harmonic
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/harmonic_oscillator.html">Harmonic
|
||||
oscillator</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/solar_system.html">Solar
|
||||
system</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html">Chaotic
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/solar_system.html">Solar system</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html">Chaotic
|
||||
systems and Lyapunov exponents</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/stiff_systems.html">Stiff
|
||||
systems</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/special_topics.html">Special
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/stiff_systems.html">Stiff systems</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/special_topics.html">Special
|
||||
topics</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/using_cuda_and_thrust.html">Using
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/using_cuda_and_thrust.html">Using
|
||||
Cuda and Thrust</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/all_examples.html">All
|
||||
examples</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/tutorial/references.html">References</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/all_examples.html">All examples</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/tutorial/references.html">References</a></span></dt>
|
||||
</dl></dd>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail.html">Odeint in
|
||||
detail</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail.html">odeint in detail</a></span></dt>
|
||||
<dd><dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail/steppers.html">Steppers</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail/generation_functions.html">Generation
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail/steppers.html">Steppers</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail/generation_functions.html">Generation
|
||||
functions</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail/integrate_functions.html">Integrate
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail/integrate_functions.html">Integrate
|
||||
functions</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html">State
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html">State
|
||||
types, algebras and operations</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail/using_boost__ref.html">Using
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail/using_boost__ref.html">Using
|
||||
boost::ref</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/odeint_in_detail/using_boost__range.html">Using
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/odeint_in_detail/using_boost__range.html">Using
|
||||
boost::range</a></span></dt>
|
||||
</dl></dd>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts.html">Concepts</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts.html">Concepts</a></span></dt>
|
||||
<dd><dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/system.html">System</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/symplectic_system.html">Symplectic
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/system.html">System</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/symplectic_system.html">Symplectic
|
||||
System</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/simple_symplectic_system.html">Simple
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/simple_symplectic_system.html">Simple
|
||||
Symplectic System</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/implicit_system.html">Implicit
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/implicit_system.html">Implicit
|
||||
System</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/observer.html">Observer</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/stepper.html">Stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/error_stepper.html">Error
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/observer.html">Observer</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/stepper.html">Stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/error_stepper.html">Error Stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/controlled_stepper.html">Controlled
|
||||
Stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/controlled_stepper.html">Controlled
|
||||
Stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/dense_output_stepper.html">Dense
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/dense_output_stepper.html">Dense
|
||||
Output Stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/state_algebra_operations.html">State
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/state_algebra_operations.html">State
|
||||
Algebra Operations</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/concepts/state_wrapper.html">State
|
||||
Wrapper</a></span></dt>
|
||||
</dl></dd>
|
||||
<dt><span class="section"><a href="odeint/reference.html">Reference</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts.html">Old Concepts</a></span></dt>
|
||||
<dd><dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts/basic_stepper.html">Basic
|
||||
stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts/error_stepper.html">Error
|
||||
stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts/controlled_stepper.html">Controlled
|
||||
stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts/dense_ouput_stepper.html">Dense
|
||||
ouput stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts/size_adjusting_stepper.html">Size
|
||||
adjusting stepper</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_concepts/compositestepper.html">CompositeStepper</a></span></dt>
|
||||
</dl></dd>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_reference.html">Old Reference</a></span></dt>
|
||||
<dd><dl>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_reference/stepper_classes.html">Stepper
|
||||
classes</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_reference/integration_functions.html">Integration
|
||||
functions</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_reference/algebras.html">Algebras</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_reference/operations.html">Operations</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_sandbox_numeric_odeint/old_reference/resizing.html">Resizing</a></span></dt>
|
||||
<dt><span class="section"><a href="boost_numeric_odeint/concepts/state_wrapper.html">State Wrapper</a></span></dt>
|
||||
</dl></dd>
|
||||
</dl>
|
||||
</div>
|
||||
</div>
|
||||
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
|
||||
<td align="left"><p><small>Last revised: April 16, 2012 at 06:44:03 GMT</small></p></td>
|
||||
<td align="left"><p><small>Last revised: May 21, 2012 at 21:01:33 GMT</small></p></td>
|
||||
<td align="right"><div class="copyright-footer"></div></td>
|
||||
</tr></table>
|
||||
<hr>
|
||||
<div class="spirit-nav"><a accesskey="n" href="boost_sandbox_numeric_odeint/getting_started.html"><img src="images/next.png" alt="Next"></a></div>
|
||||
<div class="spirit-nav"><a accesskey="n" href="boost_numeric_odeint/getting_started.html"><img src="images/next.png" alt="Next"></a></div>
|
||||
</body>
|
||||
</html>
|
||||
|
||||
@ -4,16 +4,15 @@
|
||||
<title>Reference</title>
|
||||
<link rel="stylesheet" href="../boostbook.css" type="text/css">
|
||||
<meta name="generator" content="DocBook XSL Stylesheets V1.75.2">
|
||||
<link rel="home" href="../index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
|
||||
<link rel="up" href="../index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
|
||||
<link rel="prev" href="../boost_sandbox_numeric_odeint/concepts/state_wrapper.html" title="State Wrapper">
|
||||
<link rel="next" href="../boost_sandbox_numeric_odeint/old_concepts.html" title="Old Concepts">
|
||||
<link rel="home" href="../index.html" title="Chapter 1. Boost.Numeric.Odeint">
|
||||
<link rel="up" href="../index.html" title="Chapter 1. Boost.Numeric.Odeint">
|
||||
<link rel="prev" href="../boost_numeric_odeint/concepts/state_wrapper.html" title="State Wrapper">
|
||||
</head>
|
||||
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
|
||||
<table cellpadding="2" width="100%"><tr><td valign="top"></td></tr></table>
|
||||
<hr>
|
||||
<div class="spirit-nav">
|
||||
<a accesskey="p" href="../boost_sandbox_numeric_odeint/concepts/state_wrapper.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a><a accesskey="n" href="../boost_sandbox_numeric_odeint/old_concepts.html"><img src="../images/next.png" alt="Next"></a>
|
||||
<a accesskey="p" href="../boost_numeric_odeint/concepts/state_wrapper.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a>
|
||||
</div>
|
||||
<div class="section"><div class="titlepage"><div><div><h2 class="title" style="clear: both">
|
||||
<a name="odeint.reference"></a>Reference</h2></div></div></div></div>
|
||||
@ -28,7 +27,7 @@
|
||||
</tr></table>
|
||||
<hr>
|
||||
<div class="spirit-nav">
|
||||
<a accesskey="p" href="../boost_sandbox_numeric_odeint/concepts/state_wrapper.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a><a accesskey="n" href="../boost_sandbox_numeric_odeint/old_concepts.html"><img src="../images/next.png" alt="Next"></a>
|
||||
<a accesskey="p" href="../boost_numeric_odeint/concepts/state_wrapper.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a>
|
||||
</div>
|
||||
</body>
|
||||
</html>
|
||||
|
||||
@ -4,7 +4,7 @@ This concept specifies the interface a controlled stepper has to fulfill to be u
|
||||
|
||||
[heading Description]
|
||||
|
||||
A controlled stepper following this Controlled Stepper concept provides the possibilty to perform one step of the solution /x(t)/ of an ODE with step-size /dt/ to obtain /x(t+dt)/ with a given step-size /dt/.
|
||||
A controlled stepper following this Controlled Stepper concept provides the possibility to perform one step of the solution /x(t)/ of an ODE with step-size /dt/ to obtain /x(t+dt)/ with a given step-size /dt/.
|
||||
Depending on an error estimate of the solution the step might be rejected and a smaller step-size is suggested.
|
||||
|
||||
[heading Notation]
|
||||
@ -23,7 +23,7 @@ Depending on an error estimate of the solution the step might be rejected and a
|
||||
|
||||
[table
|
||||
[[Name] [Expression] [Type] [Semantics]]
|
||||
[[Do step] [`stepper.try_step( sys , x , t , dt )`] [`controlled_step_result`] [Tries one step of step size `dt`. If the step was succesfull, `success` is returned, the resulting state is written to `x`, the new time is stored in `t` and `dt` now contains a new (possibly larger) step-size for the next step. If the error was too big, `rejected` is returned and the results are neglected - `x` and `t` are unchanged and `dt` now contains a reduced step-size to be used for the next try.] ]
|
||||
[[Do step] [`stepper.try_step( sys , x , t , dt )`] [`controlled_step_result`] [Tries one step of step size `dt`. If the step was successful, `success` is returned, the resulting state is written to `x`, the new time is stored in `t` and `dt` now contains a new (possibly larger) step-size for the next step. If the error was too big, `rejected` is returned and the results are neglected - `x` and `t` are unchanged and `dt` now contains a reduced step-size to be used for the next try.] ]
|
||||
[[Do step with reference] [`stepper.try_step( boost::ref(sys) , x , t , dt )`] [`void`] [Same as above with `System` as reference] ]
|
||||
]
|
||||
|
||||
|
||||
@ -1,9 +1,9 @@
|
||||
[section Dense Output Stepper]
|
||||
|
||||
This concept specifies the interface a dense ouput stepper has to fulfill to be used within __integrate_functions.
|
||||
This concept specifies the interface a dense output stepper has to fulfill to be used within __integrate_functions.
|
||||
|
||||
[heading Description]
|
||||
A dense ouput stepper following this Dense Output Stepper concept provides the possibility to perform a single step of the solution /x(t)/ of an ODE to obtain /x(t+dt)/.
|
||||
A dense output stepper following this Dense Output Stepper concept provides the possibility to perform a single step of the solution /x(t)/ of an ODE to obtain /x(t+dt)/.
|
||||
The step-size `dt` might be adjusted automatically due to error control.
|
||||
Dense output steppers also can interpolate the solution to calculate the state /x(t')/ at any point /t <= t' <= t+dt/.
|
||||
|
||||
|
||||
@ -21,7 +21,7 @@ Provided to an `integrate` function, observers are called at each time step, at
|
||||
|
||||
[table
|
||||
[[Name] [Expression] [Type] [Semantics]]
|
||||
[[Perform the observation] [`obs( x , t )`] [`void`] [Calls the observer which can do some analyzation of the state /x/ at time /t/.] ]
|
||||
[[Perform the observation] [`obs( x , t )`] [`void`] [Calls the observer which can do some analysis or output of the state /x/ at time /t/.] ]
|
||||
]
|
||||
|
||||
[endsect]
|
||||
@ -1,6 +1,6 @@
|
||||
[section State Algebra Operations]
|
||||
|
||||
[note The following does not apply to implicit steppers like implicit_euler or rosenbrock 4 as there the `state_type` can not be changed from `ublas::vector` and no algebra/operations are used.]
|
||||
[note The following does not apply to implicit steppers like implicit_euler or Rosenbrock 4 as there the `state_type` can not be changed from `ublas::vector` and no algebra/operations are used.]
|
||||
|
||||
[heading Description]
|
||||
|
||||
@ -91,7 +91,7 @@ In the following we list the existing `Algebra`/`Operations` configurations that
|
||||
[[`State`] [`Algebra`] [`Operations`] [Remarks]]
|
||||
[[Anything supporting __boost_range, like `std::vector`, `std::list`, `boost::array`,... based on a `value_type` that supports operators +,* (typically `double`)] [`range_algebra`] [`default_operations`] [Standard implementation, applicable for most typical situations.]]
|
||||
[[`boost::array` based on a `value_type` that supports operators +,*] [`array_algebra`] [`default_operations`] [Special implementation for boost::array with better performance than `range_algebra`]]
|
||||
[[Anything that defines operators + within itself and * with scalar (Mathematically spoken, anyhting that is a vector space).] [`vector_space_algebra`] [`default_operations`] [For the use of __controlled_stepper, the template `vector_space_reduce` has to be instantiated.]]
|
||||
[[Anything that defines operators + within itself and * with scalar (Mathematically spoken, anything that is a vector space).] [`vector_space_algebra`] [`default_operations`] [For the use of __controlled_stepper, the template `vector_space_reduce` has to be instantiated.]]
|
||||
[[`thrust::device_vector`, `thrust::host_vector`] [`thrust_algebra`] [`thrust_operations`] [For running odeint on CUDA devices by using __thrust]]
|
||||
[[`boost::array` or anything which allocates the elements in a C-like manner] [`vector_space_algebra`] [`mkl_operations`] [Using the __intel_mkl in odeint for maximum performance. Currently, only the RK4 stepper is supported.]]
|
||||
]
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
[section Stepper]
|
||||
|
||||
This concepts specifies the interface a simple stepper has to fullfill to be used within the __integrate_functions.
|
||||
This concepts specifies the interface a simple stepper has to fulfill to be used within the __integrate_functions.
|
||||
|
||||
[heading Description]
|
||||
|
||||
|
||||
@ -5,7 +5,7 @@
|
||||
The System concept models the algorithmic implementation of the rhs. of the ODE ['x' = f(x,t)].
|
||||
The only requirement for this concept is that it should be callable with a specific parameter syntax (see below).
|
||||
A System is typically implemented as a function or a functor.
|
||||
Systems fullfilling this concept are required by all Runge-Kutta steppers as well as the Bulirsch-Stoer steppers.
|
||||
Systems fulfilling this concept are required by all Runge-Kutta steppers as well as the Bulirsch-Stoer steppers.
|
||||
However, symplectic and implicit steppers work with other system concepts, see __symplectic_system and __implicit_system.
|
||||
|
||||
[heading Notation]
|
||||
|
||||
@ -1,270 +0,0 @@
|
||||
[section Old 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]
|
||||
|
||||
|
||||
[section Dense ouput stepper]
|
||||
|
||||
[endsect]
|
||||
|
||||
[section Size adjusting stepper]
|
||||
|
||||
[endsect]
|
||||
|
||||
[section CompositeStepper]
|
||||
|
||||
[endsect]
|
||||
|
||||
see the wiki
|
||||
|
||||
[endsect]
|
||||
@ -1,4 +1,4 @@
|
||||
[section Odeint in detail]
|
||||
[section odeint in detail]
|
||||
|
||||
[include details_steppers.qbk]
|
||||
|
||||
|
||||
@ -1,10 +1,10 @@
|
||||
[section Using boost::range]
|
||||
|
||||
Most steppers in odeint also accept the state give as a range. A range is sequence of values modelled by a range concept. See __boost_range for an overview over existing concepts and examples of ranges. This means that the `state_type` of the stepper must not necessarily used to call the `do_step` method.
|
||||
Most steppers in odeint also accept the state give as a range. A range is sequence of values modeled by a range concept. See __boost_range for an overview over existing concepts and examples of ranges. This means that the `state_type` of the stepper must not necessarily used to call the `do_step` method.
|
||||
|
||||
One use case for __boost_range in odeint has been shown in __tut_chaotic_system where the state consists of two parts: one for the original system and one for the perturbations. The ranges are used to initialize (solve) only the system part where the perturbation part is not touched, that is a range consisting only of the system part is used. After that the complete state including the perturbations is solved.
|
||||
|
||||
Another use case is a system consisting of coupled units where you want to initialize each unit separately with the ODE of the uncoupled unit. An example is a chain of coupled van-der-Pol-oscillators which are initialized uniformily from the uncoupled van-der-Pol-oscillator. Then you can use __boost_range to solve only one individual oscillator in the chain.
|
||||
Another use case is a system consisting of coupled units where you want to initialize each unit separately with the ODE of the uncoupled unit. An example is a chain of coupled van-der-Pol-oscillators which are initialized uniformly from the uncoupled van-der-Pol-oscillator. Then you can use __boost_range to solve only one individual oscillator in the chain.
|
||||
|
||||
In short, you can __boost_range to use one state within two system functions which expect states with different sizes.
|
||||
|
||||
|
||||
@ -7,9 +7,9 @@ In odeint all system functions and observers are passed by value. For example, i
|
||||
rk4.do_step( sys , x , t , dt ); // pass sys by value
|
||||
``
|
||||
|
||||
This behaviour is suitable for most systems, especially if your system only does not contain any data or only a few parameters. However, in some cases you might pass some large amount of data with you system function and passing them by value is not desired since the data will be copied.
|
||||
This behavior is suitable for most systems, especially if your system only does not contain any data or only a few parameters. However, in some cases you might pass some large amount of data with you system function and passing them by value is not desired since the data will be copied.
|
||||
|
||||
In such cases you can easily use `boost::ref` (and its relative `boost::cref`) which passes its argument by reference (or const reference). odeint will unpack the arguments and no copying at all of you system will take place:
|
||||
In such cases you can easily use `boost::ref` (and its relative `boost::cref`) which passes its argument by reference (or constant reference). odeint will unpack the arguments and no copying at all of you system will take place:
|
||||
|
||||
``
|
||||
rk4.do_step( boost::ref( sys ) , x , t , dt ); // pass sys as references
|
||||
|
||||
@ -24,7 +24,7 @@ This is all to use the `make_controlled` mechanism. Now you can use your control
|
||||
|
||||
For the dense_output_stepper everything works similar. Here you have to specialize `boost::numeric::odeint::get_dense_output` and `boost::numeric::odeint::dense_output_factory`. These two classes have the same syntax as their relatives `get_controller` and `controller_factory`.
|
||||
|
||||
Of course, all controllers and dense-output steppers in odeint can be used with these mechanisms. In the table below you will find, which steppers is contructed from `make_controlled` or `make_dense_output` if applied on a stepper from odeint:
|
||||
Of course, all controllers and dense-output steppers in odeint can be used with these mechanisms. In the table below you will find, which steppers is constructed from `make_controlled` or `make_dense_output` if applied on a stepper from odeint:
|
||||
|
||||
[include make_controlled_table.qbk]
|
||||
[include make_dense_output_table.qbk]
|
||||
|
||||
@ -8,7 +8,7 @@ Depending on the abilities of the stepper, the integrate functions make use of s
|
||||
|
||||
[heading Equidistant observer calls]
|
||||
|
||||
If observer calls at equdistant time intervals /dt/ are needed, the `integrate_const` function should be used.
|
||||
If observer calls at equidistant time intervals /dt/ are needed, the `integrate_const` function should be used.
|
||||
|
||||
`integrate_const( stepper , system , x0 , t0 , t1 , dt )`
|
||||
|
||||
@ -83,7 +83,7 @@ Dense output is used to obtain the states /x(t)/ at the time points from the seq
|
||||
|
||||
If exactly `n` steps with fixed step size `dt` should be executed the `integrate_n_steps` function should be used.
|
||||
|
||||
[warning This function works only with steppers fullfilling the __stepper or __error_stepper concept! Providing a __controlled_stepper or __dense_out_stepper results in a compilationerror.]
|
||||
[warning This function works only with steppers fulfilling the __stepper or __error_stepper concept! Providing a __controlled_stepper or __dense_out_stepper results in a compilation error.]
|
||||
|
||||
`integrate_n_steps( stepper , system , x0 , t0 , dt , n )`
|
||||
|
||||
@ -91,14 +91,14 @@ If exactly `n` steps with fixed step size `dt` should be executed the `integrate
|
||||
|
||||
Integrates the ODE given by `system` with subsequent steps from `stepper` starting at ['x[sub 0]] and ['t[sub 0]].
|
||||
If provided, `observer` is called after every step and at the beginning with `t0`.
|
||||
This function works only with steppers fullfilling the __stepper or __error_stepper concept!
|
||||
This function works only with steppers fulfilling the __stepper or __error_stepper concept!
|
||||
It performs exactly `n` steps with fixed step size `dt` ending at `t0 + n*dt`.
|
||||
The approximate result for ['x( t[sub 0] + n dt )] is stored in `x0`.
|
||||
This function returns the end time `t0 + n*dt`.
|
||||
|
||||
[heading Convenience integrate function]
|
||||
|
||||
Additionally to the sophisticated integrate function above odeint also provides a simple `integrate` routine which uses a dense ouput stepper based on `runge_kutta_dopri5` with standard error bounds ['10[super -6]] for the steps.
|
||||
Additionally to the sophisticated integrate function above odeint also provides a simple `integrate` routine which uses a dense output stepper based on `runge_kutta_dopri5` with standard error bounds ['10[super -6]] for the steps.
|
||||
|
||||
`integrate( system , x0 , t0 , t1 , dt )`
|
||||
|
||||
|
||||
@ -6,7 +6,7 @@ This is realized by giving the user full control over the state type and the
|
||||
mathematical operations for this state type.
|
||||
Technically, this is done by introducing three concepts: StateType, Algebra,
|
||||
Operations.
|
||||
Most of the steppers in odeint expect three classes types fullfilling these
|
||||
Most of the steppers in odeint expect three classes types fulfilling these
|
||||
concepts as template parameters.
|
||||
Note that these concepts are not fully independent of each other but rather a
|
||||
valid combination must be provided in order to make the steppers work.
|
||||
@ -282,7 +282,7 @@ private :
|
||||
};
|
||||
``
|
||||
A similar class exists for the `const` version of the iterator.
|
||||
Then we have a function returning the end iterator (similarily for `const` again):
|
||||
Then we have a function returning the end iterator (similarly for `const` again):
|
||||
``
|
||||
gsl_vector_iterator end_iterator( gsl_vector *x )
|
||||
{
|
||||
@ -310,7 +310,7 @@ Again with similar definitions for the `const` versions.
|
||||
This eventually makes odeint work with gsl vectors as state types.
|
||||
The full code for these bindings is found in [github_link
|
||||
boost/numeric/odeint/external/gsl/gsl_wrapper.hpp gsl_wrapper.hpp].
|
||||
It might look rather complicated but keep in mind that gsl is a pre compiled C
|
||||
It might look rather complicated but keep in mind that gsl is a pre-compiled C
|
||||
library so it is quite an achievement that it's possible to get it working at
|
||||
all in the first place.
|
||||
[endsect]
|
||||
@ -337,7 +337,7 @@ state_type:
|
||||
[[Assign scalar multiplication] [`x *= a`] [`state_type`] [Performs in-place multiplication of vector x with scalar a.]]
|
||||
]
|
||||
|
||||
Definign these operators makes your state type work with any basic Runge-Kutta
|
||||
Defining these operators makes your state type work with any basic Runge-Kutta
|
||||
stepper.
|
||||
However, if you want to use step-size control, some more functionality is
|
||||
required.
|
||||
@ -351,7 +351,7 @@ So for controlled steppers the following things have to be implemented:
|
||||
|
||||
[table
|
||||
[[Name] [Expression] [Type] [Semantics]]
|
||||
[[Division] [`x / y`] [`state_type`] [Calculates the elementwise division 'x/y']]
|
||||
[[Division] [`x / y`] [`state_type`] [Calculates the element-wise division 'x/y']]
|
||||
[[Absolute value] [`abs( x )`] [`state_type`] [Element wise absolute value]]
|
||||
[[Reduce] [`vector_space_reduce_impl< state_type >::reduce( state , operator , init )`] [`value_type`]
|
||||
[Performs the operation `operator` for subsequently each element of `state` and returns the aggregate value.
|
||||
@ -370,7 +370,7 @@ So for controlled steppers the following things have to be implemented:
|
||||
[section Boost.Ublas]
|
||||
As an example for the employment of the `vector_space_algebra` we will adopt
|
||||
`ublas::vector` from __ublas to work as a state type in odeint.
|
||||
This is particularily easy because `ublas::vector` supports vector-vector
|
||||
This is particularly easy because `ublas::vector` supports vector-vector
|
||||
addition and scalar-vector multiplication described above as well as `boost::size`.
|
||||
It also has a resize member function so all that has to be done in this case
|
||||
is to declare resizability:
|
||||
@ -428,7 +428,7 @@ Again, note that the two last steps were only required if you want to use
|
||||
controlled steppers.
|
||||
For simple steppers definition of the simple `+=` and `*=` operators are
|
||||
sufficient.
|
||||
Having defined such a point type, we can easily perform the integration on a lorenz
|
||||
Having defined such a point type, we can easily perform the integration on a Lorenz
|
||||
system by using the `vector_space_algebra` again:
|
||||
|
||||
[point3D_main]
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
|
||||
[import ../examples/stepper_details.cpp]
|
||||
|
||||
Solving ordinary differential equation numerically is ususally done iteratively, that is a given state of an ordinary differential equation is iterated forward ['x(t) -> x(t+dt) -> x(t+2dt)]. The steppers in odeint perform one single step. The most general stepper type is described by the __stepper concept. The stepper concepts of odeint are described in detail in section __concepts, here we briefly present the mathematical and numerical details of the steppers. The __stepper has two versions of the `do_step` method, one with an in-place transform of the currrent state and one with an out-of-place transform:
|
||||
Solving ordinary differential equation numerically is usually done iteratively, that is a given state of an ordinary differential equation is iterated forward ['x(t) -> x(t+dt) -> x(t+2dt)]. The steppers in odeint perform one single step. The most general stepper type is described by the __stepper concept. The stepper concepts of odeint are described in detail in section __concepts, here we briefly present the mathematical and numerical details of the steppers. The __stepper has two versions of the `do_step` method, one with an in-place transform of the current state and one with an out-of-place transform:
|
||||
|
||||
`do_step( sys , inout , t , dt )`
|
||||
|
||||
@ -14,7 +14,7 @@ The first parameter is always the system function - a function describing the OD
|
||||
|
||||
Up to now, we have nothing said about the system function. This function depends on the stepper. For the explicit Runge-Kutta steppers this function can be a simple callable object hence a simple (global) C-function or a functor. The parameter syntax is `sys( x , dxdt , t )` and it is assumed that it calculates ['dx/dt = f(x,t)].
|
||||
|
||||
Other types of system function represent Hamiltonian systems or system which also compute the Jacobian needed in implicit steppers. For informations which stepper uses which system function see the stepper table below. It might be possible, that odeint will introduce new system types in near future. Since the system function is strongly related to the stepper type, such an introduction of a new stepper might result in a new type of system function.
|
||||
Other types of system function represent Hamiltonian systems or system which also compute the Jacobian needed in implicit steppers. For information which stepper uses which system function see the stepper table below. It might be possible, that odeint will introduce new system types in near future. Since the system function is strongly related to the stepper type, such an introduction of a new stepper might result in a new type of system function.
|
||||
|
||||
[section Explicit steppers]
|
||||
|
||||
@ -34,7 +34,7 @@ A special class of the explicit steppers are the FSAL (first-same-as-last) stepp
|
||||
|
||||
`do_step( sys , in , dxdtin , out , dxdtout , t , dt )`
|
||||
|
||||
This method also fills the derivative at time ['t+dt] into `dxdtout`. Of course, the performance gain of such FSAL steppers only appears when combining with intergrate error estimation, like in the Runge-Kutta-Dopri5 stepper. The FSAL-trick is sometimes also referred as the Fehlberg-Trick. An example how the FSAL steppers can be used is
|
||||
This method also fills the derivative at time ['t+dt] into `dxdtout`. Of course, the performance gain of such FSAL steppers only appears when combining with integrate error estimation, like in the Runge-Kutta-Dopri5 stepper. The FSAL-trick is sometimes also referred as the Fehlberg-Trick. An example how the FSAL steppers can be used is
|
||||
|
||||
[fsal_stepper_detail_example]
|
||||
|
||||
@ -45,9 +45,9 @@ This method also fills the derivative at time ['t+dt] into `dxdtout`. Of course,
|
||||
|
||||
[section Symplectic solvers]
|
||||
|
||||
As mentioned above symplectic solvers are used for Hamiltonian systems. Symplectic solvers conserve the phase space volume exactly and if the Hamiltonian system is energy conservative they also conserve the energy approximately. A special class of symplectic systems are separable systems which can be written in the form ['dqdt/dt = f1(p)], ['dpdt/dt = f2(q)], where ['(q,p)] are the state of system. The space of ['(q,p)] is sometimes refered as the phase space and ['q] and ['p] are said the be the phase space variables. Symplectic systems in this special form occur widely in nature. For example the complete classical mechanics as written down by Newton, Lagrange and Hamilton can be forumulated in this framework. Of course, the separability of the system depends on the specific choice of coordinates.
|
||||
As mentioned above symplectic solvers are used for Hamiltonian systems. Symplectic solvers conserve the phase space volume exactly and if the Hamiltonian system is energy conservative they also conserve the energy approximately. A special class of symplectic systems are separable systems which can be written in the form ['dqdt/dt = f1(p)], ['dpdt/dt = f2(q)], where ['(q,p)] are the state of system. The space of ['(q,p)] is sometimes referred as the phase space and ['q] and ['p] are said the be the phase space variables. Symplectic systems in this special form occur widely in nature. For example the complete classical mechanics as written down by Newton, Lagrange and Hamilton can be formulated in this framework. Of course, the separability of the system depends on the specific choice of coordinates.
|
||||
|
||||
Integrable symplectic systems can be solved by odeint by means of the symplectic_euler stepper and a symplectic Runge-Kutta-Nystrom method of sixth-order. These stepper assume that the system is autonomous, hence the time will not explicitly occur. Furhter they fullfil in principle the default Stepper concept, but they expect the system to be a pair of callable objects. The first entry of this pair calculates ['f1(p)] while the second calculates ['f2(q)]. The syntax is `sys.first(p,dqdt)` and `sys.second(q,dpdt)`, where the first and second part can be again simple C-functions of functors. An example is the harmonic oscillator:
|
||||
Integrable symplectic systems can be solved by odeint by means of the symplectic_euler stepper and a symplectic Runge-Kutta-Nystrom method of sixth-order. These stepper assume that the system is autonomous, hence the time will not explicitly occur. Further they fulfill in principle the default Stepper concept, but they expect the system to be a pair of callable objects. The first entry of this pair calculates ['f1(p)] while the second calculates ['f2(q)]. The syntax is `sys.first(p,dqdt)` and `sys.second(q,dpdt)`, where the first and second part can be again simple C-functions of functors. An example is the harmonic oscillator:
|
||||
|
||||
[symplectic_stepper_detail_system_function]
|
||||
|
||||
@ -61,13 +61,13 @@ If you like to represent the system with one class you can easily bind two publi
|
||||
|
||||
[symplectic_stepper_detail_system_class_example]
|
||||
|
||||
Many Hamiltonian system can be written as ['dq/dt=p], ['dp/dt=f(q)] which is computationally much easier then the full separable system. Very often, it is also possible to transform the original equations of motion to bring the system in this simplified form. This kind of system can be used in the symplectic solvers, by simply passing ['f(p)] to the `do_step` method, again ['f(p)] will be represented by a simple C-function or a functor. Here, the above example of the harmonic oscaillator can be written as
|
||||
Many Hamiltonian system can be written as ['dq/dt=p], ['dp/dt=f(q)] which is computationally much easier then the full separable system. Very often, it is also possible to transform the original equations of motion to bring the system in this simplified form. This kind of system can be used in the symplectic solvers, by simply passing ['f(p)] to the `do_step` method, again ['f(p)] will be represented by a simple C-function or a functor. Here, the above example of the harmonic oscillator can be written as
|
||||
|
||||
[simplified_symplectic_stepper_example]
|
||||
|
||||
In this example the function `harm_osc_f1` is exactly the same function as in the above examples.
|
||||
|
||||
Note, that the state of the ODE must not be contructed explicitly via `pair< vector_type , vector_type > x`. One can also use a combination of `make_pair` and `ref`. Furthermore, a convenience version of `do_step` exists which takes q and p without combining them into a pair:
|
||||
Note, that the state of the ODE must not be constructed explicitly via `pair< vector_type , vector_type > x`. One can also use a combination of `make_pair` and `ref`. Furthermore, a convenience version of `do_step` exists which takes q and p without combining them into a pair:
|
||||
|
||||
[symplectic_stepper_detail_ref_usage]
|
||||
|
||||
@ -77,7 +77,7 @@ Note, that the state of the ODE must not be contructed explicitly via `pair< vec
|
||||
|
||||
[caution This section is not up-to-date.]
|
||||
|
||||
For some kind of systems the stability properties of the classical Runge-Kutta are not sufficient, especially if the system is said to be stiff. A stiff system possesses two or more time scales of very different order. Solvers for stiff systems are usually implicit, meaning that they solve equations like ['x(t+dt) = x(t) + dt * f(x(t+1))]. This particular scheme is the implicit euler method. Implicit methods usually solve the system of equations by a root finding algorithm like the Newton method and therefore need to know the Jacobian of the system ['J[subl ij] = df[subl i] / dx[subl j]].
|
||||
For some kind of systems the stability properties of the classical Runge-Kutta are not sufficient, especially if the system is said to be stiff. A stiff system possesses two or more time scales of very different order. Solvers for stiff systems are usually implicit, meaning that they solve equations like ['x(t+dt) = x(t) + dt * f(x(t+1))]. This particular scheme is the implicit Euler method. Implicit methods usually solve the system of equations by a root finding algorithm like the Newton method and therefore need to know the Jacobian of the system ['J[subl ij] = df[subl i] / dx[subl j]].
|
||||
|
||||
For implicit solvers the system is again a pair, where the first component computes ['f(x,t)] and the second the Jacobian. The syntax is `sys.first( x , dxdt , t )` and `sys.second( x , J , t )`. For the implicit solver the `state_type` is `ublas::vector` and the Jacobian is represented by `ublas::matrix`.
|
||||
|
||||
@ -85,9 +85,9 @@ For implicit solvers the system is again a pair, where the first component compu
|
||||
|
||||
[section Multistep methods]
|
||||
|
||||
Another large class of solvers are multi-step method. They save a small part of the history of the solution and compute the next step with the help of this history. Since multistep methods know a part of their history they do not need to compute the system function very often, usually it is only computed once. This makes multistep methods preferable if a call of the system function is expensive. Examples are ODEs defined on networks, where the computation of the interaction is usually where expensive (and might be of order O(N^2)).
|
||||
Another large class of solvers are multi-step method. They save a small part of the history of the solution and compute the next step with the help of this history. Since multi-step methods know a part of their history they do not need to compute the system function very often, usually it is only computed once. This makes multi-step methods preferable if a call of the system function is expensive. Examples are ODEs defined on networks, where the computation of the interaction is usually where expensive (and might be of order O(N^2)).
|
||||
|
||||
Multistep methods differ from the normal steppers. They safe a part of their history and this part has to be explicitly calculated and initialized. In the following example an Adams-Bashforth-stepper with a history of 5 steps is instantiated and initialized;
|
||||
Multi-step methods differ from the normal steppers. They safe a part of their history and this part has to be explicitly calculated and initialized. In the following example an Adams-Bashforth-stepper with a history of 5 steps is instantiated and initialized;
|
||||
|
||||
[multistep_detail_example]
|
||||
|
||||
@ -95,25 +95,25 @@ The initialization uses a fourth-order Runge-Kutta stepper and after the call of
|
||||
|
||||
[multistep_detail_own_stepper_initialization]
|
||||
|
||||
Many multistep methods are also explicit steppers, hence the parameter of `do_step` method do not differ from the explicit steppers.
|
||||
Many multi-step methods are also explicit steppers, hence the parameter of `do_step` method do not differ from the explicit steppers.
|
||||
|
||||
[caution The multistep methods have some internal variables which depend on the explicit solution. Hence you can not exchange the system of the state between two consecutive calls of `do_step` since then the internal variable do not correspond with the ODE and the current solution. Of course, if you use the integrate functions this will be taken into account. See the __using_steppers section for more details.]
|
||||
[caution The multi-step methods have some internal variables which depend on the explicit solution. Hence you can not exchange the system of the state between two consecutive calls of `do_step` since then the internal variable do not correspond with the ODE and the current solution. Of course, if you use the integrate functions this will be taken into account. See the __using_steppers section for more details.]
|
||||
|
||||
|
||||
[endsect]
|
||||
|
||||
[section Controlled steppers]
|
||||
|
||||
Many of the above introduced steppers possess the possibility to use adaptive stepsize control. Adaptive step size integration works in principle as follows:
|
||||
Many of the above introduced steppers possess the possibility to use adaptive step-size control. Adaptive step size integration works in principle as follows:
|
||||
|
||||
# The error of one step is calculated. This is usually done by performing two steps with different orders. The difference between these two steps is then used as a measure for the error. Stepper which can calculate the error are __error_stepper and they form a own class with an separate concept.
|
||||
# This error is compared against some predefined error tolerances. Are the tolerance violated the step is reject and the stepsize is decreases. Otherwise the step is accepted and possibly the stepsize is increased.
|
||||
# This error is compared against some predefined error tolerances. Are the tolerance violated the step is reject and the step-size is decreases. Otherwise the step is accepted and possibly the step-size is increased.
|
||||
|
||||
The class of controlled steppers has its own concept in odeint - the __controlled_stepper concept. They are usually constructed from the underlying error steppers. An example is the controller for the explicit Runge-Kutta steppers. The Runge-Kutta steppers enter the controller as a template argument. Additionally one can pass the Runge-Kutta stepper to the contructor, but this step is not neccessary; the stepper is default-constructed if possible.
|
||||
The class of controlled steppers has its own concept in odeint - the __controlled_stepper concept. They are usually constructed from the underlying error steppers. An example is the controller for the explicit Runge-Kutta steppers. The Runge-Kutta steppers enter the controller as a template argument. Additionally one can pass the Runge-Kutta stepper to the constructor, but this step is not necessary; the stepper is default-constructed if possible.
|
||||
|
||||
Different step size controlling mechanism exist. They all have in common that they somehow compare predefined error tolerance against the error and that they might reject or accept a step. If a step is rejected the step size is usually decreased and the step is made again. Then the procedure of checking the error tolerances and accepting or rejecting a step is made again and repeated until the step is accepted. The procedure is implemented in the integration functions.
|
||||
|
||||
A classical way to decide wether a step is rejected or accepted is
|
||||
A classical way to decide whether a step is rejected or accepted is
|
||||
|
||||
['val = || | err[subl i] | / ( __epsilon[subl abs] + __epsilon[subl rel] * ( a[subl x] | x[subl i] | + a[subl dxdt] | | dxdt[subl i] | )|| ]
|
||||
|
||||
@ -131,7 +131,7 @@ Here, ['O[subl S]] and ['O[subl E]] are the order of the stepper and the error s
|
||||
|
||||
[include controlled_stepper_table.qbk]
|
||||
|
||||
To ease to generation of the controlled stepper generation functions exist which take the absolute and relative error tolerances and a predefined error stepper and contruct from this knowledge an appropirate controlled stepper. The generation functions are explained in detail in XYZ.
|
||||
To ease to generation of the controlled stepper generation functions exist which take the absolute and relative error tolerances and a predefined error stepper and construct from this knowledge an appropriate controlled stepper. The generation functions are explained in detail in __generation_functions.
|
||||
|
||||
[endsect]
|
||||
|
||||
@ -141,13 +141,13 @@ A fourth class of stepper exists which are the so called dense output steppers.
|
||||
|
||||
[dense_output_detail_example]
|
||||
|
||||
Dense output stepper have their own concept. The main difference is that they control the state by them-self. If you call `do_step`, only the ODE is passed as argument. Furhtermore `do_step` return the last time interval, hence you interpolate the solution between these two times points. Another difference is that they must be initialized with `initialize`, otherwise the internal state of the stepper is default contructed which might produce funny errors or bugs.
|
||||
Dense output stepper have their own concept. The main difference is that they control the state by them-self. If you call `do_step`, only the ODE is passed as argument. Furthermore `do_step` return the last time interval, hence you interpolate the solution between these two times points. Another difference is that they must be initialized with `initialize`, otherwise the internal state of the stepper is default constructed which might produce funny errors or bugs.
|
||||
|
||||
The construction of the dense output stepper looks a little bit nasty, since in the case of the `dense_output_runge_kutta` stepper a controlled stepper and an error stepper have to be nested. To simplify the generation of the dense output stepper generation functions exist:
|
||||
|
||||
[dense_output_detail_generation1]
|
||||
|
||||
Of course, this statement is also lengthly; it demonstrates how `make_dense_output` can be used with the `result_of` protocoll. The parameters to `make_dense_output` are the absolute error tolerance, the relative error tolerance and the stepper. This explicitly assumes that the underlying stepper is a controlled stepper and that this stepper has an absolute and a relative error tolerance. For details about the generation functions see __generation_functions. Of course, the generation functions have been designed for easy use with the integrate functions:
|
||||
Of course, this statement is also lengthy; it demonstrates how `make_dense_output` can be used with the `result_of` protocol. The parameters to `make_dense_output` are the absolute error tolerance, the relative error tolerance and the stepper. This explicitly assumes that the underlying stepper is a controlled stepper and that this stepper has an absolute and a relative error tolerance. For details about the generation functions see __generation_functions. Of course, the generation functions have been designed for easy use with the integrate functions:
|
||||
|
||||
[dense_output_detail_generation2]
|
||||
|
||||
@ -155,7 +155,7 @@ Of course, this statement is also lengthly; it demonstrates how `make_dense_outp
|
||||
|
||||
[section Using steppers]
|
||||
|
||||
This section contains some general informations about the usage of the steppers in odeint.
|
||||
This section contains some general information about the usage of the steppers in odeint.
|
||||
|
||||
[* Steppers are copied by value]
|
||||
|
||||
@ -165,11 +165,11 @@ The stepper in odeint are always copied by values. They are copied for the creat
|
||||
|
||||
[caution Some of the features described in this section are not yet implemented]
|
||||
|
||||
Some steppers require to store some informations about the state of the ODE between two steps. Examples are the multistep method which store a part of the solution during the evolution of the ODE, or the FSAL steppers which store the last derivative at time ['t+dt], to be used in the next step. In both cases the steppers expect that consecutive calls of `do_step` are from the same solution and the same ODE. In this case it is absolutely neccessary that you call `do_step` with the same system function and the same state, see also the examples for the FSAL steppers above.
|
||||
Some steppers require to store some information about the state of the ODE between two steps. Examples are the multi-step method which store a part of the solution during the evolution of the ODE, or the FSAL steppers which store the last derivative at time ['t+dt], to be used in the next step. In both cases the steppers expect that consecutive calls of `do_step` are from the same solution and the same ODE. In this case it is absolutely necessary that you call `do_step` with the same system function and the same state, see also the examples for the FSAL steppers above.
|
||||
|
||||
Stepper with an internal state support two additional methods: `reset` which resets the state and `initialize` which initializes the internal state. The parameters of `initialize` depend on the specific stepper. For example the Adams-Bashforth-Moulton stepper provides two initialize methods: `initialize( system , inout , t , dt )` which initializes the internal states with the help of the Runge-Kutta 4 stepper, and `initialize( stepper , system , inout , t , dt )` which initializes with the help of `stepper`. For the case of the FSAL steppers, `initialize` is `initialize( sys , in , t )` which simply calculates the r.h.s. of the ODE and assigns its value to the internal derivative.
|
||||
|
||||
All these stepper have in common, that they initially fill their internal state by themself. Hence you are not required to call initialize. See how this works for the Adams-Bashforth-Moulton stepper: in the example we instantiate a fourth order Adams-Bashforth-Moulton stepper, meaning that it will store 4 internal derivatives of the solution at times `(t-dt,t-2dt,t-3dt,t-4dt)`.
|
||||
All these stepper have in common, that they initially fill their internal state by themselves. Hence you are not required to call initialize. See how this works for the Adams-Bashforth-Moulton stepper: in the example we instantiate a fourth order Adams-Bashforth-Moulton stepper, meaning that it will store 4 internal derivatives of the solution at times `(t-dt,t-2dt,t-3dt,t-4dt)`.
|
||||
|
||||
``
|
||||
adams_bashforth_moulton< 4 , state_type > stepper;
|
||||
@ -185,11 +185,11 @@ In the stepper table at the bottom of this page one can see which stepper have a
|
||||
|
||||
[* Stepper might be resizable]
|
||||
|
||||
Nearly all steppers in odeint need to store some intermediate results of the type `state_type` or `deriv_type`. To do this the size of these temporaries needs to be adjusted. So, most steppers in odeint provide an additional template parameter which controls the size adjustment of the internal variables - the resizer. In detail odeint provides three policy classes (resizers) `always_resizer`, `initially_resizer`, and `never_resizer`. Furthermore, all stepper have a method `adjust_size` which takes a parameter representing a state type and which manually adjusts the size of the internal variables. Before performing the actual resizing odeint always checks if the sizes of the state and the interal variable differ and only resizes if they are different.
|
||||
Nearly all steppers in odeint need to store some intermediate results of the type `state_type` or `deriv_type`. To do this the size of these temporaries needs to be adjusted. So, most steppers in odeint provide an additional template parameter which controls the size adjustment of the internal variables - the resizer. In detail odeint provides three policy classes (resizers) `always_resizer`, `initially_resizer`, and `never_resizer`. Furthermore, all stepper have a method `adjust_size` which takes a parameter representing a state type and which manually adjusts the size of the internal variables. Before performing the actual resizing odeint always checks if the sizes of the state and the internal variable differ and only resizes if they are different.
|
||||
|
||||
By default the resizing parameter is `initially_resizer`, meaning that the first call to `do_step` performs the resizing. Of course, if you have changed the size of your system and your state you have to call `adjust_size` by hand. The second resizer is the `always_resizer` which tries to resize the internal variables at every call of `do_step`. Typical use cases for this kind of resizer are self expanding lattics like shown in the tutorial or partial differential equations with an adaptive grid. The third class of resizer is the `never_resizer` which means that the internal variables are never adjusted automatically and always have the adjust by hand.
|
||||
By default the resizing parameter is `initially_resizer`, meaning that the first call to `do_step` performs the resizing. Of course, if you have changed the size of your system and your state you have to call `adjust_size` by hand. The second resizer is the `always_resizer` which tries to resize the internal variables at every call of `do_step`. Typical use cases for this kind of resizer are self expanding lattices like shown in the tutorial or partial differential equations with an adaptive grid. The third class of resizer is the `never_resizer` which means that the internal variables are never adjusted automatically and always have the adjust by hand.
|
||||
|
||||
There is a second mechanism which influences the resizing and which controls if a state type is at least resizeable - a metafunction `is_resizeable`. This metafunction returns a static boolean value if any type is resizable. For example it will return `true` for `std::vector< T >` but `false` for `tr1::array< T >`. By default and for unknown types `is_resizeable` returns `false`, so if you have your own type you need to specialize this metafunction. For more details on the resizing mechanism see the section __adapt_state_types.
|
||||
There is a second mechanism which influences the resizing and which controls if a state type is at least resizeable - a meta-function `is_resizeable`. This meta-function returns a static boolean value if any type is resizable. For example it will return `true` for `std::vector< T >` but `false` for `tr1::array< T >`. By default and for unknown types `is_resizeable` returns `false`, so if you have your own type you need to specialize this meta-function. For more details on the resizing mechanism see the section __adapt_state_types.
|
||||
|
||||
|
||||
|
||||
@ -198,7 +198,7 @@ There is a second mechanism which influences the resizing and which controls if
|
||||
odeint provides a quite large number of different steppers such that the user is left with the question of which stepper fits his needs. Our personal recommendations are:
|
||||
|
||||
* `runge_kutta_dopri5` is maybe the best default stepper. It has step size control as well as dense-output functionality. Simple create a dense-output stepper by `make_dense_output( 1.0e-6 , 1.0e-5 , runge_kutta_dopri5< state_type >() )`.
|
||||
* `runge_kutta4` is a good stepper for constant step sizes. It is widely used and very well known. If you need to create artifical time series this stepper should be the first choice.
|
||||
* `runge_kutta4` is a good stepper for constant step sizes. It is widely used and very well known. If you need to create artificial time series this stepper should be the first choice.
|
||||
* 'runge_kutta_fehlberg78' is similar to the 'runge_kutta4' with the advantage that it has higher precision. It can also be used with step size control.
|
||||
* `adams_bashforth_moulton` is very well suited for ODEs where the r.h.s. is expensive (in terms of computation time). It will calculate the system function only once during each step.
|
||||
|
||||
@ -215,13 +215,13 @@ odeint provides a quite large number of different steppers such that the user is
|
||||
|
||||
[import ../examples/stochastic_euler.cpp]
|
||||
|
||||
Of course, one can write own steppers which are fully compatible with odeint. They only have to fullfil one or several of the stepper __concepts of odeint.
|
||||
Of course, one can write own steppers which are fully compatible with odeint. They only have to fulfill one or several of the stepper __concepts of odeint.
|
||||
|
||||
We will illustrate how to write your own stepper with the example of the stochastic euler method. This method is suited to solve stochastic differential equations (SDEs). A SDE has the form
|
||||
We will illustrate how to write your own stepper with the example of the stochastic Euler method. This method is suited to solve stochastic differential equations (SDEs). A SDE has the form
|
||||
|
||||
['dx/dt = f(x) + g(x) __xi(t)]
|
||||
|
||||
where ['__xi] is Gaussian white noise with zero mean and a standard deviation ['__sigma(t)]. ['f(x)] is said to be the deterministic part while [' g(x) __xi] is the noisy part. In case ['g(x)] is independent of ['x] the SDE is said to have additive noise. It is not possible to solve SDE with the classical solvers for ODEs since the noisy part of the SDE has to be scaled differently then the deterministic part with respect to the time step. BUt there exist many solvers for SDEs. A classical and easy method is the stochastic euler solver. It works by iterating
|
||||
where ['__xi] is Gaussian white noise with zero mean and a standard deviation ['__sigma(t)]. ['f(x)] is said to be the deterministic part while [' g(x) __xi] is the noisy part. In case ['g(x)] is independent of ['x] the SDE is said to have additive noise. It is not possible to solve SDE with the classical solvers for ODEs since the noisy part of the SDE has to be scaled differently then the deterministic part with respect to the time step. But there exist many solvers for SDEs. A classical and easy method is the stochastic Euler solver. It works by iterating
|
||||
|
||||
['x(t+__Delta t) = x(t) + __Delta t f(x(t)) + __Delta t[super 1/2] g(x) __xi(t)]
|
||||
|
||||
@ -231,13 +231,13 @@ Now we will implement this method. We will call the stepper `stochastic_euler`.
|
||||
|
||||
[stochastic_euler_class_definition]
|
||||
|
||||
The types are needed in order to fulfill the stepper concept. As internal state and deriv type we use simple arrays in the stochastic euler, they are needed for the temporaries. The stepper has the order one which is returned from the `order()` function.
|
||||
The types are needed in order to fulfill the stepper concept. As internal state and deriv type we use simple arrays in the stochastic Euler, they are needed for the temporaries. The stepper has the order one which is returned from the `order()` function.
|
||||
|
||||
The system functions needs to calculate the deterministic and the stochastic part of our stochastic differential equation. So it might be suitable that the system function itself is a pair of functions, one for computing the deterministic and one for the stochastic part. The first element of the pair simply computes the deterministic part while the second the stochastic one. Then, the second part also needs to calculate the random numbers in order to simulate the stochastic process. We can now implement the `do_step` method
|
||||
|
||||
[stochastic_euler_do_step]
|
||||
|
||||
This is all. It is quite simple and the stochastic euler stepper implement here is quite general. Of course it can be enhanced, for example
|
||||
This is all. It is quite simple and the stochastic Euler stepper implement here is quite general. Of course it can be enhanced, for example
|
||||
|
||||
* use of operations and algebras as well as the resizing mechanism for maximal flexibility and portability
|
||||
* use of `boost::ref` for the system functions
|
||||
@ -248,7 +248,7 @@ Now, lets look how we use the new stepper. A nice example is the Ornstein-Uhlenb
|
||||
|
||||
['dx/dt = - x + __xi]
|
||||
|
||||
where __xi is Gaussian white noise with standard deviation ['__sigma]. Implementing the Ornstein-Uhlenbeck process is quite simple. We need two functions or funtors - one for the deterministic and one for the stochastic part:
|
||||
where __xi is Gaussian white noise with standard deviation ['__sigma]. Implementing the Ornstein-Uhlenbeck process is quite simple. We need two functions or functors - one for the deterministic and one for the stochastic part:
|
||||
|
||||
[stochastic_euler_ornstein_uhlenbeck_def]
|
||||
|
||||
|
||||
@ -8,11 +8,11 @@
|
||||
[[[github_link libs/numeric/odeint/examples/fpu.cpp fpu.cpp]] [The Fermi-Pasta-Ulam (FPU) example shows how odeint can be used to integrate lattice systems.]]
|
||||
[[[github_link libs/numeric/odeint/examples/phase_oscillator_ensemble.cpp phase_oscillator_ensemble.cpp]] [The phase oscillator ensemble example shows how globally coupled oscillators can be analyzed and how statistical measures can be computed during integration.]]
|
||||
[[[github_link libs/numeric/odeint/examples/harmonic_oscillator_units.cpp harmonic_oscillator_units.cpp]] [This examples shows how __boost_units can be used with odeint.]]
|
||||
[[[github_link libs/numeric/odeint/examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp]] [The 2D phase oscillator example shows how a two-dimenstional lattice can used with odeint and how matrix types can be used as state types in odeint.]]
|
||||
[[[github_link libs/numeric/odeint/examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp]] [The 2D phase oscillator example shows how a two-dimensional lattice can used with odeint and how matrix types can be used as state types in odeint.]]
|
||||
[[[github_link libs/numeric/odeint/examples/lorenz_gmpxx.cpp lorenz_gmpxx.cpp]] [This examples integrates the Lorenz system by means of an arbitrary precision type.]]
|
||||
[[[github_link libs/numeric/odeint/examples/thrust/phase_oscillator_ensemble.cu phase_oscillator_ensemble.cu]] [The Thrust phase oscillator ensemble example shows how globally coupled oscillators can be analyzed with Thrust and CUDA, employing the power of modern graphic devices.]]
|
||||
[[[github_link libs/numeric/odeint/examples/thrust/phase_oscillator_chain.cu phase_oscillator_chain.cu]] [The Thrust phase oscillator chain example shows how chains of nearest neighbor coupled oscillators can be integrated with Thrust and odeint.]]
|
||||
[[[github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu lorenz_parameters.cu]] [The Lorenz paramaters examples show how ensembles of ordinary differential equations can be solved by means of Thrust to study the dependence of an ODE on some parameters.]]
|
||||
[[[github_link libs/numeric/odeint/examples/mtl/gauss_packet.cpp gauss_packet.cpp]] [The MTL-Gauss-packet example shows how the mtl can be easily used with odeint.]]
|
||||
[[[github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu lorenz_parameters.cu]] [The Lorenz parameters examples show how ensembles of ordinary differential equations can be solved by means of Thrust to study the dependence of an ODE on some parameters.]]
|
||||
[[[github_link libs/numeric/odeint/examples/mtl/gauss_packet.cpp gauss_packet.cpp]] [The MTL-Gauss-packet example shows how the MTL can be easily used with odeint.]]
|
||||
[[[github_link libs/numeric/odeint/examples/stochastic_euler.cpp stochastic_euler.cpp]] [Implementation of a custom stepper - the stochastic euler - for solving stochastic differential equations.]]
|
||||
]
|
||||
|
||||
@ -2,9 +2,9 @@
|
||||
|
||||
[section Overview]
|
||||
|
||||
[caution Boost.Odeint is not an official boost library!]
|
||||
[caution Boost.Numeric.Odeint is not an official boost library!]
|
||||
|
||||
Odeint is a library for solving initial value problems (IVP) of ordinary differential equations. Mathematically, these problems are formulated as follows: ['x'(t) = f(x,t)], ['x(0) = x0]. ['x] and ['f] can be vectors and the solution is some function ['x(t)] fullfilling both equations above. In the following we will refer to ['x'(t)] also `dxdt` which is also our notation for the derivative in the source code.
|
||||
odeint is a library for solving initial value problems (IVP) of ordinary differential equations. Mathematically, these problems are formulated as follows: ['x'(t) = f(x,t)], ['x(0) = x0]. ['x] and ['f] can be vectors and the solution is some function ['x(t)] fulfilling both equations above. In the following we will refer to ['x'(t)] also `dxdt` which is also our notation for the derivative in the source code.
|
||||
|
||||
Numerical approximations for the solution ['x(t)] are calculated iteratively. The easiest algorithm is the Euler-Scheme, where starting at ['x(0)] one finds ['x(dt) = x(0) + dt f(x(0),0)]. Now one can use ['x(dt)] and obtain ['x(2dt)] in a similar way and so on. The Euler method is of order 1, that means the error at each step is ['~ dt[super 2]]. This is, of course, not very satisfying, which is why the Euler method is rarely used for real life problems and serves just as illustrative example.
|
||||
|
||||
@ -24,7 +24,7 @@ In odeint, the following algorithms are implemented:
|
||||
|
||||
[include stepper_table.qbk]
|
||||
|
||||
Ordinary differential equations occur nearly everywhere in natural sciences. For example, the whole Newtonian mechanics are described by second order differential equations. Be sure, you will find them in every discipline. They also occur if partial differential equations (PDEs) are discretized in one coordinate. Then, a system of coupled ordinary differential occurs, sometimes also refered as lattices ODEs.
|
||||
Ordinary differential equations occur nearly everywhere in natural sciences. For example, the whole Newtonian mechanics are described by second order differential equations. Be sure, you will find them in every discipline. They also occur if partial differential equations (PDEs) are discretized in one coordinate. Then, a system of coupled ordinary differential occurs, sometimes also referred as lattices ODEs.
|
||||
|
||||
[endsect]
|
||||
|
||||
@ -32,7 +32,7 @@ Ordinary differential equations occur nearly everywhere in natural sciences. For
|
||||
|
||||
[section Usage, Compilation, Headers]
|
||||
|
||||
Odeint is completely header-only, meaning that you do not link against any library. It can be include by
|
||||
odeint is completely header-only, meaning that you do not link against any library. It can be include by
|
||||
|
||||
``
|
||||
#include <boost/numeric/odeint.hpp>
|
||||
@ -52,17 +52,17 @@ Imaging, you want to numerically integrate a harmonic oscillator with friction.
|
||||
[import ../examples/harmonic_oscillator.cpp]
|
||||
[rhs_function]
|
||||
|
||||
Here we chose `vector<double>` as the state type, but others are also possible, for example `tr1::array<double,2>`. Odeint is designed in such a way that you can easily use your own state types. Next, we define the ODE which is in this case a simple function. The parameter signature of this function is crucial: the integration methods will always call them in the form `f(x, dxdt, t)`. So, even if there is no explicit time dependence, one has to define `t` as a function parameter.
|
||||
Here we chose `vector<double>` as the state type, but others are also possible, for example `tr1::array<double,2>`. odeint is designed in such a way that you can easily use your own state types. Next, we define the ODE which is in this case a simple function. The parameter signature of this function is crucial: the integration methods will always call them in the form `f(x, dxdt, t)`. So, even if there is no explicit time dependence, one has to define `t` as a function parameter.
|
||||
|
||||
Now, we have to define the initial state from which the integration should start:
|
||||
|
||||
[state_initialization]
|
||||
|
||||
For the integration itself we'll use the [funcref boost::numeric::odeint::integrate integrate] function, which is a convenient way to get quick results. It is based on the error-controlled [classref boost::numeric::odeint::runge_kutta_rk5_ck runge_kutta_rk5_ck] stepper (5th order) and uses adaptive stepsize.
|
||||
For the integration itself we'll use the [funcref boost::numeric::odeint::integrate integrate] function, which is a convenient way to get quick results. It is based on the error-controlled [classref boost::numeric::odeint::runge_kutta_rk5_ck runge_kutta_rk5_ck] stepper (5th order) and uses adaptive step-size.
|
||||
|
||||
[integration]
|
||||
|
||||
The integrate function expects as parameters the rhs of the ode as defined above, the initial state `x`, the start-and end-time of the integration as well as the initial time step. Note, that [funcref boost::numeric::odeint::integrate integrate] uses an adaptive stepsize during the integration steps so the time points will not be equally spaced. The integration returns the number of steps that were applied.
|
||||
The integrate function expects as parameters the rhs of the ode as defined above, the initial state `x`, the start-and end-time of the integration as well as the initial time step. Note, that [funcref boost::numeric::odeint::integrate integrate] uses an adaptive step-size during the integration steps so the time points will not be equally spaced. The integration returns the number of steps that were applied.
|
||||
|
||||
It is, of course, also possible to implement the ode system as a class. The rhs must then be implemented as a functor having defined the ()-operator:
|
||||
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
[library boost.sandbox.numeric.odeint
|
||||
[library Boost.Numeric.Odeint
|
||||
[quickbook 1.5]
|
||||
[id odeint]
|
||||
[dirname odeint]
|
||||
@ -69,7 +69,7 @@
|
||||
[def __tut_chaotic_system
|
||||
[link boost_sandbox_numeric_odeint.tutorial.chaotic_systems_and_lyapunov_exponents Chaotic system]]
|
||||
[def __using_steppers
|
||||
[link boost_sandbox_numeric_odeint.odeint_in_detail.steppers.using_steppers Using steppers]]
|
||||
[link boost_sandbox_numeric_odeint.odeint_in_detail.steppers.using_steppers Using steppers]]
|
||||
[def __generation_functions
|
||||
[link boost_sandbox_numeric_odeint.odeint_in_detail.generation_functions Generation functions]]
|
||||
[def __adapt_state_types
|
||||
@ -107,7 +107,7 @@
|
||||
[template subl[x]'''<subscript>'''__space[x]'''</subscript>''']
|
||||
|
||||
[template github_link[url text]'''<ulink url="https://github.com/headmyshoulder/odeint-v2/tree/master/'''[url]'''" target="_blank">'''[text]'''</ulink>''']
|
||||
# [template github_link[url text]'''<ulink url="'''[url]'''" target="_blank">'''[text]'''</ulink>''']
|
||||
[/ [template github_link[url text]'''<ulink url="'''[url]'''" target="_blank">'''[text]'''</ulink>''']]
|
||||
|
||||
[include getting_started.qbk]
|
||||
|
||||
@ -117,11 +117,7 @@
|
||||
|
||||
[include concepts.qbk]
|
||||
|
||||
[xinclude reference.xml]
|
||||
|
||||
[include concepts_old.qbk]
|
||||
|
||||
[include reference_old.qbk]
|
||||
[/ [xinclude reference.xml] ]
|
||||
|
||||
|
||||
[/
|
||||
|
||||
@ -50,7 +50,7 @@ Algebras supporting __boost_range
|
||||
[[symplectic_rkn_sb3a_mclachlan]]
|
||||
]
|
||||
|
||||
[table Algebra supporting Boost.Range
|
||||
[table Algebras supporting Boost.Range
|
||||
[[algebra]]
|
||||
[[range_algebra]]
|
||||
[[thrust_algebra]]
|
||||
|
||||
@ -1,11 +0,0 @@
|
||||
[section Integration functions]
|
||||
|
||||
[section Constant step-size functions]
|
||||
|
||||
[endsect]
|
||||
|
||||
[section Adaptive step-size functions]
|
||||
|
||||
[endsect]
|
||||
|
||||
[endsect]
|
||||
@ -1,31 +0,0 @@
|
||||
[section Old Reference]
|
||||
|
||||
[include reference_steppers.qbk]
|
||||
|
||||
[include reference_integrate_functions.qbk]
|
||||
|
||||
[section Algebras]
|
||||
|
||||
[endsect]
|
||||
|
||||
[section Operations]
|
||||
|
||||
[endsect]
|
||||
|
||||
[section Resizing]
|
||||
|
||||
[endsect]
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
[endsect]
|
||||
@ -1,118 +0,0 @@
|
||||
[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–Kutta_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–Karp_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]
|
||||
|
||||
@ -4,7 +4,7 @@
|
||||
[[Modified Midpoint] [`modified_midpoint`] [__stepper] [__system] [configurable (2)] [No] [No] [No] [Used in Bulirsch-Stoer implementation]]
|
||||
[[Runge-Kutta 4] [`runge_kutta4`] [__stepper] [__system] [4] [No] [No] [No] [The classical Runge Kutta scheme, good general scheme without error control]]
|
||||
[[Cash-Karp] [`runge_kutta_cash_karp54`] [__error_stepper] [__system] [5] [Yes (4)] [No] [No] [Good general scheme with error estimation, to be used in controlled_error_stepper]]
|
||||
[[Dormand-Prince 5] [`runge_kutta_dopri5`] [__error_stepper] [__system] [5] [Yes (4)] [Yes] [Yes] [Standard method with error control and dense ouput, to be used in controlled_error_stepper and in dense_output_controlled_explicit_fsal.]]
|
||||
[[Dormand-Prince 5] [`runge_kutta_dopri5`] [__error_stepper] [__system] [5] [Yes (4)] [Yes] [Yes] [Standard method with error control and dense output, to be used in controlled_error_stepper and in dense_output_controlled_explicit_fsal.]]
|
||||
[[Fehlberg 78] [`runge_kutta_fehlberg78`] [__error_stepper] [__system] [8] [Yes (7)] [No] [No] [Good high order method with error estimation, to be used in controlled_error_stepper.]]
|
||||
|
||||
[[Adams Bashforth] [`adams_bashforth`] [__stepper] [__system] [configurable] [No] [No] [Yes] [Multistep method]]
|
||||
@ -12,15 +12,15 @@
|
||||
[[Adams Bashforth Moulton] [`adams_bashforth_moulton`] [__stepper] [__system] [configurable] [No] [No] [Yes] [Combined multistep method]]
|
||||
|
||||
[[Controlled Runge Kutta] [`controlled_runge_kutta`] [__controlled_stepper] [__system] [depends] [Yes] [No] [depends] [Error control for __error_stepper. Requires an __error_stepper from above. Order depends on the given ErrorStepper]]
|
||||
[[Dense Output Runge Kutta] [`dense_output_runge_kutta`] [__dense_output_stepper] [__system] [depends] [No] [Yes] [Yes] [Dense ouput for __stepper and __error_stepper from above if they provide dense ouput functionality (like `euler` and `runge_kutta_dopri5`). Order depends on the given stepper.]]
|
||||
[[Dense Output Runge Kutta] [`dense_output_runge_kutta`] [__dense_output_stepper] [__system] [depends] [No] [Yes] [Yes] [Dense output for __stepper and __error_stepper from above if they provide dense output functionality (like `euler` and `runge_kutta_dopri5`). Order depends on the given stepper.]]
|
||||
|
||||
[[Bulirsch-Stoer] [`bulirsch_stoer`] [__controlled_stepper] [__system] [variable] [Yes] [No] [No] [Stepper with step size and order control. Very good if high precision is required.]]
|
||||
[[Bulirsch-Stoer Dense Output] [`bulirsch_stoer_dense_out`] [__dense_output_stepper] [__system] [variable] [Yes] [Yes] [No] [Stepper with step size and order control as well as dense ouput. Very good if high precision and dense ouput is required.]]
|
||||
[[Bulirsch-Stoer Dense Output] [`bulirsch_stoer_dense_out`] [__dense_output_stepper] [__system] [variable] [Yes] [Yes] [No] [Stepper with step size and order control as well as dense output. Very good if high precision and dense output is required.]]
|
||||
|
||||
[[Implicit Euler] [`implicit_euler`] [__stepper] [__implicit_system] [1] [No] [No] [No] [Basic implicit routine. Requires the Jacobian. Works only with __ublas vectors as state types.]]
|
||||
[[Rosenbrock 4] [`rosenbrock4`] [__error_stepper] [__implicit_system] [4] [Yes] [Yes] [No] [Good for stiff systems. Works only with __ublas vectors as state types.]]
|
||||
[[Controlled Rosenbrock 4] [`rosenbrock4_controller`] [__controlled_stepper] [__implicit_system] [4] [Yes] [Yes] [No] [Rosenbrock 4 with error control. Works only with __ublas vectors as state types.]]
|
||||
[[Dense Ouput Rosenbrock 4] [`rosenbrock4_dense_ouput`] [__dense_output_stepper] [__implicit_system] [4] [Yes] [Yes] [No] [Controlled Rosenbrock 4 with dense output. Works only with __ublas vectors as state types.]]
|
||||
[[Dense Output Rosenbrock 4] [`rosenbrock4_dense_output`] [__dense_output_stepper] [__implicit_system] [4] [Yes] [Yes] [No] [Controlled Rosenbrock 4 with dense output. Works only with __ublas vectors as state types.]]
|
||||
|
||||
[[Symplectic Euler] [`symplectic_euler`] [__stepper] [__symplectic_system __simple_symplectic_system] [1] [No] [No] [No] [Basic symplectic solver for separable Hamiltonian system]]
|
||||
[[Symplectic RKN McLachlan] [`symplectic_rkn_sb3a_mclachlan`] [__stepper] [__symplectic_system __simple_symplectic_system] [6] [No] [No] [No] [Symplectic solver for separable Hamiltonian system with order 6]]
|
||||
|
||||
@ -25,7 +25,7 @@ The following table gives an overview over all examples.
|
||||
|
||||
[section References]
|
||||
|
||||
[*General informations about numerical integration of ordinary differential equations:]
|
||||
[*General information 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,13 +2,13 @@
|
||||
|
||||
[import ../examples/chaotic_system.cpp]
|
||||
|
||||
Odeint can easily be used to investigate the properties of chaotic deterministic systems. In mathematical terms chaotic refers to an exponential growth of perturbations. In order to observe this exponential growth one usually solves the equations for the tangential dynamics which is again an ordinary differential equation. These equations are linear but time dependent and can be obtained via
|
||||
odeint can easily be used to investigate the properties of chaotic deterministic systems. In mathematical terms chaotic refers to an exponential growth of perturbations. In order to observe this exponential growth one usually solves the equations for the tangential dynamics which is again an ordinary differential equation. These equations are linear but time dependent and can be obtained via
|
||||
|
||||
['d __delta x / dt = J(x) __delta x]
|
||||
|
||||
where ['J] is the Jacobian of the system under consideration. ['__delta x] cam also be interpreted as a perturbation of the original system. In principle ['n] of these perturbations exit, they form a hypercube and evolve in the time. The Lyapunov exponents are then defined as logarithmic growth rates of the perturbations. If one Lyapunov exponent is larger then zero nearby trajectories divergy exponentially hence they are chaotic. If the largest Lyapunov exponent is zero one is usually faced with periodic motion. In the case of a largest Lyapunov exponent smaller then zero the converges the a fixed point. More informations about Lyapunov exponents and nonlinear dynamical systems can be found in many textbooks, see for example .
|
||||
where ['J] is the Jacobian of the system under consideration. ['__delta x] cam also be interpreted as a perturbation of the original system. In principle ['n] of these perturbations exit, they form a hypercube and evolve in the time. The Lyapunov exponents are then defined as logarithmic growth rates of the perturbations. If one Lyapunov exponent is larger then zero nearby trajectories diverge exponentially hence they are chaotic. If the largest Lyapunov exponent is zero one is usually faced with periodic motion. In the case of a largest Lyapunov exponent smaller then zero the converges the a fixed point. More information's about Lyapunov exponents and nonlinear dynamical systems can be found in many textbooks, see for example .
|
||||
|
||||
To calculate the Lyapunov exponents numerically one usually solves the equations of motion for ['n] perturbations und orthonormalized them every ['k] steps. The Lyapunov exponent is the average the logarithm of the stretching factor of each perturbation.
|
||||
To calculate the Lyapunov exponents numerically one usually solves the equations of motion for ['n] perturbations and orthonormalizes them every ['k] steps. The Lyapunov exponent is the average the logarithm of the stretching factor of each perturbation.
|
||||
|
||||
To demonstrate how one can use odeint to determine the Lyapunov exponents we choose the Lorenz system. It is one of the most studied dynamical systems in the nonlinear dynamics community. For the standard parameters is possesses a strange attractor with non-integer dimension. The Lyapunov exponents take values of approximately 0.9, 0 and -12.
|
||||
|
||||
@ -28,7 +28,7 @@ void lorenz( const lorenz_state_type &x , lorenz_state_type &dxdt , double t )
|
||||
dxdt[2] = -b * x[2] + x[0] * x[1];
|
||||
}
|
||||
``
|
||||
We need also need to integrate the set of the perturbations. This is done in parallel to the original system, hence within one system function. Of course, we want to use the above defintion of the Lorenz system, hence the definition of the system function including the Lorenz system itself and the perturbation could look like:
|
||||
We need also need to integrate the set of the perturbations. This is done in parallel to the original system, hence within one system function. Of course, we want to use the above definition of the Lorenz system, hence the definition of the system function including the Lorenz system itself and the perturbation could look like:
|
||||
|
||||
``
|
||||
const size_t n = 3;
|
||||
@ -92,7 +92,7 @@ state_type x;
|
||||
explicit_rk4< state_type > rk4;
|
||||
integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 );
|
||||
``
|
||||
The problem is now, that `x` is the full state containing also the perturbations and `integrate_n_steps` does not know that it should only use 3 element. In detail, odeint and its steppers determine the length of the system under consideration be determing the length of the state. In the classical solvers from the problem was solved by pointer to the state and an appropriate length, something similar to
|
||||
The problem is now, that `x` is the full state containing also the perturbations and `integrate_n_steps` does not know that it should only use 3 element. In detail, odeint and its steppers determine the length of the system under consideration be determining the length of the state. In the classical solvers from the problem was solved by pointer to the state and an appropriate length, something similar to
|
||||
|
||||
``
|
||||
void lorenz( double* x , double *dxdt , double t, void* params )
|
||||
@ -114,7 +114,7 @@ This is in principle all. Now, we only have to call `integrate_n_steps` with an
|
||||
|
||||
Having integrated a sufficient number of transients steps we are now able to calculate the Lyapunov exponents:
|
||||
|
||||
# We initialize the perturbations. They are stored linearly behind the state of the Lorenz system. The perturbation are initialized such that they fullfill
|
||||
# We initialize the perturbations. They are stored linearly behind the state of the Lorenz system. The perturbation are initialized such that they fulfill
|
||||
[' < p [subl i] , p[subl j] > = __delta [subl ij] ]
|
||||
where ['<,>] is the classical scalar product and ['__delta [subl ij]] is the Kronecker symbol.
|
||||
# Integrate 100 steps of the full system with perturbations
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
|
||||
[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 fullfils some requirements defined below.
|
||||
First of all, you have to specify the data type 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 fulfills some requirements defined below.
|
||||
|
||||
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:
|
||||
|
||||
@ -20,7 +20,7 @@ odeint can deal with instances of such classes instead of pure functions which a
|
||||
|
||||
[section Stepper Types]
|
||||
|
||||
Numerical integration works iteratively, that means you start at a state ['x(t)] and perform 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[super q]] that means the error in ['x] made by such a step is ['~dt[super q+1]]. odeint provides several steppers of different orders from which you can choose:
|
||||
Numerical integration works iteratively, that means you start at a state ['x(t)] and perform a time-step of length ['dt] to obtain the approximate state ['x(t+dt)]. There exist many different methods to perform such a time-step each of which has a certain order ['q]. If the order of a method is ['q] than it is accurate up to term ['~dt[super q]] that means the error in ['x] made by such a step is ['~dt[super q+1]]. odeint provides several steppers of different orders from which you can choose:
|
||||
|
||||
[include stepper_table.qbk]
|
||||
|
||||
@ -30,11 +30,11 @@ Some of steppers in the table above are special: Some need the Jacobian of the O
|
||||
|
||||
[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. 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:
|
||||
The basic stepper just performs one time-step 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) because 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. 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:
|
||||
|
||||
[define_const_stepper]
|
||||
|
||||
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). Each stepper defines a `do_step` method which can used directly. So, you write down the above example as
|
||||
This call integrates the system defined by `harmonic_oscillator` using the RK4 method from ['t=0] to ['10] with a step-size ['dt=0.01] and the initial condition given in `x`. The result, ['x(t=10)] is stored in `x` (in-place). Each stepper defines a `do_step` method which can used directly. So, you write down the above example as
|
||||
|
||||
[integrate_const_loop]
|
||||
|
||||
@ -42,7 +42,7 @@ This call integrates the system defined by `harmonic_oscillator` using the RK4 m
|
||||
|
||||
[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 `explicit_error_rk54_ck? -- a 5th order Runge-Kutta method with 4th order error estimation and coefficients introduced by Cash-Karp.
|
||||
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 `explicit_error_rk54_ck? -- a 5th order Runge-Kutta method with 4th order error estimation and coefficients introduced by Cash-Karp.
|
||||
|
||||
[define_adapt_stepper]
|
||||
|
||||
@ -52,7 +52,7 @@ Given the error stepper, one still needs an instance that checks the error and a
|
||||
|
||||
As above, this integrates the system defined by `harmonic_oscillator` using an adaptive step size method based on the Runge-Kutta Cash-Karp 54 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).
|
||||
|
||||
In the above example an error stepper is nested in a controlled stepper. This is a nice technique; however one drawback is that one always needs to define both steppers. Of course, one could also write the instantiation of the controlled stepper into the call of the integrate function but a complete knowledge of the underlying stepper types is still neccessary. Another point is, that the error tolerances for the step size control are not easily included into the controlled stepper. Both issues can be solved by using `make_controlled`:
|
||||
In the above example an error stepper is nested in a controlled stepper. This is a nice technique; however one drawback is that one always needs to define both steppers. Of course, one could also write the instantiation of the controlled stepper into the call of the integrate function but a complete knowledge of the underlying stepper types is still necessary. Another point is, that the error tolerances for the step size control are not easily included into the controlled stepper. Both issues can be solved by using `make_controlled`:
|
||||
|
||||
[integrate_adapt_make_controlled]
|
||||
|
||||
@ -66,7 +66,7 @@ For the Runge Kutta controller the error made during one step is compared with [
|
||||
|
||||
When using `make_controlled` the parameter ['a[sub x]] and ['a[sub dxdt]] are used with their standard values of 1.
|
||||
|
||||
In the tables below, one can find all steppers which are working with `make_controlled` and `make_dense_output` which is the analogon for the dense output steppers.
|
||||
In the tables below, one can find all steppers which are working with `make_controlled` and `make_dense_output` which is the analog for the dense output steppers.
|
||||
|
||||
[include make_controlled_table.qbk]
|
||||
|
||||
|
||||
@ -27,7 +27,7 @@ equations of motion
|
||||
|
||||
['dp[subl i] / dt = -dH / dq[subl i]]
|
||||
|
||||
In time independent Hamiltonian system the energy and the phase space volume are conserved and special integration methods have to be applied in order to ensure these conservation laws. The odeint library provides classes for separable Hamiltonian systems, which can be written in the form ['H = __Sigma p[subl i][super 2] / (2m[subl i]) + H[subl q](q)], where ['H[subl q](q)] only depends on the coordinates. 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 ['dq[subl i] / dt = p[subl i]] / m[subl i] , ['dp[subl i] / dt = f( q[subl i] )].
|
||||
In time independent Hamiltonian system the energy and the phase space volume are conserved and special integration methods have to be applied in order to ensure these conservation laws. The odeint library provides classes for separable Hamiltonian systems, which can be written in the form ['H = __Sigma p[subl i][super 2] / (2m[subl i]) + H[subl q](q)], where ['H[subl q](q)] only depends on the coordinates. Although 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 ['dq[subl i] / dt = p[subl i]] / m[subl i] , ['dp[subl i] / dt = f( q[subl i] )].
|
||||
|
||||
[endsect]
|
||||
|
||||
@ -45,7 +45,7 @@ The next step is to define a container type storing the values of ['q] and ['p]
|
||||
[import ../examples/solar_system.cpp]
|
||||
[container_type_definition]
|
||||
|
||||
The `container_type` is different from the state type of the ODE. The state type of the ode is simply a `pair< container_type , container_type >` since it needs the informations about the coordinates and the momenta.
|
||||
The `container_type` is different from the state type of the ODE. The state type of the ode is simply a `pair< container_type , container_type >` since it needs the information about the coordinates and the momenta.
|
||||
|
||||
As system function we have to provide ['f(p)] and ['f(q)]:
|
||||
|
||||
@ -53,7 +53,7 @@ As system function we have to provide ['f(p)] and ['f(q)]:
|
||||
|
||||
[momentum_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.
|
||||
In general a three body-system is chaotic, hence we can not expect that arbitrary 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.
|
||||
|
||||
As mentioned above, we need to use some special integrators in order to conserve phase space volume and energy. There is a well known family of such integrators, the so-called Runge-Kutta-Nystroem solvers, which we apply here:
|
||||
|
||||
|
||||
@ -1,7 +1,5 @@
|
||||
[section Special topics]
|
||||
|
||||
explain boost::ref
|
||||
|
||||
[section Complex state types]
|
||||
|
||||
[import ../examples/stuart_landau.cpp]
|
||||
@ -14,7 +12,7 @@ where ['__Psi] and ['i] is a complex variable. Such systems can easily be repres
|
||||
|
||||
[stuart_landau_system_function]
|
||||
|
||||
Of courst, one can also use classical functions to implement the state function. In this cast the Stuart-Landau oscillator looks like
|
||||
Of course, one can also use classical functions to implement the state function. In this cast the Stuart-Landau oscillator looks like
|
||||
|
||||
[stuart_landau_system_function_alternative]
|
||||
|
||||
@ -39,7 +37,7 @@ odeint can also be used to solved ordinary differential equations defined on lat
|
||||
|
||||
Remarkably, the Fermi-Pasta-Ulam system was the first numerical experiment which has been implemented on a computer in 1953. It was studied at Los Alamos on one of the first computer (a MANIAC I) and it triggered a whole new tree of mathematical and physical science.
|
||||
|
||||
The FPU is solved again by a symplectic solver, but in our case we can speed up the computation because the ['q] components trivially reduce to ['dq[subl i] / dt = p[subl i]]. odeint is capable of doing this peformance improvement. All you have to do is to call the symplectic solver with an state function for the ['p] components. Here, is how this function looks like
|
||||
The FPU is solved again by a symplectic solver, but in our case we can speed up the computation because the ['q] components trivially reduce to ['dq[subl i] / dt = p[subl i]]. odeint is capable of doing this performance improvement. All you have to do is to call the symplectic solver with an state function for the ['p] components. Here, is how this function looks like
|
||||
|
||||
[fpu_system_function]
|
||||
|
||||
@ -65,7 +63,7 @@ Another import high dimensional system of coupled ordinary differential equation
|
||||
|
||||
[' d__phi[subl k] / dt = __omega[subl k] + __epsilon / N __Sigma[subl j] sin( __phi[subl j] - __phi[subl k] )]
|
||||
|
||||
The natural frequencies ['__omega[subl i]] of each oscialltor follow some distribution and ['__epsilon] is the coupling strength. We choose here a Lorentzian distribution for ['__omega[subl i]]. Interestingly a phase transition can be observed if the coupling strenght exceeds a critical value. Above this value synchronization sets in and some of the oscillators oscillate with the same frequency despite their different natural frequencies. The transition is also called Kuramoto transition. Its behavior can be analyzed by employing the mean field of the phase
|
||||
The natural frequencies ['__omega[subl i]] of each oscillator follow some distribution and ['__epsilon] is the coupling strength. We choose here a Lorentzian distribution for ['__omega[subl i]]. Interestingly a phase transition can be observed if the coupling strength exceeds a critical value. Above this value synchronization sets in and some of the oscillators oscillate with the same frequency despite their different natural frequencies. The transition is also called Kuramoto transition. Its behavior can be analyzed by employing the mean field of the phase
|
||||
|
||||
['Z = K e[super i __Theta] = 1 / N __Sigma[subl k]e[super i __phi[subl k]]]
|
||||
|
||||
@ -95,7 +93,7 @@ To illustrate how odeint works with __boost_units we use the harmonic oscillator
|
||||
|
||||
[units_define_basic_quantities]
|
||||
|
||||
Note, that the `state_type` and the `deriv_type` are now a compile-time fusion sequences. They have to be explicitly defined. Next, we define the ordinary differential equation which is completety equivalent to the example in the first section of this tutorial
|
||||
Note, that the `state_type` and the `deriv_type` are now a compile-time fusion sequences. They have to be explicitly defined. Next, we define the ordinary differential equation which is completely equivalent to the example in the first section of this tutorial
|
||||
|
||||
[units_define_ode]
|
||||
|
||||
@ -121,7 +119,7 @@ The full cpp file for this example can be found here [github_link libs/numeric/o
|
||||
|
||||
odeint works well with a variety of different state types. It is not restricted to pure vector-wise types, like `vector< double >`, `array< double , N >`, `fusion::vector< double , double >`, etc. but also works with types having a different topology then simple vectors. Here, we show how odeint can be used with matrices as states type, in the next section we will show how can be used to solve ODEs defined on complex networks.
|
||||
|
||||
By default, odeint can be used with `ublas::matrix< T >` as state type for matrices. A simple example is a two-dimensional lattice of coupled phase oscillators. We like phas oscillators, they are extremly easy and might serve for different demonstration purposes. Other matrix types like `mtl::dense_matrix` or blitz arrays and matrices can used as well but need some kind of activation in order to work with odeint. This activation is described in following sections,
|
||||
By default, odeint can be used with `ublas::matrix< T >` as state type for matrices. A simple example is a two-dimensional lattice of coupled phase oscillators. We like phase oscillators, they are extremely easy and might serve for different demonstration purposes. Other matrix types like `mtl::dense_matrix` or blitz arrays and matrices can used as well but need some kind of activation in order to work with odeint. This activation is described in following sections,
|
||||
|
||||
The definition of the system is
|
||||
|
||||
@ -155,7 +153,7 @@ Ginzburg-Landau
|
||||
|
||||
[import ../examples/gmpxx/lorenz_gmpxx.cpp]
|
||||
|
||||
Besides the classical floating point number like `float`, `double`, `complex< double >` you can also use arbitrary precision types, like the types from [@http://gmplib.org/ gmp] and [@http://www.mpfr.org/ mpfr]. But you have to be carful about instantiating any numbers.
|
||||
Besides the classical floating point number like `float`, `double`, `complex< double >` you can also use arbitrary precision types, like the types from [@http://gmplib.org/ gmp] and [@http://www.mpfr.org/ mpfr]. But you have to be careful about instantiating any numbers.
|
||||
|
||||
For gmp types you have to set the default precision before any number is instantiated. This can be easily done by calling `mpf_set_default_prec( precision )` as the first function in your main program. Secondly, you can not use any global constant variables since they will not be set with the default precision you have already set.
|
||||
|
||||
@ -178,11 +176,11 @@ The full example can be found at [github_link libs/numeric/odeint/examples/loren
|
||||
|
||||
[import ../examples/resizing_lattice.cpp]
|
||||
|
||||
Odeint supports changes of the state size during integration if a state_type is used which can be resized, like `std::vector`.
|
||||
odeint supports changes of the state size during integration if a state_type is used which can be resized, like `std::vector`.
|
||||
The adjustment of the state type's size has to be done from outside and the stepper has to be instantiated with `always_resizer` as the template argument for the `resizer_type`.
|
||||
In this configuration, the stepper checks for changes in the state size and adjust it's internal storage accordingly.
|
||||
|
||||
We exemplarily show this for a Hamiltonian system of nonlinear, disordered oscillators with nonlinear nearest neighbour coupling.
|
||||
We exemplary show this for a Hamiltonian system of nonlinear, disordered oscillators with nonlinear nearest neighbor coupling.
|
||||
|
||||
The system function is implemented in terms of a class that also provides functions for calculating the energy.
|
||||
Note, that this class stores the random potential internally which is not resized, but rather a start index is kept which should be changed whenever the states' size change.
|
||||
|
||||
@ -2,13 +2,11 @@
|
||||
|
||||
[import ../examples/stiff_system.cpp]
|
||||
|
||||
An important class of ordinary differential equations are so called stiff system which are characterized by two or more time scales of different order. Examples of such systems are found in chemical systems where reaction rates of individual subreaction might differ over large ranges.
|
||||
An important class of ordinary differential equations are so called stiff system which are characterized by two or more time scales of different order. Examples of such systems are found in chemical systems where reaction rates of individual sub-reaction might differ over large ranges.
|
||||
|
||||
CHECK:
|
||||
['d S[subl 1] / dt = - 101 S[subl 2] - 100 S[subl 1]]
|
||||
|
||||
['d S[subl 1] / dt = - __gamma[subl 1] S[subl 1] S[subl 2]]
|
||||
|
||||
['d S[subl 2] / dt = - __gamma[subl 1] S[subl 1] S[subl 2]]
|
||||
['d S[subl 2] / dt = S[subl 1]]
|
||||
|
||||
|
||||
To solve stiff systems numerically the Jacobian
|
||||
|
||||
@ -14,7 +14,7 @@ Typical applications for CUDA and odeint are large systems, like lattices or dis
|
||||
|
||||
The first example is the phase oscillator ensemble from the previous section. It has a phase transition at ['__epsilon = 2] in the limit of infinite numbers of oscillators ['N]. In the case of finite ['N] this transition is smeared out but still clearly visible.
|
||||
|
||||
__thrust and Cuda are perfectly suited for such kinds of problems where one needs a large number of particles (oscillators). We start by defining the state type which is a `thrust::device_vector`. The content of this vector lives on the GPU. If you are not familar with this we recommend reading the ['Getting started] section on the __thrust website.
|
||||
__thrust and Cuda are perfectly suited for such kinds of problems where one needs a large number of particles (oscillators). We start by defining the state type which is a `thrust::device_vector`. The content of this vector lives on the GPU. If you are not familiar with this we recommend reading the ['Getting started] section on the __thrust website.
|
||||
|
||||
[thrust_phase_ensemble_state_type]
|
||||
|
||||
@ -84,12 +84,12 @@ The next example is a large, one-dimensional chain of nearest-neighbor coupled p
|
||||
|
||||
['d __phi[subl k] / dt = __omega[subl k] + sin( __phi[subl k+1] - __phi[subl k] ) + sin( __phi[subl k] - __phi[subl k-1])]
|
||||
|
||||
In principle we can use all the techniques from the previous phase oscillator ensemble example, but we have to take special care about the coupling of the oscillators. To efficiently implement the coupling you can use a very elegant way employing Thrust's permutation iterator. A permuation iterator behaves like a normal iterator on a vector but it does not iterate along the usual order of the elements.
|
||||
In principle we can use all the techniques from the previous phase oscillator ensemble example, but we have to take special care about the coupling of the oscillators. To efficiently implement the coupling you can use a very elegant way employing Thrust's permutation iterator. A permutation iterator behaves like a normal iterator on a vector but it does not iterate along the usual order of the elements.
|
||||
It rather iterates along some permutation of the elements defined by some index map. To realize the nearest neighbor coupling we create one permutation iterator which travels one step behind a usual iterator and another permutation iterator which travels one step in front. The full system class is:
|
||||
|
||||
[thrust_phase_chain_system]
|
||||
|
||||
Note, how easy you can obtain the value for the left and right neighboring oscillator in the system functor using the permutation iterators. But, the call of the `thrust::for_each` function looks relativ complicated. Every term of the r.h.s. of the ODE is resembled by one iterator packed in exactly the same way as it is unpacked in the functor above.
|
||||
Note, how easy you can obtain the value for the left and right neighboring oscillator in the system functor using the permutation iterators. But, the call of the `thrust::for_each` function looks relatively complicated. Every term of the r.h.s. of the ODE is resembled by one iterator packed in exactly the same way as it is unpacked in the functor above.
|
||||
|
||||
Now we put everything together. We create random initial conditions and decreasing frequencies such that we should get synchronization. We copy the frequencies and the initial conditions onto the device and finally initialize and perform the integration. As result we simply write out the current state, hence the phase of each oscillator.
|
||||
|
||||
@ -103,10 +103,10 @@ The full example can be found at [github_link libs/numeric/odeint/examples/thrus
|
||||
|
||||
[import ../examples/thrust/lorenz_parameters.cu]
|
||||
|
||||
Another important use case for __thrust and Cuda are parameter studies of relatively small systems. Consider for example the three-dimensional Lorenz system from the chaotic systems example in the previous section which has three parameters. If you want to study the behaviour of this system for different parameters you usually have to integrate the system for many parameter values. Using thrust and odeint you can do this integration in parallel, hence you integrate a whole ensemble of Lorenz systems where each individual realization has a different parameter value.
|
||||
Another important use case for __thrust and Cuda are parameter studies of relatively small systems. Consider for example the three-dimensional Lorenz system from the chaotic systems example in the previous section which has three parameters. If you want to study the behavior of this system for different parameters you usually have to integrate the system for many parameter values. Using thrust and odeint you can do this integration in parallel, hence you integrate a whole ensemble of Lorenz systems where each individual realization has a different parameter value.
|
||||
[/ The Lorenz system is dissipative, such that you can assume that different initial conditions will lead to the same attractor so . For Hamiltonian systems this is not the case. Here it might be interesting to study a range of initial conditions to quantify different regions in the phase space.]
|
||||
|
||||
In the following we will show how you can use __thrust to integrate the above mentioned ensemble of Lorenz systems. We will vary only the parameter ['__beta] but it is straightfoward to vary other parameters or even two or all three parameters. Furthermore, we will use the largest Lyapunov exponent to quantify the behavior of the system (chaoticity).
|
||||
In the following we will show how you can use __thrust to integrate the above mentioned ensemble of Lorenz systems. We will vary only the parameter ['__beta] but it is straightforward to vary other parameters or even two or all three parameters. Furthermore, we will use the largest Lyapunov exponent to quantify the behavior of the system (chaoticity).
|
||||
|
||||
We start be defining the range of the parameters we want to study. Of course, the state_type is again a `thrust::device_vector< value_type >`.
|
||||
|
||||
@ -116,9 +116,9 @@ The next thing we have to implement is the Lorenz system without perturbations.
|
||||
|
||||
[thrust_lorenz_parameters_define_simple_system]
|
||||
|
||||
As `state_type` a `thrust::device_vector` or a __boost_range of a `device_vector` is used. The length of the state is ['3N] where ['N] is the number of systems. The system is encoded into this vector such that all ['x] components come first, then every ['y] compoents and finally every ['z] components. Implementing the device function is then a simple task, you only have to decompose the tuple originating from the zip iterators.
|
||||
As `state_type` a `thrust::device_vector` or a __boost_range of a `device_vector` is used. The length of the state is ['3N] where ['N] is the number of systems. The system is encoded into this vector such that all ['x] components come first, then every ['y] components and finally every ['z] components. Implementing the device function is then a simple task, you only have to decompose the tuple originating from the zip iterators.
|
||||
|
||||
Besides the system without perturbations we furthermore need to calculate the system including linearized equations governing the time evolution of small perturbations. Using the method from above this is straightforward, with a small difficulty that thrust's tuples have a maximal arity of 10. But this is only a small problem since we can create a zip iterator packed with zip iterators. So the top level zip iterator contains one zip iterator for the state, one normal iterator for the parameter, and one zip iterator for the derivative. Accessing the elements of this tuple in the system function is then straightforward, you unpack the tuple with `thrust::get<>()`. We will not show the code here, it is to large. It can be found [github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu here] and is easy to understand.
|
||||
Besides the system without perturbations we furthermore need to calculate the system including linearized equations governing the time evolution of small perturbations. Using the method from above this is straightforward, with a small difficulty that Thrust's tuples have a maximal arity of 10. But this is only a small problem since we can create a zip iterator packed with zip iterators. So the top level zip iterator contains one zip iterator for the state, one normal iterator for the parameter, and one zip iterator for the derivative. Accessing the elements of this tuple in the system function is then straightforward, you unpack the tuple with `thrust::get<>()`. We will not show the code here, it is to large. It can be found [github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu here] and is easy to understand.
|
||||
|
||||
Furthermore, we need an observer which determines the norm of the perturbations, normalizes them and averages the logarithm of the norm. The device functor which is used within this observer is defined
|
||||
|
||||
@ -126,7 +126,7 @@ Furthermore, we need an observer which determines the norm of the perturbations,
|
||||
|
||||
Note, that this functor manipulates the state, i.e. the perturbations.
|
||||
|
||||
Now we complete the whole code to calculate the Lyapunov exponents. First, we have to define a state vector. This vector contains ['6N] entries, the state ['x,y,z] and its perturbations ['dx,dy,dz]. We initialize them such that ['x=y=z=10], ['dx=1], and ['dy=dz=0]. We define a stepper type, a controlled Runge-Kutta Dormand-Prince 5 stepper. We start with some integration to overcome the transient behavior. For this, we do not involve the pertubation and run the algorithm only on the state ['x,y,z] without any observer. Note, how __boost_range is used for partial integration of the state vector without pertubations (the first half of the whole state). After the transient, the full system with perturbations is integrated and the Lyapunov exponents are calculated and written to `stdout`.
|
||||
Now we complete the whole code to calculate the Lyapunov exponents. First, we have to define a state vector. This vector contains ['6N] entries, the state ['x,y,z] and its perturbations ['dx,dy,dz]. We initialize them such that ['x=y=z=10], ['dx=1], and ['dy=dz=0]. We define a stepper type, a controlled Runge-Kutta Dormand-Prince 5 stepper. We start with some integration to overcome the transient behavior. For this, we do not involve the perturbation and run the algorithm only on the state ['x,y,z] without any observer. Note, how __boost_range is used for partial integration of the state vector without perturbations (the first half of the whole state). After the transient, the full system with perturbations is integrated and the Lyapunov exponents are calculated and written to `stdout`.
|
||||
|
||||
[thrust_lorenz_parameters_integration]
|
||||
|
||||
|
||||
@ -45,7 +45,7 @@ exe rt_generic_rk4_lorenz
|
||||
;
|
||||
|
||||
exe gsl_rk4_lorenz
|
||||
: gsl_rk4_lorenz.cpp libgsl libgslcblas
|
||||
: gsl_rk4_lorenz.cpp libgslcblas libgsl
|
||||
;
|
||||
|
||||
exe odeint_rk4_phase_lattice
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user