vector_derint.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2020, 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 Derivative of a generic vector in src/deriv/vector_derint.h
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, class fp_t=double>
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  fp_t 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, class fp_t=double>
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  fp_t 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, class fp_t=double>
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  fp_t 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, class fp_t=double>
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  fp_t 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 Integral of a generic vector in src/deriv/vector_derint.h
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, class fp_t=double> fp_t vector_integ_trap
175  (size_t n, vec_t &v) {
176  if (n==0) {
177  O2SCL_ERR2("Tried to integrate zero-length vector in ",
178  "vector_integ_trap",exc_einval);
179  } else if (n==1) {
180  return v[0];
181  }
182  fp_t res=(v[0]+v[n-1])/2.0;
183  for(size_t i=1;i<n-1;i++) res+=v[i];
184  return res;
185  }
186 
187  /** \brief Integrate with an extended 3-point rule (extended
188  Simpson's rule)
189 
190  \note This function uses an untested hack I wrote for even n.
191  */
192  template<class vec_t, class fp_t=double> fp_t vector_integ_threept
193  (size_t n, vec_t &v) {
194  if (n==0) {
195  O2SCL_ERR2("Tried to integrate zero-length vector in ",
196  "vector_integ_threept",exc_einval);
197  } else if (n==1) {
198  return v[0];
199  } else if (n==2) {
200  return (v[0]+v[1])/2.0;
201  }
202  fp_t res=v[0]+v[n-1];
203  // If true, next terms have a prefactor of 4, otherwise
204  // the next terms have a prefactor of 2
205  bool four=true;
206  for(size_t i=1;i<n/2;i++) {
207  // Correct if n is even
208  if (i==n-i-2) {
209  if (four) res+=(v[i]+v[n-i-1])*3.5;
210  else res+=(v[i]+v[n-i-1])*2.5;
211  } else {
212  if (four) res+=(v[i]+v[n-i-1])*4.0;
213  else res+=(v[i]+v[n-i-1])*2.0;
214  }
215  four=!four;
216  }
217  return res/3.0;
218  }
219 
220  /** \brief Integrate with an extended rule for 4 or more points.
221 
222  This function falls back to the equivalent of
223  vector_integ_threept() for 3 points.
224  */
225  template<class vec_t, class fp_t=double> fp_t vector_integ_extended4
226  (size_t n, vec_t &v) {
227  if (n==0) {
228  O2SCL_ERR2("Tried to integrate zero-length vector in ",
229  "vector_integ_extended4",exc_einval);
230  } else if (n==1) {
231  return v[0];
232  } else if (n==2) {
233  fp_t two=2;
234  return (v[0]+v[1])/two;
235  } else if (n==3) {
236  fp_t three=3;
237  return (v[0]+4.0*v[1]+v[2])/three;
238  }
239  fp_t c1=5;
240  fp_t c2=13;
241  fp_t den=12;
242  fp_t res=(v[0]*c1+v[1]*c2+v[n-1]*c1+v[n-2]*c2)/den;
243  for(size_t i=2;i<n-2;i++) {
244  res+=v[i];
245  }
246  return res;
247  }
248 
249  /** \brief Integrate with Durand's rule for 4 or more points.
250 
251  This function falls back to the equivalent of
252  vector_integ_threept() for 3 points.
253  */
254  template<class vec_t, class fp_t=double> fp_t vector_integ_durand
255  (size_t n, vec_t &v) {
256  if (n==0) {
257  O2SCL_ERR2("Tried to integrate zero-length vector in ",
258  "vector_integ_durand",exc_einval);
259  } else if (n==1) {
260  return v[0];
261  } else if (n==2) {
262  return (v[0]+v[1])/2.0;
263  } else if (n==3) {
264  return (v[0]+4.0*v[1]+v[2])/3.0;
265  }
266  fp_t res=(v[0]*4.0+v[1]*11.0+v[n-1]*4.0+v[n-2]*11.0)/10.0;
267  for(size_t i=2;i<n-2;i++) {
268  res+=v[i];
269  }
270  return res;
271  }
272 
273  /** \brief Integrate with an extended rule for 8 or more points.
274 
275  This function falls back to vector_integ_extended4()
276  for less than 8 points.
277  */
278  template<class vec_t, class fp_t=double> fp_t vector_integ_extended8
279  (size_t n, vec_t &v) {
280  if (n<8) return vector_integ_extended4<vec_t,fp_t>(n,v);
281  fp_t c1=17;
282  fp_t c2=59;
283  fp_t c3=43;
284  fp_t c4=49;
285  fp_t den=48;
286  fp_t res=((v[0]+v[n-1])*c1+(v[1]+v[n-2])*c2+(v[2]+v[n-3])*c3+
287  (v[3]+v[n-4])*c4)/den;
288  for(size_t i=4;i<n-4;i++) {
289  res+=v[i];
290  }
291  return res;
292  }
293  //@}
294 
295 #ifndef DOXYGEN_NO_O2NS
296 }
297 #endif
298 
299 #endif
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::vector_deriv_fivept
void vector_deriv_fivept(size_t n, vec_t &v, vec2_t &dv)
Derivative of a vector with a five-point formula.
Definition: vector_derint.h:110
o2scl::vector_deriv_threept_tap
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
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::vector_deriv_threept
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
o2scl::vector_integ_extended8
fp_t vector_integ_extended8(size_t n, vec_t &v)
Integrate with an extended rule for 8 or more points.
Definition: vector_derint.h:279
o2scl::vector_deriv_fivept_tap
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...
Definition: vector_derint.h:135
o2scl::vector_integ_extended4
fp_t vector_integ_extended4(size_t n, vec_t &v)
Integrate with an extended rule for 4 or more points.
Definition: vector_derint.h:226
o2scl::vector_integ_threept
fp_t vector_integ_threept(size_t n, vec_t &v)
Integrate with an extended 3-point rule (extended Simpson's rule)
Definition: vector_derint.h:193
o2scl::vector_integ_trap
fp_t vector_integ_trap(size_t n, vec_t &v)
Integrate with an extended trapezoidal rule.
Definition: vector_derint.h:175
o2scl::vector_integ_durand
fp_t vector_integ_durand(size_t n, vec_t &v)
Integrate with Durand's rule for 4 or more points.
Definition: vector_derint.h:255

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