CppAD: A C++ Algorithmic Differentiation Package  20171217
romberg_one.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 $begin RombergOne$$16 spell 17 cppad.hpp 18 bool 19 const 20 Cpp 21 RombergOne 22$$ 23 24$section One DimensionalRomberg Integration$$25 mindex integrate Romberg$$
26
27
28 $head Syntax$$29 codei%# include <cppad/utility/romberg_one.hpp> 30 %$$ 31$icode%r% = RombergOne(%F%, %a%, %b%, %n%, %e%)%$$32 33 34 head Description$$
35 Returns the Romberg integration estimate
36 $latex r$$for a one dimensional integral 37 latex $38 r = \int_a^b F(x) {\bf d} x + O \left[ (b - a) / 2^{n-1} \right]^{2(p+1)} 39$$$ 40 41$head Include$$42 The file code cppad/romberg_one.hpp$$ is included by $code cppad/cppad.hpp$$43 but it can also be included separately with out the rest of 44 the code CppAD$$ routines. 45 46$head r$$47 The return value icode r$$ has prototype
48 $codei% 49 %Float% %r% 50 %$$51 It is the estimate computed by code RombergOne$$ for the integral above. 52 53$head F$$54 The object icode F$$ can be of any type, but it must support
55 the operation
56 $codei% 57 %F%(%x%) 58 %$$59 The argument icode x$$ to$icode F$$has prototype 60 codei% 61 const %Float% &%x% 62 %$$
63 The return value of $icode F$$is a icode Float$$ object 64 (see description of$cref/Float/RombergOne/Float/$$below). 65 66 head a$$
67 The argument $icode a$$has prototype 68 codei% 69 const %Float% &%a% 70 %$$ 71 It specifies the lower limit for the integration. 72 73$head b$$74 The argument icode b$$ has prototype
75 $codei% 76 const %Float% &%b% 77 %$$78 It specifies the upper limit for the integration. 79 80 head n$$ 81 The argument$icode n$$has prototype 82 codei% 83 size_t %n% 84 %$$
85 A total number of $latex 2^{n-1} + 1$$evaluations of icode%F%(%x%)%$$ 86 are used to estimate the integral. 87 88$head p$$89 The argument icode p$$ has prototype
90 $codei% 91 size_t %p% 92 %$$93 It must be less than or equal latex n$$ 94 and determines the accuracy order in the approximation for the integral 95 that is returned by$code RombergOne$$. 96 To be specific 97 latex $98 r = \int_a^b F(x) {\bf d} x + O \left[ (b - a) / 2^{n-1} \right]^{2(p+1)} 99$$$
100
101
102 $head e$$103 The argument icode e$$ has prototype 104$codei%
105  %Float% &%e%
106 %$$107 The input value of icode e$$ does not matter
108 and its output value is an approximation for the error in
109 the integral estimates; i.e.,
110 $latex $111 e \approx \left| r - \int_a^b F(x) {\bf d} x \right| 112$ $$113 114 head Float$$ 115 The type$icode Float$$must satisfy the conditions 116 for a cref NumericType$$ type.
117 The routine $cref CheckNumericType$$will generate an error message 118 if this is not the case. 119 In addition, if icode x$$ and$icode y$$are icode Float$$ objects,
120 $codei% 121 %x% < %y% 122 %$$123 returns the code bool$$ value true if$icode x$$is less than 124 icode y$$ and false otherwise.
125
126 $children% 127 example/utility/romberg_one.cpp 128 %$$129 head Example$$ 130$comment%
131  example/utility/romberg_one.cpp
132 %$$133 The file 134 cref romberg_one.cpp$$
135 contains an example and test a test of using this routine.
136 It returns true if it succeeds and false otherwise.
137
138 $head Source Code$$139 The source code for this routine is in the file 140 code cppad/romberg_one.hpp$$. 141 142$end
143 */
144
148
150
151 template <class Fun, class Float>
153  Fun &F ,
154  const Float &a ,
155  const Float &b ,
156  size_t n ,
157  size_t p ,
158  Float &e )
159 {
160  size_t ipow2 = 1;
161  size_t k, i;
162  Float pow2, sum, x;
163
164  Float zero = Float(0);
165  Float two = Float(2);
166
167  // check specifications for a NumericType
168  CheckNumericType<Float>();
169
171  n >= 2,
172  "RombergOne: n must be greater than or equal 2"
173  );
175
176  // set r[i] = trapazoidal rule with 2^i intervals in [a, b]
177  r[0] = ( F(a) + F(b) ) * (b - a) / two;
178  for(i = 1; i < n; i++)
179  { ipow2 *= 2;
180  // there must be a conversion from int to any numeric type
181  pow2 = Float(int(ipow2));
182  sum = zero;
183  for(k = 1; k < ipow2; k += 2)
184  { // start = a + (b-a)/pow2, increment = 2*(b-a)/pow2
185  x = ( (pow2 - Float(double(k))) * a + double(k) * b ) / pow2;
186  sum = sum + F(x);
187  }
188  // combine function evaluations in sum with those in T[i-1]
189  r[i] = r[i-1] / two + sum * (b - a) / pow2;
190  }
191
192  // now compute the higher order estimates
193  size_t ipow4 = 1; // order of accuract for previous estimate
194  Float pow4, pow4minus;
195  for(i = 0; i < p; i++)
196  { // compute estimate accurate to O[ step^(2*(i+1)) ]
197  // put resutls in r[n-1], r[n-2], ... , r[n-i+1]
198  ipow4 *= 4;
199  pow4 = Float(int(ipow4));
200  pow4minus = Float(ipow4-1);
201  for(k = n-1; k > i; k--)
202  r[k] = ( pow4 * r[k] - r[k-1] ) / pow4minus;
203  }
204
205  // error estimate for r[n]
206  e = r[n-1] - r[n-2];
207  if( e < zero )
208  e = - e;
209  return r[n-1];
210 }
211
212 } // END CppAD namespace
213
214 # endif