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