Prev Next

@(@\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 Gear's Ode Solvers

Syntax
# include <cppad/utility/ode_gear_control.hpp>
 xf = OdeGearControl(F, M, ti, tf, xi,
     smin, smax, sini, eabs, erel, ef , maxabs, nstep )

Purpose
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 OdeGear is a stiff multi-step method that can be used to approximate the solution to this equation. The routine OdeGearControl sets up this multi-step method and controls the error during such an approximation.

Include
The file cppad/ode_gear_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
 and the size of xf is equal to n (see description of Vector below). It is the approximation for @(@ X(tf) @)@.

Fun
The class Fun and the object F satisfy the prototype
     Fun &F
 This must support the following set of calls
     F.Ode(t, x, f)
     F.Ode_dep(t, x, f_x)
t
The argument t has prototype
     const Scalar &t
 (see description of Scalar below).

x
The argument x has prototype
     const Vector &x
 and has size N (see description of Vector below).

f
The argument f to F.Ode has prototype
     Vector &f
 On input and output, f is a vector of size N and the input values of the elements of f do not matter. On output, f is set equal to @(@ f(t, x) @)@ (see f(t, x) in Purpose ).

f_x
The argument f_x has prototype
     Vector &f_x
 On input and output, f_x is a vector of size @(@ N * N @)@ and the input values of the elements of f_x do not matter. On output, @[@ f\_x [i * n + j] = \partial_{x(j)} f_i ( t , x ) @]@

Warning
The arguments f , and f_x must have a call by reference in their prototypes; i.e., do not forget the & in the prototype for f and f_x .

M
The argument M has prototype
     size_t M
 It specifies the order of the multi-step method; i.e., the order of the approximating polynomial (after the initialization process). The argument M must greater than or equal one.

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 minimum value of @(@ T[M] - T[M-1] @)@ in a call to OdeGear will be @(@ smin @)@ except for the last two calls where it may be as small as @(@ smin / 2 @)@. The value of smin must be less than or equal smax .

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 @(@ T[M] - T[M-1] @)@ in a call to OdeGear.

sini
The argument sini has prototype
     Scalar &sini
 The value of sini is the minimum step size to use during initialization of the multi-step method; i.e., for calls to OdeGear where @(@ m < M @)@. The value of sini must be less than or equal smax (and can also be less than smin ).

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 | @]@

maxabs
The argument maxabs is optional in the call to OdeGearControl. 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 has the prototype
     size_t &nstep
 Its input value does not matter and its output value is the number of calls to OdeGear used by OdeGearControl.

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 OdeGearControl to use step sizes equal to smin for steps ending near @(@ tb @)@. In this case, the error relative to maxabs can be judged after OdeGearControl returns. If ef is to large relative to maxabs , OdeGearControl 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 file ode_gear_control.cpp contains an example and test a test of using this routine. It returns true if it succeeds 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_gear_control.hpp.
Input File: cppad/utility/ode_gear_control.hpp