diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index cc2e3b79..c8422144 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -430,13 +430,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests expected, Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b), 1e-10, - "DET Integral e^(-x^2 /2) from {0} to {1}", a, b); + "DET Integral of e^(-x^2 /2) from {0} to {1}", a, b); Assert.AreEqual( expected, Integrate.GaussKronrod((x) => Math.Exp(-x * x / 2), a, b), 1e-10, - "GK Integral e^(-x^2 /2) from {0} to {1}", a, b); + "GK Integral of e^(-x^2 /2) from {0} to {1}", a, b); + + Assert.AreEqual( + expected, + Integrate.GausLegendre((x) => Math.Exp(-x * x / 2), a, b, order: 128), + 1e-10, + "GL Integral of e^(-x^2 /2) from {0} to {1}", a, b); } // integral_(-oo)^(oo) sin(pi x) / (pi x) dx = 1 / pi integral_(-oo)^(oo) sin(x) / x dx @@ -453,13 +459,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests expected, factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), 1e-10, - "DET Integral sin(pi*x)/(pi*x) from -oo to oo"); + "DET Integral of sin(pi*x)/(pi*x) from -oo to oo"); Assert.AreEqual( expected, factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b), 1e-10, - "GK Integral sin(pi*x)/(pi*x) from -oo to oo"); + "GK Integral of sin(pi*x)/(pi*x) from -oo to oo"); + + Assert.AreEqual( + expected, + factor * Integrate.GausLegendre((x) => 1 / (1 + x * x), a, b, order: 128), + 1e-10, + "GL Integral of sin(pi*x)/(pi*x) from -oo to oo"); } // integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð @@ -473,30 +485,43 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests var expected = new Complex(r, i); var actualDET = ContourIntegrate.DoubleExponential((x) => 1 / new Complex(1, x * x), a, b); var actualGK = ContourIntegrate.GaussKronrod((x) => 1 / new Complex(1, x * x), a, b); + var actualGL = ContourIntegrate.GaussLegendre((x) => 1 / new Complex(1, x * x), a, b, order: 128); Assert.AreEqual( expected.Real, actualDET.Real, 1e-10, - "DET Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "DET Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Imaginary, actualDET.Imaginary, 1e-10, - "DET Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "DET Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Real, actualGK.Real, 1e-10, - "GK Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "GK Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Imaginary, actualGK.Imaginary, 1e-10, - "GK Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "GK Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + + Assert.AreEqual( + expected.Real, + actualGL.Real, + 1e-10, + "GL Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + + Assert.AreEqual( + expected.Imaginary, + actualGL.Imaginary, + 1e-10, + "GL Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); } } } diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index cc7944b3..6057ed6a 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -156,6 +156,82 @@ namespace MathNet.Numerics } } + /// + /// Approximation of the definite integral of an analytic smooth function by Gauss-Legendre quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. + /// + /// The analytic smooth function to integrate. + /// Where the interval starts. + /// Where the interval stops. + /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. + /// Approximation of the finite integral in the given interval. + public static double GausLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) + { + // Reference: + // Formula used for variable subsitution from + // 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140. + // 2. quadgk.m, GNU Octave + + if (intervalBegin > intervalEnd) + { + return -GaussLegendreRule.Integrate(f, intervalEnd, intervalBegin, order); + } + + // (-oo, oo) => [-1, 1] + // + // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = t / (1 - t^2) + // g'(t) = (1 + t^2) / (1 - t^2)^2 + if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) + { + Func u = (t) => + { + return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); + }; + return GaussLegendreRule.Integrate(u, -1, 1, order); + } + // [a, oo) => [0, 1] + // + // integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt + // = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 - s) + // g'(s) = 1 / (1 - s)^2 + else if (double.IsInfinity(intervalEnd)) + { + Func u = (s) => + { + return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); + }; + return GaussLegendreRule.Integrate(u, 0, 1, order); + } + // (-oo, b] => [-1, 0] + // + // integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt + // = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 + s) + // g'(s) = 1 / (1 + s)^2 + else if (double.IsInfinity(intervalBegin)) + { + Func u = (s) => + { + return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); + }; + return GaussLegendreRule.Integrate(u, -1, 0, order); + } + // [a, b] => [-1, 1] + // + // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2 + // g'(t) = 3 / 4 * (b - a) * (1 - t^2) + else + { + Func u = (t) => + { + return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t); + }; + return GaussLegendreRule.Integrate(u, -1, 1, order); + } + } + /// /// Approximation of the definite integral of an analytic smooth function by Gauss-Kronrod quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. /// @@ -163,8 +239,8 @@ namespace MathNet.Numerics /// Where the interval starts. /// Where the interval stops. /// The expected relative accuracy of the approximation. - /// The maximum number of interval splittings permitted before stopping - /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points + /// The maximum number of interval splittings permitted before stopping. + /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points. /// Approximation of the finite integral in the given interval. public static double GaussKronrod(Func f, double intervalBegin, double intervalEnd, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15) { @@ -270,6 +346,82 @@ namespace MathNet.Numerics } } + /// + /// Approximation of the definite integral of an analytic smooth complex function by double-exponential quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts. + /// Where the interval stops. + /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. + /// Approximation of the finite integral in the given interval. + public static Complex GaussLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) + { + // Reference: + // Formula used for variable subsitution from + // 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140. + // 2. quadgk.m, GNU Octave + + if (intervalBegin > intervalEnd) + { + return -GaussLegendre(f, intervalEnd, intervalBegin, order); + } + + // (-oo, oo) => [-1, 1] + // + // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = t / (1 - t^2) + // g'(t) = (1 + t^2) / (1 - t^2)^2 + if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) + { + Func u = (t) => + { + return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); + }; + return GaussLegendreRule.ContourIntegrate(u, -1, 1, order); + } + // [a, oo) => [0, 1] + // + // integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt + // = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 - s) + // g'(s) = 1 / (1 - s)^2 + else if (double.IsInfinity(intervalEnd)) + { + Func u = (s) => + { + return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); + }; + return GaussLegendreRule.ContourIntegrate(u, 0, 1, order); + } + // (-oo, b] => [-1, 0] + // + // integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt + // = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 + s) + // g'(s) = 1 / (1 + s)^2 + else if (double.IsInfinity(intervalBegin)) + { + Func u = (s) => + { + return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); + }; + return GaussLegendreRule.ContourIntegrate(u, -1, 0, order); + } + // [a, b] => [-1, 1] + // + // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2 + // g'(t) = 3 / 4 * (b - a) * (1 - t^2) + else + { + Func u = (t) => + { + return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t); + }; + return GaussLegendreRule.ContourIntegrate(u, -1, 1, order); + } + } + /// /// Approximation of the definite integral of an analytic smooth function by Gauss-Kronrod quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. /// diff --git a/src/Numerics/Integration/GaussLegendreRule.cs b/src/Numerics/Integration/GaussLegendreRule.cs index 60bed83f..fab03a4e 100644 --- a/src/Numerics/Integration/GaussLegendreRule.cs +++ b/src/Numerics/Integration/GaussLegendreRule.cs @@ -28,6 +28,7 @@ // using System; +using System.Numerics; using MathNet.Numerics.Integration.GaussRule; namespace MathNet.Numerics.Integration @@ -166,6 +167,48 @@ namespace MathNet.Numerics.Integration return a*sum; } + /// + /// Approximates a definite integral using an Nth order Gauss-Legendre rule. + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts, exclusive and finite. + /// Where the interval ends, exclusive and finite. + /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrate(Func f, double invervalBegin, double invervalEnd, int order) + { + GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.GetGaussPoint(order); + + Complex sum; + double ax; + int i; + int m = (order + 1) >> 1; + + double a = 0.5 * (invervalEnd - invervalBegin); + double b = 0.5 * (invervalEnd + invervalBegin); + + if (order.IsOdd()) + { + sum = gaussLegendrePoint.Weights[0] * f(b); + for (i = 1; i < m; i++) + { + ax = a * gaussLegendrePoint.Abscissas[i]; + sum += gaussLegendrePoint.Weights[i] * (f(b + ax) + f(b - ax)); + } + } + else + { + sum = 0.0; + for (i = 0; i < m; i++) + { + ax = a * gaussLegendrePoint.Abscissas[i]; + sum += gaussLegendrePoint.Weights[i] * (f(b + ax) + f(b - ax)); + } + } + + return a * sum; + } + /// /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. ///