inte_cauchy_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_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 - Could not reach requested accuracy
60  (occurs only if inte::err_nonconv is true)
61 
62  \note Currently \o2 supports only types \c double and
63  \c long \c double for the floating point type \c fp_t .
64 
65  This function is based on the CERNLIB routines RCAUCH and
66  DCAUCH which are documented at
67  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d104/top.html
68  */
69  template<class func_t, class fp_t=double,
70  const fp_t x[]=inte_gauss_cern_x_double,
71  const fp_t w[]=inte_gauss_cern_w_double>
72  class inte_cauchy_cern : public inte<func_t,fp_t> {
73 
74  public:
75 
77  it=&def_inte;
78  }
79 
80  /// The singularity (must be set before calling integ() or integ_err())
81  fp_t s;
82 
83  /** \brief Set the base integration object to use (default is \ref
84  inte_cauchy_cern::def_inte of type \ref inte_gauss_cern)
85  */
87  it=&i;
88  return 0;
89  }
90 
91  /** \brief Integrate function \c func from \c a to \c b
92  */
93  virtual int integ_err(func_t &func, fp_t a, fp_t b,
94  fp_t &res, fp_t &err) {
95  fp_t y1, y2, y3, y4;
96  size_t itx=0;
97 
98  err=0.0;
99 
100  if (s==a || s==b) {
101  O2SCL_CONV2_RET("Singularity on boundary in ",
102  "inte_cauchy_cern::integ_err().",
103  exc_einval,this->err_nonconv);
104  } else if ((s<a && s<b) || (s>a && s>b)) {
105  return it->integ_err(func,a,b,res,err);
106  }
107  fp_t h, b0;
108  if (2.0*s<a+b) {
109  h=it->integ(func,2.0*s-a,b);
110  b0=s-a;
111  } else {
112  h=it->integ(func,a,2.0*s-b);
113  b0=b-s;
114  }
115  fp_t c=0.005/b0;
116  fp_t bb=0.0;
117  bool loop1=true;
118  while(loop1==true) {
119  itx++;
120  fp_t s8, s16;
121  fp_t aa=bb;
122  bb=b0;
123  bool loop2=true;
124  while(loop2==true) {
125  fp_t c1=(bb+aa)/2.0;
126  fp_t c2=(bb-aa)/2.0;
127  fp_t c3=s+c1;
128  fp_t c4=s-c1;
129  s8=0.0;
130  s16=0.0;
131  for(int i=0;i<4;i++) {
132  fp_t u=c2*x[i];
133  y1=func(c3+u);
134  y2=func(c4-u);
135  y3=func(c3-u);
136  y4=func(c4+u);
137  s8+=w[i]*((y1+y2)+(y3+y4));
138  }
139  s8*=c2;
140  for(int i=4;i<12;i++) {
141  fp_t u=c2*x[i];
142  y1=func(c3+u);
143  y2=func(c4-u);
144  y3=func(c3-u);
145  y4=func(c4+u);
146  s16+=w[i]*((y1+y2)+(y3+y4));
147  }
148  s16*=c2;
149 
150  if (std::abs(s16-s8)<=this->tol_rel*(1.0+std::abs(s16))) {
151  loop2=false;
152  } else {
153  bb=c1;
154  if ((1.0+std::abs(c*c2))==1.0) {
155  loop2=false;
156  } else {
157  this->last_iter=itx;
158  O2SCL_CONV2_RET("Could not reach required accuracy in cern_",
159  "cauchy::integ()",exc_efailed,this->err_nonconv);
160  }
161  }
162  }
163  h+=s16;
164  if (bb==b0) loop1=false;
165  }
166  this->last_iter=itx;
167  res=h;
168  return o2scl::success;
169  }
170 
171  /// Default integration object
173 
174  protected:
175 
176 #ifndef DOXYGEN_INTERNAL
177 
178  /// The base integration object
180 
181 #endif
182 
183  };
184 
185 #ifndef DOXYGEN_NO_O2NS
186 }
187 #endif
188 
189 #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
int set_inte(inte< func_t, fp_t > &i)
Set the base integration object to use (default is inte_cauchy_cern::def_inte of type inte_gauss_cern...
fp_t tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:75
inte_gauss_cern< func_t, fp_t > def_inte
Default integration object.
invalid argument supplied by user
Definition: err_hnd.h:59
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:303
fp_t s
The singularity (must be set before calling integ() or integ_err())
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
inte< func_t, fp_t > * it
The base integration object.
static constexpr double inte_gauss_cern_w_double[12]
Integration weights for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in double precision...
Cauchy principal value integration (CERNLIB)
Success.
Definition: err_hnd.h:47
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.

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