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