Prev Next Index-> contents reference index search external Up-> CppAD utility OdeGearControl 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 OdeGearControl-> ode_gear_control.cpp Headings-> Syntax Purpose Include Notation xf Fun ---..t ---..x ---..f ---..f_x ---..Warning M ti tf xi smin smax sini eabs erel ef maxabs nstep Error Criteria Discussion Scalar Vector Example Theory Source Code

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.