CppAD: A C++ Algorithmic Differentiation Package  20171217
ode_evaluate.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15 /*
16 $begin ode_evaluate$$17 spell 18 Runge 19 fabs 20 retaped 21 Jacobian 22 const 23 Cpp 24 cppad 25 hpp 26 fp 27 namespace 28 exp 29$$ 30 31$section Evaluate a Function Defined in Terms of an ODE$$32 mindex ode_evaluate$$
33
34
35 $head Syntax$$36 codei%# include <cppad/speed/ode_evaluate.hpp> 37 %$$ 38$codei%ode_evaluate(%x%, %p%, %fp%)%$$39 40 head Purpose$$
41 This routine evaluates a function $latex f : \B{R}^n \rightarrow \B{R}^n$$42 defined by 43 latex $44 f(x) = y(x, 1) 45$$$ 46 where$latex y(x, t)$$solves the ordinary differential equation 47 latex $48 \begin{array}{rcl} 49 y(x, 0) & = & x 50 \\ 51 \partial_t y (x, t ) & = & g[ y(x,t) , t ] 52 \end{array} 53$$$
54 where $latex g : \B{R}^n \times \B{R} \rightarrow \B{R}^n$$55 is an unspecified function. 56 57 head Inclusion$$ 58 The template function$code ode_evaluate$$59 is defined in the code CppAD$$ namespace by including
60 the file $code cppad/speed/ode_evaluate.hpp$$61 (relative to the CppAD distribution directory). 62 63 head Float$$ 64 65$subhead Operation Sequence$$66 The type icode Float$$ must be a $cref NumericType$$. 67 The icode Float$$ 68$cref/operation sequence/glossary/Operation/Sequence/$$69 for this routine does not depend on the value of the argument icode x$$,
70 hence it does not need to be retaped for each value of $latex x$$. 71 72 subhead fabs$$ 73 If$icode y$$and icode z$$ are $icode Float$$objects, the syntax 74 codei% 75 %y% = fabs(%z%) 76 %$$ 77 must be supported. Note that it does not matter if the operation 78 sequence for$code fabs$$depends on icode z$$ because the
79 corresponding results are not actually used by $code ode_evaluate$$; 80 see code fabs$$ in$cref/Runge45/Runge45/Scalar/fabs/$$. 81 82 head x$$
83 The argument $icode x$$has prototype 84 codei% 85 const CppAD::vector<%Float%>& %x% 86 %$$ 87 It contains he argument value for which the function, 88 or its derivative, is being evaluated. 89 The value$latex n$$is determined by the size of the vector icode x$$.
90
91 $head p$$92 The argument icode p$$ has prototype 93$codei%
94  size_t %p%
95 %$$96 97 subhead p == 0$$
98 In this case a numerical method is used to solve the ode
99 and obtain an accurate approximation for $latex y(x, 1)$$. 100 This numerical method has a fixed 101 that does not depend on icode x$$. 102 103$subhead p = 1$$104 In this case an analytic solution for the partial derivative 105 latex \partial_x y(x, 1)$$ is returned.
106
107 $head fp$$108 The argument icode fp$$ has prototype 109$codei%
111 %$$112 The input value of the elements of icode fp$$ does not matter.
113
114 $subhead Function$$115 If icode p$$ is zero,$icode fp$$has size equal to latex n$$
116 and contains the value of $latex y(x, 1)$$. 117 118 subhead Gradient$$ 119 If$icode p$$is one, icode fp$$ has size equal to $icode n^2$$120 and for latex i = 0 , \ldots 1$$,$latex j = 0 , \ldots , n-1$$121 latex $122 \D{y[i]}{x[j]} (x, 1) = fp [ i \cdot n + j ] 123$$$
124
125 $children% 126 speed/example/ode_evaluate.cpp% 127 omh/ode_evaluate.omh 128 %$$129 130 head Example$$ 131 The file 132$cref ode_evaluate.cpp$$133 contains an example and test of code ode_evaluate.hpp$$.
134 It returns true if it succeeds and false otherwise.
135
136
137 $head Source Code$$138 The file 139 cref ode_evaluate.hpp$$ 140 contains the source code for this template function. 141 142$end
143 */
144 // BEGIN C++
148
150
151  template <class Float>
153  public:
154  // Given that y_i (0) = x_i,
155  // the following y_i (t) satisfy the ODE below:
156  // y_0 (t) = x[0]
157  // y_1 (t) = x[1] + x[0] * t
158  // y_2 (t) = x[2] + x[1] * t + x[0] * t^2/2
159  // y_3 (t) = x[3] + x[2] * t + x[1] * t^2/2 + x[0] * t^3 / 3!
160  // ...
161  void Ode(
162  const Float& t,
165  { size_t n = y.size();
166  f[0] = 0.;
167  for(size_t k = 1; k < n; k++)
168  f[k] = y[k-1];
169  }
170  };
171  //
172  template <class Float>
175  size_t p ,
178  typedef vector<Float> VectorFloat;
179
180  size_t n = x.size();
181  CPPAD_ASSERT_KNOWN( p == 0 || p == 1,
182  "ode_evaluate: p is not zero or one"
183  );
185  ((p==0) & (fp.size()==n)) || ((p==1) & (fp.size()==n*n)),
186  "ode_evaluate: the size of fp is not correct"
187  );
188  if( p == 0 )
189  { // function that defines the ode
191
192  // number of Runge45 steps to use
193  size_t M = 10;
194
195  // initial and final time
196  Float ti = 0.0;
197  Float tf = 1.0;
198
199  // initial value for y(x, t); i.e. y(x, 0)
200  // (is a reference to x)
201  const VectorFloat& yi = x;
202
203  // final value for y(x, t); i.e., y(x, 1)
204  // (is a reference to fp)
205  VectorFloat& yf = fp;
206
207  // Use fourth order Runge-Kutta to solve ODE
208  yf = CppAD::Runge45(F, M, ti, tf, yi);
209
210  return;
211  }
212  /* Compute derivaitve of y(x, 1) w.r.t x
213  y_0 (x, t) = x[0]
214  y_1 (x, t) = x[1] + x[0] * t
215  y_2 (x, t) = x[2] + x[1] * t + x[0] * t^2/2
216  y_3 (x, t) = x[3] + x[2] * t + x[1] * t^2/2 + x[0] * t^3 / 3!
217  ...
218  */
219  size_t i, j, k;
220  for(i = 0; i < n; i++)
221  { for(j = 0; j < n; j++)
222  fp[ i * n + j ] = 0.0;
223  }
224  size_t factorial = 1;
225  for(k = 0; k < n; k++)
226  { if( k > 1 )
227  factorial *= k;
228  for(i = k; i < n; i++)
229  { // partial w.r.t x[i-k] of x[i-k] * t^k / k!
230  j = i - k;
231  fp[ i * n + j ] += 1.0 / Float(factorial);
232  }
233  }
234  }
235 }
236 // END C++
237
238 # endif