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}} }@)@This is cppad-20221105 documentation. Here is a link to its current documentation .
An Arbitrary Order Gear Method

Syntax
# include <cppad/utility/ode_gear.hpp>
 OdeGear(F, m, n, T, X, e)

Purpose
This routine applies Gear's Method to solve an explicit set of ordinary differential equations. We are given @(@ f : \B{R} \times \B{R}^n \rightarrow \B{R}^n @)@ be a smooth function. This routine solves the following initial value problem @[@ \begin{array}{rcl} x( t_{m-1} ) & = & x^0 \\ x^\prime (t) & = & f[t , x(t)] \end{array} @]@ for the value of @(@ x( t_m ) @)@. If your set of ordinary differential equations are not stiff an explicit method may be better (perhaps Runge45 .)

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
    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 (highest power of @(@ t @)@) used to represent the function @(@ x(t) @)@ in the multi-step method. Upon return from OdeGear, the i-th component of the polynomial is defined by @[@ p_i ( t_j ) = X[ j * n + i ] @]@ for @(@ j = 0 , \ldots , m @)@ (where @(@ 0 \leq i < n @)@). The value of @(@ 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 @(@ x(t) @)@.

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

X
The argument X has the prototype
    Vector &X
 and size greater than or equal to @(@ (m+1) * n @)@. On input to OdeGear, for @(@ j = 0 , \ldots , m-1 @)@, and @(@ i = 0 , \ldots , n-1 @)@ @[@ X[ j * n + i ] = x_i ( t_j ) @]@ Upon return from OdeGear, for @(@ i = 0 , \ldots , n-1 @)@ @[@ X[ m * n + i ] \approx x_i ( t_m ) @]@

e
The vector e is an approximate error bound for the result; i.e., @[@ 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., @[@ e = O ( h^m ) @]@ where @(@ h @)@ is the maximum of @(@ t_{j+1} - t_j @)@.

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 less than operator (returns a bool object)

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.cpp 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 @(@ x_j @)@ for the value @(@ x ( t_j ) \in \B{R}^n @)@ which is not to be confused with @(@ x_i (t) \in \B{R} @)@ in the notation above. The interpolating polynomial @(@ p(t) @)@ is given by @[@ 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 @(@ p^\prime (t) @)@ is given by @[@ 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 @(@ t_\ell @)@ we have @[@ \begin{array}{rcl} 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{array} @]@ We define the vector @(@ \alpha \in \B{R}^{m+1} @)@ by @[@ \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 @[@ p^\prime ( t_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m @]@ Gear's method determines @(@ x_m @)@ by solving the following nonlinear equation @[@ 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 @(@ x_m^k @)@, by solving the following affine approximation of the equation above @[@ \begin{array}{rcl} 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{array} @]@ In order to initialize Newton's method; i.e. choose @(@ x_m^0 @)@ we define the vector @(@ \beta \in \B{R}^{m+1} @)@ by @[@ \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 @[@ 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 @(@ x_m^0 @)@: @[@ 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.
Input File: include/cppad/utility/ode_gear.hpp