deriv_eqi.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_DERIV_EQI_H
24 #define O2SCL_DERIV_EQI_H
25 
26 /** \file deriv_eqi.h
27  \brief File defining \ref o2scl::deriv_eqi
28 */
29 
30 #include <cmath>
31 
32 #include <boost/numeric/ublas/vector.hpp>
33 
34 #include <o2scl/deriv.h>
35 #include <o2scl/err_hnd.h>
36 
37 #ifndef DOXYGEN_NO_O2NS
38 namespace o2scl {
39 #endif
40 
41  /** \brief Derivatives for equally-spaced abscissas
42 
43  This is an implementation of the formulas for equally-spaced
44  abscissas as indicated below. The level of approximation is
45  specified in set_npoints(). The value of \f$ p \times h \f$
46  can be specified in \c xoff (default is zero).
47 
48  \note Uncertainties are not computed and the code
49  for second and third derivatives is unfinished.
50 
51  \note The derivatives given, for example, from the
52  five-point formula can sometimes be more accurate
53  than computing the derivative from the interpolation class.
54  This is especially true near the boundaries of the interpolated
55  region.
56 
57  \future Finish the second and third derivative formulas.
58 
59  Two-point formula (note that this is independent of p).
60  \f[
61  f^{\prime}(x_0+p h)=\frac{1}{h}\left[
62  f_{1}-f_{0} \right]
63  \f]
64  Three-point formula from Abramowitz and Stegun
65  \f[
66  f^{\prime}(x_0+p h)=\frac{1}{h}\left[
67  \frac{2p-1}{2}f_{-1}-2 p f_{0}+\frac{2p+1}{2}f_{1}\right]
68  \f]
69  Four-point formula from Abramowitz and Stegun
70  \f[
71  f^{\prime}(x_0+p h)=\frac{1}{h}\left[
72  -\frac{3 p^2-6 p+2}{6}f_{-1}
73  +\frac{3 p^2-4 p -1}{2}f_{0}
74  -\frac{3 p^2-2 p-2}{2}f_{1}
75  +\frac{3 p^2-1}{6}f_{2}
76  \right]
77  \f]
78  Five-point formula from Abramowitz and Stegun
79  \f{eqnarray*}
80  f^{\prime}(x_0+p h)&=&\frac{1}{h}\left[
81  \frac{2 p^3-3 p^2-p+1}{12}f_{-2}
82  -\frac{4 p^3-3p^2-8p+4}{6}f_{-1}
83  \right. \\ && \left.
84  +\frac{2p^3-5p}{2}f_{0}
85  -\frac{4p^3+3p^2-8p-4}{6}f_{1}
86  \right. \\ && \left.
87  +\frac{2p^3+3p^2-p-1}{12}f_{2}
88  \right]
89  \f}
90 
91  The relations above can be confined to give formulas
92  for second derivative formulas:
93  Three-point formula
94  \f[
95  f^{\prime}(x_0+p h)=\frac{1}{h^2}
96  \left[f_{-1}-2 f_0+f_1\right]
97  \f]
98  Four-point formula:
99  \f[
100  f^{\prime}(x_0+p h)=\frac{1}{2 h^2}
101  \left[\left(1-2p\right)f_{-1}-\left(1-6p\right)f_0
102  -\left(1+6p\right)f_1+\left(1+2p\right)f_2\right]
103  \f]
104  Five-point formula:
105  \f[
106  f^{\prime}(x_0+p h)=\frac{1}{4 h^2}
107  \left[\left(1-2p\right)^2f_{-2}
108  +\left(8p-16 p^2\right)f_{-1}
109  -\left(2-24 p^2\right)f_0
110  -\left(8p+16p^2\right)f_1
111  +\left(1+2p\right)^2 f_2\right]
112  \f]
113  Six-point formula:
114  \f{eqnarray*}
115  f^{\prime}(x_0+p h)&=&\frac{1}{12 h^2}\left[
116  \left(2-10p+15 p^2-6p^3\right)f_{-2}
117  +\left(3+14p-57p^2+30p^3\right)f_{-1}
118  \right. \\ && \left.
119  +\left(-8+20p+78 p^2-60p^3\right)f_0
120  +\left(-2-44p-42p^2+60p^3\right)f_1
121  \right. \\ && \left.
122  +\left(6+22p+3p^2-30p^3\right)f_2
123  +\left(-1-2p+3p^2+6p^3\right)f_3
124  \right]
125  \f}
126  Seven-point formula:
127  \f{eqnarray*}
128  f^{\prime}(x_0+p h)&=&\frac{1}{36 h^2}\left[
129  \left(4-24p+48p^2-36p^3+9p^4\right)f_{-3}
130  +\left(12+12p-162p^2+180p^3-54p^4\right)f_{-2}
131  \right. \\ && \left.
132  +\left(-15+120p+162p^2-360p^3+135p^4\right)f_{-1}
133  -4\left(8+48p-3p^2-90p^3+45p^4\right)f_0
134  \right. \\ && \left.
135  +3\left(14+32p-36p^2-60p^3+45p^4\right)f_1
136  +\left(-12-12p+54p^2+36p^3-54p^4\right)f_2
137  \right. \\ && \left.
138  +\left(1-6p^2+9p^4\right)f_3
139  \right]
140  \f}
141  */
142  template<class func_t=funct11,
144  class deriv_eqi : public deriv_base<func_t> {
145  public:
146 
147  deriv_eqi() {
148  h=1.0e-4;
149  xoff=0.0;
153  }
154 
155  /// Stepsize (Default \f$ 10^{-4} \f$ ).
156  double h;
157 
158  /// Offset (default 0.0)
159  double xoff;
160 
161  /** \brief Set the number of points to use for first derivatives
162  (default 5)
163 
164  Acceptable values are 2-5 (see above).
165  */
166  int set_npoints(int npoints) {
167  if (npoints==2) {
170  } else if (npoints==3) {
173  } else if (npoints==4) {
176  } else {
179  }
180  if (npoints<=1 || npoints>5) {
181  O2SCL_ERR("Invalid # of points in set_npoints(). Using default",
182  exc_einval);
183  }
184  return 0;
185  }
186 
187  /** \brief Set the number of points to use for second derivatives
188  (default 5)
189 
190  Acceptable values are 3-5 (see above).
191  */
192  int set_npoints2(int npoints) {
193  if (npoints==3) {
195  } else if (npoints==4) {
197  } else {
199  }
200  if (npoints<=2 || npoints>5) {
201  O2SCL_ERR("Invalid # of points in set_npoints2(). Using default",
202  exc_einval);
203  }
204  return 0;
205  }
206 
207  /** \brief Calculate the first derivative of \c func w.r.t. x
208  */
209  virtual int deriv_err(double x, func_t &func,
210  double &dfdx, double &err) {
211  double p=xoff/h;
212  dfdx=(this->*cp)(x,p,func)/h;
213  err=0.0;
214  return success;
215  }
216 
217  /** \brief Calculate the second derivative of \c func w.r.t. x
218  */
219  virtual int deriv2_err(double x, func_t &func,
220  double &dfdx, double &err) {
221  double p=xoff/h;
222  dfdx=(this->*c2p)(x,p,func)/h/h;
223  err=0.0;
224  return success;;
225  }
226 
227  /** \brief Calculate the third derivative of \c func w.r.t. x
228  */
229  virtual int deriv3_err(double x, func_t &func,
230  double &dfdx, double &err) {
231  double p=xoff/h;
232  dfdx=(this->*c3p)(x,h,p,func)/h;
233  err=0.0;
234  return success;;
235  }
236 
237  /** \brief Calculate the derivative at \c x given an array
238 
239  This calculates the derivative at \c x given a function
240  specified in an array \c y of size \c nx with equally spaced
241  abscissas. The first abscissa should be given as \c x0
242  and the distance between adjacent abscissas should be
243  given as \c dx. The value \c x need not be one of the
244  abscissas (i.e. it can lie in between an interval). The
245  appropriate offset is calculated automatically.
246  */
247  double deriv_vector(double x, double x0, double dx,
248  size_t nx, const vec_t &y) {
249  size_t ix=(size_t)((x-x0)/dx);
250  return (this->*cap)(x,x0,dx,nx,y,ix)/dx;
251  }
252 
253  /** \brief Calculate the second derivative at \c x given an array
254 
255  This calculates the second derivative at \c x given a function
256  specified in an array \c y of size \c nx with equally spaced
257  abscissas. The first abscissa should be given as \c x0
258  and the distance between adjacent abscissas should be
259  given as \c dx. The value \c x need not be one of the
260  abscissas (i.e. it can lie in between an interval). The
261  appropriate offset is calculated automatically.
262  */
263  double deriv2_vector(double x, double x0, double dx,
264  size_t nx, const vec_t &y)
265  {
266  size_t ix=(size_t)((x-x0)/dx);
267  return (this->*c2ap)(x,x0,dx,nx,y,ix)/dx;
268  }
269 
270  /** \brief Calculate the third derivative at \c x given an array
271 
272  This calculates the third derivative at \c x given a function
273  specified in an array \c y of size \c nx with equally spaced
274  abscissas. The first abscissa should be given as \c x0 and the
275  distance between adjacent abscissas should be given as \c
276  dx. The value \c x need not be one of the abscissas (i.e. it
277  can lie in between an interval). The appropriate offset is
278  calculated automatically.
279  */
280  double deriv3_vector(double x, double x0, double dx,
281  size_t nx, const vec_t &y)
282  {
283  size_t ix=(size_t)((x-x0)/dx);
284  return (this->*c3ap)(x,x0,dx,nx,y,ix)/dx;
285  }
286 
287  /** \brief Calculate the derivative of an entire array
288 
289  Right now this uses np=5.
290 
291  \todo generalize to other values of npoints.
292  */
293  int deriv_vector(size_t nv, double dx, const vec_t &y,
294  vec_t &dydx)
295  {
296  dydx[0]=(-25.0/12.0*y[0]+4.0*y[1]-3.0*y[2]+4.0/3.0*y[3]-0.25*y[4])/dx;
297  dydx[1]=(-0.25*y[0]-5.0/6.0*y[1]+1.5*y[2]-0.5*y[3]+1.0/12.0*y[4])/dx;
298  for(size_t i=2;i<nv-2;i++) {
299  dydx[i]=(1.0/12.0*y[i-2]-2.0/3.0*y[i-1]+2.0/3.0*y[i+1]-
300  1.0/12.0*y[i+2])/dx;
301  }
302  dydx[nv-2]=(-1.0/12.0*y[nv-5]+0.5*y[nv-4]-1.5*y[nv-3]+
303  5.0/6.0*y[nv-2]+0.25*y[nv-1])/dx;
304  dydx[nv-1]=(0.25*y[nv-5]-4.0/3.0*y[nv-4]+3.0*y[nv-3]-
305  4.0*y[nv-2]+25.0/12.0*y[nv-1])/dx;
306  return 0;
307  }
308 
309  /// Return string denoting type ("deriv_eqi")
310  virtual const char *type() { return "deriv_eqi"; }
311 
312 #ifndef DOXYGEN_NO_O2NS
313 
314  protected:
315 
316  /** \brief Calculate the first derivative of \c func w.r.t. x and the
317  uncertainty
318 
319  This function doesn't do anything, and isn't required for
320  this class since it computes higher-order derivatives directly.
321  */
322  virtual int deriv_err_int
323  (double x, funct11 &func, double &dfdx, double &err) {
324  return success;
325  }
326 
327  /// Two-point first derivative
328  double derivp2(double x, double p, func_t &func) {
329  return (func(x+h)-func(x));
330  }
331 
332  /// Three-point first derivative
333  double derivp3(double x, double p, func_t &func) {
334  if (p==0.0) {
335  return ((-0.5)*func(x-h)+(0.5)*func(x+h));
336  }
337  return ((p-0.5)*func(x-h)-2.0*p*func(x)+(p+0.5)*func(x+h));
338  }
339 
340  /// Four-point first derivative
341  double derivp4(double x, double p,
342  func_t &func) {
343  double p2=p*p;
344  return (-(3.0*p2-6.0*p+2.0)/6.0*func(x-h)+
345  (1.5*p2-2.0*p-0.5)*func(x)-
346  (1.5*p2-p-1.0)*func(x+h)+
347  (3.0*p2-1.0)/6.0*func(x+2.0*h));
348  }
349 
350  /// Five-point first derivative
351  double derivp5(double x, double p,
352  func_t &func) {
353  double p2=p*p, p3=p*p*p;
354  if (p==0.0) {
355  return ((1.0)/12.0*func(x-2.0*h)-
356  (4.0)/6.0*func(x-h)-
357  (-4.0)/6.0*func(x+h)+
358  (-1.0)/12.0*func(x+2.0*h));
359  }
360  return ((2.0*p3-3.0*p2-p+1.0)/12.0*func(x-2.0*h)-
361  (4.0*p3-3.0*p2-8.0*p+4.0)/6.0*func(x-h)+
362  (p3-2.5*p)*func(x)-
363  (4.0*p3+3.0*p2-8.0*p-4.0)/6.0*func(x+h)+
364  (2.0*p3+3.0*p2-p-1.0)/12.0*func(x+2.0*h));
365  }
366 
367 
368  /// Three-point first derivative for arrays
369  double deriv_vector3(double x, double x0, double dx, size_t nx,
370  const vec_t &y, size_t ix) {
371  double p;
372  if (ix>0 && ix<nx-1) {
373  p=x-(x0+ix*dx);
374  return ((p-0.5)*y[ix-1]-2.0*p*y[ix]+(p+0.5)*y[ix+1]);
375  } else if (ix==0) {
376  p=x-(x0+dx);
377  return ((p-0.5)*y[0]-2.0*p*y[1]+(p+0.5)*y[2]);
378  }
379  p=x-(x0+(nx-2)*dx);
380  return ((p-0.5)*y[nx-3]-2.0*p*y[nx-2]+(p+0.5)*y[nx-1]);
381  }
382 
383  /// Four-point first derivative for arrays
384  double deriv_vector4(double x, double x0, double dx, size_t nx,
385  const vec_t &y, size_t ix) {
386  double p, p2;
387  if (ix>0 && ix<nx-2) {
388  p=x-(x0+ix*dx);
389  p2=p*p;
390  return (-(3.0*p2-6.0*p+2.0)/6.0*y[ix-1]+
391  (1.5*p2-2.0*p-0.5)*y[ix]-
392  (1.5*p2-p-1.0)*y[ix+1]+
393  (3.0*p2-1.0)/6.0*y[ix+2]);
394  } else if (ix==0) {
395  p=x-(x0+dx);
396  p2=p*p;
397  return (-(3.0*p2-6.0*p+2.0)/6.0*y[0]+
398  (1.5*p2-2.0*p-0.5)*y[1]-
399  (1.5*p2-p-1.0)*y[2]+
400  (3.0*p2-1.0)/6.0*y[3]);
401  }
402  p=x-(x0+(nx-3)*dx);
403  p2=p*p;
404  return (-(3.0*p2-6.0*p+2.0)/6.0*y[nx-4]+
405  (1.5*p2-2.0*p-0.5)*y[nx-3]-
406  (1.5*p2-p-1.0)*y[nx-2]+
407  (3.0*p2-1.0)/6.0*y[nx-1]);
408  }
409 
410  /// Five-point first derivative for arrays
411  double deriv_vector5(double x, double x0,
412  double dx, size_t nx,
413  const vec_t &y, size_t ix) {
414  double p, p2, p3;
415  if (ix>1 && ix<nx-2) {
416  p=x-(x0+ix*dx);
417  p2=p*p, p3=p*p*p;
418  return ((2.0*p3-3.0*p2-p+1.0)/12.0*y[ix-2]-
419  (4.0*p3-3.0*p2-8.0*p+4.0)/6.0*y[ix-1]+
420  (p3-2.5*p)*y[ix]-
421  (4.0*p3+3.0*p2-8.0*p-4.0)/6.0*y[ix+1]+
422  (2.0*p3+3.0*p2-p-1.0)/12.0*y[ix+2]);
423  } else if (ix<=1) {
424  p=x-(x0+2*dx);
425  p2=p*p, p3=p*p*p;
426  return ((2.0*p3-3.0*p2-p+1.0)/12.0*y[0]-
427  (4.0*p3-3.0*p2-8.0*p+4.0)/6.0*y[1]+
428  (p3-2.5*p)*y[2]-
429  (4.0*p3+3.0*p2-8.0*p-4.0)/6.0*y[3]+
430  (2.0*p3+3.0*p2-p-1.0)/12.0*y[4]);
431  }
432  p=x-(x0+(nx-3)*dx);
433  p2=p*p, p3=p*p*p;
434  return ((2.0*p3-3.0*p2-p+1.0)/12.0*y[nx-5]-
435  (4.0*p3-3.0*p2-8.0*p+4.0)/6.0*y[nx-4]+
436  (p3-2.5*p)*y[nx-3]-
437  (4.0*p3+3.0*p2-8.0*p-4.0)/6.0*y[nx-2]+
438  (2.0*p3+3.0*p2-p-1.0)/12.0*y[nx-1]);
439  }
440 
441  /// Three-point second derivative
442  double deriv2p3(double x, double p, func_t &func) {
443  return (func(x-h)-2*func(x)+func(x+h));
444  }
445 
446  /// Four-point second derivative
447  double deriv2p4(double x, double p, func_t &func) {
448  return ((1.0-2.0*p)*func(x-h)-(1.0-6.0*p)*func(x)
449  -(1.0-6.0*p)*func(x+h)+(1.0+2.0*p)*func(x+2.0*h))/2.0;
450  }
451 
452  /// Five-point second derivative
453  double deriv2p5(double x, double p, func_t &func) {
454  return ((1.0-2.0*p)*(1.0-2.0*p)*func(x-2.0*h)
455  +(8.0*p-16.0*p*p)*func(x-h)
456  -(2.0-24.0*p*p)*func(x)
457  -(8.0*p+16.0*p*p)*func(x+h)
458  +(1.0+2.0*p)*(1.0+2.0*p)*func(x+2.0*h))/4.0;
459  }
460 
461  /// Pointer to the first derivative function
462  double (deriv_eqi::*cp)(double x, double p,
463  func_t &func);
464 
465  /// Pointer to the first derivative for arrays function
466  double (deriv_eqi::*cap)(double x, double x0,
467  double dx, size_t nx,
468  const vec_t &y, size_t ix);
469 
470  /// Pointer to the second derivative function
471  double (deriv_eqi::*c2p)(double x, double p,
472  func_t &func);
473 
474  /// Pointer to the second derivative for arrays function
475  double (deriv_eqi::*c2ap)(double x, double x0,
476  double dx, size_t nx,
477  const vec_t &y, size_t ix);
478 
479  /// Pointer to the third derivative function
480  double (deriv_eqi::*c3p)(double x, double h, double p,
481  func_t &func);
482 
483  /// Pointer to the third derivative for arrays function
484  double (deriv_eqi::*c3ap)(double x, double x0,
485  double dx, size_t nx,
486  const vec_t &y, size_t ix);
487 
488 #endif
489 
490  };
491 
492 #ifndef DOXYGEN_NO_O2NS
493 }
494 #endif
495 
496 #endif
double h
Stepsize (Default ).
Definition: deriv_eqi.h:156
double derivp5(double x, double p, func_t &func)
Five-point first derivative.
Definition: deriv_eqi.h:351
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
double(deriv_eqi::* c2p)(double x, double p, func_t &func)
Pointer to the second derivative function.
Definition: deriv_eqi.h:471
int deriv_vector(size_t nv, double dx, const vec_t &y, vec_t &dydx)
Calculate the derivative of an entire array.
Definition: deriv_eqi.h:293
double(deriv_eqi::* cap)(double x, double x0, double dx, size_t nx, const vec_t &y, size_t ix)
Pointer to the first derivative for arrays function.
Definition: deriv_eqi.h:466
std::function< double(double)> funct11
One-dimensional function typedef.
Definition: funct.h:44
double deriv_vector4(double x, double x0, double dx, size_t nx, const vec_t &y, size_t ix)
Four-point first derivative for arrays.
Definition: deriv_eqi.h:384
virtual int deriv3_err(double x, func_t &func, double &dfdx, double &err)
Calculate the third derivative of func w.r.t. x.
Definition: deriv_eqi.h:229
virtual int deriv2_err(double x, func_t &func, double &dfdx, double &err)
Calculate the second derivative of func w.r.t. x.
Definition: deriv_eqi.h:219
double deriv_vector5(double x, double x0, double dx, size_t nx, const vec_t &y, size_t ix)
Five-point first derivative for arrays.
Definition: deriv_eqi.h:411
invalid argument supplied by user
Definition: err_hnd.h:59
double deriv_vector(double x, double x0, double dx, size_t nx, const vec_t &y)
Calculate the derivative at x given an array.
Definition: deriv_eqi.h:247
double derivp4(double x, double p, func_t &func)
Four-point first derivative.
Definition: deriv_eqi.h:341
double derivp2(double x, double p, func_t &func)
Two-point first derivative.
Definition: deriv_eqi.h:328
double derivp3(double x, double p, func_t &func)
Three-point first derivative.
Definition: deriv_eqi.h:333
virtual int deriv_err(double x, func_t &func, double &dfdx, double &err)
Calculate the first derivative of func w.r.t. x.
Definition: deriv_eqi.h:209
virtual int deriv_err_int(double x, funct11 &func, double &dfdx, double &err)
Calculate the first derivative of func w.r.t. x and the uncertainty.
Definition: deriv_eqi.h:323
double deriv2p3(double x, double p, func_t &func)
Three-point second derivative.
Definition: deriv_eqi.h:442
double deriv2p5(double x, double p, func_t &func)
Five-point second derivative.
Definition: deriv_eqi.h:453
double deriv3_vector(double x, double x0, double dx, size_t nx, const vec_t &y)
Calculate the third derivative at x given an array.
Definition: deriv_eqi.h:280
double(deriv_eqi::* c3p)(double x, double h, double p, func_t &func)
Pointer to the third derivative function.
Definition: deriv_eqi.h:480
Numerical differentiation base [abstract base].
Definition: deriv.h:63
double deriv_vector3(double x, double x0, double dx, size_t nx, const vec_t &y, size_t ix)
Three-point first derivative for arrays.
Definition: deriv_eqi.h:369
double deriv2_vector(double x, double x0, double dx, size_t nx, const vec_t &y)
Calculate the second derivative at x given an array.
Definition: deriv_eqi.h:263
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int set_npoints2(int npoints)
Set the number of points to use for second derivatives (default 5)
Definition: deriv_eqi.h:192
Derivatives for equally-spaced abscissas.
Definition: deriv_eqi.h:144
double(deriv_eqi::* cp)(double x, double p, func_t &func)
Pointer to the first derivative function.
Definition: deriv_eqi.h:462
double xoff
Offset (default 0.0)
Definition: deriv_eqi.h:159
double deriv2p4(double x, double p, func_t &func)
Four-point second derivative.
Definition: deriv_eqi.h:447
int set_npoints(int npoints)
Set the number of points to use for first derivatives (default 5)
Definition: deriv_eqi.h:166
virtual const char * type()
Return string denoting type ("deriv_eqi")
Definition: deriv_eqi.h:310
double(deriv_eqi::* c3ap)(double x, double x0, double dx, size_t nx, const vec_t &y, size_t ix)
Pointer to the third derivative for arrays function.
Definition: deriv_eqi.h:484
double(deriv_eqi::* c2ap)(double x, double x0, double dx, size_t nx, const vec_t &y, size_t ix)
Pointer to the second derivative for arrays function.
Definition: deriv_eqi.h:475
Success.
Definition: err_hnd.h:47

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).