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