CppAD: A C++ Algorithmic Differentiation Package  20130102
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-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse 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 $codei%# include <cppad/ode_gear.hpp>
00038 %$$
00039 $codei%OdeGear(%F%, %m%, %n%, %T%, %X%, %e%)%$$
00040 
00041 
00042 $head Purpose$$
00043 This routine applies
00044 $cref/Gear's Method/OdeGear/Gear's Method/$$
00045 to solve an explicit set of ordinary differential equations.
00046 We are given 
00047 $latex f : \B{R} \times \B{R}^n \rightarrow \B{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 $cref 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 $icode Fun$$ 
00066 and the object $icode F$$ satisfy the prototype
00067 $codei%
00068      %Fun% &%F%
00069 %$$
00070 This must support the following set of calls
00071 $codei%
00072      %F%.Ode(%t%, %x%, %f%)
00073      %F%.Ode_dep(%t%, %x%, %f_x%)
00074 %$$
00075 
00076 $subhead t$$
00077 The argument $icode t$$ has prototype
00078 $codei%
00079      const %Scalar% &%t%
00080 %$$
00081 (see description of $cref/Scalar/OdeGear/Scalar/$$ below). 
00082 
00083 $subhead x$$
00084 The argument $icode x$$ has prototype
00085 $codei%
00086      const %Vector% &%x%
00087 %$$
00088 and has size $icode n$$
00089 (see description of $cref/Vector/OdeGear/Vector/$$ below). 
00090 
00091 $subhead f$$
00092 The argument $icode f$$ to $icode%F%.Ode%$$ has prototype
00093 $codei%
00094      %Vector% &%f%
00095 %$$
00096 On input and output, $icode f$$ is a vector of size $icode n$$
00097 and the input values of the elements of $icode f$$ do not matter.
00098 On output,
00099 $icode f$$ is set equal to $latex f(t, x)$$
00100 (see $icode f(t, x)$$ in $cref/Purpose/OdeGear/Purpose/$$). 
00101 
00102 $subhead f_x$$
00103 The argument $icode f_x$$ has prototype
00104 $codei%
00105      %Vector% &%f_x%
00106 %$$
00107 On input and output, $icode f_x$$ is a vector of size $latex n * n$$
00108 and the input values of the elements of $icode 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 $icode f$$, and $icode 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 $icode f$$ and $icode f_x$$.
00119 
00120 $head m$$
00121 The argument $icode m$$ has prototype
00122 $codei%
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 $icode n$$ has prototype
00137 $codei%
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 $icode T$$ has prototype
00145 $codei%
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 $icode X$$ has the prototype
00161 $codei%
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 $icode 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 $icode Scalar$$ must satisfy the conditions
00191 for a $cref NumericType$$ type.
00192 The routine $cref CheckNumericType$$ will generate an error message
00193 if this is not the case.
00194 In addition, the following operations must be defined for 
00195 $icode Scalar$$ objects $icode a$$ and $icode b$$:
00196 
00197 $table
00198 $bold Operation$$ $cnext $bold Description$$  $rnext
00199 $icode%a% < %b%$$ $cnext
00200      less than operator (returns a $code bool$$ object)
00201 $tend
00202 
00203 $head Vector$$
00204 The type $icode Vector$$ must be a $cref SimpleVector$$ class with
00205 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$.
00206 The routine $cref 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 $cref ode_gear.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 \B{R}^n$$ which is not to be confused
00225 with $latex x_i (t) \in \B{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 \B{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 \B{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           size_t(T.size()) >= (m+1),
00402           "OdeGear: size of T is not greater than or equal (m+1)"
00403      );
00404      CPPAD_ASSERT_KNOWN(
00405           size_t(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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines