CppAD: A C++ Algorithmic Differentiation Package 20110419
ode_gear.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_ODE_GEAR_INCLUDED
00003 # define CPPAD_ODE_GEAR_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin OdeGear$$
00018 $spell
00019         cppad.hpp
00020         Jan
00021         bool
00022         const
00023         CppAD
00024         dep
00025 $$
00026 
00027 $index OdeGear$$
00028 $index Ode, Gear$$
00029 $index Gear, Ode$$
00030 $index stiff, Ode$$
00031 $index differential, equation$$
00032 $index equation, differential$$
00033  
00034 $section An Arbitrary Order Gear Method$$
00035 
00036 $head Syntax$$
00037 $syntax%# include <cppad/ode_gear.hpp>
00038 %$$
00039 $syntax%OdeGear(%F%, %m%, %n%, %T%, %X%, %e%)%$$
00040 
00041 
00042 $head Purpose$$
00043 This routine applies
00044 $xref/OdeGear/Gear's Method/Gear's Method/$$
00045 to solve an explicit set of ordinary differential equations.
00046 We are given 
00047 $latex f : \R \times \R^n \rightarrow \R^n$$ be a smooth function.
00048 This routine solves the following initial value problem
00049 $latex \[
00050 \begin{array}{rcl}
00051         x( t_{m-1} )  & = & x^0    \\
00052         x^\prime (t)  & = & f[t , x(t)] 
00053 \end{array}
00054 \] $$
00055 for the value of $latex x( t_m )$$.
00056 If your set of  ordinary differential equations are not stiff
00057 an explicit method may be better (perhaps $xref/Runge45/$$.)
00058 
00059 $head Include$$
00060 The file $code cppad/ode_gear.hpp$$ is included by $code cppad/cppad.hpp$$
00061 but it can also be included separately with out the rest of 
00062 the $code CppAD$$ routines.
00063 
00064 $head Fun$$
00065 The class $italic Fun$$ 
00066 and the object $italic F$$ satisfy the prototype
00067 $syntax%
00068         %Fun% &%F%
00069 %$$
00070 This must support the following set of calls
00071 $syntax%
00072         %F%.Ode(%t%, %x%, %f%)
00073         %F%.Ode_dep(%t%, %x%, %f_x%)
00074 %$$
00075 
00076 $subhead t$$
00077 The argument $italic t$$ has prototype
00078 $syntax%
00079         const %Scalar% &%t%
00080 %$$
00081 (see description of $xref/OdeGear/Scalar/Scalar/$$ below). 
00082 
00083 $subhead x$$
00084 The argument $italic x$$ has prototype
00085 $syntax%
00086         const %Vector% &%x%
00087 %$$
00088 and has size $italic n$$
00089 (see description of $xref/OdeGear/Vector/Vector/$$ below). 
00090 
00091 $subhead f$$
00092 The argument $italic f$$ to $syntax%%F%.Ode%$$ has prototype
00093 $syntax%
00094         %Vector% &%f%
00095 %$$
00096 On input and output, $italic f$$ is a vector of size $italic n$$
00097 and the input values of the elements of $italic f$$ do not matter.
00098 On output,
00099 $italic f$$ is set equal to $latex f(t, x)$$
00100 (see $italic f(t, x)$$ in $xref/OdeGear/Purpose/Purpose/$$). 
00101 
00102 $subhead f_x$$
00103 The argument $italic f_x$$ has prototype
00104 $syntax%
00105         %Vector% &%f_x%
00106 %$$
00107 On input and output, $italic f_x$$ is a vector of size $latex n * n$$
00108 and the input values of the elements of $italic f_x$$ do not matter.
00109 On output, 
00110 $latex \[
00111         f\_x [i * n + j] = \partial_{x(j)} f_i ( t , x )
00112 \] $$ 
00113 
00114 $subhead Warning$$
00115 The arguments $italic f$$, and $italic f_x$$
00116 must have a call by reference in their prototypes; i.e.,
00117 do not forget the $code &$$ in the prototype for 
00118 $italic f$$ and $italic f_x$$.
00119 
00120 $head m$$
00121 The argument $italic m$$ has prototype
00122 $syntax%
00123         size_t %m%
00124 %$$
00125 It specifies the order (highest power of $latex t$$) 
00126 used to represent the function $latex x(t)$$ in the multi-step method. 
00127 Upon return from $code OdeGear$$,
00128 the $th i$$ component of the polynomial is defined by
00129 $latex \[
00130         p_i ( t_j ) = X[ j * n + i ]
00131 \] $$
00132 for $latex j = 0 , \ldots , m$$ (where $latex 0 \leq i < n$$).
00133 The value of $latex m$$ must be greater than or equal one.
00134 
00135 $head n$$
00136 The argument $italic n$$ has prototype
00137 $syntax%
00138         size_t %n%
00139 %$$
00140 It specifies the range space dimension of the 
00141 vector valued function $latex x(t)$$.
00142 
00143 $head T$$
00144 The argument $italic T$$ has prototype
00145 $syntax%
00146         const %Vector% &%T%
00147 %$$
00148 and size greater than or equal to $latex m+1$$.
00149 For $latex j = 0 , \ldots m$$, $latex T[j]$$ is the time 
00150 corresponding to time corresponding 
00151 to a previous point in the multi-step method.
00152 The value $latex T[m]$$ is the time 
00153 of the next point in the multi-step method.
00154 The array $latex T$$ must be monotone increasing; i.e.,
00155 $latex T[j] < T[j+1]$$.
00156 Above and below we often use the shorthand $latex t_j$$ for $latex T[j]$$.
00157 
00158 
00159 $head X$$
00160 The argument $italic X$$ has the prototype
00161 $syntax%
00162         %Vector% &%X%
00163 %$$
00164 and size greater than or equal to $latex (m+1) * n$$.
00165 On input to $code OdeGear$$,
00166 for $latex j = 0 , \ldots , m-1$$, and
00167 $latex i = 0 , \ldots , n-1$$ 
00168 $latex \[
00169         X[ j * n + i ] = x_i ( t_j )
00170 \] $$
00171 Upon return from $code OdeGear$$,
00172 for $latex i = 0 , \ldots , n-1$$ 
00173 $latex \[
00174         X[ m * n + i ] \approx x_i ( t_m ) 
00175 \] $$
00176 
00177 $head e$$
00178 The vector $italic e$$ is an approximate error bound for the result; i.e.,
00179 $latex \[
00180         e[i] \geq | X[ m * n + i ] - x_i ( t_m ) |
00181 \] $$
00182 The order of this approximation is one less than the order of
00183 the solution; i.e., 
00184 $latex \[
00185         e = O ( h^m )
00186 \] $$
00187 where $latex h$$ is the maximum of $latex t_{j+1} - t_j$$.
00188 
00189 $head Scalar$$
00190 The type $italic Scalar$$ must satisfy the conditions
00191 for a $xref/NumericType/$$ type.
00192 The routine $xref/CheckNumericType/$$ will generate an error message
00193 if this is not the case.
00194 In addition, the following operations must be defined for 
00195 $italic Scalar$$ objects $italic a$$ and $italic b$$:
00196 
00197 $table
00198 $bold Operation$$ $cnext $bold Description$$  $rnext
00199 $syntax%%a% < %b%$$ $cnext
00200         less than operator (returns a $code bool$$ object)
00201 $tend
00202 
00203 $head Vector$$
00204 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with
00205 $xref/SimpleVector/Elements of Specified Type/elements of type Scalar/$$.
00206 The routine $xref/CheckSimpleVector/$$ will generate an error message
00207 if this is not the case.
00208 
00209 $head Example$$
00210 $children%
00211         example/ode_gear.cpp
00212 %$$
00213 The file
00214 $xref/OdeGear.cpp/$$
00215 contains an example and test a test of using this routine.
00216 It returns true if it succeeds and false otherwise.
00217 
00218 $head Source Code$$
00219 The source code for this routine is in the file
00220 $code cppad/ode_gear.hpp$$.
00221 
00222 $head Theory$$
00223 For this discussion we use the shorthand $latex x_j$$ 
00224 for the value $latex x ( t_j ) \in \R^n$$ which is not to be confused
00225 with $latex x_i (t) \in \R$$ in the notation above.
00226 The interpolating polynomial $latex p(t)$$ is given by
00227 $latex \[
00228 p(t) = 
00229 \sum_{j=0}^m 
00230 x_j
00231 \frac{ 
00232         \prod_{i \neq j} ( t - t_i )
00233 }{
00234         \prod_{i \neq j} ( t_j - t_i ) 
00235 }
00236 \] $$
00237 The derivative $latex p^\prime (t)$$ is given by
00238 $latex \[
00239 p^\prime (t) = 
00240 \sum_{j=0}^m 
00241 x_j
00242 \frac{ 
00243         \sum_{i \neq j} \prod_{k \neq i,j} ( t - t_k )
00244 }{ 
00245         \prod_{k \neq j} ( t_j - t_k ) 
00246 }
00247 \] $$
00248 Evaluating the derivative at the point $latex t_\ell$$ we have
00249 $latex \[
00250 \begin{array}{rcl}
00251 p^\prime ( t_\ell ) & = & 
00252 x_\ell
00253 \frac{ 
00254         \sum_{i \neq \ell} \prod_{k \neq i,\ell} ( t_\ell - t_k )
00255 }{ 
00256         \prod_{k \neq \ell} ( t_\ell - t_k ) 
00257 }
00258 +
00259 \sum_{j \neq \ell}
00260 x_j
00261 \frac{ 
00262         \sum_{i \neq j} \prod_{k \neq i,j} ( t_\ell - t_k )
00263 }{ 
00264         \prod_{k \neq j} ( t_j - t_k ) 
00265 }
00266 \\
00267 & = &
00268 x_\ell
00269 \sum_{i \neq \ell} 
00270 \frac{ 1 }{ t_\ell - t_i }
00271 +
00272 \sum_{j \neq \ell}
00273 x_j
00274 \frac{ 
00275         \prod_{k \neq \ell,j} ( t_\ell - t_k )
00276 }{ 
00277         \prod_{k \neq j} ( t_j - t_k ) 
00278 }
00279 \\
00280 & = &
00281 x_\ell
00282 \sum_{k \neq \ell} ( t_\ell - t_k )^{-1}
00283 +
00284 \sum_{j \neq \ell}
00285 x_j
00286 ( t_j - t_\ell )^{-1}
00287 \prod_{k \neq \ell ,j} ( t_\ell - t_k ) / ( t_j - t_k )
00288 \end{array}
00289 \] $$
00290 We define the vector $latex \alpha \in \R^{m+1}$$ by
00291 $latex \[
00292 \alpha_j = \left\{ \begin{array}{ll}
00293 \sum_{k \neq m} ( t_m - t_k )^{-1}
00294         & {\rm if} \; j = m
00295 \\
00296 ( t_j - t_m )^{-1}
00297 \prod_{k \neq m,j} ( t_m - t_k ) / ( t_j - t_k )
00298         & {\rm otherwise}
00299 \end{array} \right.
00300 \] $$
00301 It follows that
00302 $latex \[
00303         p^\prime ( t_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m
00304 \] $$
00305 Gear's method determines $latex x_m$$ by solving the following 
00306 nonlinear equation
00307 $latex \[
00308         f( t_m , x_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m
00309 \] $$
00310 Newton's method for solving this equation determines iterates, 
00311 which we denote by $latex x_m^k$$, by solving the following affine 
00312 approximation of the equation above
00313 $latex \[
00314 \begin{array}{rcl}
00315 f( t_m , x_m^{k-1} ) + \partial_x f( t_m , x_m^{k-1} ) ( x_m^k - x_m^{k-1} )
00316 & = &
00317 \alpha_0 x_0^k + \alpha_1 x_1 + \cdots + \alpha_m x_m
00318 \\
00319 \left[ \alpha_m I - \partial_x f( t_m , x_m^{k-1} ) \right]  x_m
00320 & = &
00321 \left[ 
00322 f( t_m , x_m^{k-1} ) - \partial_x f( t_m , x_m^{k-1} ) x_m^{k-1} 
00323 - \alpha_0 x_0 - \cdots - \alpha_{m-1} x_{m-1}
00324 \right]
00325 \end{array}
00326 \] $$
00327 In order to initialize Newton's method; i.e. choose $latex x_m^0$$
00328 we define the vector $latex \beta \in \R^{m+1}$$ by
00329 $latex \[
00330 \beta_j = \left\{ \begin{array}{ll}
00331 \sum_{k \neq m-1} ( t_{m-1} - t_k )^{-1}
00332         & {\rm if} \; j = m-1
00333 \\
00334 ( t_j - t_{m-1} )^{-1}
00335 \prod_{k \neq m-1,j} ( t_{m-1} - t_k ) / ( t_j - t_k )
00336         & {\rm otherwise}
00337 \end{array} \right.
00338 \] $$
00339 It follows that
00340 $latex \[
00341         p^\prime ( t_{m-1} ) = \beta_0 x_0 + \cdots + \beta_m x_m
00342 \] $$
00343 We solve the following approximation of the equation above to determine
00344 $latex x_m^0$$:
00345 $latex \[
00346         f( t_{m-1} , x_{m-1} ) = 
00347         \beta_0 x_0 + \cdots + \beta_{m-1} x_{m-1} + \beta_m x_m^0
00348 \] $$
00349 
00350 
00351 $head Gear's Method$$
00352 C. W. Gear, 
00353 ``Simultaneous Numerical Solution of Differential-Algebraic Equations,'' 
00354 IEEE Transactions on Circuit Theory, 
00355 vol. 18, no. 1, pp. 89-95, Jan. 1971.
00356 
00357 
00358 $end
00359 --------------------------------------------------------------------------
00360 */
00361 
00362 # include <cstddef>
00363 # include <cppad/local/cppad_assert.hpp>
00364 # include <cppad/check_simple_vector.hpp>
00365 # include <cppad/check_numeric_type.hpp>
00366 # include <cppad/vector.hpp>
00367 # include <cppad/lu_factor.hpp>
00368 # include <cppad/lu_invert.hpp>
00369 
00370 namespace CppAD { // BEGIN CppAD namespace
00371 
00372 template <typename Vector, typename Fun>
00373 void OdeGear(
00374         Fun          &F  , 
00375         size_t        m  ,
00376         size_t        n  ,
00377         const Vector &T  , 
00378         Vector       &X  ,
00379         Vector       &e  ) 
00380 {
00381         // temporary indices
00382         size_t i, j, k;
00383 
00384         typedef typename Vector::value_type Scalar;
00385 
00386         // check numeric type specifications
00387         CheckNumericType<Scalar>();
00388 
00389         // check simple vector class specifications
00390         CheckSimpleVector<Scalar, Vector>();
00391 
00392         CPPAD_ASSERT_KNOWN(
00393                 m >= 1,
00394                 "OdeGear: m is less than one"
00395         );
00396         CPPAD_ASSERT_KNOWN(
00397                 n > 0,
00398                 "OdeGear: n is equal to zero"
00399         );
00400         CPPAD_ASSERT_KNOWN(
00401                 T.size() >= (m+1),
00402                 "OdeGear: size of T is not greater than or equal (m+1)"
00403         );
00404         CPPAD_ASSERT_KNOWN(
00405                 X.size() >= (m+1) * n,
00406                 "OdeGear: size of X is not greater than or equal (m+1) * n"
00407         );
00408         for(j = 0; j < m; j++) CPPAD_ASSERT_KNOWN(
00409                 T[j] < T[j+1],
00410                 "OdeGear: the array T is not monotone increasing"
00411         );
00412 
00413         // some constants
00414         Scalar zero(0);
00415         Scalar one(1);
00416 
00417         // vectors required by method
00418         Vector alpha(m + 1);
00419         Vector beta(m + 1);
00420         Vector f(n);
00421         Vector f_x(n * n);
00422         Vector x_m0(n);
00423         Vector x_m(n);
00424         Vector b(n);
00425         Vector A(n * n);
00426 
00427         // compute alpha[m] 
00428         alpha[m] = zero;
00429         for(k = 0; k < m; k++)
00430                 alpha[m] += one / (T[m] - T[k]);
00431 
00432         // compute beta[m-1]
00433         beta[m-1] = one / (T[m-1] - T[m]);
00434         for(k = 0; k < m-1; k++)
00435                 beta[m-1] += one / (T[m-1] - T[k]);
00436 
00437 
00438         // compute other components of alpha 
00439         for(j = 0; j < m; j++)
00440         {       // compute alpha[j]
00441                 alpha[j] = one / (T[j] - T[m]);
00442                 for(k = 0; k < m; k++)
00443                 {       if( k != j )
00444                         {       alpha[j] *= (T[m] - T[k]);
00445                                 alpha[j] /= (T[j] - T[k]);
00446                         }
00447                 }
00448         }
00449 
00450         // compute other components of beta 
00451         for(j = 0; j <= m; j++)
00452         {       if( j != m-1 )
00453                 {       // compute beta[j]
00454                         beta[j] = one / (T[j] - T[m-1]);
00455                         for(k = 0; k <= m; k++)
00456                         {       if( k != j && k != m-1 )
00457                                 {       beta[j] *= (T[m-1] - T[k]);
00458                                         beta[j] /= (T[j] - T[k]);
00459                                 }
00460                         }
00461                 }
00462         }
00463 
00464         // evaluate f(T[m-1], x_{m-1} )
00465         for(i = 0; i < n; i++)
00466                 x_m[i] = X[(m-1) * n + i];
00467         F.Ode(T[m-1], x_m, f);
00468 
00469         // solve for x_m^0
00470         for(i = 0; i < n; i++)
00471         {       x_m[i] =  f[i];
00472                 for(j = 0; j < m; j++)
00473                         x_m[i] -= beta[j] * X[j * n + i];
00474                 x_m[i] /= beta[m];
00475         }
00476         x_m0 = x_m;
00477 
00478         // evaluate partial w.r.t x of f(T[m], x_m^0)
00479         F.Ode_dep(T[m], x_m, f_x);
00480 
00481         // compute the matrix A = ( alpha[m] * I - f_x )
00482         for(i = 0; i < n; i++)
00483         {       for(j = 0; j < n; j++)
00484                         A[i * n + j]  = - f_x[i * n + j];
00485                 A[i * n + i] += alpha[m];
00486         }
00487 
00488         // LU factor (and overwrite) the matrix A
00489         int sign;
00490         CppAD::vector<size_t> ip(n) , jp(n);
00491         sign = LuFactor(ip, jp, A);
00492         CPPAD_ASSERT_KNOWN(
00493                 sign != 0,
00494                 "OdeGear: step size is to large"
00495         );
00496 
00497         // Iterations of Newton's method
00498         for(k = 0; k < 3; k++)
00499         {
00500                 // only evaluate f( T[m] , x_m ) keep f_x during iteration
00501                 F.Ode(T[m], x_m, f);
00502 
00503                 // b = f + f_x x_m - alpha[0] x_0 - ... - alpha[m-1] x_{m-1}
00504                 for(i = 0; i < n; i++)
00505                 {       b[i]         = f[i];
00506                         for(j = 0; j < n; j++)
00507                                 b[i]         -= f_x[i * n + j] * x_m[j];
00508                         for(j = 0; j < m; j++)
00509                                 b[i] -= alpha[j] * X[ j * n + i ];
00510                 }
00511                 LuInvert(ip, jp, A, b);
00512                 x_m = b;
00513         }
00514 
00515         // return estimate for x( t[k] ) and the estimated error bound
00516         for(i = 0; i < n; i++)
00517         {       X[m * n + i] = x_m[i];
00518                 e[i]         = x_m[i] - x_m0[i];
00519                 if( e[i] < zero )
00520                         e[i] = - e[i];
00521         }
00522 }
00523 
00524 } // End CppAD namespace 
00525 
00526 # endif