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 23 namespace QuantLib { 24 rescaleError(Real err,const Real resultAbs,const Real resultAsc)25 static Real rescaleError(Real err, 26 const Real resultAbs, 27 const Real resultAsc) { 28 err = std::fabs(err) ; 29 if (resultAsc != 0 && err != 0){ 30 Real scale = std::pow((200 * err / resultAsc), 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 relativeAccuracy() const218 Real GaussKronrodNonAdaptive::relativeAccuracy() const { 219 return relativeAccuracy_; 220 } 221 GaussKronrodNonAdaptive(Real absoluteAccuracy,Size maxEvaluations,Real relativeAccuracy)222 GaussKronrodNonAdaptive::GaussKronrodNonAdaptive(Real absoluteAccuracy, 223 Size maxEvaluations, 224 Real relativeAccuracy) 225 : Integrator(absoluteAccuracy, maxEvaluations), 226 relativeAccuracy_(relativeAccuracy) {} 227 228 Real integrate(const ext::function<Real (Real)> & f,Real a,Real b) const229 GaussKronrodNonAdaptive::integrate(const ext::function<Real (Real)>& f, 230 Real a, 231 Real b) const { 232 Real result; 233 //Size neval; 234 Real fv1[5], fv2[5], fv3[5], fv4[5]; 235 Real savfun[21]; /* array of function values which have been computed */ 236 Real res10, res21, res43, res87; /* 10, 21, 43 and 87 point results */ 237 Real err; 238 Real resAbs; /* approximation to the integral of abs(f) */ 239 Real resasc; /* approximation to the integral of abs(f-i/(b-a)) */ 240 int k ; 241 242 QL_REQUIRE(a<b, "b must be greater than a)"); 243 244 const Real halfLength = 0.5 * (b - a); 245 const Real center = 0.5 * (b + a); 246 const Real fCenter = f(center); 247 248 // Compute the integral using the 10- and 21-point formula. 249 250 res10 = 0; 251 res21 = w21b[5] * fCenter; 252 resAbs = w21b[5] * std::fabs(fCenter); 253 254 for (k = 0; k < 5; k++) { 255 Real abscissa = halfLength * x1[k]; 256 Real fval1 = f(center + abscissa); 257 Real fval2 = f(center - abscissa); 258 Real fval = fval1 + fval2; 259 res10 += w10[k] * fval; 260 res21 += w21a[k] * fval; 261 resAbs += w21a[k] * (std::fabs(fval1) + std::fabs(fval2)); 262 savfun[k] = fval; 263 fv1[k] = fval1; 264 fv2[k] = fval2; 265 } 266 267 for (k = 0; k < 5; k++) { 268 Real abscissa = halfLength * x2[k]; 269 Real fval1 = f(center + abscissa); 270 Real fval2 = f(center - abscissa); 271 Real fval = fval1 + fval2; 272 res21 += w21b[k] * fval; 273 resAbs += w21b[k] * (std::fabs(fval1) + std::fabs(fval2)); 274 savfun[k + 5] = fval; 275 fv3[k] = fval1; 276 fv4[k] = fval2; 277 } 278 279 result = res21 * halfLength; 280 resAbs *= halfLength ; 281 Real mean = 0.5 * res21; 282 resasc = w21b[5] * std::fabs(fCenter - mean); 283 284 for (k = 0; k < 5; k++) 285 resasc += (w21a[k] * (std::fabs(fv1[k] - mean) 286 + std::fabs(fv2[k] - mean)) 287 + w21b[k] * (std::fabs(fv3[k] - mean) 288 + std::fabs(fv4[k] - mean))); 289 290 err = rescaleError ((res21 - res10) * halfLength, resAbs, resasc) ; 291 resasc *= halfLength ; 292 293 // test for convergence. 294 if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(result)){ 295 setAbsoluteError(err); 296 setNumberOfEvaluations(21); 297 return result; 298 } 299 300 /* compute the integral using the 43-point formula. */ 301 302 res43 = w43b[11] * fCenter; 303 304 for (k = 0; k < 10; k++) 305 res43 += savfun[k] * w43a[k]; 306 307 for (k = 0; k < 11; k++){ 308 Real abscissa = halfLength * x3[k]; 309 Real fval = (f(center + abscissa) 310 + f(center - abscissa)); 311 res43 += fval * w43b[k]; 312 savfun[k + 10] = fval; 313 } 314 315 // test for convergence. 316 317 result = res43 * halfLength; 318 err = rescaleError ((res43 - res21) * halfLength, resAbs, resasc); 319 320 if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(result)){ 321 setAbsoluteError(err); 322 setNumberOfEvaluations(43); 323 return result; 324 } 325 326 /* compute the integral using the 87-point formula. */ 327 328 res87 = w87b[22] * fCenter; 329 330 for (k = 0; k < 21; k++) 331 res87 += savfun[k] * w87a[k]; 332 333 for (k = 0; k < 22; k++){ 334 Real abscissa = halfLength * x4[k]; 335 res87 += w87b[k] * (f(center + abscissa) 336 + f(center - abscissa)); 337 } 338 339 // test for convergence. 340 result = res87 * halfLength ; 341 err = rescaleError ((res87 - res43) * halfLength, resAbs, resasc); 342 343 setAbsoluteError(err); 344 setNumberOfEvaluations(87); 345 return result; 346 } 347 348 Real integrate(const ext::function<Real (Real)> & f,Real a,Real b) const349 GaussKronrodAdaptive::integrate(const ext::function<Real (Real)>& f, 350 Real a, 351 Real b) const { 352 return integrateRecursively(f, a, b, absoluteAccuracy()); 353 } 354 355 // weights for 7-point Gauss-Legendre integration 356 // (only 4 values out of 7 are given as they are symmetric) 357 static const Real g7w[] = { 0.417959183673469, 358 0.381830050505119, 359 0.279705391489277, 360 0.129484966168870 }; 361 // weights for 15-point Gauss-Kronrod integration 362 static const Real k15w[] = { 0.209482141084728, 363 0.204432940075298, 364 0.190350578064785, 365 0.169004726639267, 366 0.140653259715525, 367 0.104790010322250, 368 0.063092092629979, 369 0.022935322010529 }; 370 // abscissae (evaluation points) 371 // for 15-point Gauss-Kronrod integration 372 static const Real k15t[] = { 0.000000000000000, 373 0.207784955007898, 374 0.405845151377397, 375 0.586087235467691, 376 0.741531185599394, 377 0.864864423359769, 378 0.949107912342758, 379 0.991455371120813 }; 380 integrateRecursively(const ext::function<Real (Real)> & f,Real a,Real b,Real tolerance) const381 Real GaussKronrodAdaptive::integrateRecursively( 382 const ext::function<Real (Real)>& f, 383 Real a, 384 Real b, 385 Real tolerance) const { 386 387 Real halflength = (b - a) / 2; 388 Real center = (a + b) / 2; 389 390 Real g7; // will be result of G7 integral 391 Real k15; // will be result of K15 integral 392 393 Real t, fsum; // t (abscissa) and f(t) 394 Real fc = f(center); 395 g7 = fc * g7w[0]; 396 k15 = fc * k15w[0]; 397 398 // calculate g7 and half of k15 399 Integer j, j2; 400 for (j = 1, j2 = 2; j < 4; j++, j2 += 2) { 401 t = halflength * k15t[j2]; 402 fsum = f(center - t) + f(center + t); 403 g7 += fsum * g7w[j]; 404 k15 += fsum * k15w[j2]; 405 } 406 407 // calculate other half of k15 408 for (j2 = 1; j2 < 8; j2 += 2) { 409 t = halflength * k15t[j2]; 410 fsum = f(center - t) + f(center + t); 411 k15 += fsum * k15w[j2]; 412 } 413 414 // multiply by (a - b) / 2 415 g7 = halflength * g7; 416 k15 = halflength * k15; 417 418 // 15 more function evaluations have been used 419 increaseNumberOfEvaluations(15); 420 421 // error is <= k15 - g7 422 // if error is larger than tolerance then split the interval 423 // in two and integrate recursively 424 if (std::fabs(k15 - g7) < tolerance) { 425 return k15; 426 } else { 427 QL_REQUIRE(numberOfEvaluations()+30 <= 428 maxEvaluations(), 429 "maximum number of function evaluations " 430 "exceeded"); 431 return integrateRecursively(f, a, center, tolerance/2) 432 + integrateRecursively(f, center, b, tolerance/2); 433 } 434 } 435 436 GaussKronrodAdaptive(Real absoluteAccuracy,Size maxEvaluations)437 GaussKronrodAdaptive::GaussKronrodAdaptive(Real absoluteAccuracy, 438 Size maxEvaluations) 439 : Integrator(absoluteAccuracy, maxEvaluations) { 440 QL_REQUIRE(maxEvaluations >= 15, 441 "required maxEvaluations (" << maxEvaluations << 442 ") not allowed. It must be >= 15"); 443 } 444 } 445