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