CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_ODE_GEAR_INCLUDED 00003 # define CPPAD_ODE_GEAR_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 OdeGear$$ 00018 $spell 00019 cppad.hpp 00020 Jan 00021 bool 00022 const 00023 CppAD 00024 dep 00025 $$ 00026 00027 $index OdeGear$$ 00028 $index Ode, Gear$$ 00029 $index Gear, Ode$$ 00030 $index stiff, Ode$$ 00031 $index differential, equation$$ 00032 $index equation, differential$$ 00033 00034 $section An Arbitrary Order Gear Method$$ 00035 00036 $head Syntax$$ 00037 $codei%# include <cppad/ode_gear.hpp> 00038 %$$ 00039 $codei%OdeGear(%F%, %m%, %n%, %T%, %X%, %e%)%$$ 00040 00041 00042 $head Purpose$$ 00043 This routine applies 00044 $cref/Gear's Method/OdeGear/Gear's Method/$$ 00045 to solve an explicit set of ordinary differential equations. 00046 We are given 00047 $latex f : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$ be a smooth function. 00048 This routine solves the following initial value problem 00049 $latex \[ 00050 \begin{array}{rcl} 00051 x( t_{m-1} ) & = & x^0 \\ 00052 x^\prime (t) & = & f[t , x(t)] 00053 \end{array} 00054 \] $$ 00055 for the value of $latex x( t_m )$$. 00056 If your set of ordinary differential equations are not stiff 00057 an explicit method may be better (perhaps $cref Runge45$$.) 00058 00059 $head Include$$ 00060 The file $code cppad/ode_gear.hpp$$ is included by $code cppad/cppad.hpp$$ 00061 but it can also be included separately with out the rest of 00062 the $code CppAD$$ routines. 00063 00064 $head Fun$$ 00065 The class $icode Fun$$ 00066 and the object $icode F$$ satisfy the prototype 00067 $codei% 00068 %Fun% &%F% 00069 %$$ 00070 This must support the following set of calls 00071 $codei% 00072 %F%.Ode(%t%, %x%, %f%) 00073 %F%.Ode_dep(%t%, %x%, %f_x%) 00074 %$$ 00075 00076 $subhead t$$ 00077 The argument $icode t$$ has prototype 00078 $codei% 00079 const %Scalar% &%t% 00080 %$$ 00081 (see description of $cref/Scalar/OdeGear/Scalar/$$ below). 00082 00083 $subhead x$$ 00084 The argument $icode x$$ has prototype 00085 $codei% 00086 const %Vector% &%x% 00087 %$$ 00088 and has size $icode n$$ 00089 (see description of $cref/Vector/OdeGear/Vector/$$ below). 00090 00091 $subhead f$$ 00092 The argument $icode f$$ to $icode%F%.Ode%$$ has prototype 00093 $codei% 00094 %Vector% &%f% 00095 %$$ 00096 On input and output, $icode f$$ is a vector of size $icode n$$ 00097 and the input values of the elements of $icode f$$ do not matter. 00098 On output, 00099 $icode f$$ is set equal to $latex f(t, x)$$ 00100 (see $icode f(t, x)$$ in $cref/Purpose/OdeGear/Purpose/$$). 00101 00102 $subhead f_x$$ 00103 The argument $icode f_x$$ has prototype 00104 $codei% 00105 %Vector% &%f_x% 00106 %$$ 00107 On input and output, $icode f_x$$ is a vector of size $latex n * n$$ 00108 and the input values of the elements of $icode f_x$$ do not matter. 00109 On output, 00110 $latex \[ 00111 f\_x [i * n + j] = \partial_{x(j)} f_i ( t , x ) 00112 \] $$ 00113 00114 $subhead Warning$$ 00115 The arguments $icode f$$, and $icode f_x$$ 00116 must have a call by reference in their prototypes; i.e., 00117 do not forget the $code &$$ in the prototype for 00118 $icode f$$ and $icode f_x$$. 00119 00120 $head m$$ 00121 The argument $icode m$$ has prototype 00122 $codei% 00123 size_t %m% 00124 %$$ 00125 It specifies the order (highest power of $latex t$$) 00126 used to represent the function $latex x(t)$$ in the multi-step method. 00127 Upon return from $code OdeGear$$, 00128 the $th i$$ component of the polynomial is defined by 00129 $latex \[ 00130 p_i ( t_j ) = X[ j * n + i ] 00131 \] $$ 00132 for $latex j = 0 , \ldots , m$$ (where $latex 0 \leq i < n$$). 00133 The value of $latex m$$ must be greater than or equal one. 00134 00135 $head n$$ 00136 The argument $icode n$$ has prototype 00137 $codei% 00138 size_t %n% 00139 %$$ 00140 It specifies the range space dimension of the 00141 vector valued function $latex x(t)$$. 00142 00143 $head T$$ 00144 The argument $icode T$$ has prototype 00145 $codei% 00146 const %Vector% &%T% 00147 %$$ 00148 and size greater than or equal to $latex m+1$$. 00149 For $latex j = 0 , \ldots m$$, $latex T[j]$$ is the time 00150 corresponding to time corresponding 00151 to a previous point in the multi-step method. 00152 The value $latex T[m]$$ is the time 00153 of the next point in the multi-step method. 00154 The array $latex T$$ must be monotone increasing; i.e., 00155 $latex T[j] < T[j+1]$$. 00156 Above and below we often use the shorthand $latex t_j$$ for $latex T[j]$$. 00157 00158 00159 $head X$$ 00160 The argument $icode X$$ has the prototype 00161 $codei% 00162 %Vector% &%X% 00163 %$$ 00164 and size greater than or equal to $latex (m+1) * n$$. 00165 On input to $code OdeGear$$, 00166 for $latex j = 0 , \ldots , m-1$$, and 00167 $latex i = 0 , \ldots , n-1$$ 00168 $latex \[ 00169 X[ j * n + i ] = x_i ( t_j ) 00170 \] $$ 00171 Upon return from $code OdeGear$$, 00172 for $latex i = 0 , \ldots , n-1$$ 00173 $latex \[ 00174 X[ m * n + i ] \approx x_i ( t_m ) 00175 \] $$ 00176 00177 $head e$$ 00178 The vector $icode e$$ is an approximate error bound for the result; i.e., 00179 $latex \[ 00180 e[i] \geq | X[ m * n + i ] - x_i ( t_m ) | 00181 \] $$ 00182 The order of this approximation is one less than the order of 00183 the solution; i.e., 00184 $latex \[ 00185 e = O ( h^m ) 00186 \] $$ 00187 where $latex h$$ is the maximum of $latex t_{j+1} - t_j$$. 00188 00189 $head Scalar$$ 00190 The type $icode Scalar$$ must satisfy the conditions 00191 for a $cref NumericType$$ type. 00192 The routine $cref CheckNumericType$$ will generate an error message 00193 if this is not the case. 00194 In addition, the following operations must be defined for 00195 $icode Scalar$$ objects $icode a$$ and $icode b$$: 00196 00197 $table 00198 $bold Operation$$ $cnext $bold Description$$ $rnext 00199 $icode%a% < %b%$$ $cnext 00200 less than operator (returns a $code bool$$ object) 00201 $tend 00202 00203 $head Vector$$ 00204 The type $icode Vector$$ must be a $cref SimpleVector$$ class with 00205 $cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$. 00206 The routine $cref CheckSimpleVector$$ will generate an error message 00207 if this is not the case. 00208 00209 $head Example$$ 00210 $children% 00211 example/ode_gear.cpp 00212 %$$ 00213 The file 00214 $cref ode_gear.cpp$$ 00215 contains an example and test a test of using this routine. 00216 It returns true if it succeeds and false otherwise. 00217 00218 $head Source Code$$ 00219 The source code for this routine is in the file 00220 $code cppad/ode_gear.hpp$$. 00221 00222 $head Theory$$ 00223 For this discussion we use the shorthand $latex x_j$$ 00224 for the value $latex x ( t_j ) \in \B{R}^n$$ which is not to be confused 00225 with $latex x_i (t) \in \B{R}$$ in the notation above. 00226 The interpolating polynomial $latex p(t)$$ is given by 00227 $latex \[ 00228 p(t) = 00229 \sum_{j=0}^m 00230 x_j 00231 \frac{ 00232 \prod_{i \neq j} ( t - t_i ) 00233 }{ 00234 \prod_{i \neq j} ( t_j - t_i ) 00235 } 00236 \] $$ 00237 The derivative $latex p^\prime (t)$$ is given by 00238 $latex \[ 00239 p^\prime (t) = 00240 \sum_{j=0}^m 00241 x_j 00242 \frac{ 00243 \sum_{i \neq j} \prod_{k \neq i,j} ( t - t_k ) 00244 }{ 00245 \prod_{k \neq j} ( t_j - t_k ) 00246 } 00247 \] $$ 00248 Evaluating the derivative at the point $latex t_\ell$$ we have 00249 $latex \[ 00250 \begin{array}{rcl} 00251 p^\prime ( t_\ell ) & = & 00252 x_\ell 00253 \frac{ 00254 \sum_{i \neq \ell} \prod_{k \neq i,\ell} ( t_\ell - t_k ) 00255 }{ 00256 \prod_{k \neq \ell} ( t_\ell - t_k ) 00257 } 00258 + 00259 \sum_{j \neq \ell} 00260 x_j 00261 \frac{ 00262 \sum_{i \neq j} \prod_{k \neq i,j} ( t_\ell - t_k ) 00263 }{ 00264 \prod_{k \neq j} ( t_j - t_k ) 00265 } 00266 \\ 00267 & = & 00268 x_\ell 00269 \sum_{i \neq \ell} 00270 \frac{ 1 }{ t_\ell - t_i } 00271 + 00272 \sum_{j \neq \ell} 00273 x_j 00274 \frac{ 00275 \prod_{k \neq \ell,j} ( t_\ell - t_k ) 00276 }{ 00277 \prod_{k \neq j} ( t_j - t_k ) 00278 } 00279 \\ 00280 & = & 00281 x_\ell 00282 \sum_{k \neq \ell} ( t_\ell - t_k )^{-1} 00283 + 00284 \sum_{j \neq \ell} 00285 x_j 00286 ( t_j - t_\ell )^{-1} 00287 \prod_{k \neq \ell ,j} ( t_\ell - t_k ) / ( t_j - t_k ) 00288 \end{array} 00289 \] $$ 00290 We define the vector $latex \alpha \in \B{R}^{m+1}$$ by 00291 $latex \[ 00292 \alpha_j = \left\{ \begin{array}{ll} 00293 \sum_{k \neq m} ( t_m - t_k )^{-1} 00294 & {\rm if} \; j = m 00295 \\ 00296 ( t_j - t_m )^{-1} 00297 \prod_{k \neq m,j} ( t_m - t_k ) / ( t_j - t_k ) 00298 & {\rm otherwise} 00299 \end{array} \right. 00300 \] $$ 00301 It follows that 00302 $latex \[ 00303 p^\prime ( t_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m 00304 \] $$ 00305 Gear's method determines $latex x_m$$ by solving the following 00306 nonlinear equation 00307 $latex \[ 00308 f( t_m , x_m ) = \alpha_0 x_0 + \cdots + \alpha_m x_m 00309 \] $$ 00310 Newton's method for solving this equation determines iterates, 00311 which we denote by $latex x_m^k$$, by solving the following affine 00312 approximation of the equation above 00313 $latex \[ 00314 \begin{array}{rcl} 00315 f( t_m , x_m^{k-1} ) + \partial_x f( t_m , x_m^{k-1} ) ( x_m^k - x_m^{k-1} ) 00316 & = & 00317 \alpha_0 x_0^k + \alpha_1 x_1 + \cdots + \alpha_m x_m 00318 \\ 00319 \left[ \alpha_m I - \partial_x f( t_m , x_m^{k-1} ) \right] x_m 00320 & = & 00321 \left[ 00322 f( t_m , x_m^{k-1} ) - \partial_x f( t_m , x_m^{k-1} ) x_m^{k-1} 00323 - \alpha_0 x_0 - \cdots - \alpha_{m-1} x_{m-1} 00324 \right] 00325 \end{array} 00326 \] $$ 00327 In order to initialize Newton's method; i.e. choose $latex x_m^0$$ 00328 we define the vector $latex \beta \in \B{R}^{m+1}$$ by 00329 $latex \[ 00330 \beta_j = \left\{ \begin{array}{ll} 00331 \sum_{k \neq m-1} ( t_{m-1} - t_k )^{-1} 00332 & {\rm if} \; j = m-1 00333 \\ 00334 ( t_j - t_{m-1} )^{-1} 00335 \prod_{k \neq m-1,j} ( t_{m-1} - t_k ) / ( t_j - t_k ) 00336 & {\rm otherwise} 00337 \end{array} \right. 00338 \] $$ 00339 It follows that 00340 $latex \[ 00341 p^\prime ( t_{m-1} ) = \beta_0 x_0 + \cdots + \beta_m x_m 00342 \] $$ 00343 We solve the following approximation of the equation above to determine 00344 $latex x_m^0$$: 00345 $latex \[ 00346 f( t_{m-1} , x_{m-1} ) = 00347 \beta_0 x_0 + \cdots + \beta_{m-1} x_{m-1} + \beta_m x_m^0 00348 \] $$ 00349 00350 00351 $head Gear's Method$$ 00352 C. W. Gear, 00353 ``Simultaneous Numerical Solution of Differential-Algebraic Equations,'' 00354 IEEE Transactions on Circuit Theory, 00355 vol. 18, no. 1, pp. 89-95, Jan. 1971. 00356 00357 00358 $end 00359 -------------------------------------------------------------------------- 00360 */ 00361 00362 # include <cstddef> 00363 # include <cppad/local/cppad_assert.hpp> 00364 # include <cppad/check_simple_vector.hpp> 00365 # include <cppad/check_numeric_type.hpp> 00366 # include <cppad/vector.hpp> 00367 # include <cppad/lu_factor.hpp> 00368 # include <cppad/lu_invert.hpp> 00369 00370 namespace CppAD { // BEGIN CppAD namespace 00371 00372 template <typename Vector, typename Fun> 00373 void OdeGear( 00374 Fun &F , 00375 size_t m , 00376 size_t n , 00377 const Vector &T , 00378 Vector &X , 00379 Vector &e ) 00380 { 00381 // temporary indices 00382 size_t i, j, k; 00383 00384 typedef typename Vector::value_type Scalar; 00385 00386 // check numeric type specifications 00387 CheckNumericType<Scalar>(); 00388 00389 // check simple vector class specifications 00390 CheckSimpleVector<Scalar, Vector>(); 00391 00392 CPPAD_ASSERT_KNOWN( 00393 m >= 1, 00394 "OdeGear: m is less than one" 00395 ); 00396 CPPAD_ASSERT_KNOWN( 00397 n > 0, 00398 "OdeGear: n is equal to zero" 00399 ); 00400 CPPAD_ASSERT_KNOWN( 00401 size_t(T.size()) >= (m+1), 00402 "OdeGear: size of T is not greater than or equal (m+1)" 00403 ); 00404 CPPAD_ASSERT_KNOWN( 00405 size_t(X.size()) >= (m+1) * n, 00406 "OdeGear: size of X is not greater than or equal (m+1) * n" 00407 ); 00408 for(j = 0; j < m; j++) CPPAD_ASSERT_KNOWN( 00409 T[j] < T[j+1], 00410 "OdeGear: the array T is not monotone increasing" 00411 ); 00412 00413 // some constants 00414 Scalar zero(0); 00415 Scalar one(1); 00416 00417 // vectors required by method 00418 Vector alpha(m + 1); 00419 Vector beta(m + 1); 00420 Vector f(n); 00421 Vector f_x(n * n); 00422 Vector x_m0(n); 00423 Vector x_m(n); 00424 Vector b(n); 00425 Vector A(n * n); 00426 00427 // compute alpha[m] 00428 alpha[m] = zero; 00429 for(k = 0; k < m; k++) 00430 alpha[m] += one / (T[m] - T[k]); 00431 00432 // compute beta[m-1] 00433 beta[m-1] = one / (T[m-1] - T[m]); 00434 for(k = 0; k < m-1; k++) 00435 beta[m-1] += one / (T[m-1] - T[k]); 00436 00437 00438 // compute other components of alpha 00439 for(j = 0; j < m; j++) 00440 { // compute alpha[j] 00441 alpha[j] = one / (T[j] - T[m]); 00442 for(k = 0; k < m; k++) 00443 { if( k != j ) 00444 { alpha[j] *= (T[m] - T[k]); 00445 alpha[j] /= (T[j] - T[k]); 00446 } 00447 } 00448 } 00449 00450 // compute other components of beta 00451 for(j = 0; j <= m; j++) 00452 { if( j != m-1 ) 00453 { // compute beta[j] 00454 beta[j] = one / (T[j] - T[m-1]); 00455 for(k = 0; k <= m; k++) 00456 { if( k != j && k != m-1 ) 00457 { beta[j] *= (T[m-1] - T[k]); 00458 beta[j] /= (T[j] - T[k]); 00459 } 00460 } 00461 } 00462 } 00463 00464 // evaluate f(T[m-1], x_{m-1} ) 00465 for(i = 0; i < n; i++) 00466 x_m[i] = X[(m-1) * n + i]; 00467 F.Ode(T[m-1], x_m, f); 00468 00469 // solve for x_m^0 00470 for(i = 0; i < n; i++) 00471 { x_m[i] = f[i]; 00472 for(j = 0; j < m; j++) 00473 x_m[i] -= beta[j] * X[j * n + i]; 00474 x_m[i] /= beta[m]; 00475 } 00476 x_m0 = x_m; 00477 00478 // evaluate partial w.r.t x of f(T[m], x_m^0) 00479 F.Ode_dep(T[m], x_m, f_x); 00480 00481 // compute the matrix A = ( alpha[m] * I - f_x ) 00482 for(i = 0; i < n; i++) 00483 { for(j = 0; j < n; j++) 00484 A[i * n + j] = - f_x[i * n + j]; 00485 A[i * n + i] += alpha[m]; 00486 } 00487 00488 // LU factor (and overwrite) the matrix A 00489 int sign; 00490 CppAD::vector<size_t> ip(n) , jp(n); 00491 sign = LuFactor(ip, jp, A); 00492 CPPAD_ASSERT_KNOWN( 00493 sign != 0, 00494 "OdeGear: step size is to large" 00495 ); 00496 00497 // Iterations of Newton's method 00498 for(k = 0; k < 3; k++) 00499 { 00500 // only evaluate f( T[m] , x_m ) keep f_x during iteration 00501 F.Ode(T[m], x_m, f); 00502 00503 // b = f + f_x x_m - alpha[0] x_0 - ... - alpha[m-1] x_{m-1} 00504 for(i = 0; i < n; i++) 00505 { b[i] = f[i]; 00506 for(j = 0; j < n; j++) 00507 b[i] -= f_x[i * n + j] * x_m[j]; 00508 for(j = 0; j < m; j++) 00509 b[i] -= alpha[j] * X[ j * n + i ]; 00510 } 00511 LuInvert(ip, jp, A, b); 00512 x_m = b; 00513 } 00514 00515 // return estimate for x( t[k] ) and the estimated error bound 00516 for(i = 0; i < n; i++) 00517 { X[m * n + i] = x_m[i]; 00518 e[i] = x_m[i] - x_m0[i]; 00519 if( e[i] < zero ) 00520 e[i] = - e[i]; 00521 } 00522 } 00523 00524 } // End CppAD namespace 00525 00526 # endif