inte_cauchy_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_cauchy_cern.h
24  \brief File defining \ref o2scl::inte_cauchy_cern
25 */
26 #ifndef O2SCL_CERN_CAUCHY_H
27 #define O2SCL_CERN_CAUCHY_H
28 
29 #include <o2scl/inte.h>
30 #include <o2scl/inte_gauss_cern.h>
31 
32 #ifndef DOXYGEN_NO_O2NS
33 namespace o2scl {
34 #endif
35 
36  /** \brief Cauchy principal value integration (CERNLIB)
37 
38  The location of the singularity must be specified before-hand in
39  inte_cauchy_cern::s, and the singularity must not be at one of the
40  endpoints. Note that when integrating a function of the form \f$
41  \frac{f(x)}{(x-s)} \f$, the denominator \f$ (x-s) \f$ must be
42  specified in the argument \c func to integ(). This is different
43  from how the \ref inte_qawc_gsl operates.
44 
45  The method from \ref Longman58 is used for the decomposition of
46  the integral, and the resulting integrals are computed using
47  a user-specified base integration object.
48 
49  The uncertainty in the integral is not calculated, and is always
50  given as zero. The default base integration object is of type
51  \ref inte_gauss_cern. This is the CERNLIB default, but can be
52  modified by calling set_inte(). If the singularity is outside
53  the region of integration, then the result from the base
54  integration object is returned without calling the error
55  handler.
56 
57  Possible errors for integ() and integ_err():
58  - exc_einval - Singularity is on an endpoint
59  - exc_efailed - Couldn't reach requested accuracy
60  (occurs only if inte::err_nonconv is true)
61 
62  This function is based on the CERNLIB routines RCAUCH and
63  DCAUCH which are documented at
64  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d104/top.html
65  */
66  template<class func_t> class inte_cauchy_cern :
67  public inte<func_t> {
68 
69  public:
70 
72  x[0]=0.96028985649753623;
73  x[1]=0.79666647741362674;
74  x[2]=0.52553240991632899;
75  x[3]=0.18343464249564980;
76  x[4]=0.98940093499164993;
77  x[5]=0.94457502307323258;
78  x[6]=0.86563120238783175;
79  x[7]=0.75540440835500303;
80  x[8]=0.61787624440264375;
81  x[9]=0.45801677765722739;
82  x[10]=0.28160355077925891;
83  x[11]=0.95012509837637440e-1;
84 
85  w[0]=0.10122853629037626;
86  w[1]=0.22238103445337447;
87  w[2]=0.31370664587788729;
88  w[3]=0.36268378337836198;
89  w[4]=0.27152459411754095e-1;
90  w[5]=0.62253523938647893e-1;
91  w[6]=0.95158511682492785e-1;
92  w[7]=0.12462897125553387;
93  w[8]=0.14959598881657673;
94  w[9]=0.16915651939500254;
95  w[10]=0.18260341504492359;
96  w[11]=0.18945061045506850;
97 
98  it=&def_inte;
99  }
100 
101  /// The singularity (must be set before calling integ() or integ_err())
102  double s;
103 
104  /** \brief Set the base integration object to use (default is \ref
105  inte_cauchy_cern::def_inte of type \ref inte_gauss_cern)
106  */
108  it=&i;
109  return 0;
110  }
111 
112  /** \brief Integrate function \c func from \c a to \c b
113  */
114  virtual int integ_err(func_t &func, double a, double b,
115  double &res, double &err) {
116  double y1, y2, y3, y4;
117  size_t itx=0;
118 
119  err=0.0;
120 
121  if (s==a || s==b) {
122  O2SCL_CONV2_RET("Singularity on boundary in ",
123  "inte_cauchy_cern::integ_err().",
124  exc_einval,this->err_nonconv);
125  } else if ((s<a && s<b) || (s>a && s>b)) {
126  return it->integ_err(func,a,b,res,err);
127  }
128  double h, b0;
129  if (2.0*s<a+b) {
130  h=it->integ(func,2.0*s-a,b);
131  b0=s-a;
132  } else {
133  h=it->integ(func,a,2.0*s-b);
134  b0=b-s;
135  }
136  double c=0.005/b0;
137  double bb=0.0;
138  bool loop1=true;
139  while(loop1==true) {
140  itx++;
141  double s8, s16;
142  double aa=bb;
143  bb=b0;
144  bool loop2=true;
145  while(loop2==true) {
146  double c1=(bb+aa)/2.0;
147  double c2=(bb-aa)/2.0;
148  double c3=s+c1;
149  double c4=s-c1;
150  s8=0.0;
151  s16=0.0;
152  for(int i=0;i<4;i++) {
153  double u=c2*x[i];
154  y1=func(c3+u);
155  y2=func(c4-u);
156  y3=func(c3-u);
157  y4=func(c4+u);
158  s8+=w[i]*((y1+y2)+(y3+y4));
159  }
160  s8*=c2;
161  for(int i=4;i<12;i++) {
162  double u=c2*x[i];
163  y1=func(c3+u);
164  y2=func(c4-u);
165  y3=func(c3-u);
166  y4=func(c4+u);
167  s16+=w[i]*((y1+y2)+(y3+y4));
168  }
169  s16*=c2;
170 
171  if (fabs(s16-s8)<=this->tol_rel*(1.0+fabs(s16))) {
172  loop2=false;
173  } else {
174  bb=c1;
175  if ((1.0+fabs(c*c2))==1.0) {
176  loop2=false;
177  } else {
178  this->last_iter=itx;
179  O2SCL_CONV2_RET("Couldn't reach required accuracy in cern_",
180  "cauchy::integ()",exc_efailed,this->err_nonconv);
181  }
182  }
183  }
184  h+=s16;
185  if (bb==b0) loop1=false;
186  }
187  this->last_iter=itx;
188  res=h;
189  return o2scl::success;
190  }
191 
192  /// Default integration object
194 
195  protected:
196 
197 #ifndef DOXYGEN_INTERNAL
198 
199  /** \name Integration constants
200  */
201  //@{
202  double x[12], w[12];
203  //@}
204 
205  /// The base integration object
207 
208 #endif
209 
210  };
211 
212 #ifndef DOXYGEN_NO_O2NS
213 }
214 #endif
215 
216 #endif
int set_inte(inte< func_t > &i)
Set the base integration object to use (default is inte_cauchy_cern::def_inte of type inte_gauss_cern...
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.
invalid argument supplied by user
Definition: err_hnd.h:59
inte< func_t > * it
The base integration object.
#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
double s
The singularity (must be set before calling integ() or integ_err())
Base integration class [abstract base].
Definition: inte.h:44
inte_gauss_cern< func_t > def_inte
Default integration object.
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:68
Cauchy principal value integration (CERNLIB)
Gaussian quadrature (CERNLIB)
Success.
Definition: err_hnd.h:47

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