CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
ode_gear_control.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_UTILITY_ODE_GEAR_CONTROL_HPP
2 # define CPPAD_UTILITY_ODE_GEAR_CONTROL_HPP
3 
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6 
7 CppAD is distributed under multiple licenses. This distribution is under
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.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
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
367 # include <cppad/base_require.hpp>
368 
369 # include <cppad/utility/ode_gear.hpp>
370 
371 namespace CppAD { // Begin CppAD namespace
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 
501  bool advance;
502  if( m == M )
503  advance = one <= lambda || step <= one_plus * smin;
504  else advance = one <= lambda || step <= one_plus * sini;
505 
506 
507  if( advance )
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
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
void OdeGear(Fun &F, size_t m, size_t n, const Vector &T, Vector &X, Vector &e)
Definition: ode_gear.hpp:367
AD< Base > log(const AD< Base > &x)
AD< Base > exp(const AD< Base > &x)
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
Vector OdeGearControl(Fun &F, size_t M, const Scalar &ti, const Scalar &tf, const Vector &xi, const Scalar &smin, const Scalar &smax, Scalar &sini, const Vector &eabs, const Scalar &erel, Vector &ef, Vector &maxabs, size_t &nstep)