1 /++ 2 This module contains betterC compatible quadrature computation routines. 3 +/ 4 module mir.quadrature; 5 6 import mir.math.common: sqrt, exp; 7 import mir.math.constant: PI, LN2; 8 9 @safe pure nothrow: 10 11 @nogc extern(C) 12 { 13 double lgamma(double); 14 double tgamma(double); 15 } 16 17 /++ 18 Gauss-Hermite Quadrature 19 20 Params: 21 x = (out) user-allocated quadrature nodes in ascending order length of `N` 22 w = (out) user-allocated corresponding quadrature weights length of `N` 23 work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` 24 25 Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. 26 27 +/ 28 size_t gaussHermiteQuadrature(T)( 29 scope T[] x, 30 scope T[] w, 31 scope T[] work) @nogc 32 in { 33 assert(x.length == w.length); 34 if (x.length) 35 assert(work.length >= (x.length + 1) ^^ 2); 36 } 37 do { 38 enum T mu0 = sqrt(PI); 39 foreach (i; 0 .. x.length) 40 { 41 x[i] = 0; 42 w[i] = T(0.5) * i; 43 } 44 return golubWelsch!T(x, w, work, mu0, true); 45 } 46 47 /// 48 unittest 49 { 50 import mir.math.common; 51 52 auto n = 5; 53 auto x = new double[n]; 54 auto w = new double[n]; 55 auto work = new double[(n + 1) ^^ 2]; 56 57 gaussHermiteQuadrature(x, w, work); 58 59 static immutable xc = 60 [-2.02018287, 61 -0.95857246, 62 0. , 63 0.95857246, 64 2.02018287]; 65 66 static immutable wc = 67 [0.01995324, 68 0.39361932, 69 0.94530872, 70 0.39361932, 71 0.01995324]; 72 73 foreach (i; 0 .. n) 74 { 75 assert(x[i].approxEqual(xc[i])); 76 assert(w[i].approxEqual(wc[i])); 77 } 78 } 79 80 /++ 81 Gauss-Jacobi Quadrature 82 83 Params: 84 x = (out) user-allocated quadrature nodes in ascending order length of `N` 85 w = (out) user-allocated corresponding quadrature weights length of `N` 86 work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` 87 alpha = parameter '> -1' 88 beta = parameter '> -1' 89 90 Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. 91 92 +/ 93 size_t gaussJacobiQuadrature(T)( 94 scope T[] x, 95 scope T[] w, 96 scope T[] work, 97 const T alpha, 98 const T beta) @nogc 99 in { 100 assert(T.infinity > alpha && alpha > -1); 101 assert(T.infinity > beta && beta > -1); 102 assert(x.length == w.length); 103 if (x.length) 104 assert(work.length >= (x.length + 1) ^^ 2); 105 } 106 do { 107 if (x.length == 0) 108 return 0; 109 auto s = alpha + beta; 110 auto d = beta - alpha; 111 auto mu0 = exp(double(LN2) * (s + 1) + (lgamma(double(alpha + 1)) + lgamma(double(beta + 1)) - lgamma(double(s + 2)))); 112 x[0] = d / (s + 2); 113 const sd = s * d; 114 foreach (i; 1 .. x.length) 115 { 116 const m_i = T(1) / i; 117 const q = (2 + s * m_i); 118 x[i] = sd * (m_i * m_i) / (q * (2 + (s + 2) * m_i)); 119 w[i] = 4 * (1 + alpha * m_i) * (1 + beta * m_i) * (1 + s * m_i) 120 / ((2 + (s + 1) * m_i) * (2 + (s - 1) * m_i) * (q * q)); 121 } 122 return golubWelsch!T(x, w, work, mu0, alpha == beta); 123 } 124 125 /// 126 unittest 127 { 128 import mir.math.common; 129 130 auto n = 5; 131 auto x = new double[n]; 132 auto w = new double[n]; 133 auto work = new double[(n + 1) ^^ 2]; 134 135 gaussJacobiQuadrature(x, w, work, 2.3, 3.6); 136 137 static immutable xc = 138 [-0.6553677 , 139 -0.29480426, 140 0.09956621, 141 0.47584565, 142 0.78356514]; 143 144 static immutable wc = 145 [0.02262392, 146 0.19871672, 147 0.43585107, 148 0.32146619, 149 0.0615342 ]; 150 151 foreach (i; 0 .. n) 152 { 153 assert(x[i].approxEqual(xc[i])); 154 assert(w[i].approxEqual(wc[i])); 155 } 156 } 157 158 /++ 159 Gauss-Laguerre Quadrature 160 161 Params: 162 x = (out) user-allocated quadrature nodes in ascending order length of `N` 163 w = (out) user-allocated corresponding quadrature weights length of `N` 164 work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` 165 alpha = (optional) parameter '> -1' 166 167 Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. 168 169 +/ 170 size_t gaussLaguerreQuadrature(T)( 171 scope T[] x, 172 scope T[] w, 173 scope T[] work, 174 T alpha = 0) @nogc 175 in { 176 assert(T.infinity > alpha && alpha > -1); 177 assert(x.length == w.length); 178 if (x.length) 179 assert(work.length >= (x.length + 1) ^^ 2); 180 } 181 do { 182 183 auto mu0 = tgamma(double(alpha + 1)); 184 foreach (i; 0 .. x.length) 185 { 186 x[i] = 2 * i + (1 + alpha); 187 w[i] = i * (i + alpha); 188 } 189 return golubWelsch!T(x, w, work, mu0); 190 } 191 192 /// 193 unittest 194 { 195 import mir.math.common; 196 197 auto n = 5; 198 auto x = new double[n]; 199 auto w = new double[n]; 200 auto work = new double[(n + 1) ^^ 2]; 201 202 gaussLaguerreQuadrature(x, w, work); 203 204 static immutable xc = 205 [ 0.26356032, 206 1.41340306, 207 3.59642577, 208 7.08581001, 209 12.64080084]; 210 211 static immutable wc = 212 [5.21755611e-01, 213 3.98666811e-01, 214 7.59424497e-02, 215 3.61175868e-03, 216 2.33699724e-05]; 217 218 foreach (i; 0 .. n) 219 { 220 assert(x[i].approxEqual(xc[i])); 221 assert(w[i].approxEqual(wc[i])); 222 } 223 } 224 225 /++ 226 Gauss-Legendre Quadrature 227 228 Params: 229 x = (out) user-allocated quadrature nodes in ascending order length of `N` 230 w = (out) user-allocated corresponding quadrature weights length of `N` 231 work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` 232 233 Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. 234 235 +/ 236 size_t gaussLegendreQuadrature(T)( 237 scope T[] x, 238 scope T[] w, 239 scope T[] work) @nogc 240 in { 241 assert(x.length == w.length); 242 if (x.length) 243 assert(work.length >= (x.length + 1) ^^ 2); 244 } 245 do { 246 if (x.length == 0) 247 return 0; 248 enum mu0 = 2; 249 x[0] = 0; 250 foreach (i; 1 .. x.length) 251 { 252 const m_i = T(1) / i; 253 x[i] = 0; 254 w[i] = 1 / (4 - (m_i * m_i)); 255 } 256 return golubWelsch!T(x, w, work, mu0, true); 257 } 258 259 /// 260 unittest 261 { 262 import mir.math.common; 263 264 auto n = 5; 265 auto x = new double[n]; 266 auto w = new double[n]; 267 auto work = new double[(n + 1) ^^ 2]; 268 269 gaussLegendreQuadrature(x, w, work); 270 271 static immutable xc = 272 [-0.90617985, 273 -0.53846931, 274 0. , 275 0.53846931, 276 0.90617985]; 277 278 static immutable wc = 279 [0.23692689, 280 0.47862867, 281 0.56888889, 282 0.47862867, 283 0.23692689]; 284 285 foreach (i; 0 .. n) 286 { 287 assert(x[i].approxEqual(xc[i])); 288 assert(w[i].approxEqual(wc[i])); 289 } 290 } 291 292 /++ 293 Gauss-Lobatto Quadrature 294 295 Params: 296 x = (out) user-allocated quadrature nodes in ascending order length of `N` 297 w = (out) user-allocated corresponding quadrature weights length of `N` 298 work = (temp) user-allocated workspace length of greate or equal to `N ^^ 2` 299 alpha = (optional) parameter '> -1' 300 beta = (optional) parameter '> -1' 301 302 Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. 303 304 +/ 305 size_t gaussLobattoQuadrature(T)( 306 scope T[] x, 307 scope T[] w, 308 scope T[] work, 309 const T alpha = 0, 310 const T beta = 0, 311 ) @nogc 312 in { 313 assert(x.length >= 2); 314 assert(x.length == w.length); 315 assert(work.length >= x.length ^^ 2); 316 } 317 do { 318 auto ret = gaussJacobiQuadrature!T(x[1 .. $ - 1], w[1 .. $ - 1], work, alpha + 1, beta + 1); 319 x[0] = -1; 320 x[$ - 1] = +1; 321 auto n = x.length; 322 w[0] = w[$ - 1] = T(2) / (n * (n - 1)); 323 foreach (i; 1 .. x.length - 1) 324 w[i] /= 1 - x[i] * x[i]; 325 return ret; 326 } 327 328 /// 329 unittest 330 { 331 import mir.math.common; 332 333 auto n = 7; 334 auto x = new double[n]; 335 auto w = new double[n]; 336 auto work = new double[n ^^ 2]; 337 338 gaussLobattoQuadrature(x, w, work); 339 340 static immutable xc = 341 [-1, 342 -sqrt(5.0 / 11 + 2.0 / 11 * sqrt(5.0 / 3)), 343 -sqrt(5.0 / 11 - 2.0 / 11 * sqrt(5.0 / 3)), 344 0. , 345 +sqrt(5.0 / 11 - 2.0 / 11 * sqrt(5.0 / 3)), 346 +sqrt(5.0 / 11 + 2.0 / 11 * sqrt(5.0 / 3)), 347 1]; 348 349 static immutable wc = 350 [1.0 / 21, 351 (124.0 - 7 * sqrt(15.0)) / 350, 352 (124.0 + 7 * sqrt(15.0)) / 350, 353 256.0 / 525, 354 (124.0 + 7 * sqrt(15.0)) / 350, 355 (124.0 - 7 * sqrt(15.0)) / 350, 356 1.0 / 21]; 357 358 foreach (i; 0 .. n) 359 assert(x[i].approxEqual(xc[i])); 360 foreach (i; 0 .. n) 361 assert(w[i].approxEqual(wc[i])); 362 } 363 364 // The Golub-Welsch algorithm 365 private size_t golubWelsch(T)( 366 scope T[] alpha_x, 367 scope T[] beta_w, 368 scope T[] work, 369 double mu0, 370 bool symmetrize = false) @nogc 371 in { 372 assert(alpha_x.length == beta_w.length); 373 if (alpha_x.length) 374 assert(work.length >= (alpha_x.length + 1) ^^ 2); 375 foreach (ref b; beta_w[1 .. $]) 376 assert (T.infinity > b && b > 0); 377 } 378 do { 379 pragma(inline, false); 380 auto n = alpha_x.length; 381 if (n == 0) 382 return n; 383 foreach (ref b; beta_w[1 .. n]) 384 b = b.sqrt; 385 auto nq = n * n; 386 import mir.ndslice.slice: sliced; 387 import mir.ndslice.topology: canonical; 388 import mir.lapack: steqr; 389 auto z = work[0 .. nq].sliced(n, n); 390 auto info = steqr('I', alpha_x.sliced, beta_w[1 .. $].sliced, z.canonical, work[nq .. $].sliced); 391 foreach (i; 0 .. n) 392 { 393 auto zi0 = z[i, 0]; 394 beta_w[i] = zi0 * zi0 * mu0; 395 } 396 if (symmetrize) 397 { 398 auto h = n / 2; 399 alias x = alpha_x; 400 alias w = beta_w; 401 foreach (i; 0 .. h) 402 { 403 x[i] = -(x[n - (i + 1)] = T(0.5) * (x[n - (i + 1)] - x[i])); 404 w[i] = +(w[n - (i + 1)] = T(0.5) * (w[n - (i + 1)] + w[i])); 405 } 406 if (n % 2) 407 { 408 x[h] = 0; 409 } 410 } 411 return info; 412 }