Code-Eli  0.3.6
simpson.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef eli_mutil_quad_simpson_hpp
14 #define eli_mutil_quad_simpson_hpp
15 
16 #include <limits>
17 
18 #include "eli/code_eli.hpp"
19 
20 namespace eli
21 {
22  namespace mutil
23  {
24  namespace quad
25  {
26  template<typename data__>
27  class simpson
28  {
29  public:
31  {
34 
36  {
37  tolerance=std::sqrt(std::numeric_limits<data__>::epsilon());
38  recursion_depth=0;
39  function_count=0;
40  max_depth=30;
41  tol_factor=1.25;
42  error_factor=100;
43  coarse_value=std::numeric_limits<data__>::quiet_NaN();
44  fine_value=std::numeric_limits<data__>::quiet_NaN();
45  approximate_error=std::numeric_limits<data__>::quiet_NaN();
46  }
47  };
48 
49  public:
50  simpson() {};
51  simpson(const simpson<data__> &) {};
52  ~simpson() {};
53 
54  simpson<data__> & operator=(const simpson<data__> &) {return *this;};
55 
56  size_t order() const {return 4;}
57 
68  template<typename yit__>
69  data__ operator()(const data__ &dx, yit__ yb, yit__ ye) const
70  {
71  yit__ ym1, y(yb);
72 
73  // short circuit special cases
74  if (yb==ye)
75  return static_cast<data__>(0);
76  ++y;
77  if (y==ye)
78  return (*yb)/dx;
79  ym1=y;
80  ++y;
81  if (y==ye)
82  return ((*ym1)-(*yb))/dx;
83 
84  // cycle through the points until last odd point
85  data__ rtnval((*yb));
86  bool even_pts(true);
87 
88  while (y!=ye)
89  {
90  rtnval+=4*(*ym1)+2*(*y);
91  ym1=y;
92  ++y;
93  if (y==ye)
94  {
95  // remove extra last odd point
96  rtnval-=(*ym1);
97  even_pts=false;
98  break;
99  }
100  ym1=y;
101  ++y;
102  }
103 
104  if (even_pts) // have even number of points
105  {
106  yit__ ym2;
107 
108  // (re)set the iterators
109  --y;
110  --ym1;
111  ym2=ym1; --ym2;
112 
113  // remove extra last odd point
114  rtnval-=(*ym1);
115 
116  // add the last even segment
117  rtnval+=5*(*y)/4+2*(*ym1)-(*ym2)/4;
118  }
119 
120  // multiply by factor
121  rtnval*=dx/static_cast<data__>(3);
122 
123  return rtnval;
124  }
125 
137  template<typename xit__, typename yit__>
138  data__ operator()(xit__ x, yit__ yb, yit__ ye) const
139  {
140  data__ rtnval(0), delta;
141  xit__ xm1(x), xp1;
142  yit__ ym1(yb), y(yb), yp1(yb);
143 
144  // set the initial iterators
145  ++x; xp1=x;
146  ++y; yp1=y;
147 
148  // short circuit special cases
149  if (yb==ye)
150  return static_cast<data__>(0);
151  if (y==ye)
152  return static_cast<data__>(0);
153  ++xp1;
154  ++yp1;
155  if (yp1==ye)
156  return ((*y)-(*yb))/((*x)-(*xm1));
157 
158  // cycle through the points
159  bool even_pts(true);
160  for (; (yp1!=ye) && (y!=ye); ++xm1, ++x, ++xp1, ++ym1, ++y, ++yp1)
161  {
162  delta=((*x)-(*xm1))/((*xp1)-(*x));
163  rtnval+=(((*xp1)-(*x))/6)*(1+1/delta)*((delta*(2-delta))*(*yp1)+(1+delta)*(1+delta)*(*y)+(2*delta-1)*(*ym1));
164 
165  ++xm1; ++x; ++xp1;
166  ++ym1; ++y; ++yp1;
167  if (yp1==ye)
168  {
169  even_pts = false;
170  break;
171  }
172  }
173 
174  // test if even number of point
175  if (even_pts)
176  {
177  // reset the iterators for last segment
178  --xm1; --x; --xp1;
179  --ym1; --y; --yp1;
180 
181  // add last segment
182  delta=((*x)-(*xm1))/((*xp1)-(*x));
183  rtnval+=(((*xp1)-(*x))/6)*(((2+3*delta)/(1+delta))*(*yp1)+(3+1/delta)*(*y)-(1/(delta*(1+delta)))*(*ym1));
184  }
185 
186  return rtnval;
187  }
188 
199  template<typename f__>
200  data__ operator()(f__ fun, const data__ &x0, const data__ &x1) const
201  {
203  return operator()(fun, x0, x1, ap);
204  }
205 
218  template<typename f__>
219  data__ operator()(f__ fun, const data__ &x0, const data__ &x1, const data__ &tol, const size_t &max_depth) const
220  {
222 
223  // set adaption parameters
224  ap.tolerance=tol;
225  ap.max_depth=max_depth;
226 
227  // do the real work
228  return operator()(fun, x0, x1, ap);
229  }
230 
242  template<typename f__>
243  data__ operator()(f__ fun, const data__ &x0, const data__ &x1, typename simpson<data__>::adaptive_params &ap) const
244  {
245  data__ xc[3], fc[3];
246 
247  // create the first level evaluation of integral
248  xc[0]=x0;
249  fc[0]=fun(xc[0]);
250  xc[1]=(x0+x1)/2;
251  fc[1]=fun(xc[1]);
252  xc[2]=x1;
253  fc[2]=fun(xc[2]);
254  ap.coarse_value=simple(x0, x1, fc);
255  ap.recursion_depth=0;
256  ap.function_count=3;
257 
258  // recursively subdivide the interval to get refined integral value
259  if (ap.max_depth>0)
260  {
261  internal_recurse(fun, xc, fc, ap);
262  }
263  else
264  {
265  ap.fine_value=ap.coarse_value;
266  return ap.fine_value;
267  }
268 
269  // extrapolate final value of integral
270  return ((1<<this->order())*ap.fine_value-ap.coarse_value)/((1<<this->order())-1);
271  }
272 
273  private:
274  data__ simple(const data__ &a, const data__ &b, const data__ f[]) const
275  {
276  return((b-a)/6)*(f[0]+4*f[1]+f[2]);
277  }
278 
279  template<typename f__>
280  void internal_recurse(f__ &fun, const data__ xc[], const data__ fc[], typename simpson<data__>::adaptive_params &ap) const
281  {
282  data__ xf[5], ff[5], sf1, sf2;
283 
284  ++ap.recursion_depth;
285  // copy over and set the fine values
286  // TODO: Could generalize the recursive integration if could abstract away
287  // this portion of code. Could do this with a base class that had a
288  // template parameter that is the number of points in the stencil.
289  // That way would not have to dynamically allocate any memory.
290  xf[0]=xc[0];
291  xf[1]=(xc[0]+xc[1])/2;
292  xf[2]=xc[1];
293  xf[3]=(xc[1]+xc[2])/2;
294  xf[4]=xc[2];
295  ff[0]=fc[0];
296  ff[1]=fun(xf[1]);
297  ff[2]=fc[1];
298  ff[3]=fun(xf[3]);
299  ff[4]=fc[2];
300  ap.function_count+=2;
301 
302  // calculate the simpson rules
303  sf1=simple(xf[0], xf[2], &(ff[0]));
304  sf2=simple(xf[2], xf[4], &(ff[2]));
305  ap.fine_value=sf1+sf2;
306 
307  // check to see if error tolerance met
309  if ( (ap.approximate_error>ap.tolerance) && (ap.recursion_depth<ap.max_depth) )
310  {
311  // set up the adaption parameters based on the ones passed in
312  typename simpson<data__>::adaptive_params ap1(ap), ap2(ap);
313  ap1.function_count=0;
314  ap1.tolerance/=ap1.tol_factor;
315  ap1.coarse_value=sf1;
316  ap2.function_count=0;
317  ap2.tolerance/=ap2.tol_factor;
318  ap2.coarse_value=sf2;
319 
320  // evaluate both sides of the interval
321  internal_recurse(fun, &(xf[0]), &(ff[0]), ap1);
322  internal_recurse(fun, &(xf[2]), &(ff[2]), ap2);
323 
324  // recover the adaption information
325  ap.function_count+=ap1.function_count+ap2.function_count;
326  ap.fine_value=ap1.fine_value+ap2.fine_value;
327  ap.coarse_value=ap1.coarse_value+ap2.coarse_value;
328  ap.recursion_depth=std::max(ap1.recursion_depth, ap2.recursion_depth);
329  ap.approximate_error=ap1.approximate_error+ap2.approximate_error;
330  }
331  }
332  };
333  }
334  }
335 }
336 #endif
data__ operator()(f__ fun, const data__ &x0, const data__ &x1) const
Definition: simpson.hpp:200
Definition: math.hpp:20
data__ operator()(xit__ x, yit__ yb, yit__ ye) const
Definition: simpson.hpp:138
data__ coarse_value
Definition: simpson.hpp:33
simpson(const simpson< data__ > &)
Definition: simpson.hpp:51
size_t order() const
Definition: simpson.hpp:56
simpson()
Definition: simpson.hpp:50
data__ tolerance
Definition: simpson.hpp:33
adaptive_params()
Definition: simpson.hpp:35
size_t function_count
Definition: simpson.hpp:32
data__ tol_factor
Definition: simpson.hpp:33
data__ operator()(f__ fun, const data__ &x0, const data__ &x1, typename simpson< data__ >::adaptive_params &ap) const
Definition: simpson.hpp:243
void internal_recurse(f__ &fun, const data__ xc[], const data__ fc[], typename simpson< data__ >::adaptive_params &ap) const
Definition: simpson.hpp:280
data__ simple(const data__ &a, const data__ &b, const data__ f[]) const
Definition: simpson.hpp:274
data__ operator()(f__ fun, const data__ &x0, const data__ &x1, const data__ &tol, const size_t &max_depth) const
Definition: simpson.hpp:219
size_t recursion_depth
Definition: simpson.hpp:32
data__ operator()(const data__ &dx, yit__ yb, yit__ ye) const
Definition: simpson.hpp:69
~simpson()
Definition: simpson.hpp:52
simpson< data__ > & operator=(const simpson< data__ > &)
Definition: simpson.hpp:54
data__ approximate_error
Definition: simpson.hpp:33
data__ error_factor
Definition: simpson.hpp:33
size_t max_depth
Definition: simpson.hpp:32
Definition: simpson.hpp:27
data__ fine_value
Definition: simpson.hpp:33