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