CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_ROSEN_34_INCLUDED 00003 # define CPPAD_ROSEN_34_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 Rosen34$$ 00018 $spell 00019 cppad.hpp 00020 bool 00021 xf 00022 templated 00023 const 00024 Rosenbrock 00025 CppAD 00026 xi 00027 ti 00028 tf 00029 Karp 00030 Rosen 00031 Shampine 00032 ind 00033 dep 00034 $$ 00035 00036 $index Rosen34$$ 00037 $index ODE, Rosenbrock$$ 00038 $index Rosenbrock, ODE$$ 00039 $index solve, ODE$$ 00040 $index stiff, ODE$$ 00041 $index differential, equation$$ 00042 $index equation, differential$$ 00043 00044 $section A 3rd and 4th Order Rosenbrock ODE Solver$$ 00045 00046 $head Syntax$$ 00047 $syntax%# include <cppad/rosen_34.hpp> 00048 %$$ 00049 $syntax%%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%) 00050 %$$ 00051 $syntax%%xf% = Rosen34(%F%, %M%, %ti%, %tf%, %xi%, %e%) 00052 %$$ 00053 00054 00055 $head Description$$ 00056 This is an embedded 3rd and 4th order Rosenbrock ODE solver 00057 (see Section 16.6 of $xref/Bib/Numerical Recipes/Numerical Recipes/$$ 00058 for a description of Rosenbrock ODE solvers). 00059 In particular, we use the formulas taken from page 100 of 00060 $xref/Bib/Shampine, L.F./Shampine, L.F./$$ 00061 (except that the fraction 98/108 has been correction to be 97/108). 00062 $pre 00063 00064 $$ 00065 We use $latex n$$ for the size of the vector $italic xi$$. 00066 Let $latex \R$$ denote the real numbers 00067 and let $latex F : \R \times \R^n \rightarrow \R^n$$ be a smooth function. 00068 The return value $italic xf$$ contains a 5th order 00069 approximation for the value $latex X(tf)$$ where 00070 $latex X : [ti , tf] \rightarrow \R^n$$ is defined by 00071 the following initial value problem: 00072 $latex \[ 00073 \begin{array}{rcl} 00074 X(ti) & = & xi \\ 00075 X'(t) & = & F[t , X(t)] 00076 \end{array} 00077 \] $$ 00078 If your set of ordinary differential equations are not stiff 00079 an explicit method may be better (perhaps $xref/Runge45/$$.) 00080 00081 $head Include$$ 00082 The file $code cppad/rosen_34.hpp$$ is included by $code cppad/cppad.hpp$$ 00083 but it can also be included separately with out the rest of 00084 the $code CppAD$$ routines. 00085 00086 $head xf$$ 00087 The return value $italic xf$$ has the prototype 00088 $syntax% 00089 %Vector% %xf% 00090 %$$ 00091 and the size of $italic xf$$ is equal to $italic n$$ 00092 (see description of $xref/Rosen34/Vector/Vector/$$ below). 00093 $latex \[ 00094 X(tf) = xf + O( h^5 ) 00095 \] $$ 00096 where $latex h = (tf - ti) / M$$ is the step size. 00097 If $italic xf$$ contains not a number $cref/nan/$$, 00098 see the discussion of $cref/f/Rosen34/Fun/Nan/$$. 00099 00100 $head Fun$$ 00101 The class $italic Fun$$ 00102 and the object $italic F$$ satisfy the prototype 00103 $syntax% 00104 %Fun% &%F% 00105 %$$ 00106 This must support the following set of calls 00107 $syntax% 00108 %F%.Ode(%t%, %x%, %f%) 00109 %F%.Ode_ind(%t%, %x%, %f_t%) 00110 %F%.Ode_dep(%t%, %x%, %f_x%) 00111 %$$ 00112 00113 $subhead t$$ 00114 In all three cases, 00115 the argument $italic t$$ has prototype 00116 $syntax% 00117 const %Scalar% &%t% 00118 %$$ 00119 (see description of $xref/Rosen34/Scalar/Scalar/$$ below). 00120 00121 $subhead x$$ 00122 In all three cases, 00123 the argument $italic x$$ has prototype 00124 $syntax% 00125 const %Vector% &%x% 00126 %$$ 00127 and has size $italic n$$ 00128 (see description of $xref/Rosen34/Vector/Vector/$$ below). 00129 00130 $subhead f$$ 00131 The argument $italic f$$ to $syntax%%F%.Ode%$$ has prototype 00132 $syntax% 00133 %Vector% &%f% 00134 %$$ 00135 On input and output, $italic f$$ is a vector of size $italic n$$ 00136 and the input values of the elements of $italic f$$ do not matter. 00137 On output, 00138 $italic f$$ is set equal to $latex F(t, x)$$ 00139 (see $italic F(t, x)$$ in $xref/Rosen34/Description/Description/$$). 00140 00141 $subhead f_t$$ 00142 The argument $italic f_t$$ to $syntax%%F%.Ode_ind%$$ has prototype 00143 $syntax% 00144 %Vector% &%f_t% 00145 %$$ 00146 On input and output, $italic f_t$$ is a vector of size $italic n$$ 00147 and the input values of the elements of $italic f_t$$ do not matter. 00148 On output, the $th i$$ element of 00149 $italic f_t$$ is set equal to $latex \partial_t F_i (t, x)$$ 00150 (see $italic F(t, x)$$ in $xref/Rosen34/Description/Description/$$). 00151 00152 $subhead f_x$$ 00153 The argument $italic f_x$$ to $syntax%%F%.Ode_dep%$$ has prototype 00154 $syntax% 00155 %Vector% &%f_x% 00156 %$$ 00157 On input and output, $italic f_x$$ is a vector of size $syntax%%n%*%n%$$ 00158 and the input values of the elements of $italic f_x$$ do not matter. 00159 On output, the [$syntax%%i%*%n%+%j%$$] element of 00160 $italic f_x$$ is set equal to $latex \partial_{x(j)} F_i (t, x)$$ 00161 (see $italic F(t, x)$$ in $xref/Rosen34/Description/Description/$$). 00162 00163 $subhead Nan$$ 00164 If any of the elements of $italic f$$, $italic f_t$$, or $italic f_x$$ 00165 have the value not a number $code nan$$, 00166 the routine $code Rosen34$$ returns with all the 00167 elements of $italic xf$$ and $italic e$$ equal to $code nan$$. 00168 00169 $subhead Warning$$ 00170 The arguments $italic f$$, $italic f_t$$, and $italic f_x$$ 00171 must have a call by reference in their prototypes; i.e., 00172 do not forget the $code &$$ in the prototype for 00173 $italic f$$, $italic f_t$$ and $italic f_x$$. 00174 00175 $subhead Optimization$$ 00176 Every call of the form 00177 $syntax% 00178 %F%.Ode_ind(%t%, %x%, %f_t%) 00179 %$$ 00180 is directly followed by a call of the form 00181 $syntax% 00182 %F%.Ode_dep(%t%, %x%, %f_x%) 00183 %$$ 00184 where the arguments $italic t$$ and $italic x$$ have not changed between calls. 00185 In many cases it is faster to compute the values of $italic f_t$$ 00186 and $italic f_x$$ together and then pass them back one at a time. 00187 00188 $head M$$ 00189 The argument $italic M$$ has prototype 00190 $syntax% 00191 size_t %M% 00192 %$$ 00193 It specifies the number of steps 00194 to use when solving the differential equation. 00195 This must be greater than or equal one. 00196 The step size is given by $latex h = (tf - ti) / M$$, thus 00197 the larger $italic M$$, the more accurate the 00198 return value $italic xf$$ is as an approximation 00199 for $latex X(tf)$$. 00200 00201 $head ti$$ 00202 The argument $italic ti$$ has prototype 00203 $syntax% 00204 const %Scalar% &%ti% 00205 %$$ 00206 (see description of $xref/Rosen34/Scalar/Scalar/$$ below). 00207 It specifies the initial time for $italic t$$ in the 00208 differential equation; i.e., 00209 the time corresponding to the value $italic xi$$. 00210 00211 $head tf$$ 00212 The argument $italic tf$$ has prototype 00213 $syntax% 00214 const %Scalar% &%tf% 00215 %$$ 00216 It specifies the final time for $italic t$$ in the 00217 differential equation; i.e., 00218 the time corresponding to the value $italic xf$$. 00219 00220 $head xi$$ 00221 The argument $italic xi$$ has the prototype 00222 $syntax% 00223 const %Vector% &%xi% 00224 %$$ 00225 and the size of $italic xi$$ is equal to $italic n$$. 00226 It specifies the value of $latex X(ti)$$ 00227 00228 $head e$$ 00229 The argument $italic e$$ is optional and has the prototype 00230 $syntax% 00231 %Vector% &%e% 00232 %$$ 00233 If $italic e$$ is present, 00234 the size of $italic e$$ must be equal to $italic n$$. 00235 The input value of the elements of $italic e$$ does not matter. 00236 On output 00237 it contains an element by element 00238 estimated bound for the absolute value of the error in $italic xf$$ 00239 $latex \[ 00240 e = O( h^4 ) 00241 \] $$ 00242 where $latex h = (tf - ti) / M$$ is the step size. 00243 00244 $head Scalar$$ 00245 The type $italic Scalar$$ must satisfy the conditions 00246 for a $xref/NumericType/$$ type. 00247 The routine $xref/CheckNumericType/$$ will generate an error message 00248 if this is not the case. 00249 In addition, the following operations must be defined for 00250 $italic Scalar$$ objects $italic a$$ and $italic b$$: 00251 00252 $table 00253 $bold Operation$$ $cnext $bold Description$$ $rnext 00254 $syntax%%a% < %b%$$ $cnext 00255 less than operator (returns a $code bool$$ object) 00256 $tend 00257 00258 $head Vector$$ 00259 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with 00260 $xref/SimpleVector/Elements of Specified Type/elements of type Scalar/$$. 00261 The routine $xref/CheckSimpleVector/$$ will generate an error message 00262 if this is not the case. 00263 00264 $head Example$$ 00265 $children% 00266 example/rosen_34.cpp 00267 %$$ 00268 The file 00269 $xref/Rosen34.cpp/$$ 00270 contains an example and test a test of using this routine. 00271 It returns true if it succeeds and false otherwise. 00272 00273 $head Source Code$$ 00274 The source code for this routine is in the file 00275 $code cppad/rosen_34.hpp$$. 00276 00277 $end 00278 -------------------------------------------------------------------------- 00279 */ 00280 00281 # include <cstddef> 00282 # include <cppad/local/cppad_assert.hpp> 00283 # include <cppad/check_simple_vector.hpp> 00284 # include <cppad/check_numeric_type.hpp> 00285 # include <cppad/vector.hpp> 00286 # include <cppad/lu_factor.hpp> 00287 # include <cppad/lu_invert.hpp> 00288 00289 namespace CppAD { // BEGIN CppAD namespace 00290 00291 template <typename Scalar, typename Vector, typename Fun> 00292 Vector Rosen34( 00293 Fun &F , 00294 size_t M , 00295 const Scalar &ti , 00296 const Scalar &tf , 00297 const Vector &xi ) 00298 { Vector e( xi.size() ); 00299 return Rosen34(F, M, ti, tf, xi, e); 00300 } 00301 00302 template <typename Scalar, typename Vector, typename Fun> 00303 Vector Rosen34( 00304 Fun &F , 00305 size_t M , 00306 const Scalar &ti , 00307 const Scalar &tf , 00308 const Vector &xi , 00309 Vector &e ) 00310 { 00311 // check numeric type specifications 00312 CheckNumericType<Scalar>(); 00313 00314 // check simple vector class specifications 00315 CheckSimpleVector<Scalar, Vector>(); 00316 00317 // Parameters for Shampine's Rosenbrock method 00318 // are static to avoid recalculation on each call and 00319 // do not use Vector to avoid possible memory leak 00320 static Scalar a[3] = { 00321 Scalar(0), 00322 Scalar(1), 00323 Scalar(3) / Scalar(5) 00324 }; 00325 static Scalar b[2 * 2] = { 00326 Scalar(1), 00327 Scalar(0), 00328 Scalar(24) / Scalar(25), 00329 Scalar(3) / Scalar(25) 00330 }; 00331 static Scalar ct[4] = { 00332 Scalar(1) / Scalar(2), 00333 - Scalar(3) / Scalar(2), 00334 Scalar(121) / Scalar(50), 00335 Scalar(29) / Scalar(250) 00336 }; 00337 static Scalar cg[3 * 3] = { 00338 - Scalar(4), 00339 Scalar(0), 00340 Scalar(0), 00341 Scalar(186) / Scalar(25), 00342 Scalar(6) / Scalar(5), 00343 Scalar(0), 00344 - Scalar(56) / Scalar(125), 00345 - Scalar(27) / Scalar(125), 00346 - Scalar(1) / Scalar(5) 00347 }; 00348 static Scalar d3[3] = { 00349 Scalar(97) / Scalar(108), 00350 Scalar(11) / Scalar(72), 00351 Scalar(25) / Scalar(216) 00352 }; 00353 static Scalar d4[4] = { 00354 Scalar(19) / Scalar(18), 00355 Scalar(1) / Scalar(4), 00356 Scalar(25) / Scalar(216), 00357 Scalar(125) / Scalar(216) 00358 }; 00359 CPPAD_ASSERT_KNOWN( 00360 M >= 1, 00361 "Error in Rosen34: the number of steps is less than one" 00362 ); 00363 CPPAD_ASSERT_KNOWN( 00364 e.size() == xi.size(), 00365 "Error in Rosen34: size of e not equal to size of xi" 00366 ); 00367 size_t i, j, k, l, m; // indices 00368 00369 size_t n = xi.size(); // number of components in X(t) 00370 Scalar ns = Scalar(double(M)); // number of steps as Scalar object 00371 Scalar h = (tf - ti) / ns; // step size 00372 Scalar zero = Scalar(0); // some constants 00373 Scalar one = Scalar(1); 00374 Scalar two = Scalar(2); 00375 00376 // permutation vectors needed for LU factorization routine 00377 CppAD::vector<size_t> ip(n), jp(n); 00378 00379 // vectors used to store values returned by F 00380 Vector E(n * n), Eg(n), f_t(n); 00381 Vector g(n * 3), x3(n), x4(n), xf(n), ftmp(n), xtmp(n), nan_vec(n); 00382 00383 // initialize e = 0, nan_vec = nan 00384 for(i = 0; i < n; i++) 00385 { e[i] = zero; 00386 nan_vec[i] = nan(zero); 00387 } 00388 00389 xf = xi; // initialize solution 00390 for(m = 0; m < M; m++) 00391 { // time at beginning of this interval 00392 Scalar t = ti * (Scalar(int(M - m)) / ns) 00393 + tf * (Scalar(int(m)) / ns); 00394 00395 // value of x at beginning of this interval 00396 x3 = x4 = xf; 00397 00398 // evaluate partial derivatives at beginning of this interval 00399 F.Ode_ind(t, xf, f_t); 00400 F.Ode_dep(t, xf, E); // E = f_x 00401 if( hasnan(f_t) || hasnan(E) ) 00402 { e = nan_vec; 00403 return nan_vec; 00404 } 00405 00406 // E = I - f_x * h / 2 00407 for(i = 0; i < n; i++) 00408 { for(j = 0; j < n; j++) 00409 E[i * n + j] = - E[i * n + j] * h / two; 00410 E[i * n + i] += one; 00411 } 00412 00413 // LU factor the matrix E 00414 int sign = LuFactor(ip, jp, E); 00415 CPPAD_ASSERT_KNOWN( 00416 sign != 0, 00417 "Error in Rosen34: I - f_x * h / 2 not invertible" 00418 ); 00419 00420 // loop over integration steps 00421 for(k = 0; k < 3; k++) 00422 { // set location for next function evaluation 00423 xtmp = xf; 00424 for(l = 0; l < k; l++) 00425 { // loop over previous function evaluations 00426 Scalar bkl = b[(k-1)*2 + l]; 00427 for(i = 0; i < n; i++) 00428 { // loop over elements of x 00429 xtmp[i] += bkl * g[i*3 + l] * h; 00430 } 00431 } 00432 // ftmp = F(t + a[k] * h, xtmp) 00433 F.Ode(t + a[k] * h, xtmp, ftmp); 00434 if( hasnan(ftmp) ) 00435 { e = nan_vec; 00436 return nan_vec; 00437 } 00438 00439 // Form Eg for this integration step 00440 for(i = 0; i < n; i++) 00441 Eg[i] = ftmp[i] + ct[k] * f_t[i] * h; 00442 for(l = 0; l < k; l++) 00443 { for(i = 0; i < n; i++) 00444 Eg[i] += cg[(k-1)*3 + l] * g[i*3 + l]; 00445 } 00446 00447 // Solve the equation E * g = Eg 00448 LuInvert(ip, jp, E, Eg); 00449 00450 // save solution and advance x3, x4 00451 for(i = 0; i < n; i++) 00452 { g[i*3 + k] = Eg[i]; 00453 x3[i] += h * d3[k] * Eg[i]; 00454 x4[i] += h * d4[k] * Eg[i]; 00455 } 00456 } 00457 // Form Eg for last update to x4 only 00458 for(i = 0; i < n; i++) 00459 Eg[i] = ftmp[i] + ct[3] * f_t[i] * h; 00460 for(l = 0; l < 3; l++) 00461 { for(i = 0; i < n; i++) 00462 Eg[i] += cg[2*3 + l] * g[i*3 + l]; 00463 } 00464 00465 // Solve the equation E * g = Eg 00466 LuInvert(ip, jp, E, Eg); 00467 00468 // advance x4 and accumulate error bound 00469 for(i = 0; i < n; i++) 00470 { x4[i] += h * d4[3] * Eg[i]; 00471 00472 // cant use abs because cppad.hpp may not be included 00473 Scalar diff = x4[i] - x3[i]; 00474 if( diff < zero ) 00475 e[i] -= diff; 00476 else e[i] += diff; 00477 } 00478 00479 // advance xf for this step using x4 00480 xf = x4; 00481 } 00482 return xf; 00483 } 00484 00485 } // End CppAD namespace 00486 00487 # endif