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

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