CppAD: A C++ Algorithmic Differentiation Package 20110419
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-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 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 $syntax%# include <cppad/ode_gear_control.hpp>
00050 %$$
00051 $syntax%%xf% = OdeGearControl(%F%, %M%, %ti%, %tf%, %xi%,
00052         %smin%, %smax%, %sini%, %eabs%, %erel%, %ef% , %maxabs%, %nstep% )%$$
00053 
00054 
00055 $head Purpose$$
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 We define $latex X : [ti , tf] \rightarrow \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 $xref/OdeGearControl/Scalar/Scalar/$$ and
00079 $xref/OdeGearControl/Vector/Vector/$$ are documented below.
00080 
00081 $head xf$$
00082 The return value $italic xf$$ has the prototype
00083 $syntax%
00084         %Vector% %xf%
00085 %$$
00086 and the size of $italic xf$$ is equal to $italic n$$
00087 (see description of $xref/OdeGear/Vector/Vector/$$ below).
00088 It is the approximation for $latex X(tf)$$.
00089 
00090 $head Fun$$
00091 The class $italic Fun$$ 
00092 and the object $italic F$$ satisfy the prototype
00093 $syntax%
00094         %Fun% &%F%
00095 %$$
00096 This must support the following set of calls
00097 $syntax%
00098         %F%.Ode(%t%, %x%, %f%)
00099         %F%.Ode_dep(%t%, %x%, %f_x%)
00100 %$$
00101 
00102 $subhead t$$
00103 The argument $italic t$$ has prototype
00104 $syntax%
00105         const %Scalar% &%t%
00106 %$$
00107 (see description of $xref/OdeGear/Scalar/Scalar/$$ below). 
00108 
00109 $subhead x$$
00110 The argument $italic x$$ has prototype
00111 $syntax%
00112         const %Vector% &%x%
00113 %$$
00114 and has size $italic N$$
00115 (see description of $xref/OdeGear/Vector/Vector/$$ below). 
00116 
00117 $subhead f$$
00118 The argument $italic f$$ to $syntax%%F%.Ode%$$ has prototype
00119 $syntax%
00120         %Vector% &%f%
00121 %$$
00122 On input and output, $italic f$$ is a vector of size $italic N$$
00123 and the input values of the elements of $italic f$$ do not matter.
00124 On output,
00125 $italic f$$ is set equal to $latex f(t, x)$$
00126 (see $italic f(t, x)$$ in $xref/OdeGear/Purpose/Purpose/$$). 
00127 
00128 $subhead f_x$$
00129 The argument $italic f_x$$ has prototype
00130 $syntax%
00131         %Vector% &%f_x%
00132 %$$
00133 On input and output, $italic f_x$$ is a vector of size $latex N * N$$
00134 and the input values of the elements of $italic 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 $italic f$$, and $italic 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 $italic f$$ and $italic f_x$$.
00145 
00146 $head M$$
00147 The argument $italic M$$ has prototype
00148 $syntax%
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 $italic M$$ must greater than or equal one.
00155 
00156 $head ti$$
00157 The argument $italic ti$$ has prototype
00158 $syntax%
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 $italic tf$$ has prototype
00166 $syntax%
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 $italic xi$$ has prototype
00174 $syntax%
00175         const %Vector% &%xi%
00176 %$$
00177 and size $italic n$$.
00178 It specifies value of $latex X(ti)$$.
00179 
00180 $head smin$$
00181 The argument $italic smin$$ has prototype
00182 $syntax%
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 $italic smin$$ must be less than or equal $italic smax$$.
00189 
00190 $head smax$$
00191 The argument $italic smax$$ has prototype
00192 $syntax%
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 $italic sini$$ has prototype
00201 $syntax%
00202         %Scalar% &%sini%
00203 %$$
00204 The value of $italic 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 $italic sini$$ must be less than or equal $italic smax$$
00208 (and can also be less than $italic smin$$).
00209 
00210 $head eabs$$
00211 The argument $italic eabs$$ has prototype
00212 $syntax%
00213         const %Vector% &%eabs%
00214 %$$
00215 and size $italic n$$.
00216 Each of the elements of $italic eabs$$ must be 
00217 greater than or equal zero.
00218 It specifies a bound for the absolute
00219 error in the return value $italic xf$$ as an approximation for $latex X(tf)$$.
00220 (see the 
00221 $xref/OdeGearControl/Error Criteria Discussion/error criteria discussion/$$ 
00222 below). 
00223 
00224 $head erel$$
00225 The argument $italic erel$$ has prototype
00226 $syntax%
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 $italic xf$$ as an approximation for $latex X(tf)$$
00232 (see the 
00233 $xref/OdeGearControl/Error Criteria Discussion/error criteria discussion/$$ 
00234 below). 
00235 
00236 $head ef$$
00237 The argument value $italic ef$$ has prototype
00238 $syntax%
00239         %Vector% &%ef%
00240 %$$
00241 and size $italic 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 $italic xf$$; i.e.,
00246 $latex \[
00247         ef_i > | X( tf )_i - xf_i |
00248 \] $$
00249 
00250 $head maxabs$$
00251 The argument $italic maxabs$$ is optional in the call to $code OdeGearControl$$.
00252 If it is present, it has the prototype
00253 $syntax%
00254         %Vector% &%maxabs%
00255 %$$
00256 and size $italic 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 $italic nstep$$ has the prototype
00269 $syntax%
00270         %size_t% &%nstep%
00271 %$$
00272 Its input value does not matter and its output value
00273 is the number of calls to $xref/OdeGear/$$
00274 used by $code OdeGearControl$$.
00275 
00276 $head Error Criteria Discussion$$
00277 The relative error criteria $italic erel$$ and
00278 absolute error criteria $italic 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 $italic ta$$ is the initial step time,
00285 and $italic 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 $xref/OdeGearControl/smin/smin/$$
00297 for steps ending near $latex tb$$.
00298 In this case, the error relative to $italic maxabs$$ can be judged after
00299 $code OdeGearControl$$ returns.
00300 If $italic ef$$ is to large relative to $italic maxabs$$, 
00301 $code OdeGearControl$$ can be called again 
00302 with a smaller value of $italic smin$$.
00303 
00304 $head Scalar$$
00305 The type $italic Scalar$$ must satisfy the conditions
00306 for a $xref/NumericType/$$ type.
00307 The routine $xref/CheckNumericType/$$ will generate an error message
00308 if this is not the case.
00309 In addition, the following operations must be defined for 
00310 $italic Scalar$$ objects $italic a$$ and $italic b$$:
00311 
00312 $table
00313 $bold Operation$$ $cnext $bold Description$$  $rnext
00314 $syntax%%a% <= %b%$$ $cnext
00315         returns true (false) if $italic a$$ is less than or equal 
00316         (greater than) $italic b$$.
00317 $rnext
00318 $syntax%%a% == %b%$$ $cnext
00319         returns true (false) if $italic a$$ is equal to $italic b$$.
00320 $rnext
00321 $syntax%log(%a%)%$$ $cnext
00322         returns a $italic Scalar$$ equal to the logarithm of $italic a$$
00323 $rnext
00324 $syntax%exp(%a%)%$$ $cnext
00325         returns a $italic Scalar$$ equal to the exponential of $italic a$$
00326 $tend
00327 
00328 
00329 $head Vector$$
00330 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with
00331 $xref/SimpleVector/Elements of Specified Type/elements of type Scalar/$$.
00332 The routine $xref/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 $xref/OdeGearControl.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/std_math_unary.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 = 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                 eabs.size() == n,
00414                 "Error in OdeGearControl: size of eabs is not equal to n"
00415         );
00416         CPPAD_ASSERT_KNOWN(
00417                 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