1 /** D port of intde1.c from Takuya OOURA's DE-Quadrature package 2 ($(LINK http://www.kurims.kyoto-u.ac.jp/~ooura/)). 3 4 Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.intde.intde1; 9 10 11 import std.math; 12 13 14 15 16 /** I = integral of f(x) over (a,b) 17 18 Params: 19 f = integrand f(x) 20 a = lower limit of integration 21 b = upper limit of integration 22 eps = relative error requested 23 i = approximation to the integral 24 err = estimate of the absolute error 25 26 Remarks: 27 <pre> 28 function 29 f(x) needs to be analytic over (a,b). 30 relative error 31 eps is relative error requested excluding 32 cancellation of significant digits. 33 i.e. eps means : (absolute error) / 34 (integral_a^b |f(x)| dx). 35 eps does not mean : (absolute error) / I. 36 error message 37 err >= 0 : normal termination. 38 err < 0 : abnormal termination (m >= mmax). 39 i.e. convergent error is detected : 40 1. f(x) or (d/dx)^n f(x) has 41 discontinuous points or sharp 42 peaks over (a,b). 43 you must divide the interval 44 (a,b) at this points. 45 2. relative error of f(x) is 46 greater than eps. 47 3. f(x) has oscillatory factor 48 and frequency of the oscillation 49 is very high. 50 </pre> 51 */ 52 void intde(Real, Func)(Func f, Real a, Real b, Real eps, Real *i, Real *err) 53 { 54 /* ---- adjustable parameter ---- */ 55 enum MMAX = 256; 56 enum : Real 57 { 58 EFS = 0.1, 59 HOFF = 8.5 60 } 61 /* ------------------------------ */ 62 int m; 63 Real epsln, epsh, h0, ehp, ehm, epst, ba, ir, h, iback, 64 irback, t, ep, em, xw, xa, wg, fa, fb, errt, errh, errd; 65 66 epsln = 1 - log(EFS * eps); 67 epsh = sqrt(EFS * eps); 68 h0 = HOFF / epsln; 69 ehp = exp(h0); 70 ehm = 1 / ehp; 71 epst = exp(-ehm * epsln); 72 ba = b - a; 73 ir = f((a + b) * 0.5) * (ba * 0.25); 74 *i = ir * PI; 75 *err = fabs(*i) * epst; 76 h = 2 * h0; 77 m = 1; 78 do { 79 iback = *i; 80 irback = ir; 81 t = h * 0.5; 82 do { 83 em = exp(t); 84 ep = PI_2 * em; 85 em = PI_2 / em; 86 do { 87 xw = 1 / (1 + exp(ep - em)); 88 xa = ba * xw; 89 wg = xa * (1 - xw); 90 fa = f(a + xa) * wg; 91 fb = f(b - xa) * wg; 92 ir += fa + fb; 93 *i += (fa + fb) * (ep + em); 94 errt = (fabs(fa) + fabs(fb)) * (ep + em); 95 if (m == 1) *err += errt * epst; 96 ep *= ehp; 97 em *= ehm; 98 } while (errt > *err || xw > epsh); 99 t += h; 100 } while (t < h0); 101 if (m == 1) { 102 errh = (*err / epst) * epsh * h0; 103 errd = 1 + 2 * errh; 104 } else { 105 errd = h * (fabs(*i - 2 * iback) + 4 * fabs(ir - 2 * irback)); 106 } 107 h *= 0.5; 108 m *= 2; 109 } while (errd > errh && m < MMAX); 110 *i *= h; 111 if (errd > errh) { 112 *err = -errd * m; 113 } else { 114 *err = errh * epsh * m / (2 * EFS); 115 } 116 } 117 118 119 unittest 120 { 121 real f1(real x) { return 1.0/sqrt(x); } 122 real f2(real x) { return sqrt(4 - x*x); } 123 124 real i, err; 125 real eps = 1e-15; 126 127 intde(&f1, 0.0L, 1.0L, eps, &i, &err); 128 assert (approxEqual(i, 2.0, eps)); 129 130 intde(&f2, 0.0L, 2.0L, eps, &i, &err); 131 assert (approxEqual(i, PI, eps)); 132 } 133 134 135 136 137 /** I = integral of f(x) over (a,infinity), 138 f(x) has not oscillatory factor. 139 140 Params: 141 f = integrand f(x) 142 a = lower limit of integration 143 eps = relative error requested 144 i = approximation to the integral 145 err = estimate of the absolute error 146 147 Remarks: 148 <pre> 149 function 150 f(x) needs to be analytic over (a,infinity). 151 relative error 152 eps is relative error requested excluding 153 cancellation of significant digits. 154 i.e. eps means : (absolute error) / 155 (integral_a^infinity |f(x)| dx). 156 eps does not mean : (absolute error) / I. 157 error message 158 err >= 0 : normal termination. 159 err < 0 : abnormal termination (m >= mmax). 160 i.e. convergent error is detected : 161 1. f(x) or (d/dx)^n f(x) has 162 discontinuous points or sharp 163 peaks over (a,infinity). 164 you must divide the interval 165 (a,infinity) at this points. 166 2. relative error of f(x) is 167 greater than eps. 168 3. f(x) has oscillatory factor 169 and decay of f(x) is very slow 170 as x -> infinity. 171 </pre> 172 */ 173 void intdei(Real, Func)(Func f, Real a, Real eps, Real *i, Real *err) 174 { 175 /* ---- adjustable parameter ---- */ 176 enum MMAX = 256; 177 enum : Real 178 { 179 EFS = 0.1, 180 HOFF = 11.0 181 } 182 /* ------------------------------ */ 183 int m; 184 Real epsln, epsh, h0, ehp, ehm, epst, ir, h, iback, irback, 185 t, ep, em, xp, xm, fp, fm, errt, errh, errd; 186 187 epsln = 1 - log(EFS * eps); 188 epsh = sqrt(EFS * eps); 189 h0 = HOFF / epsln; 190 ehp = exp(h0); 191 ehm = 1 / ehp; 192 epst = exp(-ehm * epsln); 193 ir = f(a + 1); 194 *i = ir * PI_2; 195 *err = fabs(*i) * epst; 196 h = 2 * h0; 197 m = 1; 198 do { 199 iback = *i; 200 irback = ir; 201 t = h * 0.5; 202 do { 203 em = exp(t); 204 ep = PI_4 * em; 205 em = PI_4 / em; 206 do { 207 xp = exp(ep - em); 208 xm = 1 / xp; 209 fp = f(a + xp) * xp; 210 fm = f(a + xm) * xm; 211 ir += fp + fm; 212 *i += (fp + fm) * (ep + em); 213 errt = (fabs(fp) + fabs(fm)) * (ep + em); 214 if (m == 1) *err += errt * epst; 215 ep *= ehp; 216 em *= ehm; 217 } while (errt > *err || xm > epsh); 218 t += h; 219 } while (t < h0); 220 if (m == 1) { 221 errh = (*err / epst) * epsh * h0; 222 errd = 1 + 2 * errh; 223 } else { 224 errd = h * (fabs(*i - 2 * iback) + 4 * fabs(ir - 2 * irback)); 225 } 226 h *= 0.5; 227 m *= 2; 228 } while (errd > errh && m < MMAX); 229 *i *= h; 230 if (errd > errh) { 231 *err = -errd * m; 232 } else { 233 *err = errh * epsh * m / (2 * EFS); 234 } 235 } 236 237 238 unittest 239 { 240 real f3(real x) { return 1.0/(1 + x*x); } 241 real f4(real x) { return exp(-x)/sqrt(x); } 242 243 real i, err; 244 real eps = 1e-15; 245 246 intdei(&f3, 0.0L, eps, &i, &err); 247 assert (approxEqual(i, PI_2, eps)); 248 249 intdei(&f4, 0.0L, eps, &i, &err); 250 assert (approxEqual(i, sqrt(PI), eps)); 251 } 252 253 254 255 256 /** I = integral of f(x) over (a,infinity), 257 f(x) has oscillatory factor : 258 --- 259 f(x) = g(x) * sin(omega * x + theta) as x -> infinity 260 --- 261 262 Params: 263 f = integrand f(x) 264 a = lower limit of integration 265 omega = frequency of oscillation 266 eps = relative error requested 267 i = approximation to the integral 268 err = estimate of the absolute error 269 270 Remarks: 271 <pre> 272 function 273 f(x) needs to be analytic over (a,infinity). 274 relative error 275 eps is relative error requested excluding 276 cancellation of significant digits. 277 i.e. eps means : (absolute error) / 278 (integral_a^R |f(x)| dx). 279 eps does not mean : (absolute error) / I. 280 error message 281 err >= 0 : normal termination. 282 err < 0 : abnormal termination (m >= mmax). 283 i.e. convergent error is detected : 284 1. f(x) or (d/dx)^n f(x) has 285 discontinuous points or sharp 286 peaks over (a,infinity). 287 you must divide the interval 288 (a,infinity) at this points. 289 2. relative error of f(x) is 290 greater than eps. 291 </pre> 292 */ 293 void intdeo(Real, Func)(Func f, Real a, Real omega, Real eps, 294 Real *i, Real *err) 295 { 296 /* ---- adjustable parameter ---- */ 297 enum 298 { 299 MMAX = 256, 300 LMAX = 5 301 } 302 303 enum : Real 304 { 305 EFS = 0.1, 306 ENOFF = 0.40, 307 PQOFF = 2.9, 308 PPOFF = -0.72 309 } 310 /* ------------------------------ */ 311 int n, m, l, k; 312 Real epsln, epsh, frq4, per2, pp, pq, ehp, ehm, ir, h, iback, 313 irback, t, ep, em, tk, xw, wg, xa, fp, fm, errh, tn, errd; 314 315 epsln = 1 - log(EFS * eps); 316 epsh = sqrt(EFS * eps); 317 n = cast(int) (ENOFF * epsln); 318 frq4 = fabs(omega) * M_2_PI; 319 per2 = PI / fabs(omega); 320 pq = PQOFF / epsln; 321 pp = PPOFF - log(pq * pq * frq4); 322 ehp = exp(2 * pq); 323 ehm = 1 / ehp; 324 xw = exp(pp - PI_2); 325 *i = f(a + sqrt(xw * (per2 * 0.5))); 326 ir = *i * xw; 327 *i *= per2 * 0.5; 328 *err = fabs(*i); 329 h = 2; 330 m = 1; 331 do { 332 iback = *i; 333 irback = ir; 334 t = h * 0.5; 335 do { 336 em = exp(2 * pq * t); 337 ep = PI_4 * em; 338 em = PI_4 / em; 339 tk = t; 340 do { 341 xw = exp(pp - ep - em); 342 wg = sqrt(frq4 * xw + tk * tk); 343 xa = xw / (tk + wg); 344 wg = (pq * xw * (ep - em) + xa) / wg; 345 fm = f(a + xa); 346 fp = f(a + xa + per2 * tk); 347 ir += (fp + fm) * xw; 348 fm *= wg; 349 fp *= per2 - wg; 350 *i += fp + fm; 351 if (m == 1) *err += fabs(fp) + fabs(fm); 352 ep *= ehp; 353 em *= ehm; 354 tk += 1; 355 } while (ep < epsln); 356 if (m == 1) { 357 errh = *err * epsh; 358 *err *= eps; 359 } 360 tn = tk; 361 while (fabs(fm) > *err) { 362 xw = exp(pp - ep - em); 363 xa = xw / tk * 0.5; 364 wg = xa * (1 / tk + 2 * pq * (ep - em)); 365 fm = f(a + xa); 366 ir += fm * xw; 367 fm *= wg; 368 *i += fm; 369 ep *= ehp; 370 em *= ehm; 371 tk += 1; 372 } 373 fm = f(a + per2 * tn); 374 em = per2 * fm; 375 *i += em; 376 if (fabs(fp) > *err || fabs(em) > *err) { 377 l = 0; 378 for (;;) { 379 l++; 380 tn += n; 381 em = fm; 382 fm = f(a + per2 * tn); 383 xa = fm; 384 ep = fm; 385 em += fm; 386 xw = 1; 387 wg = 1; 388 for (k = 1; k <= n - 1; k++) { 389 xw = xw * (n + 1 - k) / k; 390 wg += xw; 391 fp = f(a + per2 * (tn - k)); 392 xa += fp; 393 ep += fp * wg; 394 em += fp * xw; 395 } 396 wg = per2 * n / (wg * n + xw); 397 em = wg * fabs(em); 398 if (em <= *err || l >= LMAX) break; 399 *i += per2 * xa; 400 } 401 *i += wg * ep; 402 if (em > *err) *err = em; 403 } 404 t += h; 405 } while (t < 1); 406 if (m == 1) { 407 errd = 1 + 2 * errh; 408 } else { 409 errd = h * (fabs(*i - 2 * iback) + pq * fabs(ir - 2 * irback)); 410 } 411 h *= 0.5; 412 m *= 2; 413 } while (errd > errh && m < MMAX); 414 *i *= h; 415 if (errd > errh) { 416 *err = -errd; 417 } else { 418 *err *= m * 0.5; 419 } 420 } 421 422 423 unittest 424 { 425 real f5(real x) { return sin(x)/x; } 426 real f6(real x) { return cos(x)/sqrt(x); } 427 428 real i, err; 429 real eps = 1e-15; 430 431 intdeo(&f5, 0.0L, 1.0L, eps, &i, &err); 432 assert (approxEqual(i, PI_2, eps)); 433 434 intdeo(&f6, 0.0L, 1.0L, eps, &i, &err); 435 assert (approxEqual(i, sqrt(PI_2), eps)); 436 }