QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
gaussianquadratures.cpp
Go to the documentation of this file.
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#include <ql/utilities/null.hpp>
29
30namespace QuantLib {
31
33 Size n,
34 const GaussianOrthogonalPolynomial& orthPoly)
35 : x_(n), w_(n) {
36
37 // set-up matrix to compute the roots and the weights
38 Array e(n-1);
39
40 Size i;
41 for (i=1; i < n; ++i) {
42 x_[i] = orthPoly.alpha(i);
43 e[i-1] = std::sqrt(orthPoly.beta(i));
44 }
45 x_[0] = orthPoly.alpha(0);
46
48 x_, e,
51
52 x_ = tqr.eigenvalues();
53 const Matrix& ev = tqr.eigenvectors();
54
55 Real mu_0 = orthPoly.mu_0();
56 for (i=0; i<n; ++i) {
57 w_[i] = mu_0*ev[0][i]*ev[0][i] / orthPoly.w(x_[i]);
58 }
59 }
60
61
62 namespace detail {
63 template <class Integration>
65 Size n)
66 : Integrator(Null<Real>(), n),
67 integration_(ext::make_shared<Integration>(n))
68 { }
69
70 template <class Integration>
72 const ext::function<Real (Real)>& f, Real a, Real b) const {
73
74 const Real c1 = 0.5*(b-a);
75 const Real c2 = 0.5*(a+b);
76
77 return c1*integration_->operator()(
78 [c1, c2, f](Real x) { return f(c1*x+c2);});
79 }
80
84 }
85
87 switch(order) {
88 case(6):
90 break;
91 case(7):
93 break;
94 case(12):
96 break;
97 case(20):
99 break;
100 default:
101 QL_FAIL("order " << order << " not supported");
102 }
103 }
104
105
106 // Abscissas and Weights from Abramowitz and Stegun
107
108 /* order 6 */
109 const Real TabulatedGaussLegendre::x6[3] = { 0.238619186083197,
110 0.661209386466265,
111 0.932469514203152 };
112
113 const Real TabulatedGaussLegendre::w6[3] = { 0.467913934572691,
114 0.360761573048139,
115 0.171324492379170 };
116
118
119 /* order 7 */
120 const Real TabulatedGaussLegendre::x7[4] = { 0.000000000000000,
121 0.405845151377397,
122 0.741531185599394,
123 0.949107912342759 };
124
125 const Real TabulatedGaussLegendre::w7[4] = { 0.417959183673469,
126 0.381830050505119,
127 0.279705391489277,
128 0.129484966168870 };
129
131
132 /* order 12 */
133 const Real TabulatedGaussLegendre::x12[6] = { 0.125233408511469,
134 0.367831498998180,
135 0.587317954286617,
136 0.769902674194305,
137 0.904117256370475,
138 0.981560634246719 };
139
140 const Real TabulatedGaussLegendre::w12[6] = { 0.249147045813403,
141 0.233492536538355,
142 0.203167426723066,
143 0.160078328543346,
144 0.106939325995318,
145 0.047175336386512 };
146
148
149 /* order 20 */
150 const Real TabulatedGaussLegendre::x20[10] = { 0.076526521133497,
151 0.227785851141645,
152 0.373706088715420,
153 0.510867001950827,
154 0.636053680726515,
155 0.746331906460151,
156 0.839116971822219,
157 0.912234428251326,
158 0.963971927277914,
159 0.993128599185095 };
160
161 const Real TabulatedGaussLegendre::w20[10] = { 0.152753387130726,
162 0.149172986472604,
163 0.142096109318382,
164 0.131688638449177,
165 0.118194531961518,
166 0.101930119817240,
167 0.083276741576704,
168 0.062672048334109,
169 0.040601429800387,
170 0.017614007139152 };
171
173
174}
1-D array used in linear algebra.
Definition: array.hpp:52
orthogonal polynomial for Gaussian quadratures
virtual Real alpha(Size i) const =0
virtual Real beta(Size i) const =0
virtual Real w(Real x) const =0
GaussianQuadrature(Size n, const GaussianOrthogonalPolynomial &p)
Matrix used in linear algebra.
Definition: matrix.hpp:41
template class providing a null value for a given type.
Definition: null.hpp:76
tridiag. QR eigen decomposition with explicite shift aka Wilkinson
Real integrate(const ext::function< Real(Real)> &f, Real a, Real b) const override
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
ext::function< Real(Real)> b
Integral of a 1-dimensional function using the Gauss quadratures.
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
null values
Eigenvalues/eigenvectors of a real symmetric matrix.
tridiag. QR eigen decompositions with implicit shift
std::uint64_t x_