CppAD: A C++ Algorithmic Differentiation Package  20130102
runge_45.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_RUNGE_45_INCLUDED
00003 # define CPPAD_RUNGE_45_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin Runge45$$
00018 $spell
00019      std
00020      fabs
00021      cppad.hpp
00022      bool
00023      xf
00024      templated
00025      const
00026      Runge-Kutta
00027      CppAD
00028      xi
00029      ti
00030      tf
00031      Karp
00032 $$
00033 
00034 $index Runge45$$
00035 $index ODE, Runge-Kutta$$
00036 $index Runge, ODE$$
00037 $index Kutta, ODE$$
00038 $index solve, ODE$$
00039 $index differential, equation$$
00040 $index equation, differential$$
00041  
00042 $section An Embedded 4th and 5th Order Runge-Kutta ODE Solver$$
00043 
00044 $head Syntax$$
00045 $codei%# include <cppad/runge_45.hpp>
00046 %$$
00047 $icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%)
00048 %$$
00049 $icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%, %e%)
00050 %$$
00051 
00052 
00053 $head Purpose$$
00054 This is an implementation of the
00055 Cash-Karp embedded 4th and 5th order Runge-Kutta ODE solver 
00056 described in Section 16.2 of $cref/Numerical Recipes/Bib/Numerical Recipes/$$.
00057 We use $latex n$$ for the size of the vector $icode xi$$.
00058 Let $latex \B{R}$$ denote the real numbers
00059 and let $latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ 
00060 be a smooth function.
00061 The return value $icode xf$$ contains a 5th order
00062 approximation for the value $latex X(tf)$$ where 
00063 $latex X : [ti , tf] \rightarrow \B{R}^n$$ is defined by 
00064 the following initial value problem:
00065 $latex \[
00066 \begin{array}{rcl}
00067      X(ti)  & = & xi    \\
00068      X'(t)  & = & F[t , X(t)] 
00069 \end{array}
00070 \] $$
00071 If your set of ordinary differential equations
00072 are stiff, an implicit method may be better
00073 (perhaps $cref Rosen34$$.)
00074 
00075 $head Operation Sequence$$
00076 The $cref/operation sequence/glossary/Operation/Sequence/$$ for $icode Runge$$
00077 does not depend on any of its $icode Scalar$$ input values provided that
00078 the operation sequence for
00079 $codei%
00080      %F%.Ode(%t%, %x%, %f%)
00081 %$$
00082 does not on any of its $icode Scalar$$ inputs (see below).
00083 
00084 $head Include$$
00085 The file $code cppad/runge_45.hpp$$ is included by $code cppad/cppad.hpp$$
00086 but it can also be included separately with out the rest of 
00087 the $code CppAD$$ routines.
00088 
00089 $head xf$$
00090 The return value $icode xf$$ has the prototype
00091 $codei%
00092      %Vector% %xf%
00093 %$$
00094 and the size of $icode xf$$ is equal to $icode n$$ 
00095 (see description of $cref/Vector/Runge45/Vector/$$ below).
00096 $latex \[
00097      X(tf) = xf + O( h^6 )
00098 \] $$
00099 where $latex h = (tf - ti) / M$$ is the step size.
00100 If $icode xf$$ contains not a number $cref nan$$,
00101 see the discussion for $cref/f/Runge45/Fun/f/$$.
00102 
00103 $head Fun$$
00104 The class $icode Fun$$ 
00105 and the object $icode F$$ satisfy the prototype
00106 $codei%
00107      %Fun% &%F%
00108 %$$
00109 The object $icode F$$ (and the class $icode Fun$$)
00110 must have a member function named $code Ode$$ 
00111 that supports the syntax
00112 $codei%
00113      %F%.Ode(%t%, %x%, %f%)
00114 %$$
00115 
00116 $subhead t$$
00117 The argument $icode t$$ to $icode%F%.Ode%$$ has prototype
00118 $codei%
00119      const %Scalar% &%t%
00120 %$$
00121 (see description of $cref/Scalar/Runge45/Scalar/$$ below). 
00122 
00123 $subhead x$$
00124 The argument $icode x$$ to $icode%F%.Ode%$$ has prototype
00125 $codei%
00126      const %Vector% &%x%
00127 %$$
00128 and has size $icode n$$
00129 (see description of $cref/Vector/Runge45/Vector/$$ below). 
00130 
00131 $subhead f$$
00132 The argument $icode f$$ to $icode%F%.Ode%$$ has prototype
00133 $codei%
00134      %Vector% &%f%
00135 %$$
00136 On input and output, $icode f$$ is a vector of size $icode n$$
00137 and the input values of the elements of $icode f$$ do not matter.
00138 On output,
00139 $icode f$$ is set equal to $latex F(t, x)$$ in the differential equation.
00140 If any of the elements of $icode f$$ have the value not a number $code nan$$
00141 the routine $code Runge45$$ returns with all the
00142 elements of $icode xf$$ and $icode e$$ equal to $code nan$$.
00143 
00144 $subhead Warning$$
00145 The argument $icode f$$ to $icode%F%.Ode%$$
00146 must have a call by reference in its prototype; i.e.,
00147 do not forget the $code &$$ in the prototype for $icode f$$.
00148 
00149 $head M$$
00150 The argument $icode M$$ has prototype
00151 $codei%
00152      size_t %M%
00153 %$$
00154 It specifies the number of steps
00155 to use when solving the differential equation.
00156 This must be greater than or equal one.
00157 The step size is given by $latex h = (tf - ti) / M$$, thus
00158 the larger $icode M$$, the more accurate the
00159 return value $icode xf$$ is as an approximation
00160 for $latex X(tf)$$.
00161 
00162 $head ti$$
00163 The argument $icode ti$$ has prototype
00164 $codei%
00165      const %Scalar% &%ti%
00166 %$$
00167 (see description of $cref/Scalar/Runge45/Scalar/$$ below). 
00168 It specifies the initial time for $icode t$$ in the 
00169 differential equation; i.e., 
00170 the time corresponding to the value $icode xi$$.
00171 
00172 $head tf$$
00173 The argument $icode tf$$ has prototype
00174 $codei%
00175      const %Scalar% &%tf%
00176 %$$
00177 It specifies the final time for $icode t$$ in the 
00178 differential equation; i.e., 
00179 the time corresponding to the value $icode xf$$.
00180 
00181 $head xi$$
00182 The argument $icode xi$$ has the prototype
00183 $codei%
00184      const %Vector% &%xi%
00185 %$$
00186 and the size of $icode xi$$ is equal to $icode n$$.
00187 It specifies the value of $latex X(ti)$$
00188 
00189 $head e$$
00190 The argument $icode e$$ is optional and has the prototype
00191 $codei%
00192      %Vector% &%e%
00193 %$$
00194 If $icode e$$ is present,
00195 the size of $icode e$$ must be equal to $icode n$$.
00196 The input value of the elements of $icode e$$ does not matter.
00197 On output
00198 it contains an element by element
00199 estimated bound for the absolute value of the error in $icode xf$$
00200 $latex \[
00201      e = O( h^5 )
00202 \] $$
00203 where $latex h = (tf - ti) / M$$ is the step size.
00204 If on output, $icode e$$ contains not a number $code nan$$,
00205 see the discussion for $cref/f/Runge45/Fun/f/$$.
00206 
00207 $head Scalar$$
00208 The type $icode Scalar$$ must satisfy the conditions
00209 for a $cref NumericType$$ type.
00210 The routine $cref CheckNumericType$$ will generate an error message
00211 if this is not the case.
00212 
00213 $subhead fabs$$
00214 In addition, the following function must be defined for 
00215 $icode Scalar$$ objects $icode a$$ and $icode b$$
00216 $codei%
00217      %a% = fabs(%b%)
00218 %$$
00219 Note that this operation is only used for computing $icode e$$; hence
00220 the operation sequence for $icode xf$$ can still be independent of
00221 the arguments to $code Runge45$$ even if 
00222 $codei%
00223      fabs(%b%) = std::max(-%b%, %b%)
00224 %$$. 
00225 
00226 $head Vector$$
00227 The type $icode Vector$$ must be a $cref SimpleVector$$ class with
00228 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$.
00229 The routine $cref CheckSimpleVector$$ will generate an error message
00230 if this is not the case.
00231 
00232 $head Parallel Mode$$
00233 $index parallel, Runge45$$
00234 $index Runge45, parallel$$
00235 For each set of types
00236 $cref/Scalar/Runge45/Scalar/$$, 
00237 $cref/Vector/Runge45/Vector/$$, and
00238 $cref/Fun/Runge45/Fun/$$, 
00239 the first call to $code Runge45$$
00240 must not be $cref/parallel/ta_in_parallel/$$ execution mode.
00241 
00242 
00243 $head Example$$
00244 $children%
00245      example/runge45_1.cpp%
00246      example/runge45_2.cpp
00247 %$$
00248 The file
00249 $cref runge45_1.cpp$$
00250 contains a simple example and test of $code Runge45$$.
00251 It returns true if it succeeds and false otherwise.
00252 $pre
00253 
00254 $$
00255 The file
00256 $cref runge45_2.cpp$$ contains an example using $code Runge45$$
00257 in the context of algorithmic differentiation.
00258 It also returns true if it succeeds and false otherwise.
00259 
00260 $head Source Code$$
00261 The source code for this routine is in the file
00262 $code cppad/runge_45.hpp$$.
00263 
00264 $end
00265 --------------------------------------------------------------------------
00266 */
00267 # include <cstddef>
00268 # include <cppad/local/cppad_assert.hpp>
00269 # include <cppad/check_simple_vector.hpp>
00270 # include <cppad/check_numeric_type.hpp>
00271 # include <cppad/nan.hpp>
00272 
00273 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
00274 # include <cppad/thread_alloc.hpp>
00275 
00276 namespace CppAD { // BEGIN CppAD namespace
00277 
00278 template <typename Scalar, typename Vector, typename Fun>
00279 Vector Runge45(
00280      Fun           &F , 
00281      size_t         M , 
00282      const Scalar &ti , 
00283      const Scalar &tf , 
00284      const Vector &xi )
00285 {    Vector e( xi.size() );
00286      return Runge45(F, M, ti, tf, xi, e);
00287 }
00288 
00289 template <typename Scalar, typename Vector, typename Fun>
00290 Vector Runge45(
00291      Fun           &F , 
00292      size_t         M , 
00293      const Scalar &ti , 
00294      const Scalar &tf , 
00295      const Vector &xi ,
00296      Vector       &e )
00297 {
00298      CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
00299 
00300      // check numeric type specifications
00301      CheckNumericType<Scalar>();
00302 
00303      // check simple vector class specifications
00304      CheckSimpleVector<Scalar, Vector>();
00305 
00306      // Cash-Karp parameters for embedded Runge-Kutta method
00307      // are static to avoid recalculation on each call and 
00308      // do not use Vector to avoid possible memory leak
00309      static Scalar a[6] = {
00310           Scalar(0),
00311           Scalar(1) / Scalar(5),
00312           Scalar(3) / Scalar(10),
00313           Scalar(3) / Scalar(5),
00314           Scalar(1),
00315           Scalar(7) / Scalar(8)
00316      };
00317      static Scalar b[5 * 5] = {
00318           Scalar(1) / Scalar(5),
00319           Scalar(0),
00320           Scalar(0),
00321           Scalar(0),
00322           Scalar(0),
00323           
00324           Scalar(3) / Scalar(40),
00325           Scalar(9) / Scalar(40),
00326           Scalar(0),
00327           Scalar(0),
00328           Scalar(0),
00329           
00330           Scalar(3) / Scalar(10),
00331           -Scalar(9) / Scalar(10),
00332           Scalar(6) / Scalar(5),
00333           Scalar(0),
00334           Scalar(0),
00335           
00336           -Scalar(11) / Scalar(54),
00337           Scalar(5) / Scalar(2),
00338           -Scalar(70) / Scalar(27),
00339           Scalar(35) / Scalar(27),
00340           Scalar(0),
00341           
00342           Scalar(1631) / Scalar(55296),
00343           Scalar(175) / Scalar(512),
00344           Scalar(575) / Scalar(13824),
00345           Scalar(44275) / Scalar(110592),
00346           Scalar(253) / Scalar(4096)
00347      };
00348      static Scalar c4[6] = {
00349           Scalar(2825) / Scalar(27648),
00350           Scalar(0),
00351           Scalar(18575) / Scalar(48384),
00352           Scalar(13525) / Scalar(55296),
00353           Scalar(277) / Scalar(14336),
00354           Scalar(1) / Scalar(4),
00355      };
00356      static Scalar c5[6] = {
00357           Scalar(37) / Scalar(378),
00358           Scalar(0),
00359           Scalar(250) / Scalar(621),
00360           Scalar(125) / Scalar(594),
00361           Scalar(0),
00362           Scalar(512) / Scalar(1771)
00363      };
00364 
00365      CPPAD_ASSERT_KNOWN(
00366           M >= 1,
00367           "Error in Runge45: the number of steps is less than one"
00368      );
00369      CPPAD_ASSERT_KNOWN(
00370           e.size() == xi.size(),
00371           "Error in Runge45: size of e not equal to size of xi"
00372      );
00373      size_t i, j, k, m;              // indices
00374 
00375      size_t  n = xi.size();           // number of components in X(t)
00376      Scalar  ns = Scalar(int(M));     // number of steps as Scalar object
00377      Scalar  h = (tf - ti) / ns;      // step size 
00378      Scalar  zero_or_nan = Scalar(0); // zero (nan if Ode returns has a nan)
00379      for(i = 0; i < n; i++)           // initialize e = 0
00380           e[i] = zero_or_nan;
00381 
00382      // vectors used to store values returned by F
00383      Vector fh(6 * n), xtmp(n), ftmp(n), x4(n), x5(n), xf(n);
00384 
00385      xf = xi;           // initialize solution
00386      for(m = 0; m < M; m++)
00387      {    // time at beginning of this interval
00388           // (convert to int to avoid MS compiler warning)
00389           Scalar t = ti * (Scalar(int(M - m)) / ns) 
00390                    + tf * (Scalar(int(m)) / ns);
00391 
00392           // loop over integration steps
00393           x4 = x5 = xf;   // start x4 and x5 at same point for each step
00394           for(j = 0; j < 6; j++)
00395           {    // loop over function evaluations for this step
00396                xtmp = xf;  // location for next function evaluation
00397                for(k = 0; k < j; k++)
00398                {    // loop over previous function evaluations
00399                     Scalar bjk = b[ (j-1) * 5 + k ];
00400                     for(i = 0; i < n; i++)
00401                     {    // loop over elements of x
00402                          xtmp[i] += bjk * fh[i * 6 + k];
00403                     }
00404                }
00405                // ftmp = F(t + a[j] * h, xtmp)
00406                F.Ode(t + a[j] * h, xtmp, ftmp); 
00407 
00408                // if ftmp has a nan, set zero_or_nan to nan
00409                for(i = 0; i < n; i++)
00410                     zero_or_nan *= ftmp[i];
00411 
00412                for(i = 0; i < n; i++)
00413                {    // loop over elements of x
00414                     Scalar fhi    = ftmp[i] * h;
00415                     fh[i * 6 + j] = fhi;
00416                     x4[i]        += c4[j] * fhi;
00417                     x5[i]        += c5[j] * fhi;
00418                     x5[i]        += zero_or_nan;
00419                }
00420           }
00421           // accumulate error bound
00422           for(i = 0; i < n; i++)
00423           {    // cant use abs because cppad.hpp may not be included
00424                Scalar diff = x5[i] - x4[i];
00425                e[i] += fabs(diff);
00426                e[i] += zero_or_nan;
00427           }
00428 
00429           // advance xf for this step using x5
00430           xf = x5;
00431      }
00432      return xf;
00433 }
00434 
00435 } // End CppAD namespace 
00436 
00437 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines