CppAD: A C++ Algorithmic Differentiation Package  20130918
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-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 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 $codei%# include <cppad/rosen_34.hpp>
00048 %$$
00049 $icode%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%)
00050 %$$
00051 $icode%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 $cref/Numerical Recipes/Bib/Numerical Recipes/$$
00058 for a description of Rosenbrock ODE solvers).
00059 In particular, we use the formulas taken from page 100 of
00060 $cref/Shampine, L.F./Bib/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 $icode xi$$.
00066 Let $latex \B{R}$$ denote the real numbers
00067 and let $latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ be a smooth function.
00068 The return value $icode xf$$ contains a 5th order
00069 approximation for the value $latex X(tf)$$ where 
00070 $latex X : [ti , tf] \rightarrow \B{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 $cref 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 $icode xf$$ has the prototype
00088 $codei%
00089      %Vector% %xf%
00090 %$$
00091 and the size of $icode xf$$ is equal to $icode n$$ 
00092 (see description of $cref/Vector/Rosen34/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 $icode xf$$ contains not a number $cref nan$$,
00098 see the discussion of $cref/f/Rosen34/Fun/Nan/$$.
00099 
00100 $head Fun$$
00101 The class $icode Fun$$ 
00102 and the object $icode F$$ satisfy the prototype
00103 $codei%
00104      %Fun% &%F%
00105 %$$
00106 This must support the following set of calls
00107 $codei%
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 $icode t$$ has prototype
00116 $codei%
00117      const %Scalar% &%t%
00118 %$$
00119 (see description of $cref/Scalar/Rosen34/Scalar/$$ below). 
00120 
00121 $subhead x$$
00122 In all three cases,
00123 the argument $icode x$$ has prototype
00124 $codei%
00125      const %Vector% &%x%
00126 %$$
00127 and has size $icode n$$
00128 (see description of $cref/Vector/Rosen34/Vector/$$ below). 
00129 
00130 $subhead f$$
00131 The argument $icode f$$ to $icode%F%.Ode%$$ has prototype
00132 $codei%
00133      %Vector% &%f%
00134 %$$
00135 On input and output, $icode f$$ is a vector of size $icode n$$
00136 and the input values of the elements of $icode f$$ do not matter.
00137 On output,
00138 $icode f$$ is set equal to $latex F(t, x)$$
00139 (see $icode F(t, x)$$ in $cref/Description/Rosen34/Description/$$). 
00140 
00141 $subhead f_t$$
00142 The argument $icode f_t$$ to $icode%F%.Ode_ind%$$ has prototype
00143 $codei%
00144      %Vector% &%f_t%
00145 %$$
00146 On input and output, $icode f_t$$ is a vector of size $icode n$$
00147 and the input values of the elements of $icode f_t$$ do not matter.
00148 On output, the $th i$$ element of
00149 $icode f_t$$ is set equal to $latex \partial_t F_i (t, x)$$ 
00150 (see $icode F(t, x)$$ in $cref/Description/Rosen34/Description/$$). 
00151 
00152 $subhead f_x$$
00153 The argument $icode f_x$$ to $icode%F%.Ode_dep%$$ has prototype
00154 $codei%
00155      %Vector% &%f_x%
00156 %$$
00157 On input and output, $icode f_x$$ is a vector of size $icode%n%*%n%$$
00158 and the input values of the elements of $icode f_x$$ do not matter.
00159 On output, the [$icode%i%*%n%+%j%$$] element of
00160 $icode f_x$$ is set equal to $latex \partial_{x(j)} F_i (t, x)$$ 
00161 (see $icode F(t, x)$$ in $cref/Description/Rosen34/Description/$$). 
00162 
00163 $subhead Nan$$
00164 If any of the elements of $icode f$$, $icode f_t$$, or $icode f_x$$
00165 have the value not a number $code nan$$,
00166 the routine $code Rosen34$$ returns with all the
00167 elements of $icode xf$$ and $icode e$$ equal to $code nan$$.
00168 
00169 $subhead Warning$$
00170 The arguments $icode f$$, $icode f_t$$, and $icode 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 $icode f$$, $icode f_t$$ and $icode f_x$$.
00174 
00175 $subhead Optimization$$
00176 Every call of the form 
00177 $codei%
00178      %F%.Ode_ind(%t%, %x%, %f_t%)
00179 %$$
00180 is directly followed by a call of the form 
00181 $codei%
00182      %F%.Ode_dep(%t%, %x%, %f_x%)
00183 %$$
00184 where the arguments $icode t$$ and $icode x$$ have not changed between calls.
00185 In many cases it is faster to compute the values of $icode f_t$$
00186 and $icode f_x$$ together and then pass them back one at a time.
00187 
00188 $head M$$
00189 The argument $icode M$$ has prototype
00190 $codei%
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 $icode M$$, the more accurate the
00198 return value $icode xf$$ is as an approximation
00199 for $latex X(tf)$$.
00200 
00201 $head ti$$
00202 The argument $icode ti$$ has prototype
00203 $codei%
00204      const %Scalar% &%ti%
00205 %$$
00206 (see description of $cref/Scalar/Rosen34/Scalar/$$ below). 
00207 It specifies the initial time for $icode t$$ in the 
00208 differential equation; i.e., 
00209 the time corresponding to the value $icode xi$$.
00210 
00211 $head tf$$
00212 The argument $icode tf$$ has prototype
00213 $codei%
00214      const %Scalar% &%tf%
00215 %$$
00216 It specifies the final time for $icode t$$ in the 
00217 differential equation; i.e., 
00218 the time corresponding to the value $icode xf$$.
00219 
00220 $head xi$$
00221 The argument $icode xi$$ has the prototype
00222 $codei%
00223      const %Vector% &%xi%
00224 %$$
00225 and the size of $icode xi$$ is equal to $icode n$$.
00226 It specifies the value of $latex X(ti)$$
00227 
00228 $head e$$
00229 The argument $icode e$$ is optional and has the prototype
00230 $codei%
00231      %Vector% &%e%
00232 %$$
00233 If $icode e$$ is present,
00234 the size of $icode e$$ must be equal to $icode n$$.
00235 The input value of the elements of $icode 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 $icode 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 $icode Scalar$$ must satisfy the conditions
00246 for a $cref NumericType$$ type.
00247 The routine $cref CheckNumericType$$ will generate an error message
00248 if this is not the case.
00249 In addition, the following operations must be defined for 
00250 $icode Scalar$$ objects $icode a$$ and $icode b$$:
00251 
00252 $table
00253 $bold Operation$$ $cnext $bold Description$$  $rnext
00254 $icode%a% < %b%$$ $cnext
00255      less than operator (returns a $code bool$$ object)
00256 $tend
00257 
00258 $head Vector$$
00259 The type $icode Vector$$ must be a $cref SimpleVector$$ class with
00260 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$.
00261 The routine $cref CheckSimpleVector$$ will generate an error message
00262 if this is not the case.
00263 
00264 $head Parallel Mode$$
00265 $index parallel, Rosen34$$
00266 $index Rosen34, parallel$$
00267 For each set of types
00268 $cref/Scalar/Rosen34/Scalar/$$, 
00269 $cref/Vector/Rosen34/Vector/$$, and
00270 $cref/Fun/Rosen34/Fun/$$, 
00271 the first call to $code Rosen34$$
00272 must not be $cref/parallel/ta_in_parallel/$$ execution mode.
00273 
00274 $head Example$$
00275 $children%
00276      example/rosen_34.cpp
00277 %$$
00278 The file
00279 $cref rosen_34.cpp$$
00280 contains an example and test a test of using this routine.
00281 It returns true if it succeeds and false otherwise.
00282 
00283 $head Source Code$$
00284 The source code for this routine is in the file
00285 $code cppad/rosen_34.hpp$$.
00286 
00287 $end
00288 --------------------------------------------------------------------------
00289 */
00290 
00291 # include <cstddef>
00292 # include <cppad/local/cppad_assert.hpp>
00293 # include <cppad/check_simple_vector.hpp>
00294 # include <cppad/check_numeric_type.hpp>
00295 # include <cppad/vector.hpp>
00296 # include <cppad/lu_factor.hpp>
00297 # include <cppad/lu_invert.hpp>
00298 
00299 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
00300 # include <cppad/thread_alloc.hpp>
00301 
00302 namespace CppAD { // BEGIN CppAD namespace
00303 
00304 template <typename Scalar, typename Vector, typename Fun>
00305 Vector Rosen34(
00306      Fun           &F , 
00307      size_t         M , 
00308      const Scalar &ti , 
00309      const Scalar &tf , 
00310      const Vector &xi )
00311 {    Vector e( xi.size() );
00312      return Rosen34(F, M, ti, tf, xi, e);
00313 }
00314 
00315 template <typename Scalar, typename Vector, typename Fun>
00316 Vector Rosen34(
00317      Fun           &F , 
00318      size_t         M , 
00319      const Scalar &ti , 
00320      const Scalar &tf , 
00321      const Vector &xi ,
00322      Vector       &e )
00323 {
00324      CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
00325 
00326      // check numeric type specifications
00327      CheckNumericType<Scalar>();
00328 
00329      // check simple vector class specifications
00330      CheckSimpleVector<Scalar, Vector>();
00331 
00332      // Parameters for Shampine's Rosenbrock method
00333      // are static to avoid recalculation on each call and 
00334      // do not use Vector to avoid possible memory leak
00335      static Scalar a[3] = {
00336           Scalar(0),
00337           Scalar(1),
00338           Scalar(3)   / Scalar(5)
00339      };
00340      static Scalar b[2 * 2] = {
00341           Scalar(1),
00342           Scalar(0),
00343           Scalar(24)  / Scalar(25),
00344           Scalar(3)   / Scalar(25)
00345      };
00346      static Scalar ct[4] = {
00347           Scalar(1)   / Scalar(2),
00348           - Scalar(3) / Scalar(2),
00349           Scalar(121) / Scalar(50),
00350           Scalar(29)  / Scalar(250)
00351      };
00352      static Scalar cg[3 * 3] = {
00353           - Scalar(4),
00354           Scalar(0),
00355           Scalar(0),
00356           Scalar(186) / Scalar(25),
00357           Scalar(6)   / Scalar(5),
00358           Scalar(0),
00359           - Scalar(56) / Scalar(125),
00360           - Scalar(27) / Scalar(125),
00361           - Scalar(1)  / Scalar(5)
00362      };
00363      static Scalar d3[3] = {
00364           Scalar(97) / Scalar(108),
00365           Scalar(11) / Scalar(72),
00366           Scalar(25) / Scalar(216)
00367      };
00368      static Scalar d4[4] = {
00369           Scalar(19)  / Scalar(18),
00370           Scalar(1)   / Scalar(4),
00371           Scalar(25)  / Scalar(216),
00372           Scalar(125) / Scalar(216)
00373      };
00374      CPPAD_ASSERT_KNOWN(
00375           M >= 1,
00376           "Error in Rosen34: the number of steps is less than one"
00377      );
00378      CPPAD_ASSERT_KNOWN(
00379           e.size() == xi.size(),
00380           "Error in Rosen34: size of e not equal to size of xi"
00381      );
00382      size_t i, j, k, l, m;             // indices
00383 
00384      size_t  n    = xi.size();         // number of components in X(t)
00385      Scalar  ns   = Scalar(double(M)); // number of steps as Scalar object
00386      Scalar  h    = (tf - ti) / ns;    // step size 
00387      Scalar  zero = Scalar(0);         // some constants
00388      Scalar  one  = Scalar(1);
00389      Scalar  two  = Scalar(2);
00390 
00391      // permutation vectors needed for LU factorization routine
00392      CppAD::vector<size_t> ip(n), jp(n);
00393 
00394      // vectors used to store values returned by F
00395      Vector E(n * n), Eg(n), f_t(n);
00396      Vector g(n * 3), x3(n), x4(n), xf(n), ftmp(n), xtmp(n), nan_vec(n);
00397 
00398      // initialize e = 0, nan_vec = nan
00399      for(i = 0; i < n; i++)
00400      {    e[i]       = zero;
00401           nan_vec[i] = nan(zero);
00402      }
00403 
00404      xf = xi;           // initialize solution
00405      for(m = 0; m < M; m++)
00406      {    // time at beginning of this interval
00407           Scalar t = ti * (Scalar(int(M - m)) / ns) 
00408                    + tf * (Scalar(int(m)) / ns);
00409 
00410           // value of x at beginning of this interval
00411           x3 = x4 = xf;
00412 
00413           // evaluate partial derivatives at beginning of this interval
00414           F.Ode_ind(t, xf, f_t);
00415           F.Ode_dep(t, xf, E);    // E = f_x
00416           if( hasnan(f_t) || hasnan(E) )
00417           {    e = nan_vec;
00418                return nan_vec;
00419           }
00420 
00421           // E = I - f_x * h / 2
00422           for(i = 0; i < n; i++)
00423           {    for(j = 0; j < n; j++)
00424                     E[i * n + j] = - E[i * n + j] * h / two;
00425                E[i * n + i] += one;
00426           }
00427 
00428           // LU factor the matrix E
00429 # ifndef NDEBUG
00430           int sign = LuFactor(ip, jp, E);
00431 # else
00432           LuFactor(ip, jp, E);
00433 # endif
00434           CPPAD_ASSERT_KNOWN(
00435                sign != 0,
00436                "Error in Rosen34: I - f_x * h / 2 not invertible"
00437           );
00438 
00439           // loop over integration steps
00440           for(k = 0; k < 3; k++)
00441           {    // set location for next function evaluation
00442                xtmp = xf; 
00443                for(l = 0; l < k; l++)
00444                {    // loop over previous function evaluations
00445                     Scalar bkl = b[(k-1)*2 + l];
00446                     for(i = 0; i < n; i++)
00447                     {    // loop over elements of x
00448                          xtmp[i] += bkl * g[i*3 + l] * h;
00449                     }
00450                }
00451                // ftmp = F(t + a[k] * h, xtmp)
00452                F.Ode(t + a[k] * h, xtmp, ftmp); 
00453                if( hasnan(ftmp) )
00454                {    e = nan_vec;
00455                     return nan_vec;
00456                }
00457 
00458                // Form Eg for this integration step
00459                for(i = 0; i < n; i++)
00460                     Eg[i] = ftmp[i] + ct[k] * f_t[i] * h;
00461                for(l = 0; l < k; l++)
00462                {    for(i = 0; i < n; i++)
00463                          Eg[i] += cg[(k-1)*3 + l] * g[i*3 + l];
00464                }
00465 
00466                // Solve the equation E * g = Eg
00467                LuInvert(ip, jp, E, Eg);
00468 
00469                // save solution and advance x3, x4
00470                for(i = 0; i < n; i++)
00471                {    g[i*3 + k]  = Eg[i];
00472                     x3[i]      += h * d3[k] * Eg[i];
00473                     x4[i]      += h * d4[k] * Eg[i];
00474                }
00475           }
00476           // Form Eg for last update to x4 only
00477           for(i = 0; i < n; i++)
00478                Eg[i] = ftmp[i] + ct[3] * f_t[i] * h;
00479           for(l = 0; l < 3; l++)
00480           {    for(i = 0; i < n; i++)
00481                     Eg[i] += cg[2*3 + l] * g[i*3 + l];
00482           }
00483 
00484           // Solve the equation E * g = Eg
00485           LuInvert(ip, jp, E, Eg);
00486 
00487           // advance x4 and accumulate error bound
00488           for(i = 0; i < n; i++)
00489           {    x4[i] += h * d4[3] * Eg[i];
00490 
00491                // cant use abs because cppad.hpp may not be included
00492                Scalar diff = x4[i] - x3[i];
00493                if( diff < zero )
00494                     e[i] -= diff;
00495                else e[i] += diff;
00496           }
00497 
00498           // advance xf for this step using x4
00499           xf = x4;
00500      }
00501      return xf;
00502 }
00503 
00504 } // End CppAD namespace 
00505 
00506 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines