CppAD: A C++ Algorithmic Differentiation Package  20171217
runge_45.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_UTILITY_RUNGE_45_HPP
2 # define CPPAD_UTILITY_RUNGE_45_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 Runge45$$17 spell 18 std 19 fabs 20 cppad.hpp 21 bool 22 xf 23 templated 24 const 25 Runge-Kutta 26 CppAD 27 xi 28 ti 29 tf 30 Karp 31$$ 32 33 34$section An Embedded 4th and 5th Order Runge-Kutta ODE Solver$$35 mindex Runge45 Runge Kutta solve differential equation$$
36
37 $head Syntax$$38 codei%# include <cppad/utility/runge_45.hpp> 39 %$$ 40$icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%)
41 %$$42 icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%, %e%) 43 %$$
44
45
46 $head Purpose$$47 This is an implementation of the 48 Cash-Karp embedded 4th and 5th order Runge-Kutta ODE solver 49 described in Section 16.2 of cref/Numerical Recipes/Bib/Numerical Recipes/$$. 50 We use$latex n$$for the size of the vector icode xi$$.
51 Let $latex \B{R}$$denote the real numbers 52 and let latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ 53 be a smooth function. 54 The return value$icode xf$$contains a 5th order 55 approximation for the value latex X(tf)$$ where
56 $latex X : [ti , tf] \rightarrow \B{R}^n$$is defined by 57 the following initial value problem: 58 latex $59 \begin{array}{rcl} 60 X(ti) & = & xi \\ 61 X'(t) & = & F[t , X(t)] 62 \end{array} 63$$$ 64 If your set of ordinary differential equations 65 are stiff, an implicit method may be better 66 (perhaps$cref Rosen34$$.) 67 68 head Operation Sequence$$
69 The $cref/operation sequence/glossary/Operation/Sequence/$$for icode Runge$$ 70 does not depend on any of its$icode Scalar$$input values provided that 71 the operation sequence for 72 codei% 73 %F%.Ode(%t%, %x%, %f%) 74 %$$
75 does not on any of its $icode Scalar$$inputs (see below). 76 77 head Include$$ 78 The file$code cppad/runge_45.hpp$$is included by code cppad/cppad.hpp$$
79 but it can also be included separately with out the rest of
80 the $code CppAD$$routines. 81 82 head xf$$ 83 The return value$icode xf$$has the prototype 84 codei% 85 %Vector% %xf% 86 %$$
87 and the size of $icode xf$$is equal to icode n$$ 88 (see description of$cref/Vector/Runge45/Vector/$$below). 89 latex $90 X(tf) = xf + O( h^6 ) 91$$$
92 where $latex h = (tf - ti) / M$$is the step size. 93 If icode xf$$ contains not a number$cref nan$$, 94 see the discussion for cref/f/Runge45/Fun/f/$$.
95
96 $head Fun$$97 The class icode Fun$$ 98 and the object$icode F$$satisfy the prototype 99 codei% 100 %Fun% &%F% 101 %$$
102 The object $icode F$$(and the class icode Fun$$) 103 must have a member function named$code Ode$$104 that supports the syntax 105 codei% 106 %F%.Ode(%t%, %x%, %f%) 107 %$$
108
109 $subhead t$$110 The argument icode t$$ to$icode%F%.Ode%$$has prototype 111 codei% 112 const %Scalar% &%t% 113 %$$
114 (see description of $cref/Scalar/Runge45/Scalar/$$below). 115 116 subhead x$$ 117 The argument$icode x$$to icode%F%.Ode%$$ has prototype
118 $codei% 119 const %Vector% &%x% 120 %$$121 and has size icode n$$ 122 (see description of$cref/Vector/Runge45/Vector/$$below). 123 124 subhead f$$
125 The argument $icode f$$to icode%F%.Ode%$$ has prototype 126$codei%
127  %Vector% &%f%
128 %$$129 On input and output, icode f$$ is a vector of size $icode n$$130 and the input values of the elements of icode f$$ do not matter. 131 On output, 132$icode f$$is set equal to latex F(t, x)$$ in the differential equation.
133 If any of the elements of $icode f$$have the value not a number code nan$$ 134 the routine$code Runge45$$returns with all the 135 elements of icode xf$$ and $icode e$$equal to code nan$$. 136 137$subhead Warning$$138 The argument icode f$$ to $icode%F%.Ode%$$139 must have a call by reference in its prototype; i.e., 140 do not forget the code &$$ in the prototype for$icode f$$. 141 142 head M$$
143 The argument $icode M$$has prototype 144 codei% 145 size_t %M% 146 %$$ 147 It specifies the number of steps 148 to use when solving the differential equation. 149 This must be greater than or equal one. 150 The step size is given by$latex h = (tf - ti) / M$$, thus 151 the larger icode M$$, the more accurate the
152 return value $icode xf$$is as an approximation 153 for latex X(tf)$$. 154 155$head ti$$156 The argument icode ti$$ has prototype
157 $codei% 158 const %Scalar% &%ti% 159 %$$160 (see description of cref/Scalar/Runge45/Scalar/$$ below). 161 It specifies the initial time for$icode t$$in the 162 differential equation; i.e., 163 the time corresponding to the value icode xi$$.
164
165 $head tf$$166 The argument icode tf$$ has prototype 167$codei%
168  const %Scalar% &%tf%
169 %$$170 It specifies the final time for icode t$$ in the
171 differential equation; i.e.,
172 the time corresponding to the value $icode xf$$. 173 174 head xi$$ 175 The argument$icode xi$$has the prototype 176 codei% 177 const %Vector% &%xi% 178 %$$
179 and the size of $icode xi$$is equal to icode n$$. 180 It specifies the value of$latex X(ti)$$181 182 head e$$
183 The argument $icode e$$is optional and has the prototype 184 codei% 185 %Vector% &%e% 186 %$$ 187 If$icode e$$is present, 188 the size of icode e$$ must be equal to $icode n$$. 189 The input value of the elements of icode e$$ does not matter. 190 On output 191 it contains an element by element 192 estimated bound for the absolute value of the error in$icode xf$$193 latex $194 e = O( h^5 ) 195$$$
196 where $latex h = (tf - ti) / M$$is the step size. 197 If on output, icode e$$ contains not a number$code nan$$, 198 see the discussion for cref/f/Runge45/Fun/f/$$.
199
200 $head Scalar$$201 The type icode Scalar$$ must satisfy the conditions 202 for a$cref NumericType$$type. 203 The routine cref CheckNumericType$$ will generate an error message
204 if this is not the case.
205
206 $subhead fabs$$207 In addition, the following function must be defined for 208 icode Scalar$$ objects$icode a$$and icode b$$
209 $codei% 210 %a% = fabs(%b%) 211 %$$212 Note that this operation is only used for computing icode e$$; hence 213 the operation sequence for$icode xf$$can still be independent of 214 the arguments to code Runge45$$ even if
215 $codei% 216 fabs(%b%) = std::max(-%b%, %b%) 217 %$$. 218 219 head Vector$$ 220 The type$icode Vector$$must be a cref SimpleVector$$ class with
221 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$. 222 The routine cref CheckSimpleVector$$ will generate an error message 223 if this is not the case. 224 225$head Parallel Mode$$226 For each set of types 227 cref/Scalar/Runge45/Scalar/$$,
228 $cref/Vector/Runge45/Vector/$$, and 229 cref/Fun/Runge45/Fun/$$, 230 the first call to$code Runge45$$231 must not be cref/parallel/ta_in_parallel/$$ execution mode.
232
233
234 $head Example$$235 children% 236 example/utility/runge45_1.cpp% 237 example/general/runge45_2.cpp 238 %$$ 239 The file 240$cref runge45_1.cpp$$241 contains a simple example and test of code Runge45$$.
242 It returns true if it succeeds and false otherwise.
243 $pre 244 245 $$246 The file 247 cref runge45_2.cpp$$ contains an example using$code Runge45$$248 in the context of algorithmic differentiation. 249 It also returns true if it succeeds and false otherwise. 250 251 head Source Code$$
252 The source code for this routine is in the file
253 $code cppad/runge_45.hpp$$. 254 255$end
256 --------------------------------------------------------------------------
257 */
258 # include <cstddef>
259 # include <cppad/core/cppad_assert.hpp>
262 # include <cppad/utility/nan.hpp>
263
264 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
266
267 namespace CppAD { // BEGIN CppAD namespace
268
269 template <typename Scalar, typename Vector, typename Fun>
270 Vector Runge45(
271  Fun &F ,
272  size_t M ,
273  const Scalar &ti ,
274  const Scalar &tf ,
275  const Vector &xi )
276 { Vector e( xi.size() );
277  return Runge45(F, M, ti, tf, xi, e);
278 }
279
280 template <typename Scalar, typename Vector, typename Fun>
281 Vector Runge45(
282  Fun &F ,
283  size_t M ,
284  const Scalar &ti ,
285  const Scalar &tf ,
286  const Vector &xi ,
287  Vector &e )
288 {
290
291  // check numeric type specifications
292  CheckNumericType<Scalar>();
293
294  // check simple vector class specifications
295  CheckSimpleVector<Scalar, Vector>();
296
297  // Cash-Karp parameters for embedded Runge-Kutta method
298  // are static to avoid recalculation on each call and
299  // do not use Vector to avoid possible memory leak
300  static Scalar a[6] = {
301  Scalar(0),
302  Scalar(1) / Scalar(5),
303  Scalar(3) / Scalar(10),
304  Scalar(3) / Scalar(5),
305  Scalar(1),
306  Scalar(7) / Scalar(8)
307  };
308  static Scalar b[5 * 5] = {
309  Scalar(1) / Scalar(5),
310  Scalar(0),
311  Scalar(0),
312  Scalar(0),
313  Scalar(0),
314
315  Scalar(3) / Scalar(40),
316  Scalar(9) / Scalar(40),
317  Scalar(0),
318  Scalar(0),
319  Scalar(0),
320
321  Scalar(3) / Scalar(10),
322  -Scalar(9) / Scalar(10),
323  Scalar(6) / Scalar(5),
324  Scalar(0),
325  Scalar(0),
326
327  -Scalar(11) / Scalar(54),
328  Scalar(5) / Scalar(2),
329  -Scalar(70) / Scalar(27),
330  Scalar(35) / Scalar(27),
331  Scalar(0),
332
333  Scalar(1631) / Scalar(55296),
334  Scalar(175) / Scalar(512),
335  Scalar(575) / Scalar(13824),
336  Scalar(44275) / Scalar(110592),
337  Scalar(253) / Scalar(4096)
338  };
339  static Scalar c4[6] = {
340  Scalar(2825) / Scalar(27648),
341  Scalar(0),
342  Scalar(18575) / Scalar(48384),
343  Scalar(13525) / Scalar(55296),
344  Scalar(277) / Scalar(14336),
345  Scalar(1) / Scalar(4),
346  };
347  static Scalar c5[6] = {
348  Scalar(37) / Scalar(378),
349  Scalar(0),
350  Scalar(250) / Scalar(621),
351  Scalar(125) / Scalar(594),
352  Scalar(0),
353  Scalar(512) / Scalar(1771)
354  };
355
357  M >= 1,
358  "Error in Runge45: the number of steps is less than one"
359  );
361  e.size() == xi.size(),
362  "Error in Runge45: size of e not equal to size of xi"
363  );
364  size_t i, j, k, m; // indices
365
366  size_t n = xi.size(); // number of components in X(t)
367  Scalar ns = Scalar(int(M)); // number of steps as Scalar object
368  Scalar h = (tf - ti) / ns; // step size
369  Scalar zero_or_nan = Scalar(0); // zero (nan if Ode returns has a nan)
370  for(i = 0; i < n; i++) // initialize e = 0
371  e[i] = zero_or_nan;
372
373  // vectors used to store values returned by F
374  Vector fh(6 * n), xtmp(n), ftmp(n), x4(n), x5(n), xf(n);
375
376  xf = xi; // initialize solution
377  for(m = 0; m < M; m++)
378  { // time at beginning of this interval
379  // (convert to int to avoid MS compiler warning)
380  Scalar t = ti * (Scalar(int(M - m)) / ns)
381  + tf * (Scalar(int(m)) / ns);
382
383  // loop over integration steps
384  x4 = x5 = xf; // start x4 and x5 at same point for each step
385  for(j = 0; j < 6; j++)
386  { // loop over function evaluations for this step
387  xtmp = xf; // location for next function evaluation
388  for(k = 0; k < j; k++)
389  { // loop over previous function evaluations
390  Scalar bjk = b[ (j-1) * 5 + k ];
391  for(i = 0; i < n; i++)
392  { // loop over elements of x
393  xtmp[i] += bjk * fh[i * 6 + k];
394  }
395  }
396  // ftmp = F(t + a[j] * h, xtmp)
397  F.Ode(t + a[j] * h, xtmp, ftmp);
398
399  // if ftmp has a nan, set zero_or_nan to nan
400  for(i = 0; i < n; i++)
401  zero_or_nan *= ftmp[i];
402
403  for(i = 0; i < n; i++)
404  { // loop over elements of x
405  Scalar fhi = ftmp[i] * h;
406  fh[i * 6 + j] = fhi;
407  x4[i] += c4[j] * fhi;
408  x5[i] += c5[j] * fhi;
409  x5[i] += zero_or_nan;
410  }
411  }
412  // accumulate error bound
413  for(i = 0; i < n; i++)
414  { // cant use abs because cppad.hpp may not be included
415  Scalar diff = x5[i] - x4[i];
416  e[i] += fabs(diff);
417  e[i] += zero_or_nan;
418  }
419
420  // advance xf for this step using x5
421  xf = x5;
422  }
423  return xf;
424 }
425
426 } // End CppAD namespace
427
428 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
#define CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
Check that the first call to a routine is not during parallel execution mode.
Vector Runge45(Fun &F, size_t M, const Scalar &ti, const Scalar &tf, const Vector &xi)
Definition: runge_45.hpp:270
File used to define the CppAD multi-threading allocator class.
std::complex< double > fabs(const std::complex< double > &x)