diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index e5e55ca2..1f600cfa 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -27,9 +27,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; using MathNet.Numerics.Integration; using NUnit.Framework; +using System; +using System.Numerics; namespace MathNet.Numerics.UnitTests.IntegrationTests { @@ -59,7 +60,17 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests { return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y)); } - + + /// + /// Test Function: f(x,y) = 1 / (1 + x^2) + /// + /// First input value. + /// Function result. + private static double TargetFunctionC(double x) + { + return 1 / (1 + x * x); + } + /// /// Test Function Start point. /// @@ -80,6 +91,16 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// private const double StopB = 1; + /// + /// Test Function Start point. + /// + private const double StartC = double.NegativeInfinity; + + /// + /// Test Function Stop point. + /// + private const double StopC = double.PositiveInfinity; + /// /// Target area square. /// @@ -90,6 +111,11 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// private const double TargetAreaB = 11.7078776759298776163; + /// + /// Target area. + /// + private const double TargetAreaC = Constants.Pi; + /// /// Test Integrate facade for simple use cases. /// @@ -119,6 +145,18 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests TargetAreaB, 1e-10, "Rectangle, Gauss-Legendre Order 22"); + + Assert.AreEqual( + TargetAreaC, + Integrate.DoubleExponential(TargetFunctionC, StartC, StopC), + 1e-5, + "Integral by substitution"); + + Assert.AreEqual( + TargetAreaC, + Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10), + 1e-10, + "Integral by substitution, Target 1e-10"); } /// @@ -326,7 +364,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests { Assert.AreEqual( expected, - Integrate.OnOpenInterval((x) => Math.Exp(-x * x / 2), a, b), + Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b), 1e-10, "Integral e^(-x^2 /2) from {0} to {1}", a, b); } @@ -343,9 +381,33 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests { Assert.AreEqual( expected, - factor * Integrate.OnOpenInterval((x) => 1 / (1 + x * x), a, b), + factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), 1e-10, "Integral sin(pi*x)/(pi*x) from -oo to oo"); } + + // integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð + [TestCase(double.NegativeInfinity, double.PositiveInfinity, 2.2214414690791831235, -2.2214414690791831235)] + // integral_(0)^(oo) 1/(1 + j x^2) dx = -1/2 (-1)^(3/4) ¥ð + [TestCase(0, double.PositiveInfinity, 1.1107207345395915618, -1.1107207345395915618)] + // integral_(-oo)^(0) 1/(1 + j x^2) dx = -1/2 (-1)^(3/4) ¥ð + [TestCase(double.NegativeInfinity, 0, 1.1107207345395915618, -1.1107207345395915618)] + public void TestContourIntegralBySubstitution(double a, double b, double r, double i) + { + var expected = new Complex(r, i); + var actual = ContourIntegrate.DoubleExponential((x) => 1 / new Complex(1, x * x), a, b); + + Assert.AreEqual( + expected.Real, + actual.Real, + 1e-10, + "Integral e^(-x^2 /2) / (1 + j e^x) from {0} to {1}", a, b); + + Assert.AreEqual( + expected.Imaginary, + actual.Imaginary, + 1e-10, + "Integral 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 461516e0..1a73d75a 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -28,6 +28,7 @@ // using System; +using System.Numerics; using MathNet.Numerics.Integration; namespace MathNet.Numerics @@ -50,15 +51,44 @@ namespace MathNet.Numerics return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError); } + /// + /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. + /// + /// The 2-dimensional analytic smooth function to integrate. + /// Where the interval starts for the first (inside) integral, exclusive and finite. + /// Where the interval ends for the first (inside) integral, exclusive and finite. + /// Where the interval starts for the second (outside) integral, exclusive and finite. + /// /// Where the interval ends for the second (outside) integral, 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 double OnRectangle(Func f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order) + { + return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order); + } + + /// + /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. + /// + /// The 2-dimensional analytic smooth function to integrate. + /// Where the interval starts for the first (inside) integral, exclusive and finite. + /// Where the interval ends for the first (inside) integral, exclusive and finite. + /// Where the interval starts for the second (outside) integral, exclusive and finite. + /// /// Where the interval ends for the second (outside) integral, exclusive and finite. + /// Approximation of the finite integral in the given interval. + public static double OnRectangle(Func f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB) + { + return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, 32); + } + /// /// Approximation of the definite integral of an analytic smooth function by substitution. 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, inclusive and finite. - /// Where the interval stops, inclusive and finite. + /// Where the interval starts. + /// Where the interval stops. /// The expected relative accuracy of the approximation. /// Approximation of the finite integral in the given interval. - public static double OnOpenInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) + public static double DoubleExponential(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) { // Reference: // Formula used for variable subsitution from @@ -67,7 +97,7 @@ namespace MathNet.Numerics if (intervalBegin > intervalEnd) { - return -OnOpenInterval(f, intervalEnd, intervalBegin, targetAbsoluteError); + return -DoubleExponential(f, intervalEnd, intervalBegin, targetAbsoluteError); } // (-oo, oo) => [-1, 1] @@ -81,7 +111,7 @@ namespace MathNet.Numerics { return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); }; - return OnClosedInterval(u, -1d, 1d, targetAbsoluteError); + return DoubleExponentialTransformation.Integrate(u, -1, 1, targetAbsoluteError); } // [a, oo) => [0, 1] // @@ -95,7 +125,7 @@ namespace MathNet.Numerics { return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); }; - return OnClosedInterval(u, 0d, 1d, targetAbsoluteError); + return DoubleExponentialTransformation.Integrate(u, 0, 1, targetAbsoluteError); } // (-oo, b] => [-1, 0] // @@ -109,7 +139,7 @@ namespace MathNet.Numerics { return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); }; - return OnClosedInterval(u, -1d, 0d, targetAbsoluteError); + return DoubleExponentialTransformation.Integrate(u, -1, 0, targetAbsoluteError); } // [a, b] => [-1, 1] // @@ -122,37 +152,90 @@ namespace MathNet.Numerics { return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t); }; - return OnClosedInterval(u, -1d, 1d, targetAbsoluteError); + return DoubleExponentialTransformation.Integrate(u, -1, 1, targetAbsoluteError); } } + } + /// + /// Numerical Contour Integration over a real variable, of a complex-valued function. + /// + public static class ContourIntegrate + { /// - /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. + /// Approximation of the definite integral of an analytic smooth complex function by substitution. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. /// - /// The 2-dimensional analytic smooth function to integrate. - /// Where the interval starts for the first (inside) integral, exclusive and finite. - /// Where the interval ends for the first (inside) integral, exclusive and finite. - /// Where the interval starts for the second (outside) integral, exclusive and finite. - /// /// Where the interval ends for the second (outside) integral, 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. + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts. + /// Where the interval stops. + /// The expected relative accuracy of the approximation. /// Approximation of the finite integral in the given interval. - public static double OnRectangle(Func f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order) + public static Complex DoubleExponential(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) { - return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order); - } + // 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 - /// - /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. - /// - /// The 2-dimensional analytic smooth function to integrate. - /// Where the interval starts for the first (inside) integral, exclusive and finite. - /// Where the interval ends for the first (inside) integral, exclusive and finite. - /// Where the interval starts for the second (outside) integral, exclusive and finite. - /// /// Where the interval ends for the second (outside) integral, exclusive and finite. - /// Approximation of the finite integral in the given interval. - public static double OnRectangle(Func f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB) - { - return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, 32); + if (intervalBegin > intervalEnd) + { + return -DoubleExponential(f, intervalEnd, intervalBegin, targetAbsoluteError); + } + + // (-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 DoubleExponentialTransformation.ContourIntegrate(u, -1, 1, targetAbsoluteError); + } + // [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 DoubleExponentialTransformation.ContourIntegrate(u, 0, 1, targetAbsoluteError); + } + // (-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 DoubleExponentialTransformation.ContourIntegrate(u, -1, 0, targetAbsoluteError); + } + // [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 DoubleExponentialTransformation.ContourIntegrate(u, -1, 1, targetAbsoluteError); + } } } } diff --git a/src/Numerics/Integration/DoubleExponentialTransformation.cs b/src/Numerics/Integration/DoubleExponentialTransformation.cs index bd9ddd83..9072743c 100644 --- a/src/Numerics/Integration/DoubleExponentialTransformation.cs +++ b/src/Numerics/Integration/DoubleExponentialTransformation.cs @@ -29,6 +29,7 @@ using System; using System.Linq; +using System.Numerics; namespace MathNet.Numerics.Integration { @@ -63,6 +64,25 @@ namespace MathNet.Numerics.Integration targetRelativeError); } + /// + /// Approximate the integral by the double exponential transformation + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts, inclusive and finite. + /// Where the interval stops, inclusive and finite. + /// The expected relative accuracy of the approximation. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrate(Func f, double intervalBegin, double intervalEnd, double targetRelativeError) + { + return NewtonCotesTrapeziumRule.ContourIntegrateAdaptiveTransformedOdd( + f, + intervalBegin, intervalEnd, + Enumerable.Range(0, NumberOfMaximumLevels).Select(EvaluateAbcissas), + Enumerable.Range(0, NumberOfMaximumLevels).Select(EvaluateWeights), + 1.0, + targetRelativeError); + } + /// /// Compute the abscissa vector for a single level. /// diff --git a/src/Numerics/Integration/NewtonCotesTrapeziumRule.cs b/src/Numerics/Integration/NewtonCotesTrapeziumRule.cs index d82e08a7..b29037b1 100644 --- a/src/Numerics/Integration/NewtonCotesTrapeziumRule.cs +++ b/src/Numerics/Integration/NewtonCotesTrapeziumRule.cs @@ -29,6 +29,7 @@ using System; using System.Collections.Generic; +using System.Numerics; using MathNet.Numerics.Properties; namespace MathNet.Numerics.Integration @@ -58,6 +59,23 @@ namespace MathNet.Numerics.Integration return (intervalEnd - intervalBegin)/2*(f(intervalBegin) + f(intervalEnd)); } + /// + /// Direct 2-point approximation of the definite integral in the provided interval by the trapezium rule. + /// + /// The analytic smooth complex function to integrate, defined on real domain. + /// Where the interval starts, inclusive and finite. + /// Where the interval stops, inclusive and finite. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrateTwoPoint(Func f, double intervalBegin, double intervalEnd) + { + if (f == null) + { + throw new ArgumentNullException(nameof(f)); + } + + return (intervalEnd - intervalBegin) / 2 * (f(intervalBegin) + f(intervalEnd)); + } + /// /// Composite N-point approximation of the definite integral in the provided interval by the trapezium rule. /// @@ -92,6 +110,40 @@ namespace MathNet.Numerics.Integration return step*sum; } + /// + /// Composite N-point approximation of the definite integral in the provided interval by the trapezium rule. + /// + /// The analytic smooth complex function to integrate, defined on real domain. + /// Where the interval starts, inclusive and finite. + /// Where the interval stops, inclusive and finite. + /// Number of composite subdivision partitions. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrateComposite(Func f, double intervalBegin, double intervalEnd, int numberOfPartitions) + { + if (f == null) + { + throw new ArgumentNullException(nameof(f)); + } + + if (numberOfPartitions <= 0) + { + throw new ArgumentOutOfRangeException(nameof(numberOfPartitions), Resources.ArgumentPositive); + } + + double step = (intervalEnd - intervalBegin) / numberOfPartitions; + + double offset = step; + Complex sum = 0.5 * (f(intervalBegin) + f(intervalEnd)); + for (int i = 0; i < numberOfPartitions - 1; i++) + { + // NOTE (ruegg, 2009-01-07): Do not combine intervalBegin and offset (numerical stability!) + sum += f(intervalBegin + offset); + offset += step; + } + + return step * sum; + } + /// /// Adaptive approximation of the definite integral in the provided interval by the trapezium rule. /// @@ -132,6 +184,47 @@ namespace MathNet.Numerics.Integration return sum; } + /// + /// Adaptive approximation of the definite integral in the provided interval by the trapezium rule. + /// + /// The analytic smooth complex function to integrate, define don real domain. + /// Where the interval starts, inclusive and finite. + /// Where the interval stops, inclusive and finite. + /// The expected accuracy of the approximation. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrateAdaptive(Func f, double intervalBegin, double intervalEnd, double targetError) + { + if (f == null) + { + throw new ArgumentNullException(nameof(f)); + } + + int numberOfPartitions = 1; + double step = intervalEnd - intervalBegin; + Complex sum = 0.5 * step * (f(intervalBegin) + f(intervalEnd)); + for (int k = 0; k < 20; k++) + { + Complex midpointsum = 0; + for (int i = 0; i < numberOfPartitions; i++) + { + midpointsum += f(intervalBegin + ((i + 0.5) * step)); + } + + midpointsum *= step; + sum = 0.5 * (sum + midpointsum); + step *= 0.5; + numberOfPartitions *= 2; + + if (sum.AlmostEqualRelative(midpointsum, targetError)) + { + break; + } + } + + return sum; + } + + /// /// Adaptive approximation of the definite integral by the trapezium rule. /// @@ -230,5 +323,105 @@ namespace MathNet.Numerics.Integration return sum*linearSlope; } } + + /// + /// Adaptive approximation of the definite integral by the trapezium rule. + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts, inclusive and finite. + /// Where the interval stops, inclusive and finite. + /// Abscissa vector per level provider. + /// Weight vector per level provider. + /// First Level Step + /// The expected relative accuracy of the approximation. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrateAdaptiveTransformedOdd( + Func f, + double intervalBegin, double intervalEnd, + IEnumerable levelAbscissas, IEnumerable levelWeights, + double levelOneStep, double targetRelativeError) + { + if (f == null) + { + throw new ArgumentNullException(nameof(f)); + } + + if (levelAbscissas == null) + { + throw new ArgumentNullException(nameof(levelAbscissas)); + } + + if (levelWeights == null) + { + throw new ArgumentNullException(nameof(levelWeights)); + } + + double linearSlope = 0.5 * (intervalEnd - intervalBegin); + double linearOffset = 0.5 * (intervalEnd + intervalBegin); + targetRelativeError /= 5 * linearSlope; + + using (var abcissasIterator = levelAbscissas.GetEnumerator()) + using (var weightsIterator = levelWeights.GetEnumerator()) + { + double step = levelOneStep; + + // First Level + abcissasIterator.MoveNext(); + weightsIterator.MoveNext(); + double[] abcissasL1 = abcissasIterator.Current; + double[] weightsL1 = weightsIterator.Current; + + Complex sum = f(linearOffset) * weightsL1[0]; + for (int i = 1; i < abcissasL1.Length; i++) + { + sum += weightsL1[i] * (f((linearSlope * abcissasL1[i]) + linearOffset) + f(-(linearSlope * abcissasL1[i]) + linearOffset)); + } + + sum *= step; + + // Additional Levels + double previousDelta = double.MaxValue; + for (int level = 1; abcissasIterator.MoveNext() && weightsIterator.MoveNext(); level++) + { + double[] abcissas = abcissasIterator.Current; + double[] weights = weightsIterator.Current; + + Complex midpointsum = 0; + for (int i = 0; i < abcissas.Length; i++) + { + midpointsum += weights[i] * (f((linearSlope * abcissas[i]) + linearOffset) + f(-(linearSlope * abcissas[i]) + linearOffset)); + } + + midpointsum *= step; + sum = 0.5 * (sum + midpointsum); + step *= 0.5; + + double delta = Complex.Abs(sum - midpointsum); + + if (level == 1) + { + previousDelta = delta; + continue; + } + + double r = Math.Log(delta) / Math.Log(previousDelta); + previousDelta = delta; + + if (r > 1.9 && r < 2.1) + { + // convergence region + delta = Math.Sqrt(delta); + } + + if (sum.Real.AlmostEqualNormRelative(midpointsum.Real, delta, targetRelativeError) + && sum.Imaginary.AlmostEqualNormRelative(midpointsum.Imaginary, delta, targetRelativeError)) + { + break; + } + } + + return sum * linearSlope; + } + } } }