vector_derint.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2018, 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 
160  /** \name Compute the integral over a generic vector
161 
162  These functions implement composite (sometimes called
163  extended) Newton-Cotes rules.
164 
165  Given a vector \c v of size \c n, these functions compute
166  the integral over the entire vector and return the result.
167 
168  These functions were principally based on the information at
169  http://mathworld.wolfram.com/Newton-CotesFormulas.html .
170  */
171  //@{
172  /** \brief Integrate with an extended trapezoidal rule.
173  */
174  template<class vec_t> double vector_integ_trap(size_t n, vec_t &v) {
175  if (n==0) {
176  O2SCL_ERR2("Tried to integrate zero-length vector in ",
177  "vector_integ_trap",exc_einval);
178  } else if (n==1) {
179  return v[0];
180  }
181  double res=(v[0]+v[n-1])/2.0;
182  for(size_t i=1;i<n-1;i++) res+=v[i];
183  return res;
184  }
185 
186  /** \brief Integrate with an extended 3-point rule (extended
187  Simpson's rule)
188 
189  \note This function uses an untested hack I wrote for even n.
190  */
191  template<class vec_t> double vector_integ_threept(size_t n, vec_t &v) {
192  if (n==0) {
193  O2SCL_ERR2("Tried to integrate zero-length vector in ",
194  "vector_integ_threept",exc_einval);
195  } else if (n==1) {
196  return v[0];
197  } else if (n==2) {
198  return (v[0]+v[1])/2.0;
199  }
200  double res=v[0]+v[n-1];
201  // If true, next terms have a prefactor of 4, otherwise
202  // the next terms have a prefactor of 2
203  bool four=true;
204  for(size_t i=1;i<n/2;i++) {
205  // Correct if n is even
206  if (i==n-i-2) {
207  if (four) res+=(v[i]+v[n-i-1])*3.5;
208  else res+=(v[i]+v[n-i-1])*2.5;
209  } else {
210  if (four) res+=(v[i]+v[n-i-1])*4.0;
211  else res+=(v[i]+v[n-i-1])*2.0;
212  }
213  four=!four;
214  }
215  return res/3.0;
216  }
217 
218  /** \brief Integrate with an extended rule for 4 or more points.
219 
220  This function falls back to the equivalent of
221  vector_integ_threept() for 3 points.
222  */
223  template<class vec_t> double vector_integ_extended4(size_t n, vec_t &v) {
224  if (n==0) {
225  O2SCL_ERR2("Tried to integrate zero-length vector in ",
226  "vector_integ_extended4",exc_einval);
227  } else if (n==1) {
228  return v[0];
229  } else if (n==2) {
230  return (v[0]+v[1])/2.0;
231  } else if (n==3) {
232  return (v[0]+4.0*v[1]+v[2])/3.0;
233  }
234  double res=(v[0]*5.0+v[1]*13.0+v[n-1]*5.0+v[n-2]*13.0)/12.0;
235  for(size_t i=2;i<n-2;i++) {
236  res+=v[i];
237  }
238  return res;
239  }
240 
241  /** \brief Integrate with Durand's rule for 4 or more points.
242 
243  This function falls back to the equivalent of
244  vector_integ_threept() for 3 points.
245  */
246  template<class vec_t> double vector_integ_durand(size_t n, vec_t &v) {
247  if (n==0) {
248  O2SCL_ERR2("Tried to integrate zero-length vector in ",
249  "vector_integ_durand",exc_einval);
250  } else if (n==1) {
251  return v[0];
252  } else if (n==2) {
253  return (v[0]+v[1])/2.0;
254  } else if (n==3) {
255  return (v[0]+4.0*v[1]+v[2])/3.0;
256  }
257  double res=(v[0]*4.0+v[1]*11.0+v[n-1]*4.0+v[n-2]*11.0)/10.0;
258  for(size_t i=2;i<n-2;i++) {
259  res+=v[i];
260  }
261  return res;
262  }
263 
264  /** \brief Integrate with an extended rule for 8 or more points.
265 
266  This function falls back to vector_integ_extended4()
267  for less than 8 points.
268  */
269  template<class vec_t> double vector_integ_extended8(size_t n, vec_t &v) {
270  if (n<8) return vector_integ_extended4(n,v);
271  double res=((v[0]+v[n-1])*17.0+(v[1]+v[n-2])*59.0+(v[2]+v[n-3])*43.0+
272  (v[3]+v[n-4])*49.0)/48.0;
273  for(size_t i=4;i<n-4;i++) {
274  res+=v[i];
275  }
276  return res;
277  }
278  //@}
279 
280 #ifndef DOXYGEN_NO_O2NS
281 }
282 #endif
283 
284 #endif
void vector_deriv_fivept(size_t n, vec_t &v, vec2_t &dv)
Derivative of a vector with a five-point formula.
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_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
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...
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
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).