CppAD: A C++ Algorithmic Differentiation Package  20171217
ode_gear.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_UTILITY_ODE_GEAR_HPP
2 # define CPPAD_UTILITY_ODE_GEAR_HPP
3
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15 /*
16 $begin OdeGear$$17 spell 18 cppad.hpp 19 Jan 20 bool 21 const 22 CppAD 23 dep 24$$ 25 26 27$section An Arbitrary Order Gear Method$$28 mindex OdeGear Ode stiff differential equation$$
29
30 $head Syntax$$31 codei%# include <cppad/utility/ode_gear.hpp> 32 %$$ 33$codei%OdeGear(%F%, %m%, %n%, %T%, %X%, %e%)%$$34 35 36 head Purpose$$
37 This routine applies
38 $cref/Gear's Method/OdeGear/Gear's Method/$$39 to solve an explicit set of ordinary differential equations. 40 We are given 41 latex f : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ be a smooth function. 42 This routine solves the following initial value problem 43$latex $44 \begin{array}{rcl} 45 x( t_{m-1} ) & = & x^0 \\ 46 x^\prime (t) & = & f[t , x(t)] 47 \end{array} 48$ $$49 for the value of latex x( t_m )$$.
50 If your set of ordinary differential equations are not stiff
51 an explicit method may be better (perhaps $cref Runge45$$.) 52 53 head Include$$ 54 The file$code cppad/ode_gear.hpp$$is included by code cppad/cppad.hpp$$
55 but it can also be included separately with out the rest of
56 the $code CppAD$$routines. 57 58 head Fun$$ 59 The class$icode Fun$$60 and the object icode F$$ satisfy the prototype
61 $codei% 62 %Fun% &%F% 63 %$$64 This must support the following set of calls 65 codei% 66 %F%.Ode(%t%, %x%, %f%) 67 %F%.Ode_dep(%t%, %x%, %f_x%) 68 %$$ 69 70$subhead t$$71 The argument icode t$$ has prototype
72 $codei% 73 const %Scalar% &%t% 74 %$$75 (see description of cref/Scalar/OdeGear/Scalar/$$ below). 76 77$subhead x$$78 The argument icode x$$ has prototype
79 $codei% 80 const %Vector% &%x% 81 %$$82 and has size icode n$$ 83 (see description of$cref/Vector/OdeGear/Vector/$$below). 84 85 subhead f$$
86 The argument $icode f$$to icode%F%.Ode%$$ has prototype 87$codei%
88  %Vector% &%f%
89 %$$90 On input and output, icode f$$ is a vector of size $icode n$$91 and the input values of the elements of icode f$$ do not matter. 92 On output, 93$icode f$$is set equal to latex f(t, x)$$
94 (see $icode f(t, x)$$in cref/Purpose/OdeGear/Purpose/$$). 95 96$subhead f_x$$97 The argument icode f_x$$ has prototype
98 $codei% 99 %Vector% &%f_x% 100 %$$101 On input and output, icode f_x$$ is a vector of size$latex n * n$$102 and the input values of the elements of icode f_x$$ do not matter.
103 On output,
104 $latex $105 f\_x [i * n + j] = \partial_{x(j)} f_i ( t , x ) 106$ $$107 108 subhead Warning$$ 109 The arguments$icode f$$, and icode f_x$$
110 must have a call by reference in their prototypes; i.e.,
111 do not forget the $code &$$in the prototype for 112 icode f$$ and$icode f_x$$. 113 114 head m$$
115 The argument $icode m$$has prototype 116 codei% 117 size_t %m% 118 %$$ 119 It specifies the order (highest power of$latex t$$) 120 used to represent the function latex x(t)$$ in the multi-step method.
121 Upon return from $code OdeGear$$, 122 the th i$$ component of the polynomial is defined by 123$latex $124 p_i ( t_j ) = X[ j * n + i ] 125$ $$126 for latex j = 0 , \ldots , m$$ (where $latex 0 \leq i < n$$). 127 The value of latex m$$ must be greater than or equal one. 128 129$head n$$130 The argument icode n$$ has prototype
131 $codei% 132 size_t %n% 133 %$$134 It specifies the range space dimension of the 135 vector valued function latex x(t)$$. 136 137$head T$$138 The argument icode T$$ has prototype
139 $codei% 140 const %Vector% &%T% 141 %$$142 and size greater than or equal to latex m+1$$. 143 For$latex j = 0 , \ldots m$$, latex T[j]$$ is the time
144 corresponding to time corresponding
145 to a previous point in the multi-step method.
146 The value $latex T[m]$$is the time 147 of the next point in the multi-step method. 148 The array latex T$$ must be monotone increasing; i.e., 149$latex T[j] < T[j+1]$$. 150 Above and below we often use the shorthand latex t_j$$ for $latex T[j]$$. 151 152 153 head X$$ 154 The argument$icode X$$has the prototype 155 codei% 156 %Vector% &%X% 157 %$$
158 and size greater than or equal to $latex (m+1) * n$$. 159 On input to code OdeGear$$, 160 for$latex j = 0 , \ldots , m-1$$, and 161 latex i = 0 , \ldots , n-1$$
162 $latex $163 X[ j * n + i ] = x_i ( t_j ) 164$ $$165 Upon return from code OdeGear$$, 166 for$latex i = 0 , \ldots , n-1$$167 latex $168 X[ m * n + i ] \approx x_i ( t_m ) 169$$$
170
171 $head e$$172 The vector icode e$$ is an approximate error bound for the result; i.e., 173$latex $174 e[i] \geq | X[ m * n + i ] - x_i ( t_m ) | 175$ $$176 The order of this approximation is one less than the order of 177 the solution; i.e., 178 latex $179 e = O ( h^m ) 180$$$
181 where $latex h$$is the maximum of latex t_{j+1} - t_j$$. 182 183$head Scalar$$184 The type icode Scalar$$ must satisfy the conditions
185 for a $cref NumericType$$type. 186 The routine cref CheckNumericType$$ will generate an error message 187 if this is not the case. 188 In addition, the following operations must be defined for 189$icode Scalar$$objects icode a$$ and $icode b$$: 190 191 table 192 bold Operation$$$cnext $bold Description$$rnext 193 icode%a% < %b%$$$cnext
194  less than operator (returns a $code bool$$object) 195 tend 196 197 head Vector$$ 198 The type$icode Vector$$must be a cref SimpleVector$$ class with
199 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$. 200 The routine cref CheckSimpleVector$$ will generate an error message 201 if this is not the case. 202 203$head Example$$204 children% 205 example/utility/ode_gear.cpp 206 %$$
207 The file
208 $cref ode_gear.cpp$$209 contains an example and test a test of using this routine. 210 It returns true if it succeeds and false otherwise. 211 212 head Source Code$$ 213 The source code for this routine is in the file 214$code cppad/ode_gear.hpp$$. 215 216 head Theory$$
217 For this discussion we use the shorthand $latex x_j$$218 for the value latex x ( t_j ) \in \B{R}^n$$ which is not to be confused 219 with$latex x_i (t) \in \B{R}$$in the notation above. 220 The interpolating polynomial latex p(t)$$ is given by
221 $latex $222 p(t) = 223 \sum_{j=0}^m 224 x_j 225 \frac{ 226 \prod_{i \neq j} ( t - t_i ) 227 }{ 228 \prod_{i \neq j} ( t_j - t_i ) 229 } 230$ $$231 The derivative latex p^\prime (t)$$ is given by 232$latex $233 p^\prime (t) = 234 \sum_{j=0}^m 235 x_j 236 \frac{ 237 \sum_{i \neq j} \prod_{k \neq i,j} ( t - t_k ) 238 }{ 239 \prod_{k \neq j} ( t_j - t_k ) 240 } 241$ $$242 Evaluating the derivative at the point latex t_\ell$$ we have
243 $latex $244 \begin{array}{rcl} 245 p^\prime ( t_\ell ) & = & 246 x_\ell 247 \frac{ 248 \sum_{i \neq \ell} \prod_{k \neq i,\ell} ( t_\ell - t_k ) 249 }{ 250 \prod_{k \neq \ell} ( t_\ell - t_k ) 251 } 252 + 253 \sum_{j \neq \ell} 254 x_j 255 \frac{ 256 \sum_{i \neq j} \prod_{k \neq i,j} ( t_\ell - t_k ) 257 }{ 258 \prod_{k \neq j} ( t_j - t_k ) 259 } 260 \\ 261 & = & 262 x_\ell 263 \sum_{i \neq \ell} 264 \frac{ 1 }{ t_\ell - t_i } 265 + 266 \sum_{j \neq \ell} 267 x_j 268 \frac{ 269 \prod_{k \neq \ell,j} ( t_\ell - t_k ) 270 }{ 271 \prod_{k \neq j} ( t_j - t_k ) 272 } 273 \\ 274 & = & 275 x_\ell 276 \sum_{k \neq \ell} ( t_\ell - t_k )^{-1} 277 + 278 \sum_{j \neq \ell} 279 x_j 280 ( t_j - t_\ell )^{-1} 281 \prod_{k \neq \ell ,j} ( t_\ell - t_k ) / ( t_j - t_k ) 282 \end{array} 283$ $$284 We define the vector latex \alpha \in \B{R}^{m+1}$$ by 285$latex $286 \alpha_j = \left\{ \begin{array}{ll} 287 \sum_{k \neq m} ( t_m - t_k )^{-1} 288 & {\rm if} \; j = m 289 \\ 290 ( t_j - t_m )^{-1} 291 \prod_{k \neq m,j} ( t_m - t_k ) / ( t_j - t_k ) 292 & {\rm otherwise} 293 \end{array} \right. 294$ $$295 It follows that 296 latex $297 p^\prime ( t_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m 298$$$
299 Gear's method determines $latex x_m$$by solving the following 300 nonlinear equation 301 latex $302 f( t_m , x_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m 303$$$ 304 Newton's method for solving this equation determines iterates, 305 which we denote by$latex x_m^k$$, by solving the following affine 306 approximation of the equation above 307 latex $308 \begin{array}{rcl} 309 f( t_m , x_m^{k-1} ) + \partial_x f( t_m , x_m^{k-1} ) ( x_m^k - x_m^{k-1} ) 310 & = & 311 \alpha_0 x_0^k + \alpha_1 x_1 + \cdots + \alpha_m x_m 312 \\ 313 \left[ \alpha_m I - \partial_x f( t_m , x_m^{k-1} ) \right] x_m 314 & = & 315 \left[ 316 f( t_m , x_m^{k-1} ) - \partial_x f( t_m , x_m^{k-1} ) x_m^{k-1} 317 - \alpha_0 x_0 - \cdots - \alpha_{m-1} x_{m-1} 318 \right] 319 \end{array} 320$$$
321 In order to initialize Newton's method; i.e. choose $latex x_m^0$$322 we define the vector latex \beta \in \B{R}^{m+1}$$ by 323$latex $324 \beta_j = \left\{ \begin{array}{ll} 325 \sum_{k \neq m-1} ( t_{m-1} - t_k )^{-1} 326 & {\rm if} \; j = m-1 327 \\ 328 ( t_j - t_{m-1} )^{-1} 329 \prod_{k \neq m-1,j} ( t_{m-1} - t_k ) / ( t_j - t_k ) 330 & {\rm otherwise} 331 \end{array} \right. 332$ $$333 It follows that 334 latex $335 p^\prime ( t_{m-1} ) = \beta_0 x_0 + \cdots + \beta_m x_m 336$$$
337 We solve the following approximation of the equation above to determine
338 $latex x_m^0$$: 339 latex $340 f( t_{m-1} , x_{m-1} ) = 341 \beta_0 x_0 + \cdots + \beta_{m-1} x_{m-1} + \beta_m x_m^0 342$$$ 343 344 345$head Gear's Method
346 C. W. Gear,
347 Simultaneous Numerical Solution of Differential-Algebraic Equations,''
348 IEEE Transactions on Circuit Theory,
349 vol. 18, no. 1, pp. 89-95, Jan. 1971.
350
351
352 \$end
353 --------------------------------------------------------------------------
354 */
355
356 # include <cstddef>
360 # include <cppad/utility/vector.hpp>
361 # include <cppad/utility/lu_factor.hpp>
362 # include <cppad/utility/lu_invert.hpp>
363
364 namespace CppAD { // BEGIN CppAD namespace
365
366 template <typename Vector, typename Fun>
367 void OdeGear(
368  Fun &F ,
369  size_t m ,
370  size_t n ,
371  const Vector &T ,
372  Vector &X ,
373  Vector &e )
374 {
375  // temporary indices
376  size_t i, j, k;
377
378  typedef typename Vector::value_type Scalar;
379
380  // check numeric type specifications
381  CheckNumericType<Scalar>();
382
383  // check simple vector class specifications
384  CheckSimpleVector<Scalar, Vector>();
385
387  m >= 1,
388  "OdeGear: m is less than one"
389  );
391  n > 0,
392  "OdeGear: n is equal to zero"
393  );
395  size_t(T.size()) >= (m+1),
396  "OdeGear: size of T is not greater than or equal (m+1)"
397  );
399  size_t(X.size()) >= (m+1) * n,
400  "OdeGear: size of X is not greater than or equal (m+1) * n"
401  );
402  for(j = 0; j < m; j++) CPPAD_ASSERT_KNOWN(
403  T[j] < T[j+1],
404  "OdeGear: the array T is not monotone increasing"
405  );
406
407  // some constants
408  Scalar zero(0);
409  Scalar one(1);
410
411  // vectors required by method
412  Vector alpha(m + 1);
413  Vector beta(m + 1);
414  Vector f(n);
415  Vector f_x(n * n);
416  Vector x_m0(n);
417  Vector x_m(n);
418  Vector b(n);
419  Vector A(n * n);
420
421  // compute alpha[m]
422  alpha[m] = zero;
423  for(k = 0; k < m; k++)
424  alpha[m] += one / (T[m] - T[k]);
425
426  // compute beta[m-1]
427  beta[m-1] = one / (T[m-1] - T[m]);
428  for(k = 0; k < m-1; k++)
429  beta[m-1] += one / (T[m-1] - T[k]);
430
431
432  // compute other components of alpha
433  for(j = 0; j < m; j++)
434  { // compute alpha[j]
435  alpha[j] = one / (T[j] - T[m]);
436  for(k = 0; k < m; k++)
437  { if( k != j )
438  { alpha[j] *= (T[m] - T[k]);
439  alpha[j] /= (T[j] - T[k]);
440  }
441  }
442  }
443
444  // compute other components of beta
445  for(j = 0; j <= m; j++)
446  { if( j != m-1 )
447  { // compute beta[j]
448  beta[j] = one / (T[j] - T[m-1]);
449  for(k = 0; k <= m; k++)
450  { if( k != j && k != m-1 )
451  { beta[j] *= (T[m-1] - T[k]);
452  beta[j] /= (T[j] - T[k]);
453  }
454  }
455  }
456  }
457
458  // evaluate f(T[m-1], x_{m-1} )
459  for(i = 0; i < n; i++)
460  x_m[i] = X[(m-1) * n + i];
461  F.Ode(T[m-1], x_m, f);
462
463  // solve for x_m^0
464  for(i = 0; i < n; i++)
465  { x_m[i] = f[i];
466  for(j = 0; j < m; j++)
467  x_m[i] -= beta[j] * X[j * n + i];
468  x_m[i] /= beta[m];
469  }
470  x_m0 = x_m;
471
472  // evaluate partial w.r.t x of f(T[m], x_m^0)
473  F.Ode_dep(T[m], x_m, f_x);
474
475  // compute the matrix A = ( alpha[m] * I - f_x )
476  for(i = 0; i < n; i++)
477  { for(j = 0; j < n; j++)
478  A[i * n + j] = - f_x[i * n + j];
479  A[i * n + i] += alpha[m];
480  }
481
482  // LU factor (and overwrite) the matrix A
483  CppAD::vector<size_t> ip(n) , jp(n);
484 # ifndef NDEBUG
485  int sign =
486 # endif
487  LuFactor(ip, jp, A);
489  sign != 0,
490  "OdeGear: step size is to large"
491  );
492
493  // Iterations of Newton's method
494  for(k = 0; k < 3; k++)
495  {
496  // only evaluate f( T[m] , x_m ) keep f_x during iteration
497  F.Ode(T[m], x_m, f);
498
499  // b = f + f_x x_m - alpha[0] x_0 - ... - alpha[m-1] x_{m-1}
500  for(i = 0; i < n; i++)
501  { b[i] = f[i];
502  for(j = 0; j < n; j++)
503  b[i] -= f_x[i * n + j] * x_m[j];
504  for(j = 0; j < m; j++)
505  b[i] -= alpha[j] * X[ j * n + i ];
506  }
507  LuInvert(ip, jp, A, b);
508  x_m = b;
509  }
510
511  // return estimate for x( t[k] ) and the estimated error bound
512  for(i = 0; i < n; i++)
513  { X[m * n + i] = x_m[i];
514  e[i] = x_m[i] - x_m0[i];
515  if( e[i] < zero )
516  e[i] = - e[i];
517  }
518 }
519
520 } // End CppAD namespace
521
522 # endif