CppAD: A C++ Algorithmic Differentiation Package 20110419
rosen_34.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_ROSEN_34_INCLUDED
00003 # define CPPAD_ROSEN_34_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 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 Rosen34$$
00018 $spell
00019         cppad.hpp
00020         bool
00021         xf
00022         templated
00023         const
00024         Rosenbrock
00025         CppAD
00026         xi
00027         ti
00028         tf
00029         Karp
00030         Rosen
00031         Shampine
00032         ind
00033         dep
00034 $$
00035 
00036 $index Rosen34$$
00037 $index ODE, Rosenbrock$$
00038 $index Rosenbrock, ODE$$
00039 $index solve, ODE$$
00040 $index stiff, ODE$$
00041 $index differential, equation$$
00042 $index equation, differential$$
00043  
00044 $section A 3rd and 4th Order Rosenbrock ODE Solver$$
00045 
00046 $head Syntax$$
00047 $syntax%# include <cppad/rosen_34.hpp>
00048 %$$
00049 $syntax%%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%)
00050 %$$
00051 $syntax%%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%, %e%)
00052 %$$
00053 
00054 
00055 $head Description$$
00056 This is an embedded 3rd and 4th order Rosenbrock ODE solver 
00057 (see Section 16.6 of $xref/Bib/Numerical Recipes/Numerical Recipes/$$
00058 for a description of Rosenbrock ODE solvers).
00059 In particular, we use the formulas taken from page 100 of
00060 $xref/Bib/Shampine, L.F./Shampine, L.F./$$
00061 (except that the fraction 98/108 has been correction to be 97/108).
00062 $pre
00063 
00064 $$
00065 We use $latex n$$ for the size of the vector $italic xi$$.
00066 Let $latex \R$$ denote the real numbers
00067 and let $latex F : \R \times \R^n \rightarrow \R^n$$ be a smooth function.
00068 The return value $italic xf$$ contains a 5th order
00069 approximation for the value $latex X(tf)$$ where 
00070 $latex X : [ti , tf] \rightarrow \R^n$$ is defined by 
00071 the following initial value problem:
00072 $latex \[
00073 \begin{array}{rcl}
00074         X(ti)  & = & xi    \\
00075         X'(t)  & = & F[t , X(t)] 
00076 \end{array}
00077 \] $$
00078 If your set of  ordinary differential equations are not stiff
00079 an explicit method may be better (perhaps $xref/Runge45/$$.)
00080 
00081 $head Include$$
00082 The file $code cppad/rosen_34.hpp$$ is included by $code cppad/cppad.hpp$$
00083 but it can also be included separately with out the rest of 
00084 the $code CppAD$$ routines.
00085 
00086 $head xf$$
00087 The return value $italic xf$$ has the prototype
00088 $syntax%
00089         %Vector% %xf%
00090 %$$
00091 and the size of $italic xf$$ is equal to $italic n$$ 
00092 (see description of $xref/Rosen34/Vector/Vector/$$ below).
00093 $latex \[
00094         X(tf) = xf + O( h^5 )
00095 \] $$
00096 where $latex h = (tf - ti) / M$$ is the step size.
00097 If $italic xf$$ contains not a number $cref/nan/$$,
00098 see the discussion of $cref/f/Rosen34/Fun/Nan/$$.
00099 
00100 $head Fun$$
00101 The class $italic Fun$$ 
00102 and the object $italic F$$ satisfy the prototype
00103 $syntax%
00104         %Fun% &%F%
00105 %$$
00106 This must support the following set of calls
00107 $syntax%
00108         %F%.Ode(%t%, %x%, %f%)
00109         %F%.Ode_ind(%t%, %x%, %f_t%)
00110         %F%.Ode_dep(%t%, %x%, %f_x%)
00111 %$$
00112 
00113 $subhead t$$
00114 In all three cases, 
00115 the argument $italic t$$ has prototype
00116 $syntax%
00117         const %Scalar% &%t%
00118 %$$
00119 (see description of $xref/Rosen34/Scalar/Scalar/$$ below). 
00120 
00121 $subhead x$$
00122 In all three cases,
00123 the argument $italic x$$ has prototype
00124 $syntax%
00125         const %Vector% &%x%
00126 %$$
00127 and has size $italic n$$
00128 (see description of $xref/Rosen34/Vector/Vector/$$ below). 
00129 
00130 $subhead f$$
00131 The argument $italic f$$ to $syntax%%F%.Ode%$$ has prototype
00132 $syntax%
00133         %Vector% &%f%
00134 %$$
00135 On input and output, $italic f$$ is a vector of size $italic n$$
00136 and the input values of the elements of $italic f$$ do not matter.
00137 On output,
00138 $italic f$$ is set equal to $latex F(t, x)$$
00139 (see $italic F(t, x)$$ in $xref/Rosen34/Description/Description/$$). 
00140 
00141 $subhead f_t$$
00142 The argument $italic f_t$$ to $syntax%%F%.Ode_ind%$$ has prototype
00143 $syntax%
00144         %Vector% &%f_t%
00145 %$$
00146 On input and output, $italic f_t$$ is a vector of size $italic n$$
00147 and the input values of the elements of $italic f_t$$ do not matter.
00148 On output, the $th i$$ element of
00149 $italic f_t$$ is set equal to $latex \partial_t F_i (t, x)$$ 
00150 (see $italic F(t, x)$$ in $xref/Rosen34/Description/Description/$$). 
00151 
00152 $subhead f_x$$
00153 The argument $italic f_x$$ to $syntax%%F%.Ode_dep%$$ has prototype
00154 $syntax%
00155         %Vector% &%f_x%
00156 %$$
00157 On input and output, $italic f_x$$ is a vector of size $syntax%%n%*%n%$$
00158 and the input values of the elements of $italic f_x$$ do not matter.
00159 On output, the [$syntax%%i%*%n%+%j%$$] element of
00160 $italic f_x$$ is set equal to $latex \partial_{x(j)} F_i (t, x)$$ 
00161 (see $italic F(t, x)$$ in $xref/Rosen34/Description/Description/$$). 
00162 
00163 $subhead Nan$$
00164 If any of the elements of $italic f$$, $italic f_t$$, or $italic f_x$$
00165 have the value not a number $code nan$$,
00166 the routine $code Rosen34$$ returns with all the
00167 elements of $italic xf$$ and $italic e$$ equal to $code nan$$.
00168 
00169 $subhead Warning$$
00170 The arguments $italic f$$, $italic f_t$$, and $italic f_x$$
00171 must have a call by reference in their prototypes; i.e.,
00172 do not forget the $code &$$ in the prototype for 
00173 $italic f$$, $italic f_t$$ and $italic f_x$$.
00174 
00175 $subhead Optimization$$
00176 Every call of the form 
00177 $syntax%
00178         %F%.Ode_ind(%t%, %x%, %f_t%)
00179 %$$
00180 is directly followed by a call of the form 
00181 $syntax%
00182         %F%.Ode_dep(%t%, %x%, %f_x%)
00183 %$$
00184 where the arguments $italic t$$ and $italic x$$ have not changed between calls.
00185 In many cases it is faster to compute the values of $italic f_t$$
00186 and $italic f_x$$ together and then pass them back one at a time.
00187 
00188 $head M$$
00189 The argument $italic M$$ has prototype
00190 $syntax%
00191         size_t %M%
00192 %$$
00193 It specifies the number of steps
00194 to use when solving the differential equation.
00195 This must be greater than or equal one.
00196 The step size is given by $latex h = (tf - ti) / M$$, thus
00197 the larger $italic M$$, the more accurate the
00198 return value $italic xf$$ is as an approximation
00199 for $latex X(tf)$$.
00200 
00201 $head ti$$
00202 The argument $italic ti$$ has prototype
00203 $syntax%
00204         const %Scalar% &%ti%
00205 %$$
00206 (see description of $xref/Rosen34/Scalar/Scalar/$$ below). 
00207 It specifies the initial time for $italic t$$ in the 
00208 differential equation; i.e., 
00209 the time corresponding to the value $italic xi$$.
00210 
00211 $head tf$$
00212 The argument $italic tf$$ has prototype
00213 $syntax%
00214         const %Scalar% &%tf%
00215 %$$
00216 It specifies the final time for $italic t$$ in the 
00217 differential equation; i.e., 
00218 the time corresponding to the value $italic xf$$.
00219 
00220 $head xi$$
00221 The argument $italic xi$$ has the prototype
00222 $syntax%
00223         const %Vector% &%xi%
00224 %$$
00225 and the size of $italic xi$$ is equal to $italic n$$.
00226 It specifies the value of $latex X(ti)$$
00227 
00228 $head e$$
00229 The argument $italic e$$ is optional and has the prototype
00230 $syntax%
00231         %Vector% &%e%
00232 %$$
00233 If $italic e$$ is present,
00234 the size of $italic e$$ must be equal to $italic n$$.
00235 The input value of the elements of $italic e$$ does not matter.
00236 On output
00237 it contains an element by element
00238 estimated bound for the absolute value of the error in $italic xf$$
00239 $latex \[
00240         e = O( h^4 )
00241 \] $$
00242 where $latex h = (tf - ti) / M$$ is the step size.
00243 
00244 $head Scalar$$
00245 The type $italic Scalar$$ must satisfy the conditions
00246 for a $xref/NumericType/$$ type.
00247 The routine $xref/CheckNumericType/$$ will generate an error message
00248 if this is not the case.
00249 In addition, the following operations must be defined for 
00250 $italic Scalar$$ objects $italic a$$ and $italic b$$:
00251 
00252 $table
00253 $bold Operation$$ $cnext $bold Description$$  $rnext
00254 $syntax%%a% < %b%$$ $cnext
00255         less than operator (returns a $code bool$$ object)
00256 $tend
00257 
00258 $head Vector$$
00259 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with
00260 $xref/SimpleVector/Elements of Specified Type/elements of type Scalar/$$.
00261 The routine $xref/CheckSimpleVector/$$ will generate an error message
00262 if this is not the case.
00263 
00264 $head Example$$
00265 $children%
00266         example/rosen_34.cpp
00267 %$$
00268 The file
00269 $xref/Rosen34.cpp/$$
00270 contains an example and test a test of using this routine.
00271 It returns true if it succeeds and false otherwise.
00272 
00273 $head Source Code$$
00274 The source code for this routine is in the file
00275 $code cppad/rosen_34.hpp$$.
00276 
00277 $end
00278 --------------------------------------------------------------------------
00279 */
00280 
00281 # include <cstddef>
00282 # include <cppad/local/cppad_assert.hpp>
00283 # include <cppad/check_simple_vector.hpp>
00284 # include <cppad/check_numeric_type.hpp>
00285 # include <cppad/vector.hpp>
00286 # include <cppad/lu_factor.hpp>
00287 # include <cppad/lu_invert.hpp>
00288 
00289 namespace CppAD { // BEGIN CppAD namespace
00290 
00291 template <typename Scalar, typename Vector, typename Fun>
00292 Vector Rosen34(
00293         Fun           &F , 
00294         size_t         M , 
00295         const Scalar &ti , 
00296         const Scalar &tf , 
00297         const Vector &xi )
00298 {       Vector e( xi.size() );
00299         return Rosen34(F, M, ti, tf, xi, e);
00300 }
00301 
00302 template <typename Scalar, typename Vector, typename Fun>
00303 Vector Rosen34(
00304         Fun           &F , 
00305         size_t         M , 
00306         const Scalar &ti , 
00307         const Scalar &tf , 
00308         const Vector &xi ,
00309         Vector       &e )
00310 {
00311         // check numeric type specifications
00312         CheckNumericType<Scalar>();
00313 
00314         // check simple vector class specifications
00315         CheckSimpleVector<Scalar, Vector>();
00316 
00317         // Parameters for Shampine's Rosenbrock method
00318         // are static to avoid recalculation on each call and 
00319         // do not use Vector to avoid possible memory leak
00320         static Scalar a[3] = {
00321                 Scalar(0),
00322                 Scalar(1),
00323                 Scalar(3)   / Scalar(5)
00324         };
00325         static Scalar b[2 * 2] = {
00326                 Scalar(1),
00327                 Scalar(0),
00328                 Scalar(24)  / Scalar(25),
00329                 Scalar(3)   / Scalar(25)
00330         };
00331         static Scalar ct[4] = {
00332                 Scalar(1)   / Scalar(2),
00333                 - Scalar(3) / Scalar(2),
00334                 Scalar(121) / Scalar(50),
00335                 Scalar(29)  / Scalar(250)
00336         };
00337         static Scalar cg[3 * 3] = {
00338                 - Scalar(4),
00339                 Scalar(0),
00340                 Scalar(0),
00341                 Scalar(186) / Scalar(25),
00342                 Scalar(6)   / Scalar(5),
00343                 Scalar(0),
00344                 - Scalar(56) / Scalar(125),
00345                 - Scalar(27) / Scalar(125),
00346                 - Scalar(1)  / Scalar(5)
00347         };
00348         static Scalar d3[3] = {
00349                 Scalar(97) / Scalar(108),
00350                 Scalar(11) / Scalar(72),
00351                 Scalar(25) / Scalar(216)
00352         };
00353         static Scalar d4[4] = {
00354                 Scalar(19)  / Scalar(18),
00355                 Scalar(1)   / Scalar(4),
00356                 Scalar(25)  / Scalar(216),
00357                 Scalar(125) / Scalar(216)
00358         };
00359         CPPAD_ASSERT_KNOWN(
00360                 M >= 1,
00361                 "Error in Rosen34: the number of steps is less than one"
00362         );
00363         CPPAD_ASSERT_KNOWN(
00364                 e.size() == xi.size(),
00365                 "Error in Rosen34: size of e not equal to size of xi"
00366         );
00367         size_t i, j, k, l, m;             // indices
00368 
00369         size_t  n    = xi.size();         // number of components in X(t)
00370         Scalar  ns   = Scalar(double(M)); // number of steps as Scalar object
00371         Scalar  h    = (tf - ti) / ns;    // step size 
00372         Scalar  zero = Scalar(0);         // some constants
00373         Scalar  one  = Scalar(1);
00374         Scalar  two  = Scalar(2);
00375 
00376         // permutation vectors needed for LU factorization routine
00377         CppAD::vector<size_t> ip(n), jp(n);
00378 
00379         // vectors used to store values returned by F
00380         Vector E(n * n), Eg(n), f_t(n);
00381         Vector g(n * 3), x3(n), x4(n), xf(n), ftmp(n), xtmp(n), nan_vec(n);
00382 
00383         // initialize e = 0, nan_vec = nan
00384         for(i = 0; i < n; i++)
00385         {       e[i]       = zero;
00386                 nan_vec[i] = nan(zero);
00387         }
00388 
00389         xf = xi;           // initialize solution
00390         for(m = 0; m < M; m++)
00391         {       // time at beginning of this interval
00392                 Scalar t = ti * (Scalar(int(M - m)) / ns) 
00393                          + tf * (Scalar(int(m)) / ns);
00394 
00395                 // value of x at beginning of this interval
00396                 x3 = x4 = xf;
00397 
00398                 // evaluate partial derivatives at beginning of this interval
00399                 F.Ode_ind(t, xf, f_t);
00400                 F.Ode_dep(t, xf, E);    // E = f_x
00401                 if( hasnan(f_t) || hasnan(E) )
00402                 {       e = nan_vec;
00403                         return nan_vec;
00404                 }
00405 
00406                 // E = I - f_x * h / 2
00407                 for(i = 0; i < n; i++)
00408                 {       for(j = 0; j < n; j++)
00409                                 E[i * n + j] = - E[i * n + j] * h / two;
00410                         E[i * n + i] += one;
00411                 }
00412 
00413                 // LU factor the matrix E
00414                 int sign = LuFactor(ip, jp, E);
00415                 CPPAD_ASSERT_KNOWN(
00416                         sign != 0,
00417                         "Error in Rosen34: I - f_x * h / 2 not invertible"
00418                 );
00419 
00420                 // loop over integration steps
00421                 for(k = 0; k < 3; k++)
00422                 {       // set location for next function evaluation
00423                         xtmp = xf; 
00424                         for(l = 0; l < k; l++)
00425                         {       // loop over previous function evaluations
00426                                 Scalar bkl = b[(k-1)*2 + l];
00427                                 for(i = 0; i < n; i++)
00428                                 {       // loop over elements of x
00429                                         xtmp[i] += bkl * g[i*3 + l] * h;
00430                                 }
00431                         }
00432                         // ftmp = F(t + a[k] * h, xtmp)
00433                         F.Ode(t + a[k] * h, xtmp, ftmp); 
00434                         if( hasnan(ftmp) )
00435                         {       e = nan_vec;
00436                                 return nan_vec;
00437                         }
00438 
00439                         // Form Eg for this integration step
00440                         for(i = 0; i < n; i++)
00441                                 Eg[i] = ftmp[i] + ct[k] * f_t[i] * h;
00442                         for(l = 0; l < k; l++)
00443                         {       for(i = 0; i < n; i++)
00444                                         Eg[i] += cg[(k-1)*3 + l] * g[i*3 + l];
00445                         }
00446 
00447                         // Solve the equation E * g = Eg
00448                         LuInvert(ip, jp, E, Eg);
00449 
00450                         // save solution and advance x3, x4
00451                         for(i = 0; i < n; i++)
00452                         {       g[i*3 + k]  = Eg[i];
00453                                 x3[i]      += h * d3[k] * Eg[i];
00454                                 x4[i]      += h * d4[k] * Eg[i];
00455                         }
00456                 }
00457                 // Form Eg for last update to x4 only
00458                 for(i = 0; i < n; i++)
00459                         Eg[i] = ftmp[i] + ct[3] * f_t[i] * h;
00460                 for(l = 0; l < 3; l++)
00461                 {       for(i = 0; i < n; i++)
00462                                 Eg[i] += cg[2*3 + l] * g[i*3 + l];
00463                 }
00464 
00465                 // Solve the equation E * g = Eg
00466                 LuInvert(ip, jp, E, Eg);
00467 
00468                 // advance x4 and accumulate error bound
00469                 for(i = 0; i < n; i++)
00470                 {       x4[i] += h * d4[3] * Eg[i];
00471 
00472                         // cant use abs because cppad.hpp may not be included
00473                         Scalar diff = x4[i] - x3[i];
00474                         if( diff < zero )
00475                                 e[i] -= diff;
00476                         else    e[i] += diff;
00477                 }
00478 
00479                 // advance xf for this step using x4
00480                 xf = x4;
00481         }
00482         return xf;
00483 }
00484 
00485 } // End CppAD namespace 
00486 
00487 # endif