CppAD: A C++ Algorithmic Differentiation Package  20171217
rosen_34.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 Rosen34$$17 spell 18 cppad.hpp 19 bool 20 xf 21 templated 22 const 23 Rosenbrock 24 CppAD 25 xi 26 ti 27 tf 28 Karp 29 Rosen 30 Shampine 31 ind 32 dep 33$$ 34 35 36$section A 3rd and 4th Order Rosenbrock ODE Solver$$37 mindex Rosen34 solve stiff differential equation$$
38
39 $head Syntax$$40 codei%# include <cppad/utility/rosen_34.hpp> 41 %$$ 42$icode%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%)
43 %$$44 icode%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%, %e%) 45 %$$
46
47
48 $head Description$$49 This is an embedded 3rd and 4th order Rosenbrock ODE solver 50 (see Section 16.6 of cref/Numerical Recipes/Bib/Numerical Recipes/$$ 51 for a description of Rosenbrock ODE solvers). 52 In particular, we use the formulas taken from page 100 of 53$cref/Shampine, L.F./Bib/Shampine, L.F./$$54 (except that the fraction 98/108 has been correction to be 97/108). 55 pre 56 57$$
58 We use $latex n$$for the size of the vector icode xi$$. 59 Let$latex \B{R}$$denote the real numbers 60 and let latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ be a smooth function.
61 The return value $icode xf$$contains a 5th order 62 approximation for the value latex X(tf)$$ where 63$latex X : [ti , tf] \rightarrow \B{R}^n$$is defined by 64 the following initial value problem: 65 latex $66 \begin{array}{rcl} 67 X(ti) & = & xi \\ 68 X'(t) & = & F[t , X(t)] 69 \end{array} 70$$$
71 If your set of ordinary differential equations are not stiff
72 an explicit method may be better (perhaps $cref Runge45$$.) 73 74 head Include$$ 75 The file$code cppad/rosen_34.hpp$$is included by code cppad/cppad.hpp$$
76 but it can also be included separately with out the rest of
77 the $code CppAD$$routines. 78 79 head xf$$ 80 The return value$icode xf$$has the prototype 81 codei% 82 %Vector% %xf% 83 %$$
84 and the size of $icode xf$$is equal to icode n$$ 85 (see description of$cref/Vector/Rosen34/Vector/$$below). 86 latex $87 X(tf) = xf + O( h^5 ) 88$$$
89 where $latex h = (tf - ti) / M$$is the step size. 90 If icode xf$$ contains not a number$cref nan$$, 91 see the discussion of cref/f/Rosen34/Fun/Nan/$$.
92
93 $head Fun$$94 The class icode Fun$$ 95 and the object$icode F$$satisfy the prototype 96 codei% 97 %Fun% &%F% 98 %$$
99 This must support the following set of calls
100 $codei% 101 %F%.Ode(%t%, %x%, %f%) 102 %F%.Ode_ind(%t%, %x%, %f_t%) 103 %F%.Ode_dep(%t%, %x%, %f_x%) 104 %$$105 106 subhead t$$ 107 In all three cases, 108 the argument$icode t$$has prototype 109 codei% 110 const %Scalar% &%t% 111 %$$
112 (see description of $cref/Scalar/Rosen34/Scalar/$$below). 113 114 subhead x$$ 115 In all three cases, 116 the argument$icode x$$has prototype 117 codei% 118 const %Vector% &%x% 119 %$$
120 and has size $icode n$$121 (see description of cref/Vector/Rosen34/Vector/$$ below). 122 123$subhead f$$124 The argument icode f$$ to $icode%F%.Ode%$$has prototype 125 codei% 126 %Vector% &%f% 127 %$$ 128 On input and output,$icode f$$is a vector of size icode n$$
129 and the input values of the elements of $icode f$$do not matter. 130 On output, 131 icode f$$ is set equal to$latex F(t, x)$$132 (see icode F(t, x)$$ in $cref/Description/Rosen34/Description/$$). 133 134 subhead f_t$$ 135 The argument$icode f_t$$to icode%F%.Ode_ind%$$ has prototype
136 $codei% 137 %Vector% &%f_t% 138 %$$139 On input and output, icode f_t$$ is a vector of size$icode n$$140 and the input values of the elements of icode f_t$$ do not matter.
141 On output, the $th i$$element of 142 icode f_t$$ is set equal to$latex \partial_t F_i (t, x)$$143 (see icode F(t, x)$$ in $cref/Description/Rosen34/Description/$$). 144 145 subhead f_x$$ 146 The argument$icode f_x$$to icode%F%.Ode_dep%$$ has prototype
147 $codei% 148 %Vector% &%f_x% 149 %$$150 On input and output, icode f_x$$ is a vector of size$icode%n%*%n%$$151 and the input values of the elements of icode f_x$$ do not matter.
152 On output, the [$icode%i%*%n%+%j%$$] element of 153 icode f_x$$ is set equal to$latex \partial_{x(j)} F_i (t, x)$$154 (see icode F(t, x)$$ in $cref/Description/Rosen34/Description/$$). 155 156 subhead Nan$$ 157 If any of the elements of$icode f$$, icode f_t$$, or $icode f_x$$158 have the value not a number code nan$$, 159 the routine$code Rosen34$$returns with all the 160 elements of icode xf$$ and $icode e$$equal to code nan$$. 161 162$subhead Warning$$163 The arguments icode f$$, $icode f_t$$, and icode f_x$$ 164 must have a call by reference in their prototypes; i.e., 165 do not forget the$code &$$in the prototype for 166 icode f$$, $icode f_t$$and icode f_x$$. 167 168$subhead Optimization$$169 Every call of the form 170 codei% 171 %F%.Ode_ind(%t%, %x%, %f_t%) 172 %$$
173 is directly followed by a call of the form
174 $codei% 175 %F%.Ode_dep(%t%, %x%, %f_x%) 176 %$$177 where the arguments icode t$$ and$icode x$$have not changed between calls. 178 In many cases it is faster to compute the values of icode f_t$$
179 and $icode f_x$$together and then pass them back one at a time. 180 181 head M$$ 182 The argument$icode M$$has prototype 183 codei% 184 size_t %M% 185 %$$
186 It specifies the number of steps
187 to use when solving the differential equation.
188 This must be greater than or equal one.
189 The step size is given by $latex h = (tf - ti) / M$$, thus 190 the larger icode M$$, the more accurate the 191 return value$icode xf$$is as an approximation 192 for latex X(tf)$$.
193
194 $head ti$$195 The argument icode ti$$ has prototype 196$codei%
197  const %Scalar% &%ti%
198 %$$199 (see description of cref/Scalar/Rosen34/Scalar/$$ below).
200 It specifies the initial time for $icode t$$in the 201 differential equation; i.e., 202 the time corresponding to the value icode xi$$. 203 204$head tf$$205 The argument icode tf$$ has prototype
206 $codei% 207 const %Scalar% &%tf% 208 %$$209 It specifies the final time for icode t$$ in the 210 differential equation; i.e., 211 the time corresponding to the value$icode xf$$. 212 213 head xi$$
214 The argument $icode xi$$has the prototype 215 codei% 216 const %Vector% &%xi% 217 %$$ 218 and the size of$icode xi$$is equal to icode n$$.
219 It specifies the value of $latex X(ti)$$220 221 head e$$ 222 The argument$icode e$$is optional and has the prototype 223 codei% 224 %Vector% &%e% 225 %$$
226 If $icode e$$is present, 227 the size of icode e$$ must be equal to$icode n$$. 228 The input value of the elements of icode e$$ does not matter.
229 On output
230 it contains an element by element
231 estimated bound for the absolute value of the error in $icode xf$$232 latex $233 e = O( h^4 ) 234$$$ 235 where$latex h = (tf - ti) / M$$is the step size. 236 237 head Scalar$$
238 The type $icode Scalar$$must satisfy the conditions 239 for a cref NumericType$$ type. 240 The routine$cref CheckNumericType$$will generate an error message 241 if this is not the case. 242 In addition, the following operations must be defined for 243 icode Scalar$$ objects $icode a$$and icode b$$: 244 245$table
246 $bold Operation$$cnext bold Description$$$rnext
247 $icode%a% < %b%$$cnext 248 less than operator (returns a code bool$$ object) 249$tend
250
251 $head Vector$$252 The type icode Vector$$ must be a$cref SimpleVector$$class with 253 cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$.
254 The routine $cref CheckSimpleVector$$will generate an error message 255 if this is not the case. 256 257 head Parallel Mode$$ 258 For each set of types 259$cref/Scalar/Rosen34/Scalar/$$, 260 cref/Vector/Rosen34/Vector/$$, and
261 $cref/Fun/Rosen34/Fun/$$, 262 the first call to code Rosen34$$ 263 must not be$cref/parallel/ta_in_parallel/$$execution mode. 264 265 head Example$$
266 $children% 267 example/general/rosen_34.cpp 268 %$$269 The file 270 cref rosen_34.cpp$$ 271 contains an example and test a test of using this routine. 272 It returns true if it succeeds and false otherwise. 273 274$head Source Code$$275 The source code for this routine is in the file 276 code cppad/rosen_34.hpp$$.
277
278 \$end
279 --------------------------------------------------------------------------
280 */
281
282 # include <cstddef>
289
290 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
292
294
295 template <typename Scalar, typename Vector, typename Fun>
296 Vector Rosen34(
297  Fun &F ,
298  size_t M ,
299  const Scalar &ti ,
300  const Scalar &tf ,
301  const Vector &xi )
302 { Vector e( xi.size() );
303  return Rosen34(F, M, ti, tf, xi, e);
304 }
305
306 template <typename Scalar, typename Vector, typename Fun>
307 Vector Rosen34(
308  Fun &F ,
309  size_t M ,
310  const Scalar &ti ,
311  const Scalar &tf ,
312  const Vector &xi ,
313  Vector &e )
314 {
316
317  // check numeric type specifications
318  CheckNumericType<Scalar>();
319
320  // check simple vector class specifications
321  CheckSimpleVector<Scalar, Vector>();
322
323  // Parameters for Shampine's Rosenbrock method
324  // are static to avoid recalculation on each call and
325  // do not use Vector to avoid possible memory leak
326  static Scalar a[3] = {
327  Scalar(0),
328  Scalar(1),
329  Scalar(3) / Scalar(5)
330  };
331  static Scalar b[2 * 2] = {
332  Scalar(1),
333  Scalar(0),
334  Scalar(24) / Scalar(25),
335  Scalar(3) / Scalar(25)
336  };
337  static Scalar ct[4] = {
338  Scalar(1) / Scalar(2),
339  - Scalar(3) / Scalar(2),
340  Scalar(121) / Scalar(50),
341  Scalar(29) / Scalar(250)
342  };
343  static Scalar cg[3 * 3] = {
344  - Scalar(4),
345  Scalar(0),
346  Scalar(0),
347  Scalar(186) / Scalar(25),
348  Scalar(6) / Scalar(5),
349  Scalar(0),
350  - Scalar(56) / Scalar(125),
351  - Scalar(27) / Scalar(125),
352  - Scalar(1) / Scalar(5)
353  };
354  static Scalar d3[3] = {
355  Scalar(97) / Scalar(108),
356  Scalar(11) / Scalar(72),
357  Scalar(25) / Scalar(216)
358  };
359  static Scalar d4[4] = {
360  Scalar(19) / Scalar(18),
361  Scalar(1) / Scalar(4),
362  Scalar(25) / Scalar(216),
363  Scalar(125) / Scalar(216)
364  };
366  M >= 1,
367  "Error in Rosen34: the number of steps is less than one"
368  );
370  e.size() == xi.size(),
371  "Error in Rosen34: size of e not equal to size of xi"
372  );
373  size_t i, j, k, l, m; // indices
374
375  size_t n = xi.size(); // number of components in X(t)
376  Scalar ns = Scalar(double(M)); // number of steps as Scalar object
377  Scalar h = (tf - ti) / ns; // step size
378  Scalar zero = Scalar(0); // some constants
379  Scalar one = Scalar(1);
380  Scalar two = Scalar(2);
381
382  // permutation vectors needed for LU factorization routine
384
385  // vectors used to store values returned by F
386  Vector E(n * n), Eg(n), f_t(n);
387  Vector g(n * 3), x3(n), x4(n), xf(n), ftmp(n), xtmp(n), nan_vec(n);
388
389  // initialize e = 0, nan_vec = nan
390  for(i = 0; i < n; i++)
391  { e[i] = zero;
392  nan_vec[i] = nan(zero);
393  }
394
395  xf = xi; // initialize solution
396  for(m = 0; m < M; m++)
397  { // time at beginning of this interval
398  Scalar t = ti * (Scalar(int(M - m)) / ns)
399  + tf * (Scalar(int(m)) / ns);
400
401  // value of x at beginning of this interval
402  x3 = x4 = xf;
403
404  // evaluate partial derivatives at beginning of this interval
405  F.Ode_ind(t, xf, f_t);
406  F.Ode_dep(t, xf, E); // E = f_x
407  if( hasnan(f_t) || hasnan(E) )
408  { e = nan_vec;
409  return nan_vec;
410  }
411
412  // E = I - f_x * h / 2
413  for(i = 0; i < n; i++)
414  { for(j = 0; j < n; j++)
415  E[i * n + j] = - E[i * n + j] * h / two;
416  E[i * n + i] += one;
417  }
418
419  // LU factor the matrix E
420 # ifndef NDEBUG
421  int sign = LuFactor(ip, jp, E);
422 # else
423  LuFactor(ip, jp, E);
424 # endif
426  sign != 0,
427  "Error in Rosen34: I - f_x * h / 2 not invertible"
428  );
429
430  // loop over integration steps
431  for(k = 0; k < 3; k++)
432  { // set location for next function evaluation
433  xtmp = xf;
434  for(l = 0; l < k; l++)
435  { // loop over previous function evaluations
436  Scalar bkl = b[(k-1)*2 + l];
437  for(i = 0; i < n; i++)
438  { // loop over elements of x
439  xtmp[i] += bkl * g[i*3 + l] * h;
440  }
441  }
442  // ftmp = F(t + a[k] * h, xtmp)
443  F.Ode(t + a[k] * h, xtmp, ftmp);
444  if( hasnan(ftmp) )
445  { e = nan_vec;
446  return nan_vec;
447  }
448
449  // Form Eg for this integration step
450  for(i = 0; i < n; i++)
451  Eg[i] = ftmp[i] + ct[k] * f_t[i] * h;
452  for(l = 0; l < k; l++)
453  { for(i = 0; i < n; i++)
454  Eg[i] += cg[(k-1)*3 + l] * g[i*3 + l];
455  }
456
457  // Solve the equation E * g = Eg
458  LuInvert(ip, jp, E, Eg);
459
460  // save solution and advance x3, x4
461  for(i = 0; i < n; i++)
462  { g[i*3 + k] = Eg[i];
463  x3[i] += h * d3[k] * Eg[i];
464  x4[i] += h * d4[k] * Eg[i];
465  }
466  }
467  // Form Eg for last update to x4 only
468  for(i = 0; i < n; i++)
469  Eg[i] = ftmp[i] + ct[3] * f_t[i] * h;
470  for(l = 0; l < 3; l++)
471  { for(i = 0; i < n; i++)
472  Eg[i] += cg[2*3 + l] * g[i*3 + l];
473  }
474
475  // Solve the equation E * g = Eg
476  LuInvert(ip, jp, E, Eg);
477
478  // advance x4 and accumulate error bound
479  for(i = 0; i < n; i++)
480  { x4[i] += h * d4[3] * Eg[i];
481
482  // cant use abs because cppad.hpp may not be included
483  Scalar diff = x4[i] - x3[i];
484  if( diff < zero )
485  e[i] -= diff;
486  else e[i] += diff;
487  }
488
489  // advance xf for this step using x4
490  xf = x4;
491  }
492  return xf;
493 }
494
495 } // End CppAD namespace
496
497 # endif
Check that exp is true, if not print msg and terminate execution.
Scalar nan(const Scalar &zero)
Definition: nan.hpp:187
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
void LuInvert(const SizeVector &ip, const SizeVector &jp, const FloatVector &LU, FloatVector &B)
Definition: lu_invert.hpp:169