inte_gauss_cern.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 /** \file inte_gauss_cern.h
24  \brief File defining \ref o2scl::inte_gauss_cern
25 */
26 #ifndef O2SCL_CERN_GAUSS_H
27 #define O2SCL_CERN_GAUSS_H
28
29 #include <o2scl/inte.h>
30
31 #ifndef DOXYGEN_NO_O2NS
32 namespace o2scl {
33 #endif
34
35  /** \brief Gaussian quadrature (CERNLIB)
36
37  For any interval \f$(a,b) \f$, we define \f$g_8(a,b) \f$ and
38  \f$g_{16}(a,b) \f$ to be the 8- and 16-point Gaussian
39  quadrature approximations to
40  \f[
41  I=\int_a^b f(x)~dx
42  \f]
43  and define
44  \f[
45  r(a,b)=\frac{|g_{16}(a,b)-g_{8}(a,b)|}{1+g_{16}(a,b)}
46  \f]
47  The function integ() returns \f$G \f$ given by
48  \f[
49  G=\sum_{i=1}^{k} g_{16}(x_{i-1},x_i)
50  \f]
51  where \f$x_0=a \f$ and \f$x_k=b \f$ and the subdivision
52  points \f$x_i \f$ are given by
53  \f[
54  x_i=x_{i-1}+\lambda(B-x_{i-1})
55  \f]
56  where \f$\lambda \f$ is the first number in the sequence \f$57 1,\frac{1}{2},\frac{1}{4},... \f$ for which
58  \f[
59  r(x_{i-1},x_i)<\mathrm{eps}.
60  \f]
61  If, at any stage, the ratio
62  \f[
63  q=\left| \frac{x_i-x_{i-1}}{b-a} \right|
64  \f]
65  is so small so that \f$1+0.005 q \f$ is indistinguishable from
66  unity, then the accuracy is required is not reachable and
67  the error handler is called.
68
69  Unless there is severe cancellation, inte::tol_rel may be
70  considered as specifying a bound on the relative error of the
71  integral in the case that \f$|I|>1 \f$ and an absolute error if
72  \f$|I|<1 \f$. More precisely, if \f$k \f$ is the number of
73  subintervals from above, and if
74  \f[
75  I_{abs} = \int_a^b |f(x)|~dx
76  \f]
77  then
78  \f[
79  \frac{|G-I|}{I_{abs}+k}<\mathrm{tol_rel}
80  \f]
81  will nearly always be true when no error is returned. For
82  functions with no singualarities in the interval, the accuracy
83  will usually be higher than this.
84
85  If the desired accuracy is not achieved, the integration
86  functions will call the error handler and return the best guess,
87  unless \ref inte::err_nonconv is false, in which case the error
88  handler is not called.
89
90  This function is based on the CERNLIB routines GAUSS and
91  DGAUSS which are documented at
92  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d103/top.html
93
94  \future Allow user to change \c cst?
95  */
96  template<class func_t=funct11> class inte_gauss_cern : public inte<func_t> {
97
98  public:
99
100  inte_gauss_cern() {
101  x[0]=0.96028985649753623;
102  x[1]=0.79666647741362674;
103  x[2]=0.52553240991632899;
104  x[3]=0.18343464249564980;
105  x[4]=0.98940093499164993;
106  x[5]=0.94457502307323258;
107  x[6]=0.86563120238783175;
108  x[7]=0.75540440835500303;
109  x[8]=0.61787624440264375;
110  x[9]=0.45801677765722739;
111  x[10]=0.28160355077925891;
112  x[11]=0.95012509837637440e-1;
113
114  w[0]=0.10122853629037626;
115  w[1]=0.22238103445337447;
116  w[2]=0.31370664587788729;
117  w[3]=0.36268378337836198;
118  w[4]=0.27152459411754095e-1;
119  w[5]=0.62253523938647893e-1;
120  w[6]=0.95158511682492785e-1;
121  w[7]=0.12462897125553387;
122  w[8]=0.14959598881657673;
123  w[9]=0.16915651939500254;
124  w[10]=0.18260341504492359;
125  w[11]=0.18945061045506850;
126  }
127
128  /** \brief Integrate function \c func from \c a to \c b.
129  */
130  virtual int integ_err(func_t &func, double a, double b,
131  double &res, double &err) {
132
133  double y1, y2;
134  err=0.0;
135
136  size_t itx=0;
137
138  int i;
139  bool loop=true, loop2=false;
140  static const double cst=0.005;
141  double h=0.0;
142  if (b==a) {
143  res=0.0;
144  return o2scl::success;
145  }
146  double cnst=cst/(b-a);
147  double aa=0.0, bb=a;
148  while (loop==true || loop2==true) {
149  itx++;
150  if (loop==true) {
151  aa=bb;
152  bb=b;
153  }
154  double c1=(bb+aa)/2.0;
155  double c2=(bb-aa)/2.0;
156  double s8=0.0;
157  for(i=0;i<4;i++) {
158  double u=c2*x[i];
159  y1=func(c1+u);
160  y2=func(c1-u);
161  s8+=w[i]*(y1+y2);
162  }
163  double s16=0.0;
164  for(i=4;i<12;i++) {
165  double u=c2*x[i];
166  y1=func(c1+u);
167  y2=func(c1-u);
168  s16+=w[i]*(y1+y2);
169  }
170  s16*=c2;
171
172  loop=false;
173  loop2=false;
174  if (fabs(s16-c2*s8)<this->tol_rel*(1.0+fabs(s16))) {
175  h+=s16;
176  if (bb!=b) loop=true;
177  } else {
178  bb=c1;
179  if (1.0+cnst*fabs(c2)!=1.0) {
180  loop2=true;
181  } else {
182  this->last_iter=itx;
183  O2SCL_CONV2_RET("Failed to reach required accuracy in cern_",
184  "gauss::integ().",exc_efailed,this->err_nonconv);
185  }
186  }
187  }
188  this->last_iter=itx;
189  res=h;
190  return o2scl::success;
191  }
192
193  protected:
194
195 #ifndef DOXYGEN_INTERNAL
196
197  /** \name Integration constants (Set in the constructor)
198  */
199  //@{
200  double x[12], w[12];
201  //@}
202
203 #endif
204
205  };
206
207 #ifndef DOXYGEN_NO_O2NS
208 }
209 #endif
210
211 #endif
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
virtual int integ_err(func_t &func, double a, double b, double &res, double &err)
Integrate function func from a to b.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:63
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:81
generic failure
Definition: err_hnd.h:61
Base integration class [abstract base].
Definition: inte.h:44
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:68