CppAD: A C++ Algorithmic Differentiation Package  20171217
ode_gear_control.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
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 OdeGearControl$$17 spell 18 cppad.hpp 19 CppAD 20 xf 21 xi 22 smin 23 smax 24 eabs 25 ef 26 maxabs 27 nstep 28 tf 29 sini 30 erel 31 dep 32 const 33 tb 34 ta 35 exp 36$$ 37 38 39 40$section An Error Controller for Gear's Ode Solvers$$41 mindex OdeGearControl Gear differential equation$$
42
43 $head Syntax$$44 codei%# include <cppad/utility/ode_gear_control.hpp> 45 %$$ 46$icode%xf% = OdeGearControl(%F%, %M%, %ti%, %tf%, %xi%,
47  %smin%, %smax%, %sini%, %eabs%, %erel%, %ef% , %maxabs%, %nstep% )%$$48 49 50 head Purpose$$
51 Let $latex \B{R}$$denote the real numbers 52 and let latex f : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ be a smooth function. 53 We define$latex X : [ti , tf] \rightarrow \B{R}^n$$by 54 the following initial value problem: 55 latex $56 \begin{array}{rcl} 57 X(ti) & = & xi \\ 58 X'(t) & = & f[t , X(t)] 59 \end{array} 60$$$
61 The routine $cref OdeGear$$is a stiff multi-step method that 62 can be used to approximate the solution to this equation. 63 The routine code OdeGearControl$$ sets up this multi-step method 64 and controls the error during such an approximation. 65 66$head Include$$67 The file code cppad/ode_gear_control.hpp$$
68 is included by $code cppad/cppad.hpp$$69 but it can also be included separately with out the rest of 70 the code CppAD$$ routines. 71 72$head Notation$$73 The template parameter types cref/Scalar/OdeGearControl/Scalar/$$ and
74 $cref/Vector/OdeGearControl/Vector/$$are documented below. 75 76 head xf$$ 77 The return value$icode xf$$has the prototype 78 codei% 79 %Vector% %xf% 80 %$$
81 and the size of $icode xf$$is equal to icode n$$ 82 (see description of$cref/Vector/OdeGear/Vector/$$below). 83 It is the approximation for latex X(tf)$$.
84
85 $head Fun$$86 The class icode Fun$$ 87 and the object$icode F$$satisfy the prototype 88 codei% 89 %Fun% &%F% 90 %$$
91 This must support the following set of calls
92 $codei% 93 %F%.Ode(%t%, %x%, %f%) 94 %F%.Ode_dep(%t%, %x%, %f_x%) 95 %$$96 97 subhead t$$ 98 The argument$icode t$$has prototype 99 codei% 100 const %Scalar% &%t% 101 %$$
102 (see description of $cref/Scalar/OdeGear/Scalar/$$below). 103 104 subhead x$$ 105 The argument$icode x$$has prototype 106 codei% 107 const %Vector% &%x% 108 %$$
109 and has size $icode N$$110 (see description of cref/Vector/OdeGear/Vector/$$ below). 111 112$subhead f$$113 The argument icode f$$ to $icode%F%.Ode%$$has prototype 114 codei% 115 %Vector% &%f% 116 %$$ 117 On input and output,$icode f$$is a vector of size icode N$$
118 and the input values of the elements of $icode f$$do not matter. 119 On output, 120 icode f$$ is set equal to$latex f(t, x)$$121 (see icode f(t, x)$$ in $cref/Purpose/OdeGear/Purpose/$$). 122 123 subhead f_x$$ 124 The argument$icode f_x$$has prototype 125 codei% 126 %Vector% &%f_x% 127 %$$
128 On input and output, $icode f_x$$is a vector of size latex N * N$$ 129 and the input values of the elements of$icode f_x$$do not matter. 130 On output, 131 latex $132 f\_x [i * n + j] = \partial_{x(j)} f_i ( t , x ) 133$$$
134
135 $subhead Warning$$136 The arguments icode f$$, and$icode f_x$$137 must have a call by reference in their prototypes; i.e., 138 do not forget the code &$$ in the prototype for
139 $icode f$$and icode f_x$$. 140 141$head M$$142 The argument icode M$$ has prototype
143 $codei% 144 size_t %M% 145 %$$146 It specifies the order of the multi-step method; i.e., 147 the order of the approximating polynomial 148 (after the initialization process). 149 The argument icode M$$ must greater than or equal one. 150 151$head ti$$152 The argument icode ti$$ has prototype
153 $codei% 154 const %Scalar% &%ti% 155 %$$156 It specifies the initial time for the integration of 157 the differential equation. 158 159 head tf$$ 160 The argument$icode tf$$has prototype 161 codei% 162 const %Scalar% &%tf% 163 %$$
164 It specifies the final time for the integration of
165 the differential equation.
166
167 $head xi$$168 The argument icode xi$$ has prototype 169$codei%
170  const %Vector% &%xi%
171 %$$172 and size icode n$$.
173 It specifies value of $latex X(ti)$$. 174 175 head smin$$ 176 The argument$icode smin$$has prototype 177 codei% 178 const %Scalar% &%smin% 179 %$$
180 The minimum value of $latex T[M] - T[M-1]$$in a call to code OdeGear$$ 181 will be$latex smin$$except for the last two calls where it may be 182 as small as latex smin / 2$$.
183 The value of $icode smin$$must be less than or equal icode smax$$. 184 185$head smax$$186 The argument icode smax$$ has prototype
187 $codei% 188 const %Scalar% &%smax% 189 %$$190 It specifies the maximum step size to use during the integration; 191 i.e., the maximum value for latex T[M] - T[M-1]$$ 192 in a call to$code OdeGear$$. 193 194 head sini$$
195 The argument $icode sini$$has prototype 196 codei% 197 %Scalar% &%sini% 198 %$$ 199 The value of$icode sini$$is the minimum 200 step size to use during initialization of the multi-step method; i.e., 201 for calls to code OdeGear$$ where $latex m < M$$. 202 The value of icode sini$$ must be less than or equal$icode smax$$203 (and can also be less than icode smin$$).
204
205 $head eabs$$206 The argument icode eabs$$ has prototype 207$codei%
208  const %Vector% &%eabs%
209 %$$210 and size icode n$$.
211 Each of the elements of $icode eabs$$must be 212 greater than or equal zero. 213 It specifies a bound for the absolute 214 error in the return value icode xf$$ as an approximation for$latex X(tf)$$. 215 (see the 216 cref/error criteria discussion/OdeGearControl/Error Criteria Discussion/$$
217 below).
218
219 $head erel$$220 The argument icode erel$$ has prototype 221$codei%
222  const %Scalar% &%erel%
223 %$$224 and is greater than or equal zero. 225 It specifies a bound for the relative 226 error in the return value icode xf$$ as an approximation for $latex X(tf)$$227 (see the 228 cref/error criteria discussion/OdeGearControl/Error Criteria Discussion/$$ 229 below). 230 231$head ef$$232 The argument value icode ef$$ has prototype
233 $codei% 234 %Vector% &%ef% 235 %$$236 and size icode n$$. 237 The input value of its elements does not matter. 238 On output, 239 it contains an estimated bound for the 240 absolute error in the approximation$icode xf$$; i.e., 241 latex $242 ef_i > | X( tf )_i - xf_i | 243$$$
244
245 $head maxabs$$246 The argument icode maxabs$$ is optional in the call to$code OdeGearControl$$. 247 If it is present, it has the prototype 248 codei% 249 %Vector% &%maxabs% 250 %$$
251 and size $icode n$$. 252 The input value of its elements does not matter. 253 On output, 254 it contains an estimate for the 255 maximum absolute value of latex X(t)$$; i.e., 256$latex $257 maxabs[i] \approx \max \left\{ 258 | X( t )_i | \; : \; t \in [ti, tf] 259 \right\} 260$ $$261 262 head nstep$$
263 The argument $icode nstep$$has the prototype 264 codei% 265 %size_t% &%nstep% 266 %$$ 267 Its input value does not matter and its output value 268 is the number of calls to$cref OdeGear$$269 used by code OdeGearControl$$.
270
271 $head Error Criteria Discussion$$272 The relative error criteria icode erel$$ and 273 absolute error criteria$icode eabs$$are enforced during each step of the 274 integration of the ordinary differential equations. 275 In addition, they are inversely scaled by the step size so that 276 the total error bound is less than the sum of the error bounds. 277 To be specific, if latex \tilde{X} (t)$$ is the approximate solution
278 at time $latex t$$, 279 icode ta$$ is the initial step time, 280 and$icode tb$$is the final step time, 281 latex $282 \left| \tilde{X} (tb)_j - X (tb)_j \right| 283 \leq 284 \frac{tf - ti}{tb - ta} 285 \left[ eabs[j] + erel \; | \tilde{X} (tb)_j | \right] 286$$$
287 If $latex X(tb)_j$$is near zero for some latex tb \in [ti , tf]$$, 288 and one uses an absolute error criteria$latex eabs[j]$$of zero, 289 the error criteria above will force code OdeGearControl$$
290 to use step sizes equal to
291 $cref/smin/OdeGearControl/smin/$$292 for steps ending near latex tb$$. 293 In this case, the error relative to$icode maxabs$$can be judged after 294 code OdeGearControl$$ returns.
295 If $icode ef$$is to large relative to icode maxabs$$, 296$code OdeGearControl$$can be called again 297 with a smaller value of icode smin$$.
298
299 $head Scalar$$300 The type icode Scalar$$ must satisfy the conditions 301 for a$cref NumericType$$type. 302 The routine cref CheckNumericType$$ will generate an error message
303 if this is not the case.
304 In addition, the following operations must be defined for
305 $icode Scalar$$objects icode a$$ and$icode b$$: 306 307 table 308 bold Operation$$ $cnext$bold Description$$rnext 309 icode%a% <= %b%$$ $cnext 310 returns true (false) if$icode a$$is less than or equal 311 (greater than) icode b$$.
312 $rnext 313$icode%a% == %b%$$cnext 314 returns true (false) if icode a$$ is equal to $icode b$$. 315 rnext 316 codei%log(%a%)%$$$cnext
317  returns a $icode Scalar$$equal to the logarithm of icode a$$ 318$rnext
319 $codei%exp(%a%)%$$cnext 320 returns a icode Scalar$$ equal to the exponential of$icode a$$321 tend 322 323 324 head Vector$$
325 The type $icode Vector$$must be a cref SimpleVector$$ class with 326$cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$. 327 The routine cref CheckSimpleVector$$ will generate an error message
328 if this is not the case.
329
330 $head Example$$331 children% 332 example/utility/ode_gear_control.cpp 333 %$$ 334 The file 335$cref ode_gear_control.cpp$$336 contains an example and test a test of using this routine. 337 It returns true if it succeeds and false otherwise. 338 339 head Theory$$
340 Let $latex e(s)$$be the error as a function of the 341 step size latex s$$ and suppose that there is a constant 342$latex K$$such that latex e(s) = K s^m$$.
343 Let $latex a$$be our error bound. 344 Given the value of latex e(s)$$, a step of size$latex \lambda s$$345 would be ok provided that 346 latex $347 \begin{array}{rcl} 348 a & \geq & e( \lambda s ) (tf - ti) / ( \lambda s ) \\ 349 a & \geq & K \lambda^m s^m (tf - ti) / ( \lambda s ) \\ 350 a & \geq & \lambda^{m-1} s^{m-1} (tf - ti) e(s) / s^m \\ 351 a & \geq & \lambda^{m-1} (tf - ti) e(s) / s \\ 352 \lambda^{m-1} & \leq & \frac{a}{e(s)} \frac{s}{tf - ti} 353 \end{array} 354$$$
355 Thus if the right hand side of the last inequality is greater
356 than or equal to one, the step of size $latex s$$is ok. 357 358 head Source Code$$ 359 The source code for this routine is in the file 360$code cppad/ode_gear_control.hpp.
361
362 \$end
363 --------------------------------------------------------------------------
364 */
365
366 // link exp and log for float and double
368
370
372
373 template <class Scalar, class Vector, class Fun>
375  Fun &F ,
376  size_t M ,
377  const Scalar &ti ,
378  const Scalar &tf ,
379  const Vector &xi ,
380  const Scalar &smin ,
381  const Scalar &smax ,
382  Scalar &sini ,
383  const Vector &eabs ,
384  const Scalar &erel ,
385  Vector &ef ,
386  Vector &maxabs,
387  size_t &nstep )
388 {
389  // check simple vector class specifications
390  CheckSimpleVector<Scalar, Vector>();
391
392  // dimension of the state space
393  size_t n = size_t(xi.size());
394
396  M >= 1,
397  "Error in OdeGearControl: M is less than one"
398  );
400  smin <= smax,
401  "Error in OdeGearControl: smin is greater than smax"
402  );
404  sini <= smax,
405  "Error in OdeGearControl: sini is greater than smax"
406  );
408  size_t(eabs.size()) == n,
409  "Error in OdeGearControl: size of eabs is not equal to n"
410  );
412  size_t(maxabs.size()) == n,
413  "Error in OdeGearControl: size of maxabs is not equal to n"
414  );
415
416  // some constants
417  const Scalar zero(0);
418  const Scalar one(1);
419  const Scalar one_plus( Scalar(3) / Scalar(2) );
420  const Scalar two(2);
421  const Scalar ten(10);
422
423  // temporary indices
424  size_t i, k;
425
426  // temporary Scalars
427  Scalar step, sprevious, lambda, axi, a, root, r;
428
429  // vectors of Scalars
430  Vector T (M + 1);
431  Vector X( (M + 1) * n );
432  Vector e(n);
433  Vector xf(n);
434
435  // initial integer values
436  size_t m = 1;
437  nstep = 0;
438
439  // initialize T
440  T[0] = ti;
441
442  // initialize X, ef, maxabs
443  for(i = 0; i < n; i++)
444  for(i = 0; i < n; i++)
445  { X[i] = xi[i];
446  ef[i] = zero;
447  X[i] = xi[i];
448  if( zero <= xi[i] )
449  maxabs[i] = xi[i];
450  else maxabs[i] = - xi[i];
451
452  }
453
454  // initial step size
455  step = smin;
456
457  while( T[m-1] < tf )
458  { sprevious = step;
459
460  // check maximum
461  if( smax <= step )
462  step = smax;
463
464  // check minimum
465  if( m < M )
466  { if( step <= sini )
467  step = sini;
468  }
469  else if( step <= smin )
470  step = smin;
471
472  // check if near the end
473  if( tf <= T[m-1] + one_plus * step )
474  T[m] = tf;
475  else T[m] = T[m-1] + step;
476
477  // try using this step size
478  nstep++;
479  OdeGear(F, m, n, T, X, e);
480  step = T[m] - T[m-1];
481
482  // compute value of lambda for this step
483  lambda = Scalar(10) * sprevious / step;
484  for(i = 0; i < n; i++)
485  { axi = X[m * n + i];
486  if( axi <= zero )
487  axi = - axi;
488  a = eabs[i] + erel * axi;
489  if( e[i] > zero )
490  { if( m == 1 )
491  root = (a / e[i]) / ten;
492  else
493  { r = ( a / e[i] ) * step / (tf - ti);
494  root = exp( log(r) / Scalar(m-1) );
495  }
496  if( root <= lambda )
497  lambda = root;
498  }
499  }
500
502  if( m == M )
503  advance = one <= lambda || step <= one_plus * smin;
504  else advance = one <= lambda || step <= one_plus * sini;
505
506
508  { // accept the results of this time step
509  CPPAD_ASSERT_UNKNOWN( m <= M );
510  if( m == M )
511  { // shift for next step
512  for(k = 0; k < m; k++)
513  { T[k] = T[k+1];
514  for(i = 0; i < n; i++)
515  X[k*n + i] = X[(k+1)*n + i];
516  }
517  }
518  // update ef and maxabs
519  for(i = 0; i < n; i++)
520  { ef[i] = ef[i] + e[i];
521  axi = X[m * n + i];
522  if( axi <= zero )
523  axi = - axi;
524  if( axi > maxabs[i] )
525  maxabs[i] = axi;
526  }
527  if( m != M )
528  m++; // all we need do in this case
529  }
530
531  // new step suggested by error criteria
532  step = std::min(lambda , ten) * step / two;
533  }
534  for(i = 0; i < n; i++)
535  xf[i] = X[(m-1) * n + i];
536
537  return xf;
538 }
539
540 } // End CppAD namespace
541
542 # endif