diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index ad362cde..546a46df 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 { @@ -50,7 +51,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests } /// - /// Test Function: f(x,y) = exp(-x/5) (2 + sin(x * y)) + /// Test Function: f(x,y) = exp(-x/5) (2 + sin(2 * y)) /// /// First input value. /// Second input value. @@ -60,6 +61,56 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y)); } + /// + /// Test Function: f(x) = 1 / (1 + x^2) + /// + /// First input value. + /// Function result. + private static double TargetFunctionC(double x) + { + return 1 / (1 + x * x); + } + + /// + /// Test Function: f(x) = log(x) + /// + /// First input value. + /// Function result. + private static double TargetFunctionD(double x) + { + return Math.Log(x); + } + + /// + /// Test Function: f(x) = log^2(x) + /// + /// First input value. + /// Function result. + private static double TargetFunctionE(double x) + { + return Math.Log(x) * Math.Log(x); + } + + /// + /// Test Function: f(x) = e^(-x) cos(x) + /// + /// First input value. + /// Function result. + private static double TargetFunctionF(double x) + { + return Math.Exp(-x) * Math.Cos(x); + } + + /// + /// Test Function: f(x) = sqrt(x)/sqrt(1-x^2) + /// + /// First input value. + /// Function result. + private static double TargetFunctionG(double x) + { + return Math.Sqrt(x) / Math.Sqrt(1 - x * x); + } + /// /// Test Function Start point. /// @@ -80,6 +131,56 @@ 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; + + /// + /// Test Function Start point. + /// + private const double StartD = 0; + + /// + /// Test Function Stop point. + /// + private const double StopD = 1; + + /// + /// Test Function Start point. + /// + private const double StartE = 0; + + /// + /// Test Function Stop point. + /// + private const double StopE = 1; + + /// + /// Test Function Start point. + /// + private const double StartF = 0; + + /// + /// Test Function Stop point. + /// + private const double StopF = double.PositiveInfinity; + + /// + /// Test Function Start point. + /// + private const double StartG = 0; + + /// + /// Test Function Stop point. + /// + private const double StopG = 1; + /// /// Target area square. /// @@ -90,17 +191,45 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// private const double TargetAreaB = 11.7078776759298776163; + /// + /// Target area. + /// + private const double TargetAreaC = Constants.Pi; + + /// + /// Target area. + /// + private const double TargetAreaD = -1; + + /// + /// Target area. + /// + private const double TargetAreaE = 2; + + /// + /// Target area. + /// + private const double TargetAreaF = 0.5; + + /// + /// Target area. + /// + private const double TargetAreaG = 1.1981402347355922074; + /// /// Test Integrate facade for simple use cases. /// [Test] public void TestIntegrateFacade() { + // TargetFunctionA + // integral_(0)^(10) exp(-x/5) (2 + sin(2 x)) dx = 9.1082 + Assert.AreEqual( TargetAreaA, Integrate.OnClosedInterval(TargetFunctionA, StartA, StopA), 1e-5, - "Interval"); + "Interval, Target 1e-08"); Assert.AreEqual( TargetAreaA, @@ -108,17 +237,152 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests 1e-10, "Interval, Target 1e-10"); + // TargetFunctionB + // integral_(0)^(1) integral_(0)^(10) exp(-x/5) (2 + sin(2 y)) dx dy = 11.7079 + Assert.AreEqual( Integrate.OnRectangle(TargetFunctionB, StartA, StopA, StartB, StopB), TargetAreaB, 1e-12, - "Rectangle"); + "Rectangle, order 32"); Assert.AreEqual( Integrate.OnRectangle(TargetFunctionB, StartA, StopA, StartB, StopB, 22), TargetAreaB, 1e-10, - "Rectangle, Gauss-Legendre Order 22"); + "Rectangle, Order 22"); + + // TargetFunctionC + // integral_(-oo)^(oo) 1/(1 + x^2) dx = pi + + Assert.AreEqual( + TargetAreaC, + Integrate.DoubleExponential(TargetFunctionC, StartC, StopC), + 1e-5, + "DoubleExponential, 1/(1 + x^2)"); + + Assert.AreEqual( + TargetAreaC, + Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10), + 1e-10, + "DoubleExponential, 1/(1 + x^2)"); + + // TargetFunctionD + // integral_(0)^(1) log(x) dx = -1 + + Assert.AreEqual( + TargetAreaD, + Integrate.DoubleExponential(TargetFunctionD, StartD, StopD), + 1e-10, + "DoubleExponential, log(x)"); + + Assert.AreEqual( + TargetAreaD, + Integrate.GaussLegendre(TargetFunctionD, StartD, StopD, order: 1024), + 1e-10, + "GaussLegendre, log(x), order 1024"); + + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 15), + 1e-10, + "GaussKronrod, log(x), order 15"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 21), + 1e-10, + "GaussKronrod, log(x), order 21"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 31), + 1e-10, + "GaussKronrod, log(x), order 31"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 41), + 1e-10, + "GaussKronrod, log(x), order 41"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 51), + 1e-10, + "GaussKronrod, log(x), order 51"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 61), + 1e-10, + "GaussKronrod, log(x), order 61"); + + double error, L1; + var Q = Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, out error, out L1, 1e-10, order: 15); + Assert.AreEqual( + Math.Abs(TargetAreaD), + Math.Abs(L1), + 1e-10, + "GaussKronrod, L1"); + + // TargetFunctionE + // integral_(0)^(1) log^2(x) dx = 2 + + Assert.AreEqual( + TargetAreaE, + Integrate.DoubleExponential(TargetFunctionE, StartE, StopE), + 1e-10, + "DoubleExponential, log^2(x)"); + + Assert.AreEqual( + TargetAreaE, + Integrate.GaussLegendre(TargetFunctionE, StartE, StopE, order: 128), + 1e-5, + "GaussLegendre, log^2(x), order 128"); + + Assert.AreEqual( + TargetAreaE, + Integrate.GaussKronrod(TargetFunctionE, StartE, StopE, 1e-10, order: 15), + 1e-10, + "GaussKronrod, log^2(x), order 15"); + + // TargetFunctionF + // integral_(0)^(oo) exp(-x) cos(x) dx = 1/2 + + Assert.AreEqual( + TargetAreaF, + Integrate.DoubleExponential(TargetFunctionF, StartF, StopF), + 1e-10, + "DoubleExponential, e^(-x) cos(x)"); + + Assert.AreEqual( + TargetAreaF, + Integrate.GaussLegendre(TargetFunctionF, StartF, StopF, order: 128), + 1e-10, + "GaussLegendre, e^(-x) cos(x), order 128"); + + Assert.AreEqual( + TargetAreaF, + Integrate.GaussKronrod(TargetFunctionF, StartF, StopF, 1e-10, order: 15), + 1e-10, + "GaussKronrod, e^(-x) cos(x), order 15"); + + // TargetFunctionG + // integral_(0)^(1) sqrt(x)/sqrt(1 - x^2) dx = 1.19814 + + Assert.AreEqual( + TargetAreaG, + Integrate.DoubleExponential(TargetFunctionG, StartG, StopG), + 1e-5, + "DoubleExponential, sqrt(x)/sqrt(1 - x^2)"); + + Assert.AreEqual( + TargetAreaG, + Integrate.GaussLegendre(TargetFunctionG, StartG, StopG, order: 128), + 1e-10, + "GaussLegendre, sqrt(x)/sqrt(1 - x^2), order 128"); + + Assert.AreEqual( + TargetAreaG, + Integrate.GaussKronrod(TargetFunctionG, StartG, StopG, 1e-10, order: 15), + 1e-10, + "GaussKronrod, sqrt(x)/sqrt(1 - x^2), order 15"); } /// @@ -285,7 +549,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests for (int i = 0; i < gaussLegendre.Order; i++) { - Assert.AreEqual(gaussLegendre.GetAbscissa(i),abscissa[i]); + Assert.AreEqual(gaussLegendre.GetAbscissa(i), abscissa[i]); Assert.AreEqual(gaussLegendre.GetWeight(i), weight[i]); } } @@ -311,5 +575,132 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order); Assert.AreEqual(gaussLegendre.IntervalEnd, StopA); } + + /// + /// Gauss-Kronrod rule supports integration. + /// + /// Defines an Nth order Gauss-Kronrod rule. The order also defines the number of abscissas and weights for the rule. + [TestCase(3)] + [TestCase(4)] + [TestCase(5)] + [TestCase(6)] + [TestCase(101)] + [TestCase(201)] + public void TestGaussKronrodRuleIntegration(int order) + { + double appoximateArea = GaussKronrodRule.Integrate(TargetFunctionA, StartA, StopA, out _, out _, order: order); + double relativeError = Math.Abs(TargetAreaA - appoximateArea) / TargetAreaA; + Assert.Less(relativeError, 5e-16); + } + + // integral_(-oo)^(oo) exp(-x^2/2) dx = sqrt(2 ¥ð) + // integral_(-oo)^(0) exp(-x^2/2) dx = sqrt(¥ð/2) + // integral_(0)^(oo exp(-x^2/2) dx = sqrt(¥ð/2) + // integral_(-1)^(1) exp(-x^2/2) dx = sqrt(2 ¥ð) erf(1/sqrt(2)) + // integral_(1)^(0) exp(-x^2/2) dx = -sqrt(¥ð/2) erf(1/sqrt(2)) + [TestCase(double.NegativeInfinity, double.PositiveInfinity, Constants.Sqrt2Pi)] + [TestCase(double.NegativeInfinity, 0, Constants.SqrtPiOver2)] + [TestCase(0, double.PositiveInfinity, Constants.SqrtPiOver2)] + [TestCase(-1, 1, 1.7112487837842976063)] + [TestCase(1, 0, -0.85562439189214880317)] + public void TestIntegralOfGaussian(double a, double b, double expected) + { + Assert.AreEqual( + expected, + Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b), + 1e-10, + "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 of e^(-x^2 /2) from {0} to {1}", a, b); + + Assert.AreEqual( + expected, + Integrate.GaussLegendre((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 + // = 1 / pi integral_(oo)^(oo) 1 / (1 + t^2) dt + // = 1 + // or = 2 / pi integral_(0)^(oo) 1 / (1 + t^2) dt + // or = 2 / pi integral_(-oo)^(0) 1 / (1 + t^2) dt + [TestCase(double.NegativeInfinity, double.PositiveInfinity, 1, Constants.InvPi)] + [TestCase(0, double.PositiveInfinity, 1, Constants.TwoInvPi)] + [TestCase(double.NegativeInfinity, 0, 1, Constants.TwoInvPi)] + public void TestIntegralOfSinc(double a, double b, double expected, double factor) + { + Assert.AreEqual( + expected, + factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), + 1e-10, + "DET Integral of sin(pi*x)/(pi*x) from {0} to {1}", a, b); + + Assert.AreEqual( + expected, + factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b), + 1e-10, + "GK Integral of sin(pi*x)/(pi*x) from {0} to {1}", a, b); + + Assert.AreEqual( + expected, + factor * Integrate.GaussLegendre((x) => 1 / (1 + x * x), a, b, order: 128), + 1e-10, + "GL Integral of sin(pi*x)/(pi*x) from {0} to {1}", a, b); + } + + // integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð + // integral_(0)^(oo) 1/(1 + j x^2) dx = -1/2 (-1)^(3/4) ¥ð + // integral_(-oo)^(0) 1/(1 + j x^2) dx = -1/2 (-1)^(3/4) ¥ð + [TestCase(double.NegativeInfinity, double.PositiveInfinity, 2.2214414690791831235, -2.2214414690791831235)] + [TestCase(0, double.PositiveInfinity, 1.1107207345395915618, -1.1107207345395915618)] + [TestCase(double.NegativeInfinity, 0, 1.1107207345395915618, -1.1107207345395915618)] + public void TestContourIntegral(double a, double b, double r, double i) + { + 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 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 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 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 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 2d708fdf..b7218248 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 @@ -90,5 +91,361 @@ namespace MathNet.Numerics { return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, 32); } + + /// + /// Approximation of the definite integral of an analytic smooth 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 function to integrate. + /// 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 DoubleExponential(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) + { + // 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 -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.Integrate(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.Integrate(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.Integrate(u, -1, 0, targetAbsoluteError); + } + else + { + return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError); + } + } + + /// + /// 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 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.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. + /// + /// The analytic smooth function to integrate. + /// 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. + /// 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) + { + return GaussKronrodRule.Integrate(f, intervalBegin, intervalEnd, out _, out _, targetRelativeError: targetRelativeError, maximumDepth: maximumDepth, order: 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. + /// + /// The analytic smooth function to integrate. + /// Where the interval starts. + /// Where the interval stops. + /// The difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation + /// The L1 norm of the result, if there is a significant difference between this and the returned value, then the result is likely to be ill-conditioned. + /// 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, 21, 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, out double error, out double L1Norm, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15) + { + return GaussKronrodRule.Integrate(f, intervalBegin, intervalEnd, out error, out L1Norm, targetRelativeError: targetRelativeError, maximumDepth: maximumDepth, order: order); + } + } + + /// + /// Numerical Contour Integration of a complex-valued function over a real variable,. + /// + public static class ContourIntegrate + { + /// + /// 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. + /// The expected relative accuracy of the approximation. + /// Approximation of the finite integral in the given interval. + public static Complex DoubleExponential(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) + { + // 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 -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); + } + else + { + return DoubleExponentialTransformation.ContourIntegrate(f, intervalBegin, intervalEnd, targetAbsoluteError); + } + } + + /// + /// 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. + /// + /// 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. + /// The maximum number of interval splittings permitted before stopping + /// The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points + /// Approximation of the finite integral in the given interval. + public static Complex GaussKronrod(Func f, double intervalBegin, double intervalEnd, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15) + { + return GaussKronrodRule.ContourIntegrate(f, intervalBegin, intervalEnd, out _, out _, targetRelativeError: targetRelativeError, maximumDepth: maximumDepth, order: 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. + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts. + /// Where the interval stops. + /// The difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation + /// The L1 norm of the result, if there is a significant difference between this and the returned value, then the result is likely to be ill-conditioned. + /// 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, 21, 31, 41, 51 and 61 points + /// Approximation of the finite integral in the given interval. + public static Complex GaussKronrod(Func f, double intervalBegin, double intervalEnd, out double error, out double L1Norm, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15) + { + return GaussKronrodRule.ContourIntegrate(f, intervalBegin, intervalEnd, out error, out L1Norm, targetRelativeError: targetRelativeError, maximumDepth: maximumDepth, order: order); + } } } 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/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs new file mode 100644 index 00000000..467dc076 --- /dev/null +++ b/src/Numerics/Integration/GaussKronrodRule.cs @@ -0,0 +1,428 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2019 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// This file uses code from the Boost Project. +// Copyright John Maddock 2017. +// Copyright Nick Thompson 2017. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// https://github.com/boostorg/math/blob/develop/include/boost/math/quadrature/gauss_kronrod.hpp + +using MathNet.Numerics.Integration.GaussRule; +using System; +using System.Numerics; + +namespace MathNet.Numerics.Integration +{ + public class GaussKronrodRule + { + private readonly GaussPointPair gaussKronrodPoint; + + /// + /// Getter for the order. + /// + public int Order + { + get + { + return gaussKronrodPoint.Order; + } + } + + /// + /// Getter that returns a clone of the array containing the Kronrod abscissas. + /// + public double[] KronrodAbscissas + { + get + { + return gaussKronrodPoint.Abscissas.Clone() as double[]; + } + } + + /// + /// Getter that returns a clone of the array containing the Kronrod weights. + /// + public double[] KronrodWeights + { + get + { + return gaussKronrodPoint.Weights.Clone() as double[]; + } + } + + /// + /// Getter that returns a clone of the array containing the Gauss weights. + /// + public double[] GaussWeights + { + get + { + return gaussKronrodPoint.SecondWeights.Clone() as double[]; + } + } + + public GaussKronrodRule(int order) + { + gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); + } + + /// + /// Performs adaptive Gauss-Kronrod quadrature on function f over the range (a,b) + /// + /// The analytic smooth function to integrate + /// Where the interval starts + /// Where the interval stops + /// The difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation + /// The L1 norm of the result, if there is a significant difference between this and the returned value, then the result is likely to be ill-conditioned. + /// The maximum relative error in the result + /// The maximum number of interval splittings permitted before stopping + /// The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points + public static double Integrate(Func f, double intervalBegin, double intervalEnd, out double error, out double L1Norm, double targetRelativeError = 1E-10, int maximumDepth = 15, int order = 15) + { + // 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 (f == null) + { + throw new ArgumentNullException(nameof(f)); + } + + if (intervalBegin > intervalEnd) + { + return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); + } + + GaussPointPair gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(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 ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue)) + { + Func u = (t) => + { + return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); + }; + return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); + } + // [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 (intervalEnd > double.MaxValue) + { + Func u = (s) => + { + return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); + }; + return recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); + } + // (-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 (intervalBegin < double.MinValue) + { + Func u = (s) => + { + return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); + }; + return recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); + } + // [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 recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm, gaussKronrodPoint); + } + } + + /// + /// Performs adaptive Gauss-Kronrod quadrature on function f over the range (a,b) + /// + /// The analytic smooth complex function to integrate, defined on the real axis. + /// Where the interval starts + /// Where the interval stops + /// The difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation + /// The L1 norm of the result, if there is a significant difference between this and the returned value, then the result is likely to be ill-conditioned. + /// The maximum relative error in the result + /// The maximum number of interval splittings permitted before stopping + /// The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points + /// + public static Complex ContourIntegrate(Func f, double intervalBegin, double intervalEnd, out double error, out double L1Norm, double targetRelativeError = 1E-10, int maximumDepth = 15, int order = 15) + { + // 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 (f == null) + { + throw new ArgumentNullException(nameof(f)); + } + + if (intervalBegin > intervalEnd) + { + return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); + } + + GaussPointPair gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(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 ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue)) + { + Func u = (t) => + { + return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); + }; + return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); + } + // [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 (intervalEnd > double.MaxValue) + { + Func u = (s) => + { + return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); + }; + return contour_recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); + } + // (-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 (intervalBegin < double.MinValue) + { + Func u = (s) => + { + return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); + }; + return contour_recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); + } + // [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 contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm, gaussKronrodPoint); + } + } + + private static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPointPair gaussKronrodPoint) + { + int gauss_start = 2; + int kronrod_start = 1; + int gauss_order = (gaussKronrodPoint.Order - 1) / 2; + + double kronrod_result = 0d; + double gauss_result = 0d; + double fp, fm; + + var KAbscissa = gaussKronrodPoint.Abscissas; + var KWeights = gaussKronrodPoint.Weights; + var GWeights = gaussKronrodPoint.SecondWeights; + + if ((gauss_order & 1) == 1) + { + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_result += fp * GWeights[0]; + } + else + { + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_start = 1; + kronrod_start = 2; + } + double L1 = Math.Abs(kronrod_result); + + for (int i = gauss_start; i < KAbscissa.Length; i += 2) + { + fp = f(KAbscissa[i]); + fm = f(-KAbscissa[i]); + kronrod_result += (fp + fm) * KWeights[i]; + L1 += (Math.Abs(fp) + Math.Abs(fm)) * KWeights[i]; + gauss_result += (fp + fm) * GWeights[i / 2]; + } + for (int i = kronrod_start; i < KAbscissa.Length; i += 2) + { + fp = f(KAbscissa[i]); + fm = f(-KAbscissa[i]); + kronrod_result += (fp + fm) * KWeights[i]; + L1 += (Math.Abs(fp) + Math.Abs(fm)) * KWeights[i]; + } + pL1 = L1; + error = Math.Max(Math.Abs(kronrod_result - gauss_result), Math.Abs(kronrod_result * Precision.MachineEpsilon * 2d)); + return kronrod_result; + } + + private static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPointPair gaussKronrodPoint) + { + int gauss_start = 2; + int kronrod_start = 1; + int gauss_order = (gaussKronrodPoint.Order - 1) / 2; + + Complex kronrod_result = new Complex(); + Complex gauss_result = new Complex(); + Complex fp, fm; + + var KAbscissa = gaussKronrodPoint.Abscissas; + var KWeights = gaussKronrodPoint.Weights; + var GWeights = gaussKronrodPoint.SecondWeights; + + if (gauss_order.IsOdd()) + { + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_result += fp * GWeights[0]; + } + else + { + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_start = 1; + kronrod_start = 2; + } + double L1 = Complex.Abs(kronrod_result); + + for (int i = gauss_start; i < KAbscissa.Length; i += 2) + { + fp = f(KAbscissa[i]); + fm = f(-KAbscissa[i]); + kronrod_result += (fp + fm) * KWeights[i]; + L1 += (Complex.Abs(fp) + Complex.Abs(fm)) * KWeights[i]; + gauss_result += (fp + fm) * GWeights[i / 2]; + } + for (int i = kronrod_start; i < KAbscissa.Length; i += 2) + { + fp = f(KAbscissa[i]); + fm = f(-KAbscissa[i]); + kronrod_result += (fp + fm) * KWeights[i]; + L1 += (Complex.Abs(fp) + Complex.Abs(fm)) * KWeights[i]; + } + pL1 = L1; + error = Math.Max(Complex.Abs(kronrod_result - gauss_result), Complex.Abs(kronrod_result * Precision.MachineEpsilon * 2d)); + return kronrod_result; + } + + private static double recursive_adaptive_integrate(Func f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1, GaussPointPair gaussKronrodPoint) + { + double error_local; + double mean = (b + a) / 2; + double scale = (b - a) / 2; + + var r1 = integrate_non_adaptive_m1_1((x) => f(scale * x + mean), out error_local, out L1, gaussKronrodPoint); + var estimate = scale * r1; + + var tmp = estimate * rel_tol; + var abs_tol1 = Math.Abs(tmp); + if (abs_tol == 0) + { + abs_tol = abs_tol1; + } + + if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local)) + { + double mid = (a + b) / 2d; + double L1_local; + estimate = recursive_adaptive_integrate(f, a, mid, max_levels - 1, rel_tol, abs_tol / 2, out error, out L1, gaussKronrodPoint); + estimate += recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local, gaussKronrodPoint); + error += error_local; + L1 += L1_local; + return estimate; + } + L1 *= scale; + error = error_local; + return estimate; + } + + private static Complex contour_recursive_adaptive_integrate(Func f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1, GaussPointPair gaussKronrodPoint) + { + double error_local; + double mean = (b + a) / 2; + double scale = (b - a) / 2; + + var r1 = contour_integrate_non_adaptive_m1_1((x) => f(scale * x + mean), out error_local, out L1, gaussKronrodPoint); + var estimate = scale * r1; + + var tmp = estimate * rel_tol; + var abs_tol1 = Complex.Abs(tmp); + if (abs_tol == 0) + { + abs_tol = abs_tol1; + } + + if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local)) + { + double mid = (a + b) / 2d; + double L1_local; + estimate = contour_recursive_adaptive_integrate(f, a, mid, max_levels - 1, rel_tol, abs_tol / 2, out error, out L1, gaussKronrodPoint); + estimate += contour_recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local, gaussKronrodPoint); + error += error_local; + L1 += L1_local; + return estimate; + } + L1 *= scale; + error = error_local; + return estimate; + } + } +} 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]. /// diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs new file mode 100644 index 00000000..e1f4e93a --- /dev/null +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs @@ -0,0 +1,618 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Numerics; + +namespace MathNet.Numerics.Integration.GaussRule +{ + /// + /// Contains a method to compute the Gauss-Kronrod abscissas/weights and precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61. + /// + internal static partial class GaussKronrodPoint + { + /// + /// Precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61. + /// + internal static readonly Dictionary PreComputed = new Dictionary + { + { 15, new GaussPointPair(15, + new[] // 15-point Gauss-Kronrod Abscissa + { + 0.00000000000000000e+00, + 2.07784955007898468e-01, + 4.05845151377397167e-01, + 5.86087235467691130e-01, + 7.41531185599394440e-01, + 8.64864423359769073e-01, + 9.49107912342758525e-01, + 9.91455371120812639e-01, + }, + new[] // 15-point Gauss-Kronrod Weights + { + 2.09482141084727828e-01, + 2.04432940075298892e-01, + 1.90350578064785410e-01, + 1.69004726639267903e-01, + 1.40653259715525919e-01, + 1.04790010322250184e-01, + 6.30920926299785533e-02, + 2.29353220105292250e-02, + }, 7, + new[] // 7-point Gauss Weights + { + 4.17959183673469388e-01, + 3.81830050505118945e-01, + 2.79705391489276668e-01, + 1.29484966168869693e-01, + }) + }, + { 21, new GaussPointPair(21, + new[] // 21-point Gauss-Kronrod Abscissa + { + 0.00000000000000000e+00, + 1.48874338981631211e-01, + 2.94392862701460198e-01, + 4.33395394129247191e-01, + 5.62757134668604683e-01, + 6.79409568299024406e-01, + 7.80817726586416897e-01, + 8.65063366688984511e-01, + 9.30157491355708226e-01, + 9.73906528517171720e-01, + 9.95657163025808081e-01, + }, + new[] // 21-point Gauss-Kronrod Weights + { + 1.49445554002916906e-01, + 1.47739104901338491e-01, + 1.42775938577060081e-01, + 1.34709217311473326e-01, + 1.23491976262065851e-01, + 1.09387158802297642e-01, + 9.31254545836976055e-02, + 7.50396748109199528e-02, + 5.47558965743519960e-02, + 3.25581623079647275e-02, + 1.16946388673718743e-02, + }, 10, + new[] // 10-point Gauss Weights + { + 2.95524224714752870e-01, + 2.69266719309996355e-01, + 2.19086362515982044e-01, + 1.49451349150580593e-01, + 6.66713443086881376e-02, + }) + }, + { 31, new GaussPointPair(31, + new[] // 31-point Gauss-Kronrod Abscissa + { + 0.00000000000000000e+00, + 1.01142066918717499e-01, + 2.01194093997434522e-01, + 2.99180007153168812e-01, + 3.94151347077563370e-01, + 4.85081863640239681e-01, + 5.70972172608538848e-01, + 6.50996741297416971e-01, + 7.24417731360170047e-01, + 7.90418501442465933e-01, + 8.48206583410427216e-01, + 8.97264532344081901e-01, + 9.37273392400705904e-01, + 9.67739075679139134e-01, + 9.87992518020485428e-01, + 9.98002298693397060e-01, + }, + new[] // 31-point Gauss-Kronrod Weights + { + 1.01330007014791549e-01, + 1.00769845523875595e-01, + 9.91735987217919593e-02, + 9.66427269836236785e-02, + 9.31265981708253212e-02, + 8.85644430562117706e-02, + 8.30805028231330210e-02, + 7.68496807577203789e-02, + 6.98541213187282587e-02, + 6.20095678006706403e-02, + 5.34815246909280873e-02, + 4.45897513247648766e-02, + 3.53463607913758462e-02, + 2.54608473267153202e-02, + 1.50079473293161225e-02, + 5.37747987292334899e-03, + }, 15, + new[] // 15-point Gauss Weights + { + 2.02578241925561273e-01, + 1.98431485327111576e-01, + 1.86161000015562211e-01, + 1.66269205816993934e-01, + 1.39570677926154314e-01, + 1.07159220467171935e-01, + 7.03660474881081247e-02, + 3.07532419961172684e-02, + }) + }, + { 41, new GaussPointPair(41, + new[] // 41-point Gauss-Kronrod Abscissa + { + 0.00000000000000000e+00, + 7.65265211334973338e-02, + 1.52605465240922676e-01, + 2.27785851141645078e-01, + 3.01627868114913004e-01, + 3.73706088715419561e-01, + 4.43593175238725103e-01, + 5.10867001950827098e-01, + 5.75140446819710315e-01, + 6.36053680726515025e-01, + 6.93237656334751385e-01, + 7.46331906460150793e-01, + 7.95041428837551198e-01, + 8.39116971822218823e-01, + 8.78276811252281976e-01, + 9.12234428251325906e-01, + 9.40822633831754754e-01, + 9.63971927277913791e-01, + 9.81507877450250259e-01, + 9.93128599185094925e-01, + 9.98859031588277664e-01, + }, + new[] // 41-point Gauss-Kronrod Weights + { + 7.66007119179996564e-02, + 7.63778676720807367e-02, + 7.57044976845566747e-02, + 7.45828754004991890e-02, + 7.30306903327866675e-02, + 7.10544235534440683e-02, + 6.86486729285216193e-02, + 6.58345971336184221e-02, + 6.26532375547811680e-02, + 5.91114008806395724e-02, + 5.51951053482859947e-02, + 5.09445739237286919e-02, + 4.64348218674976747e-02, + 4.16688733279736863e-02, + 3.66001697582007980e-02, + 3.12873067770327990e-02, + 2.58821336049511588e-02, + 2.03883734612665236e-02, + 1.46261692569712530e-02, + 8.60026985564294220e-03, + 3.07358371852053150e-03, + }, 20, + new[] // 20-point Gauss Weights + { + 1.52753387130725851e-01, + 1.49172986472603747e-01, + 1.42096109318382051e-01, + 1.31688638449176627e-01, + 1.18194531961518417e-01, + 1.01930119817240435e-01, + 8.32767415767047487e-02, + 6.26720483341090636e-02, + 4.06014298003869413e-02, + 1.76140071391521183e-02, + }) + }, + { 51, new GaussPointPair(51, + new[] // 51-point Gauss-Kronrod Abscissa + { + 0.00000000000000000e+00, + 6.15444830056850789e-02, + 1.22864692610710396e-01, + 1.83718939421048892e-01, + 2.43866883720988432e-01, + 3.03089538931107830e-01, + 3.61172305809387838e-01, + 4.17885382193037749e-01, + 4.73002731445714961e-01, + 5.26325284334719183e-01, + 5.77662930241222968e-01, + 6.26810099010317413e-01, + 6.73566368473468364e-01, + 7.17766406813084388e-01, + 7.59259263037357631e-01, + 7.97873797998500059e-01, + 8.33442628760834001e-01, + 8.65847065293275595e-01, + 8.94991997878275369e-01, + 9.20747115281701562e-01, + 9.42974571228974339e-01, + 9.61614986425842512e-01, + 9.76663921459517511e-01, + 9.88035794534077248e-01, + 9.95556969790498098e-01, + 9.99262104992609834e-01, + }, + new[] // 51-point Gauss-Kronrod Weights + { + 6.15808180678329351e-02, + 6.14711898714253167e-02, + 6.11285097170530483e-02, + 6.05394553760458629e-02, + 5.97203403241740600e-02, + 5.86896800223942080e-02, + 5.74371163615678329e-02, + 5.59508112204123173e-02, + 5.42511298885454901e-02, + 5.23628858064074759e-02, + 5.02776790807156720e-02, + 4.79825371388367139e-02, + 4.55029130499217889e-02, + 4.28728450201700495e-02, + 4.00838255040323821e-02, + 3.71162714834155436e-02, + 3.40021302743293378e-02, + 3.07923001673874889e-02, + 2.74753175878517378e-02, + 2.40099456069532162e-02, + 2.04353711458828355e-02, + 1.68478177091282982e-02, + 1.32362291955716748e-02, + 9.47397338617415161e-03, + 5.56193213535671376e-03, + 1.98738389233031593e-03, + }, 25, + new[] // 25-point Gauss Weights + { + 1.23176053726715451e-01, + 1.22242442990310042e-01, + 1.19455763535784772e-01, + 1.14858259145711648e-01, + 1.08519624474263653e-01, + 1.00535949067050644e-01, + 9.10282619829636498e-02, + 8.01407003350010180e-02, + 6.80383338123569172e-02, + 5.49046959758351919e-02, + 4.09391567013063127e-02, + 2.63549866150321373e-02, + 1.13937985010262879e-02, + }) + }, + { 61, new GaussPointPair(61, + new[] // 61-point Gauss-Kronrod Abscissa + { + 0.00000000000000000e+00, + 5.14718425553176958e-02, + 1.02806937966737030e-01, + 1.53869913608583547e-01, + 2.04525116682309891e-01, + 2.54636926167889846e-01, + 3.04073202273625077e-01, + 3.52704725530878113e-01, + 4.00401254830394393e-01, + 4.47033769538089177e-01, + 4.92480467861778575e-01, + 5.36624148142019899e-01, + 5.79345235826361692e-01, + 6.20526182989242861e-01, + 6.60061064126626961e-01, + 6.97850494793315797e-01, + 7.33790062453226805e-01, + 7.67777432104826195e-01, + 7.99727835821839083e-01, + 8.29565762382768397e-01, + 8.57205233546061099e-01, + 8.82560535792052682e-01, + 9.05573307699907799e-01, + 9.26200047429274326e-01, + 9.44374444748559979e-01, + 9.60021864968307512e-01, + 9.73116322501126268e-01, + 9.83668123279747210e-01, + 9.91630996870404595e-01, + 9.96893484074649540e-01, + 9.99484410050490638e-01, + }, + new[] // 61-point Gauss-Kronrod Weights + { + 5.14947294294515676e-02, + 5.14261285374590259e-02, + 5.12215478492587722e-02, + 5.08817958987496065e-02, + 5.04059214027823468e-02, + 4.97956834270742064e-02, + 4.90554345550297789e-02, + 4.81858617570871291e-02, + 4.71855465692991539e-02, + 4.60592382710069881e-02, + 4.48148001331626632e-02, + 4.34525397013560693e-02, + 4.19698102151642461e-02, + 4.03745389515359591e-02, + 3.86789456247275930e-02, + 3.68823646518212292e-02, + 3.49793380280600241e-02, + 3.29814470574837260e-02, + 3.09072575623877625e-02, + 2.87540487650412928e-02, + 2.65099548823331016e-02, + 2.41911620780806014e-02, + 2.18280358216091923e-02, + 1.94141411939423812e-02, + 1.69208891890532726e-02, + 1.43697295070458048e-02, + 1.18230152534963417e-02, + 9.27327965951776343e-03, + 6.63070391593129217e-03, + 3.89046112709988405e-03, + 1.38901369867700762e-03, + }, 30, + new[] // 30-point Gauss Weights + { + 1.02852652893558840e-01, + 1.01762389748405505e-01, + 9.95934205867952671e-02, + 9.63687371746442596e-02, + 9.21225222377861287e-02, + 8.68997872010829798e-02, + 8.07558952294202154e-02, + 7.37559747377052063e-02, + 6.59742298821804951e-02, + 5.74931562176190665e-02, + 4.84026728305940529e-02, + 3.87991925696270496e-02, + 2.87847078833233693e-02, + 1.84664683110909591e-02, + 7.96819249616660562e-03, + }) + }, + }; + } + + /// + /// Contains a method to compute the Gauss-Kronrod abscissas/weights. + /// + internal static partial class GaussKronrodPoint + { + /// + /// Computes the Gauss-Kronrod abscissas/weights and Gauss weights. + /// + /// Defines an Nth order Gauss-Kronrod rule. The order also defines the number of abscissas and weights for the rule. + /// Required precision to compute the abscissas/weights. + /// Object containing the non-negative abscissas/weights, order. + internal static GaussPointPair Generate(int order, double eps) + { + int gaussOrder = (order - 1) / 2; + int gaussStart = gaussOrder.IsOdd() ? 0 : 1; + int kronrodStart = gaussOrder.IsOdd() ? 1 : 0; + + var gaussPoint = GaussLegendrePointFactory.GetGaussPoint(gaussOrder); + var gaussAbscissas = gaussPoint.Abscissas; + var gaussWeights = gaussPoint.Weights; + + // Calculate Kronrod polynomial in terms of Legendre polynomials + // K(x) = c0*P(0, x) + c1*P(1, x) + ... + + var c = StieltjesP(gaussOrder + 1); + + // Calculate Abscissas for Kronrod polynomial + + int r = gaussOrder.IsOdd() ? (gaussOrder - 1) / 2 + 1 : gaussOrder / 2 + 1; + var kronrodAbscissas = new double[r]; + + for (int k = 1; k <= gaussOrder + 1; k = k + 2) + { + var x0 = (1.0 - (1.0 - 1.0 / gaussOrder) / (8 * gaussOrder * gaussOrder)) * Math.Cos((k - 0.5) * Math.PI / (2.0 * gaussOrder + 1.0)); + var dx = 0d; + var j = 1; // iterations + + // Newton iterations + do + { + var E = LegendreSeries(c, x0); + dx = E.Item1 / E.Item2; + x0 = x0 - dx; + j++; + } + while (Math.Abs(dx) > eps && j < 100); + + if (Math.Abs(x0) < Precision.MachineEpsilon) x0 = 0.0; + + kronrodAbscissas[(k - 1) / 2] = x0; + } + + // Concatenate two abscissas + + var abscissas = new double[gaussAbscissas.Length + kronrodAbscissas.Length]; + gaussAbscissas.CopyTo(abscissas, 0); + kronrodAbscissas.CopyTo(abscissas, gaussAbscissas.Length); + abscissas = abscissas.OrderBy(v => v).ToArray(); + + // Calculate weights for abscissas + + var weights = new double[gaussAbscissas.Length + kronrodAbscissas.Length]; + for (int i = gaussStart; i < abscissas.Length; i += 2) + { + var x = abscissas[i]; + + var E = LegendreSeries(c, x); + var L = LegendreP(gaussOrder, x); + + var p = L.Item2; + var w2 = 2.0 / ((1.0 - x * x) * p * p); // Gauss weight + weights[i] = w2 + 2.0 / ((gaussOrder + 1.0) * p * E.Item1); + } + for (int i = kronrodStart; i < abscissas.Length; i += 2) + { + var x = abscissas[i]; + + var E = LegendreSeries(c, x); + var L = LegendreP(gaussOrder, x); + + weights[i] = 2.0 / ((gaussOrder + 1.0) * L.Item1 * E.Item2); + } + + return new GaussPointPair(order, abscissas, weights, gaussOrder, gaussWeights); + } + + /// + /// Returns coefficients of a Stieltjes polynomial in terms of Legendre polynomials. + /// + internal static double[] StieltjesP(int order) + { + // Reference: + // 1. Patterson, Thomas NL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856. + // 2. Piessens, Robert, and Maria Branders. "A note on the optimal addition of abscissas to quadrature formulas of Gauss and Lobatto type." Mathematics of Computation (1974): 135-139. + // 3. Legendre-Stieltjes Polynomials, Boost.Math + // + // Here, we are using Patterson algorithm, expanding the Stieltjes polynomial in terms of Legendre polynomials. + // + // Kronrod Polynomial K[n + 1, x] is expanded in terms of Legendre Polynomial P[n, x]. + // + // K[n + 1, x] = sum_(n=1)^r a[i] P[2 * i - 1 - q, x] + // + // where P[n, x] is the Legendre polynomial of degree n, + // [x] denotes the integer part of x, + // q = n - 2[n/2] + // r = [(n + 3)/2] + // + // The added n + 1 Kronrod abscissae is the roots of the Kronrod polynomial. + + if (order == 1) // P(1, x) + return new double[] { 0, 1 }; + else if (order == 2) // -2/5 * P(0, x) + P(2, x) + return new double[] { -0.4, 0, 1 }; + else if (order == 3) // -9/14 * P(1, x) + P(3, x) + return new double[] { 0, -0.642857142857142857142857142857, 0, 1 }; + else if (order == 4) // 14/891 * P(0, x) - 20/27 * P(2, x) + P(4, x) + return new double[] { 0.0157126823793490460157126823793, 0, -0.740740740740740740740740740741, 0, 1 }; + else if (order == 5) // 135/12584 * P(1, x) - 35/44 * P(3, x) + P(5, x) + return new double[] { 0, 0.0107279084551811824539097266370, 0, -0.795454545454545454545454545455, 0, 1 }; + + int n = order - 1; + int q = n.IsOdd() ? 1 : 0; + int r = n.IsOdd() ? (n - 1) / 2 + 2 : n / 2 + 1; + + double[] a = new double[r + 1]; + + // Calculate a[i] for i = 1, ..., r + // + // a[r] = 1; + // a[r - 1] = -a[r] * S[r, 1] / S[r - 1, 1]; + // a[r - 2] = -a[r] * S[r, 2] / S[r - 2, 2] - a[r - 1] * S[r - 1, 2] / S[r - 2, 2]; + // ... + // a[1] = -a[r] * S[r, r - 1] / S[1, r - 1] - a[r - 1] * S[r - 1, r - 1] / S[1, r - 1] - ... - a[2] * S[2, r - 1] / S[1, r - 1]; + // + // S[i, k] / S[r - k, k] = S[i - 1, k] / S[r - k, k] + // * ((n - q + 2 * (i + k - 1)) * (n + q + 2 * (k - i + 1)) * (n - 1 - q + 2 * (i - k)) * (2 * (k + i - 1) - 1 - q - n)) + // / ((n - q + 2 * (i - k)) * (2 * (k + i - 1) - q - n) * (n + 1 + q + 2 * (k - i)) * (n - 1 - q + 2 * (i + k))); + + a[r] = 1.0; + for (int k = 1; k < r; k++) + { + double ratio = 1.0; + a[r - k] = 0.0; + for (int i = r + 1 - k; i <= r; i++) + { + double numerator = (n - q + 2 * (i + k - 1)) * (n + q + 2 * (k - i + 1)) * (n - 1 - q + 2 * (i - k)) * (2 * (k + i - 1) - 1 - q - n); + double denominator = (n - q + 2 * (i - k)) * (2 * (k + i - 1) - q - n) * (n + 1 + q + 2 * (k - i)) * (n - 1 - q + 2 * (i + k)); + ratio = ratio * numerator / denominator; + a[r - k] -= a[i] * ratio; + } + } + + // K = sum c[k] P[k, x] + + double[] c = new double[2 * r - q]; + for (int i = 1; i < a.Length; i++) + { + c[2 * i - 1 - q] = a[i]; + } + + return c; + } + + /// + /// Return value and derivative of a Legendre series at given points. + /// + internal static Tuple LegendreSeries(double[] a, double x) + { + // S = a[0]*P[0, x] + ... + a[k]*P[k, x] + ... + a[n]*P[n, x] + // where P[k, x] is the Legendre polynomial of order k + // + // According to the Clenshaw algorithm, S can be written by + // S = a[0] + x*b[1, x] - 1/2 * b[2,x] + // + // b[n + 1, x] = 0 + // b[n + 2, x] = 0 + // b[k, x] = a[k] + (2k + 1)/(k + 1)*x*b[k + 1, x] - (k + 1)/(k + 2)*b[k + 2, x] + // + // Derivative of S is given by + // S' = b[1, x] + x*b'[1, x] - 1/2 * b'[2,x] + // + // b'[k, x] = (2k + 1)/(k + 1)*b[k + 1, x] + (2k + 1)/(k + 1)*x*b'[k + 1, x] - (k + 1)/(k + 2)*b'[k + 2, x] + + if (a.Length == 1) + return new Tuple(a[0], 0); + if (a.Length == 2) + return new Tuple(a[0] + a[1] * x, a[1]); + + double b0 = 0.0, b1 = 0.0, b2 = 0.0; + double p0 = 0.0, p1 = 0.0, p2 = 0.0; + + for (int k = a.Length - 1; k >= 1; k--) + { + b0 = a[k] + (2.0 * k + 1.0) / (k + 1.0) * x * b1 - (k + 1.0) / (k + 2.0) * b2; + p0 = (2.0 * k + 1.0) / (k + 1.0) * (b1 + x * p1) - (k + 1.0) / (k + 2.0) * p2; + + b2 = b1; + b1 = b0; + p2 = p1; + p1 = p0; + } + + var value = a[0] + b1 * x - 0.5 * b2; + var derivative = b1 + p1 * x - 0.5 * p2; + + return new Tuple( value, derivative ); + } + + /// + /// Return value and derivative of a Legendre polynomial of order at given points. + /// + internal static Tuple LegendreP(int order, double x) + { + // The Legendre polynomial, P[n, x], is defined by the recurrence relation: + // + // P[0, x] = 1 + // P[1, x] = x + // (n + 1) * P[n + 1, x] = (2 * n + 1) * x * P[n, x] - n * P[n - 1, x] + // + // The derivative of the Legendre polynomial, P'[n, x] is given by + // P'[0, x] = 0 + // P'[1, x] = 1 + // (n + 1) * P'[n + 1, x] = (2 * n + 1) * P[n, x] + (2 * n + 1) * x * P'[n, x] - n * P'[n - 1, x] + // = (2 * n + 1) * (P[n, x] + x * P'[n, x]) - n * P'[n - 1, x] + + if (order == 0) + return new Tuple(1.0, 0.0); + if (order == 1) + return new Tuple(x, 1.0); + + double b0 = 0.0, b1 = 1.0, b2 = 0.0; + double p0 = 0.0, p1 = 0.0, p2 = 0.0; + + for (int k = 1; k <= order; k++) + { + b0 = (2.0 * k - 1.0) / k * x * b1 - (k - 1.0) / k * b2; // L(k, x) + p0 = (2.0 * k - 1.0) / k * (b1 + x * p1) - (k - 1.0) / k * p2; // L'(k, x) + + b2 = b1; + b1 = b0; + p2 = p1; + p1 = p0; + } + + var value = b0; + var derivative = p0; + + return new Tuple(value, derivative); + } + } +} diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs new file mode 100644 index 00000000..372fb577 --- /dev/null +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs @@ -0,0 +1,34 @@ +using System; + +namespace MathNet.Numerics.Integration.GaussRule +{ + /// + /// Creates a Gauss-Kronrod point. + /// + internal static class GaussKronrodPointFactory + { + [ThreadStatic] + private static GaussPointPair gaussKronrodPoint; + + /// + /// Getter for the GaussKronrodPoint. + /// + /// Defines an Nth order Gauss-Kronrod rule. Precomputed Gauss-Kronrod abscissas/weights for orders 15, 21, 31, 41, 51, 61 are used, otherwise they're calculated on the fly. + /// Object containing the non-negative abscissas/weights, and order. + public static GaussPointPair GetGaussPoint(int order) + { + // Try to get the GaussKronrodPoint from the cached static field. + bool gaussKronrodPointIsCached = gaussKronrodPoint != null && gaussKronrodPoint.Order == order; + if (!gaussKronrodPointIsCached) + { + // Try to find the GaussKronrodPoint in the precomputed dictionary. + if (!GaussKronrodPoint.PreComputed.TryGetValue(order, out gaussKronrodPoint)) + { + gaussKronrodPoint = GaussKronrodPoint.Generate(order, 1E-10); + } + } + + return gaussKronrodPoint; + } + } +} diff --git a/src/Numerics/Integration/GaussRule/GaussPointPair.cs b/src/Numerics/Integration/GaussRule/GaussPointPair.cs new file mode 100644 index 00000000..c870474d --- /dev/null +++ b/src/Numerics/Integration/GaussRule/GaussPointPair.cs @@ -0,0 +1,40 @@ +namespace MathNet.Numerics.Integration.GaussRule +{ + /// + /// Contains two GaussPoint. + /// + internal class GaussPointPair + { + internal int Order { get; private set; } + + internal double[] Abscissas { get; private set; } + + internal double[] Weights { get; private set; } + + internal int SecondOrder { get; private set; } + + internal double[] SecondAbscissas { get; private set; } + + internal double[] SecondWeights { get; private set; } + + internal double IntervalBegin { get; private set; } + + internal double IntervalEnd { get; private set; } + + internal GaussPointPair(double intervalBegin, double intervalEnd, int order, double[] abscissas, double[] weights, int secondOrder, double[] secondAbscissas, double[] secondWeights) + { + IntervalBegin = intervalBegin; + IntervalEnd = intervalEnd; + Order = order; + Abscissas = abscissas; + Weights = weights; + SecondOrder = secondOrder; + SecondAbscissas = secondAbscissas; + SecondWeights = secondWeights; + } + + internal GaussPointPair(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondWeights) + : this(-1, 1, order, abscissas, weights, secondOrder, null, secondWeights) + { } + } +} 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; + } + } } }