CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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