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