1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2006 François du Vignaud
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#include <ql/math/integrals/kronrodintegral.hpp>
21#include <ql/types.hpp>
22
23namespace QuantLib {
24
25 static Real rescaleError(Real err,
26 const Real resultAbs,
27 const Real resultAsc) {
28 err = std::fabs(x: err) ;
29 if (resultAsc != 0 && err != 0){
30 Real scale = std::pow(x: (200 * err / resultAsc), y: 1.5) ;
31 if (scale < 1)
32 err = resultAsc * scale ;
33 else
34 err = resultAsc ;
35 }
36 if (resultAbs > QL_MIN_POSITIVE_REAL / (50 * QL_EPSILON )){
37 Real min_err = 50 * QL_EPSILON * resultAbs ;
38 if (min_err > err)
39 err = min_err ;
40 }
41 return err ;
42 }
43
44 /* Gauss-Kronrod-Patterson quadrature coefficients for use in
45 quadpack routine qng. These coefficients were calculated with
46 101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov
47 1981. */
48
49 /* x1, abscissae common to the 10-, 21-, 43- and 87-point rule */
50 static const Real x1[5] = {
51 0.973906528517171720077964012084452,
52 0.865063366688984510732096688423493,
53 0.679409568299024406234327365114874,
54 0.433395394129247190799265943165784,
55 0.148874338981631210884826001129720
56 } ;
57
58 /* w10, weights of the 10-point formula */
59 static const Real w10[5] = {
60 0.066671344308688137593568809893332,
61 0.149451349150580593145776339657697,
62 0.219086362515982043995534934228163,
63 0.269266719309996355091226921569469,
64 0.295524224714752870173892994651338
65 } ;
66
67 /* x2, abscissae common to the 21-, 43- and 87-point rule */
68 static const Real x2[5] = {
69 0.995657163025808080735527280689003,
70 0.930157491355708226001207180059508,
71 0.780817726586416897063717578345042,
72 0.562757134668604683339000099272694,
73 0.294392862701460198131126603103866
74 } ;
75
76 /* w21a, weights of the 21-point formula for abscissae x1 */
77 static const Real w21a[5] = {
78 0.032558162307964727478818972459390,
79 0.075039674810919952767043140916190,
80 0.109387158802297641899210590325805,
81 0.134709217311473325928054001771707,
82 0.147739104901338491374841515972068
83 } ;
84
85 /* w21b, weights of the 21-point formula for abscissae x2 */
86 static const Real w21b[6] = {
87 0.011694638867371874278064396062192,
88 0.054755896574351996031381300244580,
89 0.093125454583697605535065465083366,
90 0.123491976262065851077958109831074,
91 0.142775938577060080797094273138717,
92 0.149445554002916905664936468389821
93 } ;
94
95 /* x3, abscissae common to the 43- and 87-point rule */
96 static const Real x3[11] = {
97 0.999333360901932081394099323919911,
98 0.987433402908088869795961478381209,
99 0.954807934814266299257919200290473,
100 0.900148695748328293625099494069092,
101 0.825198314983114150847066732588520,
102 0.732148388989304982612354848755461,
103 0.622847970537725238641159120344323,
104 0.499479574071056499952214885499755,
105 0.364901661346580768043989548502644,
106 0.222254919776601296498260928066212,
107 0.074650617461383322043914435796506
108 } ;
109
110 /* w43a, weights of the 43-point formula for abscissae x1, x3 */
111 static const Real w43a[10] = {
112 0.016296734289666564924281974617663,
113 0.037522876120869501461613795898115,
114 0.054694902058255442147212685465005,
115 0.067355414609478086075553166302174,
116 0.073870199632393953432140695251367,
117 0.005768556059769796184184327908655,
118 0.027371890593248842081276069289151,
119 0.046560826910428830743339154433824,
120 0.061744995201442564496240336030883,
121 0.071387267268693397768559114425516
122 } ;
123
124 /* w43b, weights of the 43-point formula for abscissae x3 */
125 static const Real w43b[12] = {
126 0.001844477640212414100389106552965,
127 0.010798689585891651740465406741293,
128 0.021895363867795428102523123075149,
129 0.032597463975345689443882222526137,
130 0.042163137935191811847627924327955,
131 0.050741939600184577780189020092084,
132 0.058379395542619248375475369330206,
133 0.064746404951445885544689259517511,
134 0.069566197912356484528633315038405,
135 0.072824441471833208150939535192842,
136 0.074507751014175118273571813842889,
137 0.074722147517403005594425168280423
138 } ;
139
140 /* x4, abscissae of the 87-point rule */
141 static const Real x4[22] = {
142 0.999902977262729234490529830591582,
143 0.997989895986678745427496322365960,
144 0.992175497860687222808523352251425,
145 0.981358163572712773571916941623894,
146 0.965057623858384619128284110607926,
147 0.943167613133670596816416634507426,
148 0.915806414685507209591826430720050,
149 0.883221657771316501372117548744163,
150 0.845710748462415666605902011504855,
151 0.803557658035230982788739474980964,
152 0.757005730685495558328942793432020,
153 0.706273209787321819824094274740840,
154 0.651589466501177922534422205016736,
155 0.593223374057961088875273770349144,
156 0.531493605970831932285268948562671,
157 0.466763623042022844871966781659270,
158 0.399424847859218804732101665817923,
159 0.329874877106188288265053371824597,
160 0.258503559202161551802280975429025,
161 0.185695396568346652015917141167606,
162 0.111842213179907468172398359241362,
163 0.037352123394619870814998165437704
164 } ;
165
166 /* w87a, weights of the 87-point formula for abscissae x1, x2, x3 */
167 static const Real w87a[21] = {
168 0.008148377384149172900002878448190,
169 0.018761438201562822243935059003794,
170 0.027347451050052286161582829741283,
171 0.033677707311637930046581056957588,
172 0.036935099820427907614589586742499,
173 0.002884872430211530501334156248695,
174 0.013685946022712701888950035273128,
175 0.023280413502888311123409291030404,
176 0.030872497611713358675466394126442,
177 0.035693633639418770719351355457044,
178 0.000915283345202241360843392549948,
179 0.005399280219300471367738743391053,
180 0.010947679601118931134327826856808,
181 0.016298731696787335262665703223280,
182 0.021081568889203835112433060188190,
183 0.025370969769253827243467999831710,
184 0.029189697756475752501446154084920,
185 0.032373202467202789685788194889595,
186 0.034783098950365142750781997949596,
187 0.036412220731351787562801163687577,
188 0.037253875503047708539592001191226
189 } ;
190
191 /* w87b, weights of the 87-point formula for abscissae x4 */
192 static const Real w87b[23] = {
193 0.000274145563762072350016527092881,
194 0.001807124155057942948341311753254,
195 0.004096869282759164864458070683480,
196 0.006758290051847378699816577897424,
197 0.009549957672201646536053581325377,
198 0.012329447652244853694626639963780,
199 0.015010447346388952376697286041943,
200 0.017548967986243191099665352925900,
201 0.019938037786440888202278192730714,
202 0.022194935961012286796332102959499,
203 0.024339147126000805470360647041454,
204 0.026374505414839207241503786552615,
205 0.028286910788771200659968002987960,
206 0.030052581128092695322521110347341,
207 0.031646751371439929404586051078883,
208 0.033050413419978503290785944862689,
209 0.034255099704226061787082821046821,
210 0.035262412660156681033782717998428,
211 0.036076989622888701185500318003895,
212 0.036698604498456094498018047441094,
213 0.037120549269832576114119958413599,
214 0.037334228751935040321235449094698,
215 0.037361073762679023410321241766599
216 } ;
217
218 void GaussKronrodNonAdaptive::setRelativeAccuracy(Real relativeAccuracy) {
219 relativeAccuracy_ = relativeAccuracy;
220 }
221
222
223 Real GaussKronrodNonAdaptive::relativeAccuracy() const {
224 return relativeAccuracy_;
225 }
226
227 GaussKronrodNonAdaptive::GaussKronrodNonAdaptive(Real absoluteAccuracy,
228 Size maxEvaluations,
229 Real relativeAccuracy)
230 : Integrator(absoluteAccuracy, maxEvaluations),
231 relativeAccuracy_(relativeAccuracy) {}
232
233 Real
234 GaussKronrodNonAdaptive::integrate(const ext::function<Real (Real)>& f,
235 Real a,
236 Real b) const {
237 Real result;
238 //Size neval;
239 Real fv1[5], fv2[5], fv3[5], fv4[5];
240 Real savfun[21]; /* array of function values which have been computed */
241 Real res10, res21, res43, res87; /* 10, 21, 43 and 87 point results */
242 Real err;
243 Real resAbs; /* approximation to the integral of abs(f) */
244 Real resasc; /* approximation to the integral of abs(f-i/(b-a)) */
245 int k ;
246
247 QL_REQUIRE(a<b, "b must be greater than a)");
248
249 const Real halfLength = 0.5 * (b - a);
250 const Real center = 0.5 * (b + a);
251 const Real fCenter = f(center);
252
253 // Compute the integral using the 10- and 21-point formula.
254
255 res10 = 0;
256 res21 = w21b[5] * fCenter;
257 resAbs = w21b[5] * std::fabs(x: fCenter);
258
259 for (k = 0; k < 5; k++) {
260 Real abscissa = halfLength * x1[k];
261 Real fval1 = f(center + abscissa);
262 Real fval2 = f(center - abscissa);
263 Real fval = fval1 + fval2;
264 res10 += w10[k] * fval;
265 res21 += w21a[k] * fval;
266 resAbs += w21a[k] * (std::fabs(x: fval1) + std::fabs(x: fval2));
267 savfun[k] = fval;
268 fv1[k] = fval1;
269 fv2[k] = fval2;
270 }
271
272 for (k = 0; k < 5; k++) {
273 Real abscissa = halfLength * x2[k];
274 Real fval1 = f(center + abscissa);
275 Real fval2 = f(center - abscissa);
276 Real fval = fval1 + fval2;
277 res21 += w21b[k] * fval;
278 resAbs += w21b[k] * (std::fabs(x: fval1) + std::fabs(x: fval2));
279 savfun[k + 5] = fval;
280 fv3[k] = fval1;
281 fv4[k] = fval2;
282 }
283
284 result = res21 * halfLength;
285 resAbs *= halfLength ;
286 Real mean = 0.5 * res21;
287 resasc = w21b[5] * std::fabs(x: fCenter - mean);
288
289 for (k = 0; k < 5; k++)
290 resasc += (w21a[k] * (std::fabs(x: fv1[k] - mean)
291 + std::fabs(x: fv2[k] - mean))
292 + w21b[k] * (std::fabs(x: fv3[k] - mean)
293 + std::fabs(x: fv4[k] - mean)));
294
295 err = rescaleError (err: (res21 - res10) * halfLength, resultAbs: resAbs, resultAsc: resasc) ;
296 resasc *= halfLength ;
297
298 // test for convergence.
299 if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(x: result)){
300 setAbsoluteError(err);
301 setNumberOfEvaluations(21);
302 return result;
303 }
304
305 /* compute the integral using the 43-point formula. */
306
307 res43 = w43b[11] * fCenter;
308
309 for (k = 0; k < 10; k++)
310 res43 += savfun[k] * w43a[k];
311
312 for (k = 0; k < 11; k++){
313 Real abscissa = halfLength * x3[k];
314 Real fval = (f(center + abscissa)
315 + f(center - abscissa));
316 res43 += fval * w43b[k];
317 savfun[k + 10] = fval;
318 }
319
320 // test for convergence.
321
322 result = res43 * halfLength;
323 err = rescaleError (err: (res43 - res21) * halfLength, resultAbs: resAbs, resultAsc: resasc);
324
325 if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(x: result)){
326 setAbsoluteError(err);
327 setNumberOfEvaluations(43);
328 return result;
329 }
330
331 /* compute the integral using the 87-point formula. */
332
333 res87 = w87b[22] * fCenter;
334
335 for (k = 0; k < 21; k++)
336 res87 += savfun[k] * w87a[k];
337
338 for (k = 0; k < 22; k++){
339 Real abscissa = halfLength * x4[k];
340 res87 += w87b[k] * (f(center + abscissa)
341 + f(center - abscissa));
342 }
343
344 // test for convergence.
345 result = res87 * halfLength ;
346 err = rescaleError (err: (res87 - res43) * halfLength, resultAbs: resAbs, resultAsc: resasc);
347
348 setAbsoluteError(err);
349 setNumberOfEvaluations(87);
350 return result;
351 }
352
353 Real
354 GaussKronrodAdaptive::integrate(const ext::function<Real (Real)>& f,
355 Real a,
356 Real b) const {
357 return integrateRecursively(f, a, b, tolerance: absoluteAccuracy());
358 }
359
360 // weights for 7-point Gauss-Legendre integration
361 // (only 4 values out of 7 are given as they are symmetric)
362 static const Real g7w[] = { 0.417959183673469,
363 0.381830050505119,
364 0.279705391489277,
365 0.129484966168870 };
366 // weights for 15-point Gauss-Kronrod integration
367 static const Real k15w[] = { 0.209482141084728,
368 0.204432940075298,
369 0.190350578064785,
370 0.169004726639267,
371 0.140653259715525,
372 0.104790010322250,
373 0.063092092629979,
374 0.022935322010529 };
375 // abscissae (evaluation points)
376 // for 15-point Gauss-Kronrod integration
377 static const Real k15t[] = { 0.000000000000000,
378 0.207784955007898,
379 0.405845151377397,
380 0.586087235467691,
381 0.741531185599394,
382 0.864864423359769,
383 0.949107912342758,
384 0.991455371120813 };
385
386 Real GaussKronrodAdaptive::integrateRecursively(
387 const ext::function<Real (Real)>& f,
388 Real a,
389 Real b,
390 Real tolerance) const {
391
392 Real halflength = (b - a) / 2;
393 Real center = (a + b) / 2;
394
395 Real g7; // will be result of G7 integral
396 Real k15; // will be result of K15 integral
397
398 Real t, fsum; // t (abscissa) and f(t)
399 Real fc = f(center);
400 g7 = fc * g7w[0];
401 k15 = fc * k15w[0];
402
403 // calculate g7 and half of k15
404 Integer j, j2;
405 for (j = 1, j2 = 2; j < 4; j++, j2 += 2) {
406 t = halflength * k15t[j2];
407 fsum = f(center - t) + f(center + t);
408 g7 += fsum * g7w[j];
409 k15 += fsum * k15w[j2];
410 }
411
412 // calculate other half of k15
413 for (j2 = 1; j2 < 8; j2 += 2) {
414 t = halflength * k15t[j2];
415 fsum = f(center - t) + f(center + t);
416 k15 += fsum * k15w[j2];
417 }
418
419 // multiply by (a - b) / 2
420 g7 = halflength * g7;
421 k15 = halflength * k15;
422
423 // 15 more function evaluations have been used
424 increaseNumberOfEvaluations(increase: 15);
425
426 // error is <= k15 - g7
427 // if error is larger than tolerance then split the interval
428 // in two and integrate recursively
429 if (std::fabs(x: k15 - g7) < tolerance) {
430 return k15;
431 } else {
432 QL_REQUIRE(numberOfEvaluations()+30 <=
433 maxEvaluations(),
434 "maximum number of function evaluations "
435 "exceeded");
436 return integrateRecursively(f, a, b: center, tolerance: tolerance/2)
437 + integrateRecursively(f, a: center, b, tolerance: tolerance/2);
438 }
439 }
440
441
442 GaussKronrodAdaptive::GaussKronrodAdaptive(Real absoluteAccuracy,
443 Size maxEvaluations)
444 : Integrator(absoluteAccuracy, maxEvaluations) {
445 QL_REQUIRE(maxEvaluations >= 15,
446 "required maxEvaluations (" << maxEvaluations <<
447 ") not allowed. It must be >= 15");
448 }
449}
450

source code of quantlib/ql/math/integrals/kronrodintegral.cpp

Morty Proxy This is a proxified and sanitized view of the page, visit original site.