\(\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}} }\)

OdeGearControl

View page source

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{eqnarray} X(ti) & = & xi \\ X'(t) & = & f[t , X(t)] \end{eqnarray}

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/utility/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 . 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.

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{eqnarray} 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{eqnarray}

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 .