Prev Next

An Arbitrary Order Gear Method

Syntax
# include <cppad/ode_gear.hpp>
OdeGear(FmnTXe)

Purpose
This routine applies Gear's Method to solve an explicit set of ordinary differential equations. We are given  f : \R \times \R^n \rightarrow \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/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(txf)
     
F.Ode_dep(txf_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 OdeGear.cpp contains an example and test a test of using this routine. It returns true if it succeeds and false otherwise.

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 \R^n which is not to be confused with  x_i (t) \in \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 \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 \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: cppad/ode_gear.hpp