Prev

Next

Index-> contents reference index search external

Up-> CppAD utility OdeErrControl

CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix

utility-> ErrorHandler NearEqual speed_test SpeedTest time_test test_boolofvoid NumericType CheckNumericType SimpleVector CheckSimpleVector nan pow_int Poly LuDetAndSolve RombergOne RombergMul Runge45 Rosen34 OdeErrControl OdeGear OdeGearControl CppAD_vector thread_alloc index_sort to_string set_union sparse_rc sparse_rcv

OdeErrControl-> ode_err_control.cpp ode_err_maxabs.cpp

Headings-> Syntax Description Include Notation xf Method ---..step ---..Nan ---..order ti tf xi smin smax scur eabs erel ef maxabs nstep Error Criteria Discussion Scalar Vector Example Theory Source Code

@(@\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }@)@An Error Controller for ODE Solvers
Syntax
# include <cppad/utility/ode_err_control.hpp>
xf = OdeErrControl(method, ti, tf, xi,
     smin, smax, scur, eabs, erel, ef , maxabs, nstep )

Description
Let @(@ \B{R} @)@ denote the real numbers and let @(@ F : \B{R} \times \B{R}^n \rightarrow \B{R}^n @)@ be a smooth function. We define @(@ X : [ti , tf] \rightarrow \B{R}^n @)@ by the following initial value problem: @[@ \begin{array}{rcl} X(ti) & = & xi \\ X'(t) & = & F[t , X(t)] \end{array} @]@ The routine OdeErrControl can be used to adjust the step size used an arbitrary integration methods in order to be as fast as possible and still with in a requested error bound.

Include
The file cppad/ode_err_control.hpp is included by cppad/cppad.hpp but it can also be included separately with out the rest of the CppAD routines.

Notation
The template parameter types Scalar and Vector are documented below.

xf
The return value xf has the prototype
     Vector xf
(see description of Vector below). and the size of xf is equal to n . If xf contains not a number nan , see the discussion of step .

Method
The class Method and the object method satisfy the following syntax
     Method &method
The object method must support step and order member functions defined below:

step
The syntax
     method.step(ta, tb, xa, xb, eb)
executes one step of the integration method.

ta
The argument ta has prototype
     const Scalar &ta
It specifies the initial time for this step in the ODE integration. (see description of Scalar below).

tb
The argument tb has prototype
     const Scalar &tb
It specifies the final time for this step in the ODE integration.

xa
The argument xa has prototype
     const Vector &xa
and size n . It specifies the value of @(@ X(ta) @)@. (see description of Vector below).

xb
The argument value xb has prototype
     Vector &xb
and size n . The input value of its elements does not matter. On output, it contains the approximation for @(@ X(tb) @)@ that the method obtains.

eb
The argument value eb has prototype
     Vector &eb
and size n . The input value of its elements does not matter. On output, it contains an estimate for the error in the approximation xb . It is assumed (locally) that the error bound in this approximation nearly equal to @(@ K (tb - ta)^m @)@ where K is a fixed constant and m is the corresponding argument to CodeControl.

Nan
If any element of the vector eb or xb are not a number nan, the current step is considered to large. If this happens with the current step size equal to smin , OdeErrControl returns with xf and ef as vectors of nan.

order
If m is size_t, the object method must also support the following syntax
     m = method.order()
The return value m is the order of the error estimate; i.e., there is a constant K such that if @(@ ti \leq ta \leq tb \leq tf @)@, @[@ | eb(tb) | \leq K | tb - ta |^m @]@ where ta , tb , and eb are as in method.step(ta, tb, xa, xb, eb)

ti
The argument ti has prototype
     const Scalar &ti
It specifies the initial time for the integration of the differential equation.

tf
The argument tf has prototype
     const Scalar &tf
It specifies the final time for the integration of the differential equation.

xi
The argument xi has prototype
     const Vector &xi
and size n . It specifies value of @(@ X(ti) @)@.

smin
The argument smin has prototype
     const Scalar &smin
The step size during a call to method is defined as the corresponding value of @(@ tb - ta @)@. If @(@ tf - ti \leq smin @)@, the integration will be done in one step of size tf - ti . Otherwise, the minimum value of tb - ta will be @(@ smin @)@ except for the last two calls to method where it may be as small as @(@ smin / 2 @)@.

smax
The argument smax has prototype
     const Scalar &smax
It specifies the maximum step size to use during the integration; i.e., the maximum value for @(@ tb - ta @)@ in a call to method . The value of smax must be greater than or equal smin .

scur
The argument scur has prototype
     Scalar &scur
The value of scur is the suggested next step size, based on error criteria, to try in the next call to method . On input it corresponds to the first call to method , in this call to OdeErrControl (where @(@ ta = ti @)@). On output it corresponds to the next call to method , in a subsequent call to OdeErrControl (where ta = tf ).

eabs
The argument eabs has prototype
     const Vector &eabs
and size n . Each of the elements of eabs must be greater than or equal zero. It specifies a bound for the absolute error in the return value xf as an approximation for @(@ X(tf) @)@. (see the error criteria discussion below).

erel
The argument erel has prototype
     const Scalar &erel
and is greater than or equal zero. It specifies a bound for the relative error in the return value xf as an approximation for @(@ X(tf) @)@ (see the error criteria discussion below).

ef
The argument value ef has prototype
     Vector &ef
and size n . The input value of its elements does not matter. On output, it contains an estimated bound for the absolute error in the approximation xf ; i.e., @[@ ef_i > | X( tf )_i - xf_i | @]@ If on output ef contains not a number nan, see the discussion of step .

maxabs
The argument maxabs is optional in the call to OdeErrControl. If it is present, it has the prototype
     Vector &maxabs
and size n . The input value of its elements does not matter. On output, it contains an estimate for the maximum absolute value of @(@ X(t) @)@; i.e., @[@ maxabs[i] \approx \max \left\{ | X( t )_i | \; : \; t \in [ti, tf] \right\} @]@

nstep
The argument nstep is optional in the call to OdeErrControl. If it is present, it has the prototype
     size_t &nstep
Its input value does not matter and its output value is the number of calls to method.step used by OdeErrControl.

Error Criteria Discussion
The relative error criteria erel and absolute error criteria eabs are enforced during each step of the integration of the ordinary differential equations. In addition, they are inversely scaled by the step size so that the total error bound is less than the sum of the error bounds. To be specific, if @(@ \tilde{X} (t) @)@ is the approximate solution at time @(@ t @)@, ta is the initial step time, and tb is the final step time, @[@ \left| \tilde{X} (tb)_j - X (tb)_j \right| \leq \frac{tf - ti}{tb - ta} \left[ eabs[j] + erel \; | \tilde{X} (tb)_j | \right] @]@ If @(@ X(tb)_j @)@ is near zero for some @(@ tb \in [ti , tf] @)@, and one uses an absolute error criteria @(@ eabs[j] @)@ of zero, the error criteria above will force OdeErrControl to use step sizes equal to smin for steps ending near @(@ tb @)@. In this case, the error relative to maxabs can be judged after OdeErrControl returns. If ef is to large relative to maxabs , OdeErrControl can be called again with a smaller value of smin .

Scalar
The type Scalar must satisfy the conditions for a NumericType type. The routine CheckNumericType will generate an error message if this is not the case. In addition, the following operations must be defined for Scalar objects a and b :

Operation	Description
a <= b	returns true (false) if a is less than or equal (greater than) b .
a == b	returns true (false) if a is equal to b .
log(a)	returns a Scalar equal to the logarithm of a
exp(a)	returns a Scalar equal to the exponential of a

Vector
The type Vector must be a SimpleVector class with elements of type Scalar . The routine CheckSimpleVector will generate an error message if this is not the case.

Example
The files ode_err_control.cpp and ode_err_maxabs.cpp contain examples and tests of using this routine. They return true if they succeed and false otherwise.

Theory
Let @(@ e(s) @)@ be the error as a function of the step size @(@ s @)@ and suppose that there is a constant @(@ K @)@ such that @(@ e(s) = K s^m @)@. Let @(@ a @)@ be our error bound. Given the value of @(@ e(s) @)@, a step of size @(@ \lambda s @)@ would be ok provided that @[@ \begin{array}{rcl} a & \geq & e( \lambda s ) (tf - ti) / ( \lambda s ) \\ a & \geq & K \lambda^m s^m (tf - ti) / ( \lambda s ) \\ a & \geq & \lambda^{m-1} s^{m-1} (tf - ti) e(s) / s^m \\ a & \geq & \lambda^{m-1} (tf - ti) e(s) / s \\ \lambda^{m-1} & \leq & \frac{a}{e(s)} \frac{s}{tf - ti} \end{array} @]@ Thus if the right hand side of the last inequality is greater than or equal to one, the step of size @(@ s @)@ is ok.

Source Code
The source code for this routine is in the file cppad/ode_err_control.hpp.

---

Input File: cppad/utility/ode_err_control.hpp