1 // Written in the D programming language. 2 3 /** 4 This module provides numerical integration by double-exponential (DE) formula. 5 */ 6 module deint; 7 8 import std.algorithm; 9 import std.array; 10 import std.functional; 11 import std.math; 12 import std.typecons; 13 14 15 16 /** A struct for performing numerical integration by double-exponential (DE) formula. 17 * It is also known as "Tanh-sinh quadrature". 18 * In DE formula, the integration (1) is converted to (2). 19 * (1) int_{xa}^{xb} f(x) w(x) dx 20 * (2) int_{ta}^{tb} f(g(t)) w(g(t)) g'(t) dt 21 * 22 * This struct calculates reusable state in advance. 23 * The type of DE formula is automatically decided from the given interval of the integration. 24 * 25 * Reference: 26 * $(HTTP en.wikipedia.org/wiki/Tanh-sinh_quadrature) 27 */ 28 struct DEInt(F) 29 { 30 private static 31 F returnConst1(F x) { return F(1); } 32 33 34 /** Initialize an object for computing DE integration. 35 * Params: 36 * xa = starting value of original integration. 37 * xb = end value of original integration. 38 * weightFn = weight function. The default value of this parameter is `(F x) => F(1)`. 39 * isExpDecay = if the integration is formed as int_a^b f(x) exp(-x) dx, this value is Yes. otherwise No. 40 * trapN = division points of trapezoidal quadrature. 41 * ta = starting value of integration transformed by DE-formula. 42 * tb = starting value of integration transformed by DE-formula. 43 */ 44 this( 45 F xa, F xb, 46 scope F delegate(F) weightFn = toDelegate(&returnConst1), 47 Flag!"isExpDecay" isExpDecay = No.isExpDecay, 48 size_t trapN = 100, 49 F ta = -5, F tb = 5) 50 { 51 setParams(xa, xb, weightFn, isExpDecay, trapN, ta, tb); 52 } 53 54 55 /// ditto 56 void setParams( 57 F xa, F xb, 58 scope F delegate(F) weightFn = toDelegate(&returnConst1), 59 Flag!"isExpDecay" isExpDecay = No.isExpDecay, 60 size_t trapN = 100, 61 F ta = -5, F tb = 5) 62 in{ 63 if(isExpDecay) 64 assert(_xa != -F.infinity); 65 } 66 do { 67 if(xa > xb) { 68 this.setParams(xb, xa, weightFn, isExpDecay, trapN, ta, tb); 69 _weights = _weights.map!"cast(immutable)(-a)".array; 70 } else if(xa == -F.infinity && xb != F.infinity) { 71 //assert(!isExpDecay); 72 this.setParams(-xb, F.infinity, (F x) => weightFn(-x), isExpDecay, trapN, ta, tb); 73 _xs = _xs.map!"cast(immutable)(-a)".array; 74 } 75 76 _xa = xa; 77 _xb = xb; 78 _ta = ta; 79 _tb = tb; 80 _trapN = trapN; 81 _isExpDecay = cast(bool)isExpDecay; 82 83 if(_xs !is null) 84 return; 85 86 87 if(xa == -F.infinity && xb == F.infinity){ 88 assert(!isExpDecay); 89 90 auto params = _makeParamsImpl(ta, tb, trapN, weightFn, delegate(F t){ 91 immutable F 92 sinht = sinh(t), 93 x = sinh(PI / 2 * sinht), 94 dx = cosh(PI / 2 * sinht) * PI / 2 * cosh(t); 95 96 return cast(F[2])[x, dx]; 97 }); 98 99 100 _xs = params[0]; 101 _weights = params[1]; 102 _intType = isExpDecay ? "IIE" : "II"; 103 }else if(xb == F.infinity) { 104 auto params = _makeParamsImpl(ta, tb, trapN, weightFn, delegate(F t){ 105 if(!isExpDecay){ 106 real x = exp(PI / 2 * sinh(t)), 107 dx = x * PI / 2 * cosh(t); 108 109 return cast(F[2])[x + xa, dx]; 110 }else{ 111 real expmt = exp(-t), 112 x = exp(t - expmt), 113 dx = (1 + expmt) * x; 114 115 116 return cast(F[2])[x + xa, dx]; 117 } 118 }); 119 120 _xs = params[0]; 121 _weights = params[1]; 122 _intType = isExpDecay ? "FIE" : "FI"; 123 }else{ 124 immutable F diff2 = (xb - xa) / 2, 125 avg2 = (xb + xa) / 2; 126 127 auto params = _makeParamsImpl(ta, tb, trapN, weightFn, delegate(F t){ 128 immutable F 129 cosht = cosh(t), 130 sinht = sinh(t), 131 x = tanh(PI / 2 * sinht) * diff2 + avg2, 132 cosh2 = cosh(PI / 2 * sinht)^^2, 133 dx = PI / 2 * cosht / cosh2; 134 135 return cast(F[2])[x, dx * diff2]; 136 }); 137 138 _xs = params[0]; 139 _weights = params[1]; 140 _intType = isExpDecay ? "FFE" : "FF"; 141 } 142 } 143 144 145 146 /** Execute integration of func by DE formula. 147 */ 148 F integrate(Fn)(scope Fn func) const 149 { 150 F sum = 0; 151 foreach(i; 0 .. _xs.length) 152 sum += func(_xs[i]) * _weights[i]; 153 154 return sum; 155 } 156 157 158 @property const 159 { 160 /** Return type of integration. 161 * The First and second charactor of the return value are explain the the starting value or end value of integration is finite value ('F') or infinity 'I'. 162 * If the return value has the third character and its value is 'E', the integration is formed as int_{xa}^{xb} f(x) exp(-x) dx 163 */ 164 string type() { return _intType; } 165 166 /// Return the set value 167 F xa() { return _xa; } 168 169 /// ditto 170 F xb() { return _xb; } 171 172 /// ditto 173 F ta() { return _ta; } 174 175 /// ditto 176 F tb() { return _tb; } 177 178 /// ditto 179 size_t trapN() { return _trapN; } 180 181 /// ditto 182 bool isExpDecay() { return _isExpDecay; } 183 184 /// Return division points (computing points) of trapezoidal quadrature for DE formula. 185 immutable(F)[] xs() { return _xs; } 186 187 /// Return weights of each division points. 188 immutable(F)[] weights() { return _weights; } 189 } 190 191 192 private: 193 string _intType; 194 F _xa, _xb, _ta, _tb; 195 size_t _trapN; 196 bool _isExpDecay; 197 immutable(F)[] _xs; 198 immutable(F)[] _weights; 199 200 static 201 immutable(F)[][2] _makeParamsImpl(F ta, F tb, size_t trapN, scope F delegate(F) weightFn, scope F[2] delegate(F) fn) 202 { 203 immutable(F)[] xs, weights; 204 immutable F h = (tb - ta) / (trapN-1); 205 foreach(i; 0 .. trapN) { 206 immutable xWt = fn(i * h + ta); 207 xs ~= xWt[0]; 208 weights ~= xWt[1] * h * weightFn(xWt[0]); 209 } 210 211 return [xs, weights]; 212 } 213 } 214 215 /// 216 unittest 217 { 218 // integration on [0, 1] 219 auto int01 = DEInt!real(0, 1); 220 assert(int01.type == "FF"); 221 222 // int_0^1 x dx = 0.5 223 assert(int01.integrate((real x) => x).approxEqual(0.5)); 224 225 // int_0^1 x^^2 dx = 1/3 226 assert(int01.integrate((real x) => x^^2).approxEqual(1/3.0)); 227 228 229 // integration on [-inf, inf] 230 auto intII = DEInt!real(-real.infinity, real.infinity); 231 assert(intII.type == "II"); 232 233 // Gaussian integral 234 assert(intII.integrate((real x) => exp(-x^^2)).approxEqual(sqrt(PI))); 235 236 import std.mathspecial; 237 // integration int_1^inf x * exp(-x) dx = Gamma(2, 1) 238 auto intFI = DEInt!real(1, real.infinity, (real x) => exp(-x), Yes.isExpDecay); 239 assert(intFI.type == "FIE"); 240 241 // incomplete gamma function 242 assert(intFI.integrate((real x) => x).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2))); 243 } 244 245 // Test of \int_a^b exp(x) dx 246 unittest 247 { 248 // [a, b] 249 real[2][] dataset = [ 250 [0.0, 1.0], 251 [-1.0, 1.0], 252 [0.0, -2.0], 253 [2.0, 0.0], 254 [-10.0, 0.0], 255 [0.0, 10.0] 256 ]; 257 258 foreach(data; dataset) { 259 real intDE = DEInt!real(data[0], data[1]).integrate((real x) => exp(x)); 260 real truevalue = exp(data[1]) - exp(data[0]); 261 262 assert(approxEqual(intDE, truevalue)); 263 } 264 265 266 // reverse interval 267 foreach(data; dataset) { 268 //real intDE = intAToB_DE(data[0], data[1], 100, (real x) => exp(x)); 269 real intDE = DEInt!real(data[1], data[0]).integrate((real x) => exp(x)); 270 real truevalue = -(exp(data[1]) - exp(data[0])); 271 272 assert(approxEqual(intDE, truevalue)); 273 } 274 } 275 276 // Complementary error function: erfc(a) = 2/sqrt(pi) * int_a^inf e^{-x^2} dx 277 unittest 278 { 279 import std.mathspecial; 280 281 real[] dataset = [ 282 0.0, 283 0.01, 284 0.1, 285 1, 286 10, 287 100, 288 ]; 289 290 foreach(a; dataset) { 291 auto intDE = DEInt!real(a, real.infinity).integrate((real x) => exp(-x^^2)); 292 intDE *= M_2_SQRTPI; 293 294 auto truevalue = erfc(a); 295 assert(approxEqual(intDE, truevalue)); 296 } 297 } 298 299 // Complementary error function: erfc(a) = -2/sqrt(pi) * int_(-inf)^(-a) e^{-x^2} dx 300 unittest 301 { 302 import std.mathspecial; 303 304 real[] dataset = [ 305 0.0, 306 0.01, 307 0.1, 308 1, 309 10, 310 100, 311 ]; 312 313 foreach(a; dataset) { 314 auto intDE = DEInt!real(-real.infinity, -a).integrate((real x) => exp(-x^^2)); 315 intDE *= M_2_SQRTPI; 316 317 auto truevalue = erfc(a); 318 assert(approxEqual(intDE, truevalue)); 319 } 320 } 321 322 // int_inf^inf 1/(1+x^2) dx = pi 323 unittest 324 { 325 auto intDE1 = DEInt!real(-real.infinity, real.infinity) 326 .integrate((real x) => 1/(1 + x^^2)); 327 assert(approxEqual(intDE1, PI)); 328 } 329 330 // int exp(-x^2) = sqrt(pi) 331 unittest 332 { 333 immutable val1 = DEInt!real(0, real.infinity, (real x) => exp(-x), Yes.isExpDecay, 100) 334 .integrate((real x) => 1.0/(2*sqrt(x))); 335 336 immutable val2 = DEInt!real(real.infinity, 0, (real x) => exp(-x), Yes.isExpDecay, 100) 337 .integrate((real x) => -1.0/(2*sqrt(x))); 338 339 immutable val3 = DEInt!real(-real.infinity, 0, (real x) => exp(x), Yes.isExpDecay, 100) 340 .integrate((real x) => 1.0/(2*sqrt(-x))); 341 342 immutable val4 = DEInt!real(0, -real.infinity, (real x) => exp(x), Yes.isExpDecay, 100) 343 .integrate((real x) => -1.0/(2*sqrt(-x))); 344 345 assert(approxEqual(val1, sqrt(PI)/2)); 346 assert(approxEqual(val2, sqrt(PI)/2)); 347 assert(approxEqual(val3, sqrt(PI)/2)); 348 assert(approxEqual(val4, sqrt(PI)/2)); 349 }