# ifndef CPPAD_RUNGE_45_INCLUDED # define CPPAD_RUNGE_45_INCLUDED /* -------------------------------------------------------------------------- CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-08 Bradley M. Bell CppAD is distributed under multiple licenses. This distribution is under the terms of the Common Public License Version 1.0. A copy of this license is included in the COPYING file of this distribution. Please visit http://www.coin-or.org/CppAD/ for information on other licenses. -------------------------------------------------------------------------- */ /* $begin Runge45$$spell cppad.hpp bool xf templated const Runge-Kutta CppAD xi ti tf Karp$$$index Runge45$$index ODE, Runge-Kutta$$ $index Runge, ODE$$index Kutta, ODE$$$index solve, ODE$$index differential, equation$$ $index equation, differential$$section An Embedded 4th and 5th Order Runge-Kutta ODE Solver$$$head Syntax$$syntax%# include %$$ $syntax%%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%) %$$syntax%%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%, %e%) %$$$head Purpose$$This is an implementation of the Cash-Karp embedded 4th and 5th order Runge-Kutta ODE solver described in Section 16.2 of xref/Bib/Numerical Recipes/Numerical Recipes/$$. We use $latex n$$for the size of the vector italic xi$$. Let$latex \R$$denote the real numbers and let latex F : \R \times \R^n \rightarrow \R^n$$ be a smooth function. The return value $italic xf$$contains a 5th order approximation for the value latex X(tf)$$ where$latex X : [ti , tf] \rightarrow \R^n$$is defined by the following initial value problem: latex $\begin{array}{rcl} X(ti) & = & xi \\ X'(t) & = & F[t , X(t)] \end{array}$$$ If your set of ordinary differential equations are stiff, an implicit method may be better (perhaps $xref/Rosen34/$$.) head Include$$ The file$code cppad/runge_45.hpp$$is included by code cppad/cppad.hpp$$ but it can also be included separately with out the rest of the $code CppAD$$routines. head xf$$ The return value$italic xf$$has the prototype syntax% %Vector% %xf% %$$ and the size of $italic xf$$is equal to italic n$$ (see description of$xref/Runge45/Vector/Vector/$$below). latex $X(tf) = xf + O( h^6 )$$$ where $latex h = (tf - ti) / M$$is the step size. If italic xf$$ contains not a number$cref/nan/$$, see the discussion for cref/f/Runge45/Fun/f/$$. $head Fun$$The class italic Fun$$ and the object$italic F$$satisfy the prototype syntax% %Fun% &%F% %$$ The object $italic F$$(and the class italic Fun$$) must have a member function named$code Ode$$that supports the syntax syntax% %F%.Ode(%t%, %x%, %f%) %$$ $subhead t$$The argument italic t$$ to$syntax%%F%.Ode%$$has prototype syntax% const %Scalar% &%t% %$$ (see description of $xref/Runge45/Scalar/Scalar/$$below). subhead x$$ The argument$italic x$$to syntax%%F%.Ode%$$ has prototype $syntax% const %Vector% &%x% %$$and has size italic n$$ (see description of$xref/Runge45/Vector/Vector/$$below). subhead f$$ The argument $italic f$$to syntax%%F%.Ode%$$ has prototype$syntax% %Vector% &%f% %$$On input and output, italic f$$ is a vector of size $italic n$$and the input values of the elements of italic f$$ do not matter. On output,$italic f$$is set equal to latex F(t, x)$$ in the differential equation. If any of the elements of $italic f$$have the value not a number code nan$$ the routine$code Runge45$$returns with all the elements of italic xf$$ and $italic e$$equal to code nan$$.$subhead Warning$$The argument italic f$$ to $syntax%%F%.Ode%$$must have a call by reference in its prototype; i.e., do not forget the code &$$ in the prototype for$italic f$$. head M$$ The argument $italic M$$has prototype syntax% size_t %M% %$$ It specifies the number of steps to use when solving the differential equation. This must be greater than or equal one. The step size is given by$latex h = (tf - ti) / M$$, thus the larger italic M$$, the more accurate the return value $italic xf$$is as an approximation for latex X(tf)$$.$head ti$$The argument italic ti$$ has prototype $syntax% const %Scalar% &%ti% %$$(see description of xref/Runge45/Scalar/Scalar/$$ below). It specifies the initial time for$italic t$$in the differential equation; i.e., the time corresponding to the value italic xi$$. $head tf$$The argument italic tf$$ has prototype$syntax% const %Scalar% &%tf% %$$It specifies the final time for italic t$$ in the differential equation; i.e., the time corresponding to the value $italic xf$$. head xi$$ The argument$italic xi$$has the prototype syntax% const %Vector% &%xi% %$$ and the size of $italic xi$$is equal to italic n$$. It specifies the value of$latex X(ti)$$head e$$ The argument $italic e$$is optional and has the prototype syntax% %Vector% &%e% %$$ If$italic e$$is present, the size of italic e$$ must be equal to $italic n$$. The input value of the elements of italic e$$ does not matter. On output it contains an element by element estimated bound for the absolute value of the error in$italic xf$$latex $e = O( h^5 )$$$ where $latex h = (tf - ti) / M$$is the step size. If on output, italic e$$ contains not a number$code nan$$, see the discussion for cref/f/Runge45/Fun/f/$$. $head Scalar$$The type italic Scalar$$ must satisfy the conditions for a$xref/NumericType/$$type. The routine xref/CheckNumericType/$$ will generate an error message if this is not the case. In addition, the following operations must be defined for $italic Scalar$$objects italic a$$ and$italic b$$: table bold Operation$$ $cnext$bold Description$$rnext syntax%%a% < %b%$$ $cnext less than operator (returns a$code bool$$object) tend head Vector$$ The type $italic Vector$$must be a xref/SimpleVector/$$ class with$xref/SimpleVector/Elements of Specified Type/elements of type Scalar/$$. The routine xref/CheckSimpleVector/$$ will generate an error message if this is not the case. $head Example$$children% example/runge_45.cpp %$$ The file$xref/Runge45.cpp/$$contains an example and test a test of using this routine. It returns true if it succeeds and false otherwise. head Source Code$$ The source code for this routine is in the file $code cppad/runge_45.hpp$$.$end -------------------------------------------------------------------------- */ # include # include # include # include # include namespace CppAD { // BEGIN CppAD namespace template Vector Runge45( Fun &F , size_t M , const Scalar &ti , const Scalar &tf , const Vector &xi ) { Vector e( xi.size() ); return Runge45(F, M, ti, tf, xi, e); } template Vector Runge45( Fun &F , size_t M , const Scalar &ti , const Scalar &tf , const Vector &xi , Vector &e ) { // check numeric type specifications CheckNumericType(); // check simple vector class specifications CheckSimpleVector(); // Cash-Karp parameters for embedded Runge-Kutta method // are static to avoid recalculation on each call and // do not use Vector to avoid possible memory leak static Scalar a[6] = { Scalar(0), Scalar(1) / Scalar(5), Scalar(3) / Scalar(10), Scalar(3) / Scalar(5), Scalar(1), Scalar(7) / Scalar(8) }; static Scalar b[5 * 5] = { Scalar(1) / Scalar(5), Scalar(0), Scalar(0), Scalar(0), Scalar(0), Scalar(3) / Scalar(40), Scalar(9) / Scalar(40), Scalar(0), Scalar(0), Scalar(0), Scalar(3) / Scalar(10), -Scalar(9) / Scalar(10), Scalar(6) / Scalar(5), Scalar(0), Scalar(0), -Scalar(11) / Scalar(54), Scalar(5) / Scalar(2), -Scalar(70) / Scalar(27), Scalar(35) / Scalar(27), Scalar(0), Scalar(1631) / Scalar(55296), Scalar(175) / Scalar(512), Scalar(575) / Scalar(13824), Scalar(44275) / Scalar(110592), Scalar(253) / Scalar(4096) }; static Scalar c4[6] = { Scalar(2825) / Scalar(27648), Scalar(0), Scalar(18575) / Scalar(48384), Scalar(13525) / Scalar(55296), Scalar(277) / Scalar(14336), Scalar(1) / Scalar(4), }; static Scalar c5[6] = { Scalar(37) / Scalar(378), Scalar(0), Scalar(250) / Scalar(621), Scalar(125) / Scalar(594), Scalar(0), Scalar(512) / Scalar(1771) }; CPPAD_ASSERT_KNOWN( M >= 1, "Error in Runge45: the number of steps is less than one" ); CPPAD_ASSERT_KNOWN( e.size() == xi.size(), "Error in Runge45: size of e not equal to size of xi" ); size_t i, j, k, m; // indices size_t n = xi.size(); // number of components in X(t) Scalar ns = Scalar(int(M)); // number of steps as Scalar object Scalar h = (tf - ti) / ns; // step size Scalar z = Scalar(0); // zero for(i = 0; i < n; i++) // initialize e = 0 e[i] = z; // vectors used to store values returned by F Vector fh(6 * n), xtmp(n), ftmp(n), x4(n), x5(n), xf(n), nan_vec(n); // vector of nans for(i = 0; i < n; i++) nan_vec[i] = nan(z); xf = xi; // initialize solution for(m = 0; m < M; m++) { // time at beginning of this interval // (convert to int to avoid MS compiler warning) Scalar t = ti * (Scalar(int(M - m)) / ns) + tf * (Scalar(int(m)) / ns); // loop over integration steps x4 = x5 = xf; // start x4 and x5 at same point for each step for(j = 0; j < 6; j++) { // loop over function evaluations for this step xtmp = xf; // location for next function evaluation for(k = 0; k < j; k++) { // loop over previous function evaluations Scalar bjk = b[ (j-1) * 5 + k ]; for(i = 0; i < n; i++) { // loop over elements of x xtmp[i] += bjk * fh[i * 6 + k]; } } // ftmp = F(t + a[j] * h, xtmp) F.Ode(t + a[j] * h, xtmp, ftmp); if( hasnan(ftmp) ) { e = nan_vec; return nan_vec; } for(i = 0; i < n; i++) { // loop over elements of x Scalar fhi = ftmp[i] * h; fh[i * 6 + j] = fhi; x4[i] += c4[j] * fhi; x5[i] += c5[j] * fhi; } } // accumulate error bound for(i = 0; i < n; i++) { // cant use abs because cppad.hpp may not be included Scalar diff = x5[i] - x4[i]; if( diff < z ) e[i] -= diff; else e[i] += diff; } // advance xf for this step using x5 xf = x5; } return xf; } } // End CppAD namespace # endif