1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2005 Klaus Spanderen
5 Copyright (C) 2005 Gary Kennedy
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21/*! \file gaussianquadratures.hpp
22 \brief Integral of a 1-dimensional function using the Gauss quadratures
23*/
24
25#ifndef quantlib_gaussian_quadratures_hpp
26#define quantlib_gaussian_quadratures_hpp
27
28#include <ql/math/array.hpp>
29#include <ql/math/integrals/integral.hpp>
30#include <ql/math/integrals/gaussianorthogonalpolynomial.hpp>
31
32namespace QuantLib {
33 class GaussianOrthogonalPolynomial;
34
35 //! Integral of a 1-dimensional function using the Gauss quadratures method
36 /*! References:
37 Gauss quadratures and orthogonal polynomials
38
39 G.H. Gloub and J.H. Welsch: Calculation of Gauss quadrature rule.
40 Math. Comput. 23 (1986), 221-230
41
42 "Numerical Recipes in C", 2nd edition,
43 Press, Teukolsky, Vetterling, Flannery,
44
45 \test the correctness of the result is tested by checking it
46 against known good values.
47 */
48 class GaussianQuadrature {
49 public:
50 GaussianQuadrature(Size n,
51 const GaussianOrthogonalPolynomial& p);
52
53#if defined(__GNUC__) && (__GNUC__ >= 7)
54#pragma GCC diagnostic push
55#pragma GCC diagnostic ignored "-Wnoexcept-type"
56#endif
57
58 template <class F>
59 Real operator()(const F& f) const {
60 Real sum = 0.0;
61 for (Integer i = Integer(order())-1; i >= 0; --i) {
62 sum += w_[i] * f(x_[i]);
63 }
64 return sum;
65 }
66
67#if defined(__GNUC__) && (__GNUC__ >= 7)
68#pragma GCC diagnostic pop
69#endif
70
71 Size order() const { return x_.size(); }
72 const Array& weights() { return w_; }
73 const Array& x() { return x_; }
74
75 protected:
76 Array x_, w_;
77 };
78
79
80 //! generalized Gauss-Laguerre integration
81 /*! This class performs a 1-dimensional Gauss-Laguerre integration.
82 \f[
83 \int_{0}^{\inf} f(x) \mathrm{d}x
84 \f]
85 The weighting function is
86 \f[
87 w(x;s)=x^s \exp{-x}
88 \f]
89 and \f[ s > -1 \f]
90 */
91 class GaussLaguerreIntegration : public GaussianQuadrature {
92 public:
93 explicit GaussLaguerreIntegration(Size n, Real s = 0.0)
94 : GaussianQuadrature(n, GaussLaguerrePolynomial(s)) {}
95 };
96
97 //! generalized Gauss-Hermite integration
98 /*! This class performs a 1-dimensional Gauss-Hermite integration.
99 \f[
100 \int_{-\inf}^{\inf} f(x) \mathrm{d}x
101 \f]
102 The weighting function is
103 \f[
104 w(x;\mu)=|x|^{2\mu} \exp{-x*x}
105 \f]
106 and \f[ \mu > -0.5 \f]
107 */
108 class GaussHermiteIntegration : public GaussianQuadrature {
109 public:
110 explicit GaussHermiteIntegration(Size n, Real mu = 0.0)
111 : GaussianQuadrature(n, GaussHermitePolynomial(mu)) {}
112 };
113
114 //! Gauss-Jacobi integration
115 /*! This class performs a 1-dimensional Gauss-Jacobi integration.
116 \f[
117 \int_{-1}^{1} f(x) \mathrm{d}x
118 \f]
119 The weighting function is
120 \f[
121 w(x;\alpha,\beta)=(1-x)^\alpha (1+x)^\beta
122 \f]
123 */
124 class GaussJacobiIntegration : public GaussianQuadrature {
125 public:
126 GaussJacobiIntegration(Size n, Real alpha, Real beta)
127 : GaussianQuadrature(n, GaussJacobiPolynomial(alpha, beta)) {}
128 };
129
130 //! Gauss-Hyperbolic integration
131 /*! This class performs a 1-dimensional Gauss-Hyperbolic integration.
132 \f[
133 \int_{-\inf}^{\inf} f(x) \mathrm{d}x
134 \f]
135 The weighting function is
136 \f[
137 w(x)=1/cosh(x)
138 \f]
139 */
140 class GaussHyperbolicIntegration : public GaussianQuadrature {
141 public:
142 explicit GaussHyperbolicIntegration(Size n)
143 : GaussianQuadrature(n, GaussHyperbolicPolynomial()) {}
144 };
145
146 //! Gauss-Legendre integration
147 /*! This class performs a 1-dimensional Gauss-Legendre integration.
148 \f[
149 \int_{-1}^{1} f(x) \mathrm{d}x
150 \f]
151 The weighting function is
152 \f[
153 w(x)=1
154 \f]
155 */
156 class GaussLegendreIntegration : public GaussianQuadrature {
157 public:
158 explicit GaussLegendreIntegration(Size n)
159 : GaussianQuadrature(n, GaussJacobiPolynomial(0.0, 0.0)) {}
160 };
161
162 //! Gauss-Chebyshev integration
163 /*! This class performs a 1-dimensional Gauss-Chebyshev integration.
164 \f[
165 \int_{-1}^{1} f(x) \mathrm{d}x
166 \f]
167 The weighting function is
168 \f[
169 w(x)=(1-x^2)^{-1/2}
170 \f]
171 */
172 class GaussChebyshevIntegration : public GaussianQuadrature {
173 public:
174 explicit GaussChebyshevIntegration(Size n)
175 : GaussianQuadrature(n, GaussJacobiPolynomial(-0.5, -0.5)) {}
176 };
177
178 //! Gauss-Chebyshev integration (second kind)
179 /*! This class performs a 1-dimensional Gauss-Chebyshev integration.
180 \f[
181 \int_{-1}^{1} f(x) \mathrm{d}x
182 \f]
183 The weighting function is
184 \f[
185 w(x)=(1-x^2)^{1/2}
186 \f]
187 */
188 class GaussChebyshev2ndIntegration : public GaussianQuadrature {
189 public:
190 explicit GaussChebyshev2ndIntegration(Size n)
191 : GaussianQuadrature(n, GaussJacobiPolynomial(0.5, 0.5)) {}
192 };
193
194 //! Gauss-Gegenbauer integration
195 /*! This class performs a 1-dimensional Gauss-Gegenbauer integration.
196 \f[
197 \int_{-1}^{1} f(x) \mathrm{d}x
198 \f]
199 The weighting function is
200 \f[
201 w(x)=(1-x^2)^{\lambda-1/2}
202 \f]
203 */
204 class GaussGegenbauerIntegration : public GaussianQuadrature {
205 public:
206 GaussGegenbauerIntegration(Size n, Real lambda)
207 : GaussianQuadrature(n, GaussJacobiPolynomial(lambda-0.5, lambda-0.5))
208 {}
209 };
210
211
212 namespace detail {
213 template <class Integration>
214 class GaussianQuadratureIntegrator: public Integrator {
215 public:
216 explicit GaussianQuadratureIntegrator(Size n);
217
218 ext::shared_ptr<Integration> getIntegration() const { return integration_; }
219
220 private:
221 Real integrate(const ext::function<Real (Real)>& f,
222 Real a,
223 Real b) const override;
224
225 const ext::shared_ptr<Integration> integration_;
226 };
227 }
228
229 typedef detail::GaussianQuadratureIntegrator<GaussLegendreIntegration>
230 GaussLegendreIntegrator;
231
232 typedef detail::GaussianQuadratureIntegrator<GaussChebyshevIntegration>
233 GaussChebyshevIntegrator;
234
235 typedef detail::GaussianQuadratureIntegrator<GaussChebyshev2ndIntegration>
236 GaussChebyshev2ndIntegrator;
237
238 //! tabulated Gauss-Legendre quadratures
239 class TabulatedGaussLegendre {
240 public:
241 explicit TabulatedGaussLegendre(Size n = 20) { order(n); }
242 template <class F>
243 Real operator() (const F& f) const {
244 QL_ASSERT(w_ != nullptr, "Null weights");
245 QL_ASSERT(x_ != nullptr, "Null abscissas");
246 Size startIdx;
247 Real val;
248
249 const Size isOrderOdd = order_ & 1;
250
251 if (isOrderOdd) {
252 QL_ASSERT((n_>0), "assume at least 1 point in quadrature");
253 val = w_[0]*f(x_[0]);
254 startIdx=1;
255 } else {
256 val = 0.0;
257 startIdx=0;
258 }
259
260 for (Size i=startIdx; i<n_; ++i) {
261 val += w_[i]*f( x_[i]);
262 val += w_[i]*f(-x_[i]);
263 }
264 return val;
265 }
266
267 void order(Size);
268 Size order() const { return order_; }
269
270 private:
271 Size order_;
272
273 const Real* w_;
274 const Real* x_;
275 Size n_;
276
277 static const Real w6[3];
278 static const Real x6[3];
279 static const Size n6;
280
281 static const Real w7[4];
282 static const Real x7[4];
283 static const Size n7;
284
285 static const Real w12[6];
286 static const Real x12[6];
287 static const Size n12;
288
289 static const Real w20[10];
290 static const Real x20[10];
291 static const Size n20;
292 };
293
294}
295
296#endif
297

source code of quantlib/ql/math/integrals/gaussianquadratures.hpp