CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_RUNGE_45_INCLUDED 00003 # define CPPAD_RUNGE_45_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 Runge45$$ 00018 $spell 00019 std 00020 fabs 00021 cppad.hpp 00022 bool 00023 xf 00024 templated 00025 const 00026 Runge-Kutta 00027 CppAD 00028 xi 00029 ti 00030 tf 00031 Karp 00032 $$ 00033 00034 $index Runge45$$ 00035 $index ODE, Runge-Kutta$$ 00036 $index Runge, ODE$$ 00037 $index Kutta, ODE$$ 00038 $index solve, ODE$$ 00039 $index differential, equation$$ 00040 $index equation, differential$$ 00041 00042 $section An Embedded 4th and 5th Order Runge-Kutta ODE Solver$$ 00043 00044 $head Syntax$$ 00045 $codei%# include <cppad/runge_45.hpp> 00046 %$$ 00047 $icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%) 00048 %$$ 00049 $icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%, %e%) 00050 %$$ 00051 00052 00053 $head Purpose$$ 00054 This is an implementation of the 00055 Cash-Karp embedded 4th and 5th order Runge-Kutta ODE solver 00056 described in Section 16.2 of $cref/Numerical Recipes/Bib/Numerical Recipes/$$. 00057 We use $latex n$$ for the size of the vector $icode xi$$. 00058 Let $latex \B{R}$$ denote the real numbers 00059 and let $latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ 00060 be a smooth function. 00061 The return value $icode xf$$ contains a 5th order 00062 approximation for the value $latex X(tf)$$ where 00063 $latex X : [ti , tf] \rightarrow \B{R}^n$$ is defined by 00064 the following initial value problem: 00065 $latex \[ 00066 \begin{array}{rcl} 00067 X(ti) & = & xi \\ 00068 X'(t) & = & F[t , X(t)] 00069 \end{array} 00070 \] $$ 00071 If your set of ordinary differential equations 00072 are stiff, an implicit method may be better 00073 (perhaps $cref Rosen34$$.) 00074 00075 $head Operation Sequence$$ 00076 The $cref/operation sequence/glossary/Operation/Sequence/$$ for $icode Runge$$ 00077 does not depend on any of its $icode Scalar$$ input values provided that 00078 the operation sequence for 00079 $codei% 00080 %F%.Ode(%t%, %x%, %f%) 00081 %$$ 00082 does not on any of its $icode Scalar$$ inputs (see below). 00083 00084 $head Include$$ 00085 The file $code cppad/runge_45.hpp$$ is included by $code cppad/cppad.hpp$$ 00086 but it can also be included separately with out the rest of 00087 the $code CppAD$$ routines. 00088 00089 $head xf$$ 00090 The return value $icode xf$$ has the prototype 00091 $codei% 00092 %Vector% %xf% 00093 %$$ 00094 and the size of $icode xf$$ is equal to $icode n$$ 00095 (see description of $cref/Vector/Runge45/Vector/$$ below). 00096 $latex \[ 00097 X(tf) = xf + O( h^6 ) 00098 \] $$ 00099 where $latex h = (tf - ti) / M$$ is the step size. 00100 If $icode xf$$ contains not a number $cref nan$$, 00101 see the discussion for $cref/f/Runge45/Fun/f/$$. 00102 00103 $head Fun$$ 00104 The class $icode Fun$$ 00105 and the object $icode F$$ satisfy the prototype 00106 $codei% 00107 %Fun% &%F% 00108 %$$ 00109 The object $icode F$$ (and the class $icode Fun$$) 00110 must have a member function named $code Ode$$ 00111 that supports the syntax 00112 $codei% 00113 %F%.Ode(%t%, %x%, %f%) 00114 %$$ 00115 00116 $subhead t$$ 00117 The argument $icode t$$ to $icode%F%.Ode%$$ has prototype 00118 $codei% 00119 const %Scalar% &%t% 00120 %$$ 00121 (see description of $cref/Scalar/Runge45/Scalar/$$ below). 00122 00123 $subhead x$$ 00124 The argument $icode x$$ to $icode%F%.Ode%$$ has prototype 00125 $codei% 00126 const %Vector% &%x% 00127 %$$ 00128 and has size $icode n$$ 00129 (see description of $cref/Vector/Runge45/Vector/$$ below). 00130 00131 $subhead f$$ 00132 The argument $icode f$$ to $icode%F%.Ode%$$ has prototype 00133 $codei% 00134 %Vector% &%f% 00135 %$$ 00136 On input and output, $icode f$$ is a vector of size $icode n$$ 00137 and the input values of the elements of $icode f$$ do not matter. 00138 On output, 00139 $icode f$$ is set equal to $latex F(t, x)$$ in the differential equation. 00140 If any of the elements of $icode f$$ have the value not a number $code nan$$ 00141 the routine $code Runge45$$ returns with all the 00142 elements of $icode xf$$ and $icode e$$ equal to $code nan$$. 00143 00144 $subhead Warning$$ 00145 The argument $icode f$$ to $icode%F%.Ode%$$ 00146 must have a call by reference in its prototype; i.e., 00147 do not forget the $code &$$ in the prototype for $icode f$$. 00148 00149 $head M$$ 00150 The argument $icode M$$ has prototype 00151 $codei% 00152 size_t %M% 00153 %$$ 00154 It specifies the number of steps 00155 to use when solving the differential equation. 00156 This must be greater than or equal one. 00157 The step size is given by $latex h = (tf - ti) / M$$, thus 00158 the larger $icode M$$, the more accurate the 00159 return value $icode xf$$ is as an approximation 00160 for $latex X(tf)$$. 00161 00162 $head ti$$ 00163 The argument $icode ti$$ has prototype 00164 $codei% 00165 const %Scalar% &%ti% 00166 %$$ 00167 (see description of $cref/Scalar/Runge45/Scalar/$$ below). 00168 It specifies the initial time for $icode t$$ in the 00169 differential equation; i.e., 00170 the time corresponding to the value $icode xi$$. 00171 00172 $head tf$$ 00173 The argument $icode tf$$ has prototype 00174 $codei% 00175 const %Scalar% &%tf% 00176 %$$ 00177 It specifies the final time for $icode t$$ in the 00178 differential equation; i.e., 00179 the time corresponding to the value $icode xf$$. 00180 00181 $head xi$$ 00182 The argument $icode xi$$ has the prototype 00183 $codei% 00184 const %Vector% &%xi% 00185 %$$ 00186 and the size of $icode xi$$ is equal to $icode n$$. 00187 It specifies the value of $latex X(ti)$$ 00188 00189 $head e$$ 00190 The argument $icode e$$ is optional and has the prototype 00191 $codei% 00192 %Vector% &%e% 00193 %$$ 00194 If $icode e$$ is present, 00195 the size of $icode e$$ must be equal to $icode n$$. 00196 The input value of the elements of $icode e$$ does not matter. 00197 On output 00198 it contains an element by element 00199 estimated bound for the absolute value of the error in $icode xf$$ 00200 $latex \[ 00201 e = O( h^5 ) 00202 \] $$ 00203 where $latex h = (tf - ti) / M$$ is the step size. 00204 If on output, $icode e$$ contains not a number $code nan$$, 00205 see the discussion for $cref/f/Runge45/Fun/f/$$. 00206 00207 $head Scalar$$ 00208 The type $icode Scalar$$ must satisfy the conditions 00209 for a $cref NumericType$$ type. 00210 The routine $cref CheckNumericType$$ will generate an error message 00211 if this is not the case. 00212 00213 $subhead fabs$$ 00214 In addition, the following function must be defined for 00215 $icode Scalar$$ objects $icode a$$ and $icode b$$ 00216 $codei% 00217 %a% = fabs(%b%) 00218 %$$ 00219 Note that this operation is only used for computing $icode e$$; hence 00220 the operation sequence for $icode xf$$ can still be independent of 00221 the arguments to $code Runge45$$ even if 00222 $codei% 00223 fabs(%b%) = std::max(-%b%, %b%) 00224 %$$. 00225 00226 $head Vector$$ 00227 The type $icode Vector$$ must be a $cref SimpleVector$$ class with 00228 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$. 00229 The routine $cref CheckSimpleVector$$ will generate an error message 00230 if this is not the case. 00231 00232 $head Parallel Mode$$ 00233 $index parallel, Runge45$$ 00234 $index Runge45, parallel$$ 00235 For each set of types 00236 $cref/Scalar/Runge45/Scalar/$$, 00237 $cref/Vector/Runge45/Vector/$$, and 00238 $cref/Fun/Runge45/Fun/$$, 00239 the first call to $code Runge45$$ 00240 must not be $cref/parallel/ta_in_parallel/$$ execution mode. 00241 00242 00243 $head Example$$ 00244 $children% 00245 example/runge45_1.cpp% 00246 example/runge45_2.cpp 00247 %$$ 00248 The file 00249 $cref runge45_1.cpp$$ 00250 contains a simple example and test of $code Runge45$$. 00251 It returns true if it succeeds and false otherwise. 00252 $pre 00253 00254 $$ 00255 The file 00256 $cref runge45_2.cpp$$ contains an example using $code Runge45$$ 00257 in the context of algorithmic differentiation. 00258 It also returns true if it succeeds and false otherwise. 00259 00260 $head Source Code$$ 00261 The source code for this routine is in the file 00262 $code cppad/runge_45.hpp$$. 00263 00264 $end 00265 -------------------------------------------------------------------------- 00266 */ 00267 # include <cstddef> 00268 # include <cppad/local/cppad_assert.hpp> 00269 # include <cppad/check_simple_vector.hpp> 00270 # include <cppad/check_numeric_type.hpp> 00271 # include <cppad/nan.hpp> 00272 00273 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL 00274 # include <cppad/thread_alloc.hpp> 00275 00276 namespace CppAD { // BEGIN CppAD namespace 00277 00278 template <typename Scalar, typename Vector, typename Fun> 00279 Vector Runge45( 00280 Fun &F , 00281 size_t M , 00282 const Scalar &ti , 00283 const Scalar &tf , 00284 const Vector &xi ) 00285 { Vector e( xi.size() ); 00286 return Runge45(F, M, ti, tf, xi, e); 00287 } 00288 00289 template <typename Scalar, typename Vector, typename Fun> 00290 Vector Runge45( 00291 Fun &F , 00292 size_t M , 00293 const Scalar &ti , 00294 const Scalar &tf , 00295 const Vector &xi , 00296 Vector &e ) 00297 { 00298 CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL; 00299 00300 // check numeric type specifications 00301 CheckNumericType<Scalar>(); 00302 00303 // check simple vector class specifications 00304 CheckSimpleVector<Scalar, Vector>(); 00305 00306 // Cash-Karp parameters for embedded Runge-Kutta method 00307 // are static to avoid recalculation on each call and 00308 // do not use Vector to avoid possible memory leak 00309 static Scalar a[6] = { 00310 Scalar(0), 00311 Scalar(1) / Scalar(5), 00312 Scalar(3) / Scalar(10), 00313 Scalar(3) / Scalar(5), 00314 Scalar(1), 00315 Scalar(7) / Scalar(8) 00316 }; 00317 static Scalar b[5 * 5] = { 00318 Scalar(1) / Scalar(5), 00319 Scalar(0), 00320 Scalar(0), 00321 Scalar(0), 00322 Scalar(0), 00323 00324 Scalar(3) / Scalar(40), 00325 Scalar(9) / Scalar(40), 00326 Scalar(0), 00327 Scalar(0), 00328 Scalar(0), 00329 00330 Scalar(3) / Scalar(10), 00331 -Scalar(9) / Scalar(10), 00332 Scalar(6) / Scalar(5), 00333 Scalar(0), 00334 Scalar(0), 00335 00336 -Scalar(11) / Scalar(54), 00337 Scalar(5) / Scalar(2), 00338 -Scalar(70) / Scalar(27), 00339 Scalar(35) / Scalar(27), 00340 Scalar(0), 00341 00342 Scalar(1631) / Scalar(55296), 00343 Scalar(175) / Scalar(512), 00344 Scalar(575) / Scalar(13824), 00345 Scalar(44275) / Scalar(110592), 00346 Scalar(253) / Scalar(4096) 00347 }; 00348 static Scalar c4[6] = { 00349 Scalar(2825) / Scalar(27648), 00350 Scalar(0), 00351 Scalar(18575) / Scalar(48384), 00352 Scalar(13525) / Scalar(55296), 00353 Scalar(277) / Scalar(14336), 00354 Scalar(1) / Scalar(4), 00355 }; 00356 static Scalar c5[6] = { 00357 Scalar(37) / Scalar(378), 00358 Scalar(0), 00359 Scalar(250) / Scalar(621), 00360 Scalar(125) / Scalar(594), 00361 Scalar(0), 00362 Scalar(512) / Scalar(1771) 00363 }; 00364 00365 CPPAD_ASSERT_KNOWN( 00366 M >= 1, 00367 "Error in Runge45: the number of steps is less than one" 00368 ); 00369 CPPAD_ASSERT_KNOWN( 00370 e.size() == xi.size(), 00371 "Error in Runge45: size of e not equal to size of xi" 00372 ); 00373 size_t i, j, k, m; // indices 00374 00375 size_t n = xi.size(); // number of components in X(t) 00376 Scalar ns = Scalar(int(M)); // number of steps as Scalar object 00377 Scalar h = (tf - ti) / ns; // step size 00378 Scalar zero_or_nan = Scalar(0); // zero (nan if Ode returns has a nan) 00379 for(i = 0; i < n; i++) // initialize e = 0 00380 e[i] = zero_or_nan; 00381 00382 // vectors used to store values returned by F 00383 Vector fh(6 * n), xtmp(n), ftmp(n), x4(n), x5(n), xf(n); 00384 00385 xf = xi; // initialize solution 00386 for(m = 0; m < M; m++) 00387 { // time at beginning of this interval 00388 // (convert to int to avoid MS compiler warning) 00389 Scalar t = ti * (Scalar(int(M - m)) / ns) 00390 + tf * (Scalar(int(m)) / ns); 00391 00392 // loop over integration steps 00393 x4 = x5 = xf; // start x4 and x5 at same point for each step 00394 for(j = 0; j < 6; j++) 00395 { // loop over function evaluations for this step 00396 xtmp = xf; // location for next function evaluation 00397 for(k = 0; k < j; k++) 00398 { // loop over previous function evaluations 00399 Scalar bjk = b[ (j-1) * 5 + k ]; 00400 for(i = 0; i < n; i++) 00401 { // loop over elements of x 00402 xtmp[i] += bjk * fh[i * 6 + k]; 00403 } 00404 } 00405 // ftmp = F(t + a[j] * h, xtmp) 00406 F.Ode(t + a[j] * h, xtmp, ftmp); 00407 00408 // if ftmp has a nan, set zero_or_nan to nan 00409 for(i = 0; i < n; i++) 00410 zero_or_nan *= ftmp[i]; 00411 00412 for(i = 0; i < n; i++) 00413 { // loop over elements of x 00414 Scalar fhi = ftmp[i] * h; 00415 fh[i * 6 + j] = fhi; 00416 x4[i] += c4[j] * fhi; 00417 x5[i] += c5[j] * fhi; 00418 x5[i] += zero_or_nan; 00419 } 00420 } 00421 // accumulate error bound 00422 for(i = 0; i < n; i++) 00423 { // cant use abs because cppad.hpp may not be included 00424 Scalar diff = x5[i] - x4[i]; 00425 e[i] += fabs(diff); 00426 e[i] += zero_or_nan; 00427 } 00428 00429 // advance xf for this step using x5 00430 xf = x5; 00431 } 00432 return xf; 00433 } 00434 00435 } // End CppAD namespace 00436 00437 # endif