1 // Written in the D programming language. 2 3 /** 4 This module provides numerical integration by double-exponential (DE) formula. 5 */ 6 module deint.deint; 7 8 import std.algorithm; 9 import std.array; 10 import std.functional; 11 import std.math; 12 import std.typecons; 13 import std.traits; 14 15 16 /** 17 This template checks that a type `T` is a kind of NumInt. 18 */ 19 enum bool isNumericalIntegrate(T) 20 = is(typeof((T integ){ 21 foreach(xw; lockstep(integ.xs, integ.ws)) {} 22 })); 23 24 25 26 /** This structure stores computing points and weights for numerical integration. 27 * In this type, a numerical integration is computed as 28 * int f(x) dx = sum_i f(xs[i]) ws[i] 29 */ 30 struct NumInt(X, W) 31 { 32 /** Constructor 33 * Params: 34 * xs = computing points 35 * ws = weights for each computing points 36 */ 37 this(immutable(X)[] xs, immutable(W)[] ws) 38 in{ 39 assert(xs.length == ws.length); 40 } 41 do { 42 _xs = xs; 43 _ws = ws; 44 } 45 46 47 const @property 48 { 49 immutable(X)[] xs() { return _xs; } 50 immutable(W)[] ws() { return _ws; } 51 } 52 53 54 /** Execute integration 55 */ 56 auto integrate(Fn)(scope Fn func) const 57 { 58 alias T = Unqual!(typeof(func(_xs[0]) * _ws[0])); 59 60 T sum = 0; 61 foreach(i; 0 .. _xs.length) 62 sum += func(_xs[i]) * _ws[i]; 63 64 return sum; 65 } 66 67 68 private: 69 immutable(X)[] _xs; 70 immutable(W)[] _ws; 71 } 72 73 74 /** This function create a instance of NumInt!(F, F) for the DE formula. 75 * It is also known as "Tanh-sinh quadrature". 76 * In the DE formula, the integration (1) is converted to (2). 77 * (1) int_{xa}^{xb} f(x) dx 78 * (2) int_{ta}^{tb} f(g(t)) g'(t) dt 79 * 80 * The type of DE formula is automatically decided from the given interval of the integration. 81 * 82 * Params: 83 * xa = starting value of original integration. 84 * xb = end value of original integration. 85 * isExpDecay = if the integration is formed as int_a^b f(x) exp(-x) dx, this value is Yes. otherwise No. 86 * trapN = division points of trapezoidal quadrature. 87 * ta = starting value of integration transformed by DE-formula. 88 * tb = starting value of integration transformed by DE-formula. 89 * 90 * Reference: 91 * $(HTTP en.wikipedia.org/wiki/Tanh-sinh_quadrature) 92 */ 93 NumInt!(F, F) makeDEInt(F)( 94 F xa, F xb, 95 Flag!"isExpDecay" isExpDecay = No.isExpDecay, 96 size_t trapN = 100, 97 F ta = -5, F tb = 5) 98 { 99 NumInt!(F, F) _makeParamsImpl(scope F[2] delegate(F) fn) 100 { 101 immutable(F)[] xs; 102 immutable(F)[] ws; 103 if (!__ctfe) 104 { 105 xs.reserve(trapN); 106 ws.reserve(trapN); 107 } 108 immutable F h = (tb - ta) / (trapN-1); 109 foreach(i; 0 .. trapN) { 110 immutable xWt = fn(i * h + ta); 111 xs ~= xWt[0]; 112 ws ~= xWt[1] * h; 113 } 114 115 return NumInt!(F, F)(xs, ws); 116 } 117 118 119 if(xa > xb) { 120 auto invInt = .makeDEInt!F(xb, xa, isExpDecay, trapN, ta, tb); 121 auto newws = invInt.ws.map!"cast(immutable)(-a)".array; 122 return NumInt!(F, F)(invInt.xs, newws); 123 } else if(xa == -F.infinity && xb != F.infinity) { 124 auto tmpInt = .makeDEInt!F(-xb, F.infinity, isExpDecay, trapN, ta, tb); 125 auto newxs = tmpInt.xs.map!"cast(immutable)(-a)".array; 126 return NumInt!(F, F)(newxs, tmpInt.ws); 127 } 128 129 130 if(xa == -F.infinity && xb == F.infinity){ 131 assert(!isExpDecay); 132 133 return _makeParamsImpl(delegate(F t){ 134 immutable F 135 sinht = sinh(t), 136 x = sinh(PI / 2 * sinht), 137 dx = cosh(PI / 2 * sinht) * PI / 2 * cosh(t); 138 139 return cast(F[2])[x, dx]; 140 }); 141 }else if(xb == F.infinity) { 142 return _makeParamsImpl(delegate(F t){ 143 if(!isExpDecay){ 144 real x = exp(PI / 2 * sinh(t)), 145 dx = x * PI / 2 * cosh(t); 146 147 return cast(F[2])[x + xa, dx]; 148 }else{ 149 real expmt = exp(-t), 150 x = exp(t - expmt), 151 dx = (1 + expmt) * x; 152 153 154 return cast(F[2])[x + xa, dx]; 155 } 156 }); 157 }else{ 158 if(xa == 0 || xb == 0){ 159 return _makeParamsImpl(delegate(F t){ 160 immutable F 161 cosht = cosh(t), 162 sinht = sinh(t), 163 epsinht = exp(PI * sinht), 164 x = (xa + xb*epsinht)/(1 + epsinht), 165 cosh2 = cosh(PI / 2 * sinht)^^2, 166 dx = PI / 2 * cosht / cosh2; 167 168 return cast(F[2])[x, dx * (xb - xa)/2]; 169 }); 170 }else{ 171 immutable F diff2 = (xb - xa)/2; 172 immutable F avg2 = (xb + xa)/2; 173 174 return _makeParamsImpl(delegate(F t){ 175 immutable F 176 cosht = cosh(t), 177 sinht = sinh(t), 178 x = tanh(PI / 2 * sinht) * diff2 + avg2, 179 cosh2 = cosh(PI / 2 * sinht)^^2, 180 dx = PI / 2 * cosht / cosh2; 181 182 return cast(F[2])[x, dx * diff2]; 183 }); 184 } 185 } 186 } 187 188 189 /// 190 unittest 191 { 192 // integration on [0, 1] 193 auto int01 = makeDEInt!real(0, 1); 194 195 // int_0^1 x dx = 0.5 196 assert(int01.integrate((real x) => x).approxEqual(0.5)); 197 198 // int_0^1 x^^2 dx = 1/3 199 assert(int01.integrate((real x) => x^^2).approxEqual(1/3.0)); 200 201 // int_0^1 log(x)/sqrt(x) dx = -4 202 assert(int01.integrate((real x) => log(x)/sqrt(x)).approxEqual(-4, 1E-12, 1E-12)); 203 204 // int_{-1}^0 log(-x)/sqrt(-x) dx = -4 205 assert(makeDEInt!real(-1, 0).integrate((real x) => log(-x)/sqrt(-x)).approxEqual(-4, 1E-12, 1E-12)); 206 207 // integration on [-1, 1] 208 auto int11 = makeDEInt!real(-1, 1); 209 210 // int_{-1}^1 1/sqrt(1-x^2) dx = pi 211 assert(int11.integrate((real x) => 1/(1 + x^^2)).approxEqual(PI/2, 1E-12, 1E-12)); 212 213 // integration on [-inf, inf] 214 auto intII = makeDEInt!real(-real.infinity, real.infinity); 215 216 // Gaussian integral 217 assert(intII.integrate((real x) => exp(-x^^2)).approxEqual(sqrt(PI))); 218 219 220 import std.mathspecial; 221 // integration on [1, inf] and integrand is expected to decay exponentially 222 auto intFI = makeDEInt!real(1, real.infinity, Yes.isExpDecay); 223 224 // compute incomplete gamma function: int_1^inf x * exp(-x) dx = Gamma(2, 1) 225 assert(intFI.integrate((real x) => x * exp(-x)).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2))); 226 } 227 228 // Test of \int_a^b exp(x) dx 229 unittest 230 { 231 // [a, b] 232 real[2][] dataset = [ 233 [0.0, 1.0], 234 [-1.0, 1.0], 235 [0.0, -2.0], 236 [2.0, 0.0], 237 [-10.0, 0.0], 238 [0.0, 10.0] 239 ]; 240 241 foreach(data; dataset) { 242 real intDE = makeDEInt!real(data[0], data[1]).integrate((real x) => exp(x)); 243 real truevalue = exp(data[1]) - exp(data[0]); 244 245 assert(approxEqual(intDE, truevalue)); 246 } 247 248 249 // reverse interval 250 foreach(data; dataset) { 251 //real intDE = intAToB_DE(data[0], data[1], 100, (real x) => exp(x)); 252 real intDE = makeDEInt!real(data[1], data[0]).integrate((real x) => exp(x)); 253 real truevalue = -(exp(data[1]) - exp(data[0])); 254 255 assert(approxEqual(intDE, truevalue)); 256 } 257 } 258 259 // Complementary error function: erfc(a) = 2/sqrt(pi) * int_a^inf e^{-x^2} dx 260 unittest 261 { 262 import std.mathspecial; 263 264 real[] dataset = [ 265 0.0, 266 0.01, 267 0.1, 268 1, 269 10, 270 100, 271 ]; 272 273 foreach(a; dataset) { 274 auto intDE = makeDEInt!real(a, real.infinity).integrate((real x) => exp(-x^^2)); 275 intDE *= M_2_SQRTPI; 276 277 auto truevalue = erfc(a); 278 assert(approxEqual(intDE, truevalue)); 279 } 280 } 281 282 // Complementary error function: erfc(a) = -2/sqrt(pi) * int_(-inf)^(-a) e^{-x^2} dx 283 unittest 284 { 285 import std.mathspecial; 286 287 real[] dataset = [ 288 0.0, 289 0.01, 290 0.1, 291 1, 292 10, 293 100, 294 ]; 295 296 foreach(a; dataset) { 297 auto intDE = makeDEInt!real(-real.infinity, -a).integrate((real x) => exp(-x^^2)); 298 intDE *= M_2_SQRTPI; 299 300 auto truevalue = erfc(a); 301 assert(approxEqual(intDE, truevalue)); 302 } 303 } 304 305 // int_inf^inf 1/(1+x^2) dx = pi 306 unittest 307 { 308 auto intDE1 = makeDEInt!real(-real.infinity, real.infinity) 309 .integrate((real x) => 1/(1 + x^^2)); 310 assert(approxEqual(intDE1, PI)); 311 } 312 313 // int exp(-x^2) = sqrt(pi) 314 unittest 315 { 316 immutable val1 = makeDEInt!real(0, real.infinity, Yes.isExpDecay, 100) 317 .integrate((real x) => 1.0/(2*sqrt(x)) * exp(-x) ); 318 319 immutable val2 = makeDEInt!real(real.infinity, 0, Yes.isExpDecay, 100) 320 .integrate((real x) => -1.0/(2*sqrt(x)) * exp(-x) ); 321 322 immutable val3 = makeDEInt!real(-real.infinity, 0, Yes.isExpDecay, 100) 323 .integrate((real x) => 1.0/(2*sqrt(-x)) * exp(x) ); 324 325 immutable val4 = makeDEInt!real(0, -real.infinity, Yes.isExpDecay, 100) 326 .integrate((real x) => -1.0/(2*sqrt(-x)) * exp(x) ); 327 328 assert(approxEqual(val1, sqrt(PI)/2)); 329 assert(approxEqual(val2, sqrt(PI)/2)); 330 assert(approxEqual(val3, sqrt(PI)/2)); 331 assert(approxEqual(val4, sqrt(PI)/2)); 332 } 333 334 // int exp(-x^2) * 1i = 1i * sqrt(pi) 335 unittest 336 { 337 import std.complex : Complex; 338 339 immutable val1 = makeDEInt!real(0, real.infinity, Yes.isExpDecay) 340 .integrate((real x) => 1.0/(2*sqrt(x)) * exp(-x) * Complex!real(0, 1)); 341 342 assert(val1.re.approxEqual(0)); 343 assert(val1.im.approxEqual(sqrt(PI)/2)); 344 } 345 346 // unittest for README.md 347 unittest 348 { 349 import std.mathspecial; 350 import deint.utils; 351 352 // integrate f(x)exp(-x) on (1, inf) 353 // Now, we know that the integrand f(x)exp(-x) decay exponentially. 354 auto intFI = makeDEInt!real(1, real.infinity, Yes.isExpDecay); 355 356 // incomplete gamma function 357 assert(intFI.integrate((real x) => x * exp(-x)).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2))); 358 359 // Also, you can use `withWeight` which pre-computes and stores weight exp(-x). 360 // The `withWeight` is useful when integrands have same weights. 361 auto intFIW = intFI.withWeight((real x) => exp(-x)); 362 363 // incomplete gamma function 364 assert(intFIW.integrate((real x) => x).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2))); 365 assert(intFIW.integrate((real x) => x^^2).approxEqual(gammaIncompleteCompl(3, 1) * gamma(3))); 366 assert(intFIW.integrate((real x) => x^^3).approxEqual(gammaIncompleteCompl(4, 1) * gamma(4))); 367 } 368 369 unittest 370 { 371 // integration on [-inf, inf] 372 auto intII = makeDEInt!real(-real.infinity, real.infinity); 373 374 // Gaussian integral 375 assert(intII.integrate((real x) => exp(-x^^2)).approxEqual(sqrt(PI))); 376 } 377 378 unittest 379 { 380 // DEInt!real is a struct which computes x_k and w_k in advance. 381 auto int01 = makeDEInt!real(0, 1); 382 383 // When f(x) = x, int_0^1 x dx = 0.5 384 assert(int01.integrate((real x) => x).approxEqual(0.5)); 385 386 // `DEInt!real` is reusable. 387 // When f(x) = x^^2, int_0^1 x^^2 dx = 1/3 388 assert(int01.integrate((real x) => x^^2).approxEqual(1/3.0)); 389 } 390 391 392 393 /** This function create a instance of NumInt!(F, F) for computing fourtier-type integration by the DE formula. 394 * On the DE formula for the fourier-type integration, the integration (1) is converted to (2). 395 * (1) int_0^{inf} g(x) sin(omega * x) dx 396 * (2) sum_{k=-Nlow}^{Nhigh} g(M*phi(h*k)) sin(M*phi(h*k)/omega) M*phi'(h*k)/omega * h 397 * 398 * Params: 399 * isSine = `Yes` if the type of integration is "sine". `No` if the type is "cosine". 400 * omega = angular frequency 401 * stepH = step value of each interval 402 * Nlow = number of negative-valued computing points 403 * Nhigh = number of positive-valued computing points 404 * 405 * Reference: 406 * $(HTTP www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1040-20.pdf) 407 */ 408 NumInt!(F, F) makeDEIntFourier(F)( 409 Flag!"isSine" isSine, 410 F omega, 411 F stepH = 0.1, 412 size_t Nlow = 49, size_t Nhigh = 50) 413 { 414 immutable F beta = 1.0L/4; 415 immutable F h = stepH; 416 immutable F M = PI / h; 417 immutable F alpha = beta / sqrt(1 + M*log(1 + M)/(4*PI)); 418 419 F[2] phi(F t) { 420 F num = t; 421 F den = 1 - exp(-2*t - alpha*(1 - exp(-t)) - beta*(exp(t) - 1)); 422 F dden = (den - 1) * (-2 - alpha*exp(-t) - beta * exp(t)); 423 424 F num2 = den - num * dden; 425 F den2 = den^^2; 426 427 if(abs(t) < 1E-3) { 428 F ddden = (den-1) * (-2 - alpha*exp(-t) - beta * exp(t))^^2 429 + (den-1) * ( + alpha*exp(-t) - beta * exp(t)); 430 431 num = 1; 432 den = dden; 433 434 num2 = dden - (dden + num * ddden); 435 den2 = 2*den*dden; 436 } 437 438 return [num/den, num2/den2]; 439 } 440 441 immutable(F)[] xs; 442 immutable(F)[] ws; 443 if (!__ctfe) 444 { 445 xs.reserve(Nlow + Nhigh); 446 ws.reserve(Nlow + Nhigh); 447 } 448 449 foreach(long i; -cast(long)Nlow .. cast(long)Nhigh) { 450 F t = h * i; 451 if(isSine) { 452 auto p = phi(t); 453 xs ~= M * p[0] / omega; 454 ws ~= M * p[1] / omega * h; 455 } else { 456 auto p = phi(t - PI/(2*M)); 457 xs ~= M * p[0] / omega; 458 ws ~= M * p[1] / omega * h; 459 } 460 } 461 462 return NumInt!(F, F)(xs, ws); 463 } 464 465 unittest { 466 auto de = makeDEIntFourier!real(Yes.isSine, 1, 0.2); 467 auto dirichlet = de.integrate((real x) => sin(x)/x); 468 469 assert(dirichlet.approxEqual(PI/2)); 470 }