vector_derint.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-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_VECTOR_DERINT_H
24 #define O2SCL_VECTOR_DERINT_H
25 
26 /** \file vector_derint.h
27  \brief Derivatives of integrals of functions stored in vectors
28  with implicit fixed-size grid
29 
30  These functions differentiate or integrate a function over a
31  fixed-size grid specified in a vector. Derivatives and integrals
32  are always specified without the factor defining the grid size
33  (i.e. \f$ dx \f$), so the user must always multiply or divide the
34  result by the grid size afterwards as necessary.
35 
36  The derivative functions will call the error handler if they are
37  given empty vectors or vectors of length 1. The integration rules
38  often expect a minimum number of points, so for smaller vectors
39  they fall back onto rules which use fewer points. For empty
40  vectors they return zero, for vectors of length 1, they always
41  return the sole element of the vector, and for vectors of length
42  2, they always return the average of the two elements. If the
43  vector length is zero, the integration functions call the
44  error handler.
45 
46  More points does not always mean higher accuracy.
47 */
48 #include <o2scl/interp.h>
49 
50 #ifndef DOXYGEN_NO_O2NS
51 namespace o2scl {
52 #endif
53 
54  /** \name Compute the derivative at every point of a generic vector
55 
56  Given a vector \c v of size \c n, these functions compute
57  the derivative at every point and store the result in \c dv.
58  */
59  //@{
60  /** \brief Derivative of a vector with a three-point formula
61  */
62  template<class vec_t, class vec2_t>
63  void vector_deriv_threept(size_t n, vec_t &v, vec2_t &dv) {
64  if (n<=1) {
65  O2SCL_ERR2("Requested derivative of zero or one-element vector ",
66  "in vector_deriv_threept().",exc_einval);
67  }
68  if (n==2) {
69  double d=v[1]-v[0];
70  dv[0]=d;
71  dv[1]=d;
72  return;
73  }
74  dv[0]=-1.5*v[0]+2.0*v[1]-0.5*v[2];
75  dv[n-1]=1.5*v[n-1]-2.0*v[n-2]+0.5*v[n-3];
76  for(size_t i=1;i<n-1;i++) {
77  dv[i]=(v[i+1]-v[i-1])/2.0;
78  }
79  return;
80  }
81 
82  /** \brief Derivative of a vector with a three-point formula
83  using two-point at the edges
84  */
85  template<class vec_t, class vec2_t>
86  void vector_deriv_threept_tap(size_t n, vec_t &v, vec2_t &dv) {
87  if (n<=1) {
88  O2SCL_ERR2("Requested derivative of zero or one-element vector ",
89  "in vector_deriv_threept().",exc_einval);
90  }
91  if (n==2) {
92  double d=v[1]-v[0];
93  dv[0]=d;
94  dv[1]=d;
95  return;
96  }
97  // 2-point
98  dv[0]=v[1]-v[0];
99  dv[n-1]=v[n-1]-v[n-2];
100  // 3-point
101  for(size_t i=1;i<n-1;i++) {
102  dv[i]=(v[i+1]-v[i-1])/2.0;
103  }
104  return;
105  }
106 
107  /** \brief Derivative of a vector with a five-point formula
108  */
109  template<class vec_t, class vec2_t>
110  void vector_deriv_fivept(size_t n, vec_t &v, vec2_t &dv) {
111  if (n<=1) {
112  O2SCL_ERR2("Requested derivative of zero or one-element vector ",
113  "in vector_deriv_fivept().",exc_einval);
114  }
115  if (n==2) {
116  double d=v[1]-v[0];
117  dv[0]=d;
118  dv[1]=d;
119  return;
120  }
121  dv[0]=-25.0/12.0*v[0]+4.0*v[1]-3.0*v[2]+4.0/3.0*v[3]-0.25*v[4];
122  dv[1]=-0.25*v[0]-5.0/6.0*v[1]+1.5*v[2]-0.5*v[3]+v[4]/12.0;
123  dv[n-2]=-v[n-5]/12.0+0.5*v[n-4]-1.5*v[n-3]+5.0/6.0*v[n-2]+0.25*v[n-1];
124  dv[n-1]=0.25*v[n-5]-4.0*v[n-4]/3.0+3.0*v[n-3]-4.0*v[n-2]+25.0*v[n-1]/12.0;
125  for(size_t i=2;i<n-2;i++) {
126  dv[i]=1.0/12.0*(v[i-2]-v[i+2])+2.0/3.0*(v[i+1]-v[i-1]);
127  }
128  return;
129  }
130 
131  /** \brief Derivative of a vector with a five-point formula with
132  four- and three-point formulas used at the edges
133  */
134  template<class vec_t, class vec2_t>
135  void vector_deriv_fivept_tap(size_t n, vec_t &v, vec2_t &dv) {
136  if (n<=1) {
137  O2SCL_ERR2("Requested derivative of zero or one-element vector ",
138  "in vector_deriv_fivept().",exc_einval);
139  }
140  if (n==2) {
141  double d=v[1]-v[0];
142  dv[0]=d;
143  dv[1]=d;
144  return;
145  }
146  // 3-point
147  dv[0]=-1.5*v[0]+2.0*v[1]-0.5*v[2];
148  dv[n-1]=1.5*v[n-1]-2.0*v[n-2]+0.5*v[n-3];
149  // 4-point
150  dv[1]=-v[0]/3.0-v[1]/2.0+v[2]-v[3]/6.0;
151  dv[n-2]=v[n-4]/6.0-v[n-3]+v[n-2]/2.0+v[n-1]/3.0;
152  // 5-point
153  for(size_t i=2;i<n-2;i++) {
154  dv[i]=1.0/12.0*(v[i-2]-v[i+2])+2.0/3.0*(v[i+1]-v[i-1]);
155  }
156  return;
157  }
158 
159  /** \brief Derivative from interpolation object
160  */
161  template<class ovec_t, class vec2_t>
162  void vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv,
163  size_t interp_type=itp_cspline) {
164  ovec_t grid(n);
165  for(size_t i=0;i<n;i++) grid[i]=((double)i);
166  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
167  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(((double)i));
168  return;
169  }
170 
171  /** \brief Derivative from interpolation object
172  */
173  template<class vec_t, class vec2_t, class vec3_t>
174  void vector_deriv_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
175  size_t interp_type=itp_cspline) {
176  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
177  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
178  return;
179  }
180  //@}
181 
182  /** \name Compute the integral over a generic vector
183 
184  These functions implement composite (sometimes called
185  extended) Newton-Cotes rules.
186 
187  Given a vector \c v of size \c n, these functions compute
188  the integral over the entire vector and return the result.
189 
190  These functions were principally based on the information at
191  http://mathworld.wolfram.com/Newton-CotesFormulas.html .
192  */
193  //@{
194  /** \brief Integrate with an extended trapezoidal rule.
195  */
196  template<class vec_t> double vector_integ_trap(size_t n, vec_t &v) {
197  if (n==0) {
198  O2SCL_ERR2("Tried to integrate zero-length vector in ",
199  "vector_integ_trap",exc_einval);
200  } else if (n==1) {
201  return v[0];
202  }
203  double res=(v[0]+v[n-1])/2.0;
204  for(size_t i=1;i<n-1;i++) res+=v[i];
205  return res;
206  }
207 
208  /** \brief Integrate with an extended 3-point rule (extended
209  Simpson's rule)
210 
211  \note This function uses an untested hack I wrote for even n.
212  */
213  template<class vec_t> double vector_integ_threept(size_t n, vec_t &v) {
214  if (n==0) {
215  O2SCL_ERR2("Tried to integrate zero-length vector in ",
216  "vector_integ_threept",exc_einval);
217  } else if (n==1) {
218  return v[0];
219  } else if (n==2) {
220  return (v[0]+v[1])/2.0;
221  }
222  double res=v[0]+v[n-1];
223  // If true, next terms have a prefactor of 4, otherwise
224  // the next terms have a prefactor of 2
225  bool four=true;
226  for(size_t i=1;i<n/2;i++) {
227  // Correct if n is even
228  if (i==n-i-2) {
229  if (four) res+=(v[i]+v[n-i-1])*3.5;
230  else res+=(v[i]+v[n-i-1])*2.5;
231  } else {
232  if (four) res+=(v[i]+v[n-i-1])*4.0;
233  else res+=(v[i]+v[n-i-1])*2.0;
234  }
235  four=!four;
236  }
237  return res/3.0;
238  }
239 
240  /** \brief Integrate with an extended rule for 4 or more points.
241 
242  This function falls back to the equivalent of
243  vector_integ_threept() for 3 points.
244  */
245  template<class vec_t> double vector_integ_extended4(size_t n, vec_t &v) {
246  if (n==0) {
247  O2SCL_ERR2("Tried to integrate zero-length vector in ",
248  "vector_integ_extended4",exc_einval);
249  } else if (n==1) {
250  return v[0];
251  } else if (n==2) {
252  return (v[0]+v[1])/2.0;
253  } else if (n==3) {
254  return (v[0]+4.0*v[1]+v[2])/3.0;
255  }
256  double res=(v[0]*5.0+v[1]*13.0+v[n-1]*5.0+v[n-2]*13.0)/12.0;
257  for(size_t i=2;i<n-2;i++) {
258  res+=v[i];
259  }
260  return res;
261  }
262 
263  /** \brief Integrate with Durand's rule for 4 or more points.
264 
265  This function falls back to the equivalent of
266  vector_integ_threept() for 3 points.
267  */
268  template<class vec_t> double vector_integ_durand(size_t n, vec_t &v) {
269  if (n==0) {
270  O2SCL_ERR2("Tried to integrate zero-length vector in ",
271  "vector_integ_durand",exc_einval);
272  } else if (n==1) {
273  return v[0];
274  } else if (n==2) {
275  return (v[0]+v[1])/2.0;
276  } else if (n==3) {
277  return (v[0]+4.0*v[1]+v[2])/3.0;
278  }
279  double res=(v[0]*4.0+v[1]*11.0+v[n-1]*4.0+v[n-2]*11.0)/10.0;
280  for(size_t i=2;i<n-2;i++) {
281  res+=v[i];
282  }
283  return res;
284  }
285 
286  /** \brief Integrate with an extended rule for 8 or more points.
287 
288  This function falls back to vector_integ_extended4()
289  for less than 8 points.
290  */
291  template<class vec_t> double vector_integ_extended8(size_t n, vec_t &v) {
292  if (n<8) return vector_integ_extended4(n,v);
293  double res=((v[0]+v[n-1])*17.0+(v[1]+v[n-2])*59.0+(v[2]+v[n-3])*43.0+
294  (v[3]+v[n-4])*49.0)/48.0;
295  for(size_t i=4;i<n-4;i++) {
296  res+=v[i];
297  }
298  return res;
299  }
300 
301  /** \brief Integral from interpolation object
302  */
303  template<class ovec_t>
304  double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type) {
305  ovec_t grid(n);
306  for(size_t i=0;i<n;i++) grid[i]=((double)i);
307  interp_vec<ovec_t> oi(n,grid,v,interp_type);
308  return oi.integ(0.0,((double)(n-1)));
309  }
310  //@}
311 
312 #ifndef DOXYGEN_NO_O2NS
313 }
314 #endif
315 
316 #endif
void vector_deriv_fivept(size_t n, vec_t &v, vec2_t &dv)
Derivative of a vector with a five-point formula.
Interpolation class for pre-specified vector.
Definition: interp.h:1685
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 vector_integ_durand(size_t n, vec_t &v)
Integrate with Durand&#39;s rule for 4 or more points.
double vector_integ_extended8(size_t n, vec_t &v)
Integrate with an extended rule for 8 or more points.
void vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv, size_t interp_type=itp_cspline)
Derivative from interpolation object.
void vector_deriv_threept_tap(size_t n, vec_t &v, vec2_t &dv)
Derivative of a vector with a three-point formula using two-point at the edges.
Definition: vector_derint.h:86
invalid argument supplied by user
Definition: err_hnd.h:59
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1820
double vector_integ_trap(size_t n, vec_t &v)
Integrate with an extended trapezoidal rule.
double vector_integ_threept(size_t n, vec_t &v)
Integrate with an extended 3-point rule (extended Simpson&#39;s rule)
void vector_deriv_fivept_tap(size_t n, vec_t &v, vec2_t &dv)
Derivative of a vector with a five-point formula with four- and three-point formulas used at the edge...
Cubic spline for natural boundary conditions.
Definition: interp.h:72
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral from interpolation object.
double vector_integ_extended4(size_t n, vec_t &v)
Integrate with an extended rule for 4 or more points.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
virtual double integ(const double x1, const double x2) const
Give the value of the integral .
Definition: interp.h:1840
void vector_deriv_threept(size_t n, vec_t &v, vec2_t &dv)
Derivative of a vector with a three-point formula.
Definition: vector_derint.h:63

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