QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
gausslobattointegral.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) 2008 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file gausslabottointegral.cpp
21 \brief integral of a one-dimensional function using the adaptive
22 Gauss-Lobatto integral
23*/
24
26#include <algorithm>
27
28namespace QuantLib {
29
30 const Real GaussLobattoIntegral::alpha_ = std::sqrt(2.0/3.0);
31 const Real GaussLobattoIntegral::beta_ = 1.0/std::sqrt(5.0);
32 const Real GaussLobattoIntegral::x1_ = 0.94288241569547971906;
33 const Real GaussLobattoIntegral::x2_ = 0.64185334234578130578;
34 const Real GaussLobattoIntegral::x3_ = 0.23638319966214988028;
35
37 Real absAccuracy,
38 Real relAccuracy,
39 bool useConvergenceEstimate)
40 : Integrator(absAccuracy, maxIterations),
41 relAccuracy_(relAccuracy),
42 useConvergenceEstimate_(useConvergenceEstimate) {
43 }
44
46 const ext::function<Real (Real)>& f,
47 Real a, Real b) const {
48
50 const Real calcAbsTolerance = calculateAbsTolerance(f, a, b);
51
53 return adaptivGaussLobattoStep(f, a, b, f(a), f(b), calcAbsTolerance);
54 }
55
57 const ext::function<Real (Real)>& f,
58 Real a, Real b) const {
59
60
61 Real relTol = std::max(relAccuracy_, QL_EPSILON);
62
63 const Real m = (a+b)/2;
64 const Real h = (b-a)/2;
65 const Real y1 = f(a);
66 const Real y3 = f(m-alpha_*h);
67 const Real y5 = f(m-beta_*h);
68 const Real y7 = f(m);
69 const Real y9 = f(m+beta_*h);
70 const Real y11= f(m+alpha_*h);
71 const Real y13= f(b);
72
73 const Real f1 = f(m-x1_*h);
74 const Real f2 = f(m+x1_*h);
75 const Real f3 = f(m-x2_*h);
76 const Real f4 = f(m+x2_*h);
77 const Real f5 = f(m-x3_*h);
78 const Real f6 = f(m+x3_*h);
79
80 Real acc=h*(0.0158271919734801831*(y1+y13)
81 +0.0942738402188500455*(f1+f2)
82 +0.1550719873365853963*(y3+y11)
83 +0.1888215739601824544*(f3+f4)
84 +0.1997734052268585268*(y5+y9)
85 +0.2249264653333395270*(f5+f6)
86 +0.2426110719014077338*y7);
87
89 if (acc == 0.0 && ( f1 != 0.0 || f2 != 0.0 || f3 != 0.0
90 || f4 != 0.0 || f5 != 0.0 || f6 != 0.0)) {
91 QL_FAIL("can not calculate absolute accuracy "
92 "from relative accuracy");
93 }
94
95 Real r = 1.0;
97 const Real integral2 = (h/6)*(y1+y13+5*(y5+y9));
98 const Real integral1 = (h/1470)*(77*(y1+y13)+432*(y3+y11)+
99 625*(y5+y9)+672*y7);
100
101 if (std::fabs(integral2-acc) != 0.0)
102 r = std::fabs(integral1-acc)/std::fabs(integral2-acc);
103 if (r == 0.0 || r > 1.0)
104 r = 1.0;
105 }
106
107 if (relAccuracy_ != Null<Real>())
108 return std::min(absoluteAccuracy(), acc*relTol)/(r*QL_EPSILON);
109 else {
110 return absoluteAccuracy()/(r*QL_EPSILON);
111 }
112 }
113
115 const ext::function<Real (Real)>& f,
116 Real a, Real b, Real fa, Real fb,
117 Real acc) const {
119 "max number of iterations reached");
120
121 const Real h=(b-a)/2;
122 const Real m=(a+b)/2;
123
124 const Real mll=m-alpha_*h;
125 const Real ml =m-beta_*h;
126 const Real mr =m+beta_*h;
127 const Real mrr=m+alpha_*h;
128
129 const Real fmll= f(mll);
130 const Real fml = f(ml);
131 const Real fm = f(m);
132 const Real fmr = f(mr);
133 const Real fmrr= f(mrr);
135
136 const Real integral2=(h/6)*(fa+fb+5*(fml+fmr));
137 const Real integral1=(h/1470)*(77*(fa+fb)
138 +432*(fmll+fmrr)+625*(fml+fmr)+672*fm);
139
140 // avoid 80 bit logic on x86 cpu
141 const Real dist = acc + (integral1-integral2);
142 if(dist==acc || mll<=a || b<=mrr) {
143 QL_REQUIRE(m>a && b>m,"Interval contains no more machine number");
144 return integral1;
145 }
146 else {
147 return adaptivGaussLobattoStep(f,a,mll,fa,fmll,acc)
148 + adaptivGaussLobattoStep(f,mll,ml,fmll,fml,acc)
149 + adaptivGaussLobattoStep(f,ml,m,fml,fm,acc)
150 + adaptivGaussLobattoStep(f,m,mr,fm,fmr,acc)
151 + adaptivGaussLobattoStep(f,mr,mrr,fmr,fmrr,acc)
152 + adaptivGaussLobattoStep(f,mrr,b,fmrr,fb,acc);
153 }
154 }
155}
GaussLobattoIntegral(Size maxIterations, Real absAccuracy, Real relAccuracy=Null< Real >(), bool useConvergenceEstimate=true)
Real integrate(const ext::function< Real(Real)> &f, Real a, Real b) const override
Real adaptivGaussLobattoStep(const ext::function< Real(Real)> &f, Real a, Real b, Real fa, Real fb, Real is) const
Real calculateAbsTolerance(const ext::function< Real(Real)> &f, Real a, Real b) const
Real absoluteAccuracy() const
Definition: integral.cpp:43
Size numberOfEvaluations() const
Definition: integral.cpp:59
Size maxEvaluations() const
Definition: integral.cpp:47
void increaseNumberOfEvaluations(Size increase) const
Definition: integral.cpp:67
void setNumberOfEvaluations(Size evaluations) const
Definition: integral.cpp:63
template class providing a null value for a given type.
Definition: null.hpp:76
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#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 one-dimensional function using the adaptive Gauss-Lobatto integral
#define QL_EPSILON
Definition: qldefines.hpp:178
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
ext::shared_ptr< YieldTermStructure > r