CppAD: A C++ Algorithmic Differentiation Package  20130102
ode_err_control.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_ODE_ERR_CONTROL_INCLUDED
00003 # define CPPAD_ODE_ERR_CONTROL_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 OdeErrControl$$
00018 $spell
00019      cppad.hpp
00020      nstep
00021      maxabs
00022      exp
00023      scur
00024      CppAD
00025      xf
00026      tf
00027      xi
00028      smin
00029      smax
00030      eabs
00031      erel
00032      ef
00033      ta
00034      tb
00035      xa
00036      xb
00037      const
00038      eb
00039 $$
00040 
00041 $index OdeErrControl$$
00042 $index ODE, control error$$
00043 $index control, ODE error$$
00044 $index error, control ODE$$
00045 $index differential, ODE error control$$
00046 $index equation, ODE error control$$
00047 
00048  
00049 $section An Error Controller for ODE Solvers$$
00050 
00051 $head Syntax$$
00052 $code # include <cppad/ode_err_control.hpp>$$
00053 $pre
00054 $$
00055 $icode%xf% = OdeErrControl(%method%, %ti%, %tf%, %xi%,
00056      %smin%, %smax%, %scur%, %eabs%, %erel%, %ef% , %maxabs%, %nstep% )%$$
00057 
00058 
00059 $head Description$$
00060 Let $latex \B{R}$$ denote the real numbers
00061 and let $latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ be a smooth function.
00062 We define $latex X : [ti , tf] \rightarrow \B{R}^n$$ by 
00063 the following initial value problem:
00064 $latex \[
00065 \begin{array}{rcl}
00066      X(ti)  & = & xi    \\
00067      X'(t)  & = & F[t , X(t)] 
00068 \end{array}
00069 \] $$
00070 The routine $code OdeErrControl$$ can be used to adjust the step size
00071 used an arbitrary integration methods in order to be as fast as possible
00072 and still with in a requested error bound.
00073 
00074 $head Include$$
00075 The file $code cppad/ode_err_control.hpp$$ is included by 
00076 $code cppad/cppad.hpp$$
00077 but it can also be included separately with out the rest of 
00078 the $code CppAD$$ routines.
00079 
00080 $head Notation$$
00081 The template parameter types $cref/Scalar/OdeErrControl/Scalar/$$ and
00082 $cref/Vector/OdeErrControl/Vector/$$ are documented below.
00083 
00084 $head xf$$
00085 The return value $icode xf$$ has the prototype
00086 $codei%
00087      %Vector% %xf%
00088 %$$
00089 (see description of $cref/Vector/OdeErrControl/Vector/$$ below). 
00090 and the size of $icode xf$$ is equal to $icode n$$.
00091 If $icode xf$$ contains not a number $cref nan$$,
00092 see the discussion of $cref/step/OdeErrControl/Method/Nan/$$.
00093 
00094 $head Method$$
00095 The class $icode Method$$
00096 and the object $icode method$$ satisfy the following syntax
00097 $codei%
00098      %Method% &%method%
00099 %$$
00100 The object $icode method$$ must support $code step$$ and 
00101 $code order$$ member functions defined below:
00102 
00103 $subhead step$$
00104 The syntax
00105 $codei%
00106      %method%.step(%ta%, %tb%, %xa%, %xb%, %eb%)
00107 %$$
00108 executes one step of the integration method. 
00109 $codei%
00110 
00111 %ta%
00112 %$$
00113 The argument $icode ta$$ has prototype
00114 $codei%
00115      const %Scalar% &%ta%
00116 %$$
00117 It specifies the initial time for this step in the 
00118 ODE integration.
00119 (see description of $cref/Scalar/OdeErrControl/Scalar/$$ below). 
00120 $codei%
00121 
00122 %tb%
00123 %$$
00124 The argument $icode tb$$ has prototype
00125 $codei%
00126      const %Scalar% &%tb%
00127 %$$
00128 It specifies the final time for this step in the 
00129 ODE integration.
00130 $codei%
00131 
00132 %xa%
00133 %$$
00134 The argument $icode xa$$ has prototype
00135 $codei%
00136      const %Vector% &%xa%
00137 %$$
00138 and size $icode n$$.
00139 It specifies the value of $latex X(ta)$$.
00140 (see description of $cref/Vector/OdeErrControl/Vector/$$ below). 
00141 $codei%
00142 
00143 %xb%
00144 %$$
00145 The argument value $icode xb$$ has prototype
00146 $codei%
00147      %Vector% &%xb%
00148 %$$
00149 and size $icode n$$.
00150 The input value of its elements does not matter.
00151 On output, 
00152 it contains the approximation for $latex X(tb)$$ that the method obtains.
00153 $codei%
00154 
00155 %eb%
00156 %$$
00157 The argument value $icode eb$$ has prototype
00158 $codei%
00159      %Vector% &%eb%
00160 %$$
00161 and size $icode n$$.
00162 The input value of its elements does not matter.
00163 On output, 
00164 it contains an estimate for the error in the approximation $icode xb$$.
00165 It is assumed (locally) that the error bound in this approximation 
00166 nearly equal to $latex K (tb - ta)^m$$ 
00167 where $icode K$$ is a fixed constant and $icode m$$
00168 is the corresponding argument to $code CodeControl$$.
00169 
00170 $subhead Nan$$
00171 If any element of the vector $icode eb$$ or $icode xb$$ are
00172 not a number $code nan$$,
00173 the current step is considered to large.
00174 If this happens with the current step size equal to $icode smin$$, 
00175 $code OdeErrControl$$ returns with $icode xf$$ and $icode ef$$ as vectors
00176 of $code nan$$.
00177 
00178 $subhead order$$
00179 If $icode m$$ is $code size_t$$,
00180 the object $icode method$$ must also support the following syntax 
00181 $codei%
00182      %m% = %method%.order()
00183 %$$
00184 The return value $icode m$$ is the order of the error estimate;
00185 i.e., there is a constant K such that if $latex ti \leq ta \leq tb \leq tf$$,
00186 $latex \[
00187      | eb(tb) | \leq K | tb - ta |^m
00188 \] $$
00189 where $icode ta$$, $icode tb$$, and $icode eb$$ are as in 
00190 $icode%method%.step(%ta%, %tb%, %xa%, %xb%, %eb%)%$$
00191 
00192 
00193 $head ti$$
00194 The argument $icode ti$$ has prototype
00195 $codei%
00196      const %Scalar% &%ti%
00197 %$$
00198 It specifies the initial time for the integration of 
00199 the differential equation.
00200 
00201 
00202 $head tf$$
00203 The argument $icode tf$$ has prototype
00204 $codei%
00205      const %Scalar% &%tf%
00206 %$$
00207 It specifies the final time for the integration of 
00208 the differential equation.
00209 
00210 $head xi$$
00211 The argument $icode xi$$ has prototype
00212 $codei%
00213      const %Vector% &%xi%
00214 %$$
00215 and size $icode n$$.
00216 It specifies value of $latex X(ti)$$.
00217 
00218 $head smin$$
00219 The argument $icode smin$$ has prototype
00220 $codei%
00221      const %Scalar% &%smin%
00222 %$$
00223 The step size during a call to $icode method$$ is defined as
00224 the corresponding value of $latex tb - ta$$.
00225 If $latex tf - ti \leq smin$$,
00226 the integration will be done in one step of size $icode tf - ti$$.
00227 Otherwise,
00228 the minimum value of $icode tb - ta$$ will be $latex smin$$
00229 except for the last two calls to $icode method$$ where it may be
00230 as small as $latex smin / 2$$.
00231 
00232 $head smax$$
00233 The argument $icode smax$$ has prototype
00234 $codei%
00235      const %Scalar% &%smax%
00236 %$$
00237 It specifies the maximum step size to use during the integration; 
00238 i.e., the maximum value for $latex tb - ta$$ in a call to $icode method$$.
00239 The value of $icode smax$$ must be greater than or equal $icode smin$$.
00240 
00241 $head scur$$
00242 The argument $icode scur$$ has prototype
00243 $codei%
00244      %Scalar% &%scur%
00245 %$$
00246 The value of $icode scur$$ is the suggested next step size,
00247 based on error criteria, to try in the next call to $icode method$$.
00248 On input it corresponds to the first call to $icode method$$,
00249 in this call to $code OdeErrControl$$ (where $latex ta = ti$$).
00250 On output it corresponds to the next call to $icode method$$,
00251 in a subsequent call to $code OdeErrControl$$ (where $icode ta = tf$$).
00252 
00253 $head eabs$$
00254 The argument $icode eabs$$ has prototype
00255 $codei%
00256      const %Vector% &%eabs%
00257 %$$
00258 and size $icode n$$.
00259 Each of the elements of $icode eabs$$ must be 
00260 greater than or equal zero.
00261 It specifies a bound for the absolute
00262 error in the return value $icode xf$$ as an approximation for $latex X(tf)$$.
00263 (see the 
00264 $cref/error criteria discussion/OdeErrControl/Error Criteria Discussion/$$ 
00265 below). 
00266 
00267 $head erel$$
00268 The argument $icode erel$$ has prototype
00269 $codei%
00270      const %Scalar% &%erel%
00271 %$$
00272 and is greater than or equal zero.
00273 It specifies a bound for the relative 
00274 error in the return value $icode xf$$ as an approximation for $latex X(tf)$$
00275 (see the 
00276 $cref/error criteria discussion/OdeErrControl/Error Criteria Discussion/$$ 
00277 below). 
00278 
00279 $head ef$$
00280 The argument value $icode ef$$ has prototype
00281 $codei%
00282      %Vector% &%ef%
00283 %$$
00284 and size $icode n$$.
00285 The input value of its elements does not matter.
00286 On output, 
00287 it contains an estimated bound for the 
00288 absolute error in the approximation $icode xf$$; i.e.,
00289 $latex \[
00290      ef_i > | X( tf )_i - xf_i |
00291 \] $$
00292 If on output $icode ef$$ contains not a number $code nan$$,
00293 see the discussion of $cref/step/OdeErrControl/Method/Nan/$$.
00294 
00295 $head maxabs$$
00296 The argument $icode maxabs$$ is optional in the call to $code OdeErrControl$$.
00297 If it is present, it has the prototype
00298 $codei%
00299      %Vector% &%maxabs%
00300 %$$
00301 and size $icode n$$.
00302 The input value of its elements does not matter.
00303 On output, 
00304 it contains an estimate for the 
00305 maximum absolute value of $latex X(t)$$; i.e.,
00306 $latex \[
00307      maxabs[i] \approx \max \left\{ 
00308           | X( t )_i | \; : \;  t \in [ti, tf] 
00309      \right\}
00310 \] $$
00311 
00312 $head nstep$$
00313 The argument $icode nstep$$ is optional in the call to $code OdeErrControl$$.
00314 If it is present, it has the prototype
00315 $codei%
00316      %size_t% &%nstep%
00317 %$$
00318 Its input value does not matter and its output value
00319 is the number of calls to $icode%method%.step%$$ 
00320 used by $code OdeErrControl$$.
00321 
00322 $head Error Criteria Discussion$$
00323 The relative error criteria $icode erel$$ and
00324 absolute error criteria $icode eabs$$ are enforced during each step of the
00325 integration of the ordinary differential equations.
00326 In addition, they are inversely scaled by the step size so that
00327 the total error bound is less than the sum of the error bounds.
00328 To be specific, if $latex \tilde{X} (t)$$ is the approximate solution
00329 at time $latex t$$, 
00330 $icode ta$$ is the initial step time,
00331 and $icode tb$$ is the final step time,
00332 $latex \[
00333 \left| \tilde{X} (tb)_j  - X (tb)_j \right| 
00334 \leq 
00335 \frac{tf - ti}{tb - ta}
00336 \left[ eabs[j] + erel \;  | \tilde{X} (tb)_j | \right] 
00337 \] $$
00338 If $latex X(tb)_j$$ is near zero for some $latex tb \in [ti , tf]$$,
00339 and one uses an absolute error criteria $latex eabs[j]$$ of zero,
00340 the error criteria above will force $code OdeErrControl$$
00341 to use step sizes equal to 
00342 $cref/smin/OdeErrControl/smin/$$
00343 for steps ending near $latex tb$$.
00344 In this case, the error relative to $icode maxabs$$ can be judged after
00345 $code OdeErrControl$$ returns.
00346 If $icode ef$$ is to large relative to $icode maxabs$$, 
00347 $code OdeErrControl$$ can be called again 
00348 with a smaller value of $icode smin$$.
00349 
00350 $head Scalar$$
00351 The type $icode Scalar$$ must satisfy the conditions
00352 for a $cref NumericType$$ type.
00353 The routine $cref CheckNumericType$$ will generate an error message
00354 if this is not the case.
00355 In addition, the following operations must be defined for 
00356 $icode Scalar$$ objects $icode a$$ and $icode b$$:
00357 
00358 $table
00359 $bold Operation$$ $cnext $bold Description$$  $rnext
00360 $icode%a% <= %b%$$ $cnext
00361      returns true (false) if $icode a$$ is less than or equal 
00362      (greater than) $icode b$$.
00363 $rnext
00364 $icode%a% == %b%$$ $cnext
00365      returns true (false) if $icode a$$ is equal to $icode b$$.
00366 $rnext
00367 $codei%log(%a%)%$$ $cnext
00368      returns a $icode Scalar$$ equal to the logarithm of $icode a$$
00369 $rnext
00370 $codei%exp(%a%)%$$ $cnext
00371      returns a $icode Scalar$$ equal to the exponential of $icode a$$
00372 $tend
00373 
00374 
00375 $head Vector$$
00376 The type $icode Vector$$ must be a $cref SimpleVector$$ class with
00377 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$.
00378 The routine $cref CheckSimpleVector$$ will generate an error message
00379 if this is not the case.
00380 
00381 $head Example$$
00382 $children%
00383      example/ode_err_control.cpp%
00384      example/ode_err_maxabs.cpp
00385 %$$
00386 The files
00387 $cref ode_err_control.cpp$$
00388 and
00389 $cref ode_err_maxabs.cpp$$
00390 contain examples and tests of using this routine.
00391 They return true if they succeed and false otherwise.
00392 
00393 $head Theory$$
00394 Let $latex e(s)$$ be the error as a function of the
00395 step size $latex s$$ and suppose that there is a constant
00396 $latex K$$ such that $latex e(s) = K s^m$$.
00397 Let $latex a$$ be our error bound.
00398 Given the value of $latex e(s)$$, a step of size $latex \lambda s$$
00399 would be ok provided that
00400 $latex \[
00401 \begin{array}{rcl}
00402      a  & \geq & e( \lambda s ) (tf - ti) / ( \lambda s ) \\
00403      a  & \geq & K \lambda^m s^m (tf - ti) / ( \lambda s ) \\
00404      a  & \geq & \lambda^{m-1} s^{m-1} (tf - ti) e(s) / s^m \\
00405      a  & \geq & \lambda^{m-1} (tf - ti) e(s) / s           \\
00406      \lambda^{m-1} & \leq & \frac{a}{e(s)} \frac{s}{tf - ti}
00407 \end{array}
00408 \] $$
00409 Thus if the right hand side of the last inequality is greater 
00410 than or equal to one, the step of size $latex s$$ is ok. 
00411 
00412 $head Source Code$$
00413 The source code for this routine is in the file
00414 $code cppad/ode_err_control.hpp$$.
00415 
00416 $end
00417 --------------------------------------------------------------------------
00418 */
00419 
00420 // link exp and log for float and double
00421 # include <cppad/base_require.hpp>
00422 
00423 # include <cppad/local/cppad_assert.hpp>
00424 # include <cppad/check_simple_vector.hpp>
00425 # include <cppad/nan.hpp>
00426 
00427 namespace CppAD { // Begin CppAD namespace
00428 
00429 template <typename Scalar, typename Vector, typename Method>
00430 Vector OdeErrControl(
00431      Method          &method, 
00432      const Scalar    &ti    , 
00433      const Scalar    &tf    , 
00434      const Vector    &xi    , 
00435      const Scalar    &smin  , 
00436      const Scalar    &smax  , 
00437      Scalar          &scur  ,
00438      const Vector    &eabs  , 
00439      const Scalar    &erel  , 
00440      Vector          &ef    ,
00441      Vector          &maxabs,
00442      size_t          &nstep ) 
00443 {
00444      // check simple vector class specifications
00445      CheckSimpleVector<Scalar, Vector>();
00446 
00447      size_t n = size_t(xi.size());
00448 
00449      CPPAD_ASSERT_KNOWN(
00450           smin <= smax,
00451           "Error in OdeErrControl: smin > smax"
00452      );
00453      CPPAD_ASSERT_KNOWN(
00454           size_t(eabs.size()) == n,
00455           "Error in OdeErrControl: size of eabs is not equal to n"
00456      );
00457      CPPAD_ASSERT_KNOWN(
00458           size_t(maxabs.size()) == n,
00459           "Error in OdeErrControl: size of maxabs is not equal to n"
00460      );
00461      size_t m = method.order();
00462      CPPAD_ASSERT_KNOWN(
00463           m > 1,
00464           "Error in OdeErrControl: m is less than or equal one"
00465      );
00466 
00467      bool    ok;
00468      bool    minimum_step;
00469      size_t  i;
00470      Vector xa(n), xb(n), eb(n), nan_vec(n);
00471 
00472      // initialization
00473      Scalar zero(0);
00474      Scalar one(1);
00475      Scalar two(2);
00476      Scalar three(3);
00477      Scalar m1(m-1);
00478      Scalar ta = ti;
00479      for(i = 0; i < n; i++)
00480      {    nan_vec[i] = nan(zero);
00481           ef[i]      = zero;
00482           xa[i]      = xi[i];
00483           if( zero <= xi[i] )
00484                maxabs[i] = xi[i];
00485           else maxabs[i] = - xi[i];
00486 
00487      }  
00488      nstep = 0;
00489 
00490      Scalar tb, step, lambda, axbi, a, r, root;
00491      while( ! (ta == tf) )
00492      {    // start with value suggested by error criteria
00493           step = scur;
00494 
00495           // check maximum
00496           if( smax <= step )
00497                step = smax;
00498 
00499           // check minimum
00500           minimum_step = step <= smin;
00501           if( minimum_step )
00502                step = smin;
00503 
00504           // check if near the end
00505           if( tf <= ta + step * three / two )
00506                tb = tf;
00507           else tb = ta + step;
00508 
00509           // try using this step size
00510           nstep++;
00511           method.step(ta, tb, xa, xb, eb);
00512           step = tb - ta;
00513 
00514           // check if this steps error estimate is ok
00515           ok = ! (hasnan(xb) || hasnan(eb));
00516           if( (! ok) && minimum_step )
00517           {    ef = nan_vec;
00518                return nan_vec;
00519           }
00520 
00521           // compute value of lambda for this step
00522           lambda = Scalar(10) * scur / step;
00523           for(i = 0; i < n; i++)
00524           {    if( zero <= xb[i] )
00525                     axbi = xb[i];
00526                else axbi = - xb[i];
00527                a    = eabs[i] + erel * axbi;
00528                if( ! (eb[i] == zero) )
00529                {    r = ( a / eb[i] ) * step / (tf - ti);
00530                     root = exp( log(r) / m1 ); 
00531                     if( root <= lambda )
00532                          lambda = root;
00533                }
00534           }
00535           if( ok && ( one <= lambda || step <= smin * three / two) )
00536           {    // this step is within error limits or 
00537                // close to the minimum size
00538                ta = tb;
00539                for(i = 0; i < n; i++)
00540                {    xa[i] = xb[i];
00541                     ef[i] = ef[i] + eb[i];
00542                     if( zero <= xb[i] )
00543                          axbi = xb[i];
00544                     else axbi = - xb[i];
00545                     if( axbi > maxabs[i] )
00546                          maxabs[i] = axbi;
00547                }
00548           }
00549           if( ! ok )
00550           {    // decrease step an see if method will work this time
00551                scur = step / two;
00552           }
00553           else if( ! (ta == tf) )
00554           {    // step suggested by the error criteria is not used
00555                // on the last step because it may be very small.
00556                scur = lambda * step / two;
00557           }
00558      }
00559      return xa;
00560 }
00561 
00562 template <typename Scalar, typename Vector, typename Method>
00563 Vector OdeErrControl(
00564      Method          &method, 
00565      const Scalar    &ti    , 
00566      const Scalar    &tf    , 
00567      const Vector    &xi    , 
00568      const Scalar    &smin  , 
00569      const Scalar    &smax  , 
00570      Scalar          &scur  ,
00571      const Vector    &eabs  , 
00572      const Scalar    &erel  , 
00573      Vector          &ef    ) 
00574 {    Vector maxabs(xi.size());
00575      size_t nstep;
00576      return OdeErrControl(
00577      method, ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep
00578      );
00579 }
00580 
00581 template <typename Scalar, typename Vector, typename Method>
00582 Vector OdeErrControl(
00583      Method          &method, 
00584      const Scalar    &ti    , 
00585      const Scalar    &tf    , 
00586      const Vector    &xi    , 
00587      const Scalar    &smin  , 
00588      const Scalar    &smax  , 
00589      Scalar          &scur  ,
00590      const Vector    &eabs  , 
00591      const Scalar    &erel  , 
00592      Vector          &ef    ,
00593      Vector          &maxabs) 
00594 {    size_t nstep;
00595      return OdeErrControl(
00596      method, ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep
00597      );
00598 }
00599 
00600 } // End CppAD namespace 
00601 
00602 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines