lines 9-396 of file: include/cppad/utility/ode_gear.hpp

{xrst_begin OdeGear}
{xrst_spell
   affine
   dep
   interpolating
   iterates
   pp
   test test
}

An Arbitrary Order Gear Method
##############################

Syntax
******

   # ``include <cppad/utility/ode_gear.hpp>``

``OdeGear`` ( *F* , *m* , *n* , *T* , *X* , *e* )

Purpose
*******
This routine applies
:ref:`OdeGear@Gear's Method`
to solve an explicit set of ordinary differential equations.
We are given
:math:`f : \B{R} \times \B{R}^n \rightarrow \B{R}^n` be a smooth function.
This routine solves the following initial value problem

.. math::
   :nowrap:

   \begin{eqnarray}
      x( t_{m-1} )  & = & x^0    \\
      x^\prime (t)  & = & f[t , x(t)]
   \end{eqnarray}

for the value of :math:`x( t_m )`.
If your set of  ordinary differential equations are not stiff
an explicit method may be better (perhaps :ref:`Runge45-name` .)

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

Fun
***
The class *Fun*
and the object *F* satisfy the prototype

   *Fun* & *F*

This must support the following set of calls

| |tab| *F* . ``Ode`` ( *t* , *x* , *f* )
| |tab| *F* . ``Ode_dep`` ( *t* , *x* , *f_x* )

t
=
The argument *t* has prototype

   ``const`` *Scalar* & *t*

(see description of :ref:`OdeGear@Scalar` below).

x
=
The argument *x* has prototype

   ``const`` *Vector* & *x*

and has size *n*
(see description of :ref:`OdeGear@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 :math:`f(t, x)`
(see *f* ( *t* , *x* ) in :ref:`OdeGear@Purpose` ).

f_x
===
The argument *f_x* has prototype

   *Vector* & *f_x*

On input and output, *f_x* is a vector of size :math:`n * n`
and the input values of the elements of *f_x* do not matter.
On output,

.. math::

   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 (highest power of :math:`t`)
used to represent the function :math:`x(t)` in the multi-step method.
Upon return from ``OdeGear`` ,
the *i*-th component of the polynomial is defined by

.. math::

   p_i ( t_j ) = X[ j * n + i ]

for :math:`j = 0 , \ldots , m` (where :math:`0 \leq i < n`).
The value of :math:`m` must be greater than or equal one.

n
*
The argument *n* has prototype

   ``size_t`` *n*

It specifies the range space dimension of the
vector valued function :math:`x(t)`.

T
*
The argument *T* has prototype

   ``const`` *Vector* & *T*

and size greater than or equal to :math:`m+1`.
For :math:`j = 0 , \ldots m`, :math:`T[j]` is the time
corresponding to time corresponding
to a previous point in the multi-step method.
The value :math:`T[m]` is the time
of the next point in the multi-step method.
The array :math:`T` must be monotone increasing; i.e.,
:math:`T[j] < T[j+1]`.
Above and below we often use the shorthand :math:`t_j` for :math:`T[j]`.

X
*
The argument *X* has the prototype

   *Vector* & *X*

and size greater than or equal to :math:`(m+1) * n`.
On input to ``OdeGear`` ,
for :math:`j = 0 , \ldots , m-1`, and
:math:`i = 0 , \ldots , n-1`

.. math::

   X[ j * n + i ] = x_i ( t_j )

Upon return from ``OdeGear`` ,
for :math:`i = 0 , \ldots , n-1`

.. math::

   X[ m * n + i ] \approx x_i ( t_m )

e
*
The vector *e* is an approximate error bound for the result; i.e.,

.. math::

   e[i] \geq | X[ m * n + i ] - x_i ( t_m ) |

The order of this approximation is one less than the order of
the solution; i.e.,

.. math::

   e = O ( h^m )

where :math:`h` is the maximum of :math:`t_{j+1} - t_j`.

Scalar
******
The type *Scalar* must satisfy the conditions
for a :ref:`NumericType-name` type.
The routine :ref:`CheckNumericType-name` 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* :

.. list-table::
   :widths: auto

   * - **Operation**
     - **Description**
   * - *a* < *b*
     - less than operator (returns a ``bool`` object)

Vector
******
The type *Vector* must be a :ref:`SimpleVector-name` class with
:ref:`elements of type Scalar<SimpleVector@Elements of Specified Type>` .
The routine :ref:`CheckSimpleVector-name` will generate an error message
if this is not the case.

Example
*******
{xrst_toc_hidden
   example/utility/ode_gear.cpp
}
The file
:ref:`ode_gear.cpp-name`
contains an example and test a test of using this routine.

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

Theory
******
For this discussion we use the shorthand :math:`x_j`
for the value :math:`x ( t_j ) \in \B{R}^n` which is not to be confused
with :math:`x_i (t) \in \B{R}` in the notation above.
The interpolating polynomial :math:`p(t)` is given by

.. math::

   p(t) =
   \sum_{j=0}^m
   x_j
   \frac{
      \prod_{i \neq j} ( t - t_i )
   }{
      \prod_{i \neq j} ( t_j - t_i )
   }

The derivative :math:`p^\prime (t)` is given by

.. math::

   p^\prime (t) =
   \sum_{j=0}^m
   x_j
   \frac{
      \sum_{i \neq j} \prod_{k \neq i,j} ( t - t_k )
   }{
      \prod_{k \neq j} ( t_j - t_k )
   }

Evaluating the derivative at the point :math:`t_\ell` we have

.. math::
   :nowrap:

   \begin{eqnarray}
   p^\prime ( t_\ell ) & = &
   x_\ell
   \frac{
      \sum_{i \neq \ell} \prod_{k \neq i,\ell} ( t_\ell - t_k )
   }{
      \prod_{k \neq \ell} ( t_\ell - t_k )
   }
   +
   \sum_{j \neq \ell}
   x_j
   \frac{
      \sum_{i \neq j} \prod_{k \neq i,j} ( t_\ell - t_k )
   }{
      \prod_{k \neq j} ( t_j - t_k )
   }
   \\
   & = &
   x_\ell
   \sum_{i \neq \ell}
   \frac{ 1 }{ t_\ell - t_i }
   +
   \sum_{j \neq \ell}
   x_j
   \frac{
      \prod_{k \neq \ell,j} ( t_\ell - t_k )
   }{
      \prod_{k \neq j} ( t_j - t_k )
   }
   \\
   & = &
   x_\ell
   \sum_{k \neq \ell} ( t_\ell - t_k )^{-1}
   +
   \sum_{j \neq \ell}
   x_j
   ( t_j - t_\ell )^{-1}
   \prod_{k \neq \ell ,j} ( t_\ell - t_k ) / ( t_j - t_k )
   \end{eqnarray}

We define the vector :math:`\alpha \in \B{R}^{m+1}` by

.. math::

   \alpha_j = \left\{ \begin{array}{ll}
   \sum_{k \neq m} ( t_m - t_k )^{-1}
      & {\rm if} \; j = m
   \\
   ( t_j - t_m )^{-1}
   \prod_{k \neq m,j} ( t_m - t_k ) / ( t_j - t_k )
      & {\rm otherwise}
   \end{array} \right.

It follows that

.. math::

   p^\prime ( t_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m

Gear's method determines :math:`x_m` by solving the following
nonlinear equation

.. math::

   f( t_m , x_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m

Newton's method for solving this equation determines iterates,
which we denote by :math:`x_m^k`, by solving the following affine
approximation of the equation above

.. math::
   :nowrap:

   \begin{eqnarray}
   f( t_m , x_m^{k-1} ) + \partial_x f( t_m , x_m^{k-1} ) ( x_m^k - x_m^{k-1} )
   & = &
   \alpha_0 x_0^k + \alpha_1 x_1 + \cdots + \alpha_m x_m
   \\
   \left[ \alpha_m I - \partial_x f( t_m , x_m^{k-1} ) \right]  x_m
   & = &
   \left[
   f( t_m , x_m^{k-1} ) - \partial_x f( t_m , x_m^{k-1} ) x_m^{k-1}
   - \alpha_0 x_0 - \cdots - \alpha_{m-1} x_{m-1}
   \right]
   \end{eqnarray}

In order to initialize Newton's method; i.e. choose :math:`x_m^0`
we define the vector :math:`\beta \in \B{R}^{m+1}` by

.. math::

   \beta_j = \left\{ \begin{array}{ll}
   \sum_{k \neq m-1} ( t_{m-1} - t_k )^{-1}
      & {\rm if} \; j = m-1
   \\
   ( t_j - t_{m-1} )^{-1}
   \prod_{k \neq m-1,j} ( t_{m-1} - t_k ) / ( t_j - t_k )
      & {\rm otherwise}
   \end{array} \right.

It follows that

.. math::

   p^\prime ( t_{m-1} ) = \beta_0 x_0 + \cdots + \beta_m x_m

We solve the following approximation of the equation above to determine
:math:`x_m^0`:

.. math::

   f( t_{m-1} , x_{m-1} ) =
   \beta_0 x_0 + \cdots + \beta_{m-1} x_{m-1} + \beta_m x_m^0

Gear's Method
*************
C. W. Gear,
*Simultaneous Numerical Solution of Differential-Algebraic Equations* ,
IEEE Transactions on Circuit Theory,
vol. 18, no. 1, pp. 89-95, Jan. 1971.

{xrst_end OdeGear}
