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