CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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-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 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 $syntax%%xf% = OdeErrControl(%method%, %ti%, %tf%, %xi%, 00056 %smin%, %smax%, %scur%, %eabs%, %erel%, %ef% , %maxabs%, %nstep% )%$$ 00057 00058 00059 $head Description$$ 00060 Let $latex \R$$ denote the real numbers 00061 and let $latex F : \R \times \R^n \rightarrow \R^n$$ be a smooth function. 00062 We define $latex X : [ti , tf] \rightarrow \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 $xref/OdeErrControl/Scalar/Scalar/$$ and 00082 $xref/OdeErrControl/Vector/Vector/$$ are documented below. 00083 00084 $head xf$$ 00085 The return value $italic xf$$ has the prototype 00086 $syntax% 00087 %Vector% %xf% 00088 %$$ 00089 (see description of $xref/OdeErrControl/Vector/Vector/$$ below). 00090 and the size of $italic xf$$ is equal to $italic n$$. 00091 If $italic xf$$ contains not a number $cref/nan/$$, 00092 see the discussion of $cref/step/OdeErrControl/Method/Nan/$$. 00093 00094 $head Method$$ 00095 The class $italic Method$$ 00096 and the object $italic method$$ satisfy the following syntax 00097 $syntax% 00098 %Method% &%method% 00099 %$$ 00100 The object $italic method$$ must support $code step$$ and 00101 $code order$$ member functions defined below: 00102 00103 $subhead step$$ 00104 The syntax 00105 $syntax% 00106 %method%.step(%ta%, %tb%, %xa%, %xb%, %eb%) 00107 %$$ 00108 executes one step of the integration method. 00109 $syntax% 00110 00111 %ta% 00112 %$$ 00113 The argument $italic ta$$ has prototype 00114 $syntax% 00115 const %Scalar% &%ta% 00116 %$$ 00117 It specifies the initial time for this step in the 00118 ODE integration. 00119 (see description of $xref/OdeErrControl/Scalar/Scalar/$$ below). 00120 $syntax% 00121 00122 %tb% 00123 %$$ 00124 The argument $italic tb$$ has prototype 00125 $syntax% 00126 const %Scalar% &%tb% 00127 %$$ 00128 It specifies the final time for this step in the 00129 ODE integration. 00130 $syntax% 00131 00132 %xa% 00133 %$$ 00134 The argument $italic xa$$ has prototype 00135 $syntax% 00136 const %Vector% &%xa% 00137 %$$ 00138 and size $italic n$$. 00139 It specifies the value of $latex X(ta)$$. 00140 (see description of $xref/OdeErrControl/Vector/Vector/$$ below). 00141 $syntax% 00142 00143 %xb% 00144 %$$ 00145 The argument value $italic xb$$ has prototype 00146 $syntax% 00147 %Vector% &%xb% 00148 %$$ 00149 and size $italic 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 $syntax% 00154 00155 %eb% 00156 %$$ 00157 The argument value $italic eb$$ has prototype 00158 $syntax% 00159 %Vector% &%eb% 00160 %$$ 00161 and size $italic 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 $italic xb$$. 00165 It is assumed (locally) that the error bound in this approximation 00166 nearly equal to $latex K (tb - ta)^m$$ 00167 where $italic K$$ is a fixed constant and $italic m$$ 00168 is the corresponding argument to $code CodeControl$$. 00169 00170 $subhead Nan$$ 00171 If any element of the vector $italic eb$$ or $italic 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 $italic smin$$, 00175 $code OdeErrControl$$ returns with $italic xf$$ and $italic ef$$ as vectors 00176 of $code nan$$. 00177 00178 $subhead order$$ 00179 If $italic m$$ is $code size_t$$, 00180 the object $italic method$$ must also support the following syntax 00181 $syntax% 00182 %m% = %method%.order() 00183 %$$ 00184 The return value $italic 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 $italic ta$$, $italic tb$$, and $italic eb$$ are as in 00190 $syntax%%method%.step(%ta%, %tb%, %xa%, %xb%, %eb%)%$$ 00191 00192 00193 $head ti$$ 00194 The argument $italic ti$$ has prototype 00195 $syntax% 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 $italic tf$$ has prototype 00204 $syntax% 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 $italic xi$$ has prototype 00212 $syntax% 00213 const %Vector% &%xi% 00214 %$$ 00215 and size $italic n$$. 00216 It specifies value of $latex X(ti)$$. 00217 00218 $head smin$$ 00219 The argument $italic smin$$ has prototype 00220 $syntax% 00221 const %Scalar% &%smin% 00222 %$$ 00223 The step size during a call to $italic 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 $italic tf - ti$$. 00227 Otherwise, 00228 the minimum value of $italic tb - ta$$ will be $latex smin$$ 00229 except for the last two calls to $italic method$$ where it may be 00230 as small as $latex smin / 2$$. 00231 00232 $head smax$$ 00233 The argument $italic smax$$ has prototype 00234 $syntax% 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 $italic method$$. 00239 The value of $italic smax$$ must be greater than or equal $italic smin$$. 00240 00241 $head scur$$ 00242 The argument $italic scur$$ has prototype 00243 $syntax% 00244 %Scalar% &%scur% 00245 %$$ 00246 The value of $italic scur$$ is the suggested next step size, 00247 based on error criteria, to try in the next call to $italic method$$. 00248 On input it corresponds to the first call to $italic method$$, 00249 in this call to $code OdeErrControl$$ (where $latex ta = ti$$). 00250 On output it corresponds to the next call to $italic method$$, 00251 in a subsequent call to $code OdeErrControl$$ (where $italic ta = tf$$). 00252 00253 $head eabs$$ 00254 The argument $italic eabs$$ has prototype 00255 $syntax% 00256 const %Vector% &%eabs% 00257 %$$ 00258 and size $italic n$$. 00259 Each of the elements of $italic eabs$$ must be 00260 greater than or equal zero. 00261 It specifies a bound for the absolute 00262 error in the return value $italic xf$$ as an approximation for $latex X(tf)$$. 00263 (see the 00264 $xref/OdeErrControl/Error Criteria Discussion/error criteria discussion/$$ 00265 below). 00266 00267 $head erel$$ 00268 The argument $italic erel$$ has prototype 00269 $syntax% 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 $italic xf$$ as an approximation for $latex X(tf)$$ 00275 (see the 00276 $xref/OdeErrControl/Error Criteria Discussion/error criteria discussion/$$ 00277 below). 00278 00279 $head ef$$ 00280 The argument value $italic ef$$ has prototype 00281 $syntax% 00282 %Vector% &%ef% 00283 %$$ 00284 and size $italic 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 $italic xf$$; i.e., 00289 $latex \[ 00290 ef_i > | X( tf )_i - xf_i | 00291 \] $$ 00292 If on output $italic ef$$ contains not a number $code nan$$, 00293 see the discussion of $cref/step/OdeErrControl/Method/Nan/$$. 00294 00295 $head maxabs$$ 00296 The argument $italic maxabs$$ is optional in the call to $code OdeErrControl$$. 00297 If it is present, it has the prototype 00298 $syntax% 00299 %Vector% &%maxabs% 00300 %$$ 00301 and size $italic 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 $italic nstep$$ is optional in the call to $code OdeErrControl$$. 00314 If it is present, it has the prototype 00315 $syntax% 00316 %size_t% &%nstep% 00317 %$$ 00318 Its input value does not matter and its output value 00319 is the number of calls to $syntax%%method%.step%$$ 00320 used by $code OdeErrControl$$. 00321 00322 $head Error Criteria Discussion$$ 00323 The relative error criteria $italic erel$$ and 00324 absolute error criteria $italic 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 $italic ta$$ is the initial step time, 00331 and $italic 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 $xref/OdeErrControl/smin/smin/$$ 00343 for steps ending near $latex tb$$. 00344 In this case, the error relative to $italic maxabs$$ can be judged after 00345 $code OdeErrControl$$ returns. 00346 If $italic ef$$ is to large relative to $italic maxabs$$, 00347 $code OdeErrControl$$ can be called again 00348 with a smaller value of $italic smin$$. 00349 00350 $head Scalar$$ 00351 The type $italic Scalar$$ must satisfy the conditions 00352 for a $xref/NumericType/$$ type. 00353 The routine $xref/CheckNumericType/$$ will generate an error message 00354 if this is not the case. 00355 In addition, the following operations must be defined for 00356 $italic Scalar$$ objects $italic a$$ and $italic b$$: 00357 00358 $table 00359 $bold Operation$$ $cnext $bold Description$$ $rnext 00360 $syntax%%a% <= %b%$$ $cnext 00361 returns true (false) if $italic a$$ is less than or equal 00362 (greater than) $italic b$$. 00363 $rnext 00364 $syntax%%a% == %b%$$ $cnext 00365 returns true (false) if $italic a$$ is equal to $italic b$$. 00366 $rnext 00367 $syntax%log(%a%)%$$ $cnext 00368 returns a $italic Scalar$$ equal to the logarithm of $italic a$$ 00369 $rnext 00370 $syntax%exp(%a%)%$$ $cnext 00371 returns a $italic Scalar$$ equal to the exponential of $italic a$$ 00372 $tend 00373 00374 00375 $head Vector$$ 00376 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with 00377 $xref/SimpleVector/Elements of Specified Type/elements of type Scalar/$$. 00378 The routine $xref/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 $xref/OdeErrControl.cpp/$$ 00388 and 00389 $xref/OdeErrMaxabs.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/std_math_unary.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 = xi.size(); 00448 00449 CPPAD_ASSERT_KNOWN( 00450 smin <= smax, 00451 "Error in OdeErrControl: smin > smax" 00452 ); 00453 CPPAD_ASSERT_KNOWN( 00454 eabs.size() == n, 00455 "Error in OdeErrControl: size of eabs is not equal to n" 00456 ); 00457 CPPAD_ASSERT_KNOWN( 00458 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