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