From 92ae8e358ac69e04af8b789832061a7f805333c4 Mon Sep 17 00:00:00 2001 From: diluculo Date: Tue, 10 Sep 2019 20:24:42 +0900 Subject: [PATCH 01/13] Add Integrate.OnOpenInterval --- .paket/Paket.Restore.targets | 77 ++++++++++++------- .../IntegrationTests/IntegrationTest.cs | 38 ++++++++- src/Numerics/Integrate.cs | 72 ++++++++++++++++- 3 files changed, 153 insertions(+), 34 deletions(-) diff --git a/.paket/Paket.Restore.targets b/.paket/Paket.Restore.targets index 818b4eca..952ad42f 100644 --- a/.paket/Paket.Restore.targets +++ b/.paket/Paket.Restore.targets @@ -5,6 +5,11 @@ $(MSBuildAllProjects);$(MSBuildThisFileFullPath) + + $(MSBuildVersion) + 15.0.0 + false + true true $(MSBuildThisFileDirectory) @@ -59,7 +64,7 @@ - + true true @@ -73,37 +78,47 @@ + + + - + true $(NoWarn);NU1603;NU1604;NU1605;NU1608 + false + true - - - /usr/bin/shasum "$(PaketRestoreCacheFile)" | /usr/bin/awk '{ print $1 }' - /usr/bin/shasum "$(PaketLockFilePath)" | /usr/bin/awk '{ print $1 }' + + + + + + + $([System.IO.File]::ReadAllText('$(PaketRestoreCacheFile)')) + + + + + + + $([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `": "`)[0].Replace(`"`, ``).Replace(` `, ``)) + $([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `": "`)[1].Replace(`"`, ``).Replace(` `, ``)) + + + + + %(PaketRestoreCachedKeyValue.Value) + %(PaketRestoreCachedKeyValue.Value) - - - - - - - - - - - - - - $([System.IO.File]::ReadAllText('$(PaketRestoreCacheFile)')) - $([System.IO.File]::ReadAllText('$(PaketLockFilePath)')) + + true - false + false true @@ -116,7 +131,10 @@ + + + @@ -126,7 +144,7 @@ - + $(PaketIntermediateOutputPath)\$(MSBuildProjectFile).paket.references.cached @@ -163,6 +181,7 @@ + @@ -224,8 +243,6 @@ false - $(MSBuildVersion) - 15.8.0 @@ -252,10 +269,12 @@ - <_NuspecFiles Include="$(AdjustedNuspecOutputPath)\*.$(PackageVersion).nuspec"/> + <_NuspecFiles Include="$(AdjustedNuspecOutputPath)\*.$(PackageVersion.Split(`+`)[0]).nuspec"/> - + + + @@ -281,7 +300,7 @@ DevelopmentDependency="$(DevelopmentDependency)" BuildOutputInPackage="@(_BuildOutputInPackage)" TargetPathsToSymbols="@(_TargetPathsToSymbols)" - SymbolPackageFormat="symbols.nupkg" + SymbolPackageFormat="$(SymbolPackageFormat)" TargetFrameworks="@(_TargetFrameworks)" AssemblyName="$(AssemblyName)" PackageOutputPath="$(PackageOutputAbsolutePath)" @@ -328,7 +347,7 @@ DevelopmentDependency="$(DevelopmentDependency)" BuildOutputInPackage="@(_BuildOutputInPackage)" TargetPathsToSymbols="@(_TargetPathsToSymbols)" - SymbolPackageFormat="symbols.nupkg" + SymbolPackageFormat="$(SymbolPackageFormat)" TargetFrameworks="@(_TargetFrameworks)" AssemblyName="$(AssemblyName)" PackageOutputPath="$(PackageOutputAbsolutePath)" diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index ad362cde..e5e55ca2 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -59,7 +59,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests { return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y)); } - + /// /// Test Function Start point. /// @@ -311,5 +311,41 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order); Assert.AreEqual(gaussLegendre.IntervalEnd, StopA); } + + // integral_(-oo)^(oo) exp(-x^2/2) dx = sqrt(2 ¥ð) + [TestCase(double.NegativeInfinity, double.PositiveInfinity, Constants.Sqrt2Pi)] + // integral_(-oo)^(0) exp(-x^2/2) dx = sqrt(¥ð/2) + [TestCase(double.NegativeInfinity, 0, Constants.SqrtPiOver2)] + // integral_(0)^(oo exp(-x^2/2) dx = sqrt(¥ð/2) + [TestCase(0, double.PositiveInfinity, Constants.SqrtPiOver2)] + // integral_(-1)^(1) exp(-x^2/2) dx = sqrt(2 ¥ð) erf(1/sqrt(2)) + [TestCase(-1, 1, 1.7112487837842976063)] + // integral_(1)^(0) exp(-x^2/2) dx = -sqrt(¥ð/2) erf(1/sqrt(2)) + [TestCase(1, 0, -0.85562439189214880317)] + public void TestGaussianIntegralBySubstitution(double a, double b, double expected) + { + Assert.AreEqual( + expected, + Integrate.OnOpenInterval((x) => Math.Exp(-x * x / 2), a, b), + 1e-10, + "Integral 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 TestSincIntegralBySubstitution(double a, double b, double expected, double factor) + { + Assert.AreEqual( + expected, + factor * Integrate.OnOpenInterval((x) => 1 / (1 + x * x), a, b), + 1e-10, + "Integral sin(pi*x)/(pi*x) from -oo to oo"); + } } } diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index 2d708fdf..461516e0 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -45,21 +45,85 @@ namespace MathNet.Numerics /// 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 double OnClosedInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError) + public static double OnClosedInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) { return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError); } /// - /// Approximation of the definite integral of an analytic smooth function on a closed interval. + /// 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. + /// The expected relative accuracy of the approximation. /// Approximation of the finite integral in the given interval. - public static double OnClosedInterval(Func f, double intervalBegin, double intervalEnd) + public static double OnOpenInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) { - return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, 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 -OnOpenInterval(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 OnClosedInterval(u, -1d, 1d, 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 OnClosedInterval(u, 0d, 1d, 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 OnClosedInterval(u, -1d, 0d, 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 OnClosedInterval(u, -1d, 1d, targetAbsoluteError); + } } /// From 067ec6fed3edf75ecf7987309e14051ffc332c3f Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 11 Sep 2019 15:02:10 +0900 Subject: [PATCH 02/13] Add complex version of double-exponential integral --- .../IntegrationTests/IntegrationTest.cs | 70 ++++++- src/Numerics/Integrate.cs | 143 ++++++++++--- .../DoubleExponentialTransformation.cs | 20 ++ .../Integration/NewtonCotesTrapeziumRule.cs | 193 ++++++++++++++++++ 4 files changed, 392 insertions(+), 34 deletions(-) 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; + } + } } } From cca30a3a0e4de716016986626f110c0e6d4aad4e Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 11 Sep 2019 18:36:48 +0900 Subject: [PATCH 03/13] Add Gauss-Kronrod integration rule --- .../IntegrationTests/IntegrationTest.cs | 137 ++- src/Numerics/Integrate.cs | 72 +- src/Numerics/Integration/GaussKronrodRule.cs | 803 ++++++++++++++++++ 3 files changed, 984 insertions(+), 28 deletions(-) create mode 100644 src/Numerics/Integration/GaussKronrodRule.cs diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index 1f600cfa..cc2e3b79 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -71,6 +71,16 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests return 1 / (1 + x * x); } + /// + /// Test Function: f(x,y) = log(x) + /// + /// First input value. + /// Function result. + private static double TargetFunctionD(double x) + { + return Math.Log(x); + } + /// /// Test Function Start point. /// @@ -101,6 +111,16 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// private const double StopC = double.PositiveInfinity; + /// + /// Test Function Start point. + /// + private const double StartD = 0; + + /// + /// Test Function Stop point. + /// + private const double StopD = 1; + /// /// Target area square. /// @@ -116,6 +136,11 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// private const double TargetAreaC = Constants.Pi; + /// + /// Target area. + /// + private const double TargetAreaD = -1; + /// /// Test Integrate facade for simple use cases. /// @@ -150,13 +175,52 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests TargetAreaC, Integrate.DoubleExponential(TargetFunctionC, StartC, StopC), 1e-5, - "Integral by substitution"); + "DoubleExponential"); Assert.AreEqual( TargetAreaC, Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10), 1e-10, - "Integral by substitution, Target 1e-10"); + "DoubleExponential, Target 1e-10"); + + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 15), + 1e-10, + "GaussKronrod, Target 1e-10, order 15"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 21), + 1e-10, + "GaussKronrod, Target 1e-10, order 21"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 31), + 1e-10, + "GaussKronrod, Target 1e-10, order 31"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 41), + 1e-10, + "GaussKronrod, Target 1e-10, order 41"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 51), + 1e-10, + "GaussKronrod, Target 1e-10, order 51"); + Assert.AreEqual( + TargetAreaD, + Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 61), + 1e-10, + "GaussKronrod, Target 1e-10, 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"); } /// @@ -351,22 +415,28 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests } // integral_(-oo)^(oo) exp(-x^2/2) dx = sqrt(2 ¥ð) - [TestCase(double.NegativeInfinity, double.PositiveInfinity, Constants.Sqrt2Pi)] // integral_(-oo)^(0) exp(-x^2/2) dx = sqrt(¥ð/2) - [TestCase(double.NegativeInfinity, 0, Constants.SqrtPiOver2)] // integral_(0)^(oo exp(-x^2/2) dx = sqrt(¥ð/2) - [TestCase(0, double.PositiveInfinity, Constants.SqrtPiOver2)] // integral_(-1)^(1) exp(-x^2/2) dx = sqrt(2 ¥ð) erf(1/sqrt(2)) - [TestCase(-1, 1, 1.7112487837842976063)] // 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 TestGaussianIntegralBySubstitution(double a, double b, double expected) + public void TestIntegralOfGaussian(double a, double b, double expected) { Assert.AreEqual( - expected, - Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b), - 1e-10, - "Integral e^(-x^2 /2) from {0} to {1}", a, b); + expected, + Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b), + 1e-10, + "DET Integral e^(-x^2 /2) from {0} to {1}", a, b); + + Assert.AreEqual( + expected, + Integrate.GaussKronrod((x) => Math.Exp(-x * x / 2), a, b), + 1e-10, + "GK Integral e^(-x^2 /2) from {0} to {1}", a, b); } // integral_(-oo)^(oo) sin(pi x) / (pi x) dx = 1 / pi integral_(-oo)^(oo) sin(x) / x dx @@ -377,37 +447,56 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests [TestCase(double.NegativeInfinity, double.PositiveInfinity, 1, Constants.InvPi)] [TestCase(0, double.PositiveInfinity, 1, Constants.TwoInvPi)] [TestCase(double.NegativeInfinity, 0, 1, Constants.TwoInvPi)] - public void TestSincIntegralBySubstitution(double a, double b, double expected, double factor) + 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, - "Integral sin(pi*x)/(pi*x) from -oo to oo"); + expected, + factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), + 1e-10, + "DET Integral sin(pi*x)/(pi*x) from -oo to oo"); + + Assert.AreEqual( + expected, + factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b), + 1e-10, + "GK Integral sin(pi*x)/(pi*x) from -oo to oo"); } // 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, double.PositiveInfinity, 2.2214414690791831235, -2.2214414690791831235)] + [TestCase(0, double.PositiveInfinity, 1.1107207345395915618, -1.1107207345395915618)] [TestCase(double.NegativeInfinity, 0, 1.1107207345395915618, -1.1107207345395915618)] - public void TestContourIntegralBySubstitution(double a, double b, double r, double i) + public void TestContourIntegral(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); + 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); + + Assert.AreEqual( + expected.Real, + actualDET.Real, + 1e-10, + "DET Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + + Assert.AreEqual( + expected.Imaginary, + actualDET.Imaginary, + 1e-10, + "DET Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Real, - actual.Real, + actualGK.Real, 1e-10, - "Integral e^(-x^2 /2) / (1 + j e^x) from {0} to {1}", a, b); + "GK Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Imaginary, - actual.Imaginary, + actualGK.Imaginary, 1e-10, - "Integral e^(-x^2 /2) / (1 + j e^x) from {0} to {1}", a, b); + "GK Integral 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 1a73d75a..cc7944b3 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -79,9 +79,9 @@ namespace MathNet.Numerics { 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. + /// 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. @@ -155,15 +155,47 @@ namespace MathNet.Numerics return DoubleExponentialTransformation.Integrate(u, -1, 1, targetAbsoluteError); } } + + /// + /// 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 over a real variable, of a complex-valued function. + /// 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 substitution. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. + /// 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. @@ -237,5 +269,37 @@ namespace MathNet.Numerics return DoubleExponentialTransformation.ContourIntegrate(u, -1, 1, targetAbsoluteError); } } + + /// + /// 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/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs new file mode 100644 index 00000000..82e1bf36 --- /dev/null +++ b/src/Numerics/Integration/GaussKronrodRule.cs @@ -0,0 +1,803 @@ +// +// 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 System; +using System.Numerics; + +namespace MathNet.Numerics.Integration +{ + public static class GaussKronrodRule + { + const double epsilon = 2.2204460492503131e-016; + + /// + /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points. + /// + static int Order = 15; + + static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1) + { + int gauss_start = 2; + int kronrod_start = 1; + int gauss_order = ((int)Order - 1) / 2; + + double kronrod_result = 0d; + double gauss_result = 0d; + double fp, fm; + + var KAbscissa = KronrodAbscissa(); + var KWeights = KronrodWeights(); + var GWeights = GaussWeights(); + + 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 * epsilon * 2d)); + return kronrod_result; + } + + static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1) + { + int gauss_start = 2; + int kronrod_start = 1; + int gauss_order = ((int)Order - 1) / 2; + + Complex kronrod_result = new Complex(); + Complex gauss_result = new Complex(); + Complex fp, fm; + + var KAbscissa = KronrodAbscissa(); + var KWeights = KronrodWeights(); + var GWeights = GaussWeights(); + + 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 * epsilon * 2d)); + return kronrod_result; + } + + 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) + { + 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); + 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); + estimate += recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local); + error += error_local; + L1 += L1_local; + return estimate; + } + L1 *= scale; + error = error_local; + return estimate; + } + + 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) + { + 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); + 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); + estimate += contour_recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local); + error += error_local; + L1 += L1_local; + return estimate; + } + L1 *= scale; + error = error_local; + return estimate; + } + + /// + /// 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)); + } + + Order = order; + + if (intervalBegin > intervalEnd) + { + return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, 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); + } + // [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); + } + // (-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); + } + // [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); + } + } + + /// + /// 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)); + } + + Order = order; + + if (intervalBegin > intervalEnd) + { + return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, 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); + } + // [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); + } + // (-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); + } + // [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); + } + } + + #region Pre-computed Abscissa and weights + + static double[] KronrodAbscissa() + { + switch (Order) + { + default: + case 15: + return PrecomputedKronrodAbscissas[0]; + case 21: + return PrecomputedKronrodAbscissas[1]; + case 31: + return PrecomputedKronrodAbscissas[2]; + case 41: + return PrecomputedKronrodAbscissas[3]; + case 51: + return PrecomputedKronrodAbscissas[4]; + case 61: + return PrecomputedKronrodAbscissas[5]; + } + } + + static double[] KronrodWeights() + { + switch (Order) + { + default: + case 15: + return PrecomputedKronrodWeights[0]; + case 21: + return PrecomputedKronrodWeights[1]; + case 31: + return PrecomputedKronrodWeights[2]; + case 41: + return PrecomputedKronrodWeights[3]; + case 51: + return PrecomputedKronrodWeights[4]; + case 61: + return PrecomputedKronrodWeights[5]; + } + } + + static double[] GaussWeights() + { + switch (Order) + { + default: + case 15: + return PrecomputedGaussWeights[0]; + case 21: + return PrecomputedGaussWeights[1]; + case 31: + return PrecomputedGaussWeights[2]; + case 41: + return PrecomputedGaussWeights[3]; + case 51: + return PrecomputedGaussWeights[4]; + case 61: + return PrecomputedGaussWeights[5]; + } + } + + /// + /// precomputed abscissa vector per order 15, 21, 31, 41, 51 and 61 + /// + static readonly double[][] PrecomputedKronrodAbscissas = + { + new[] // 15-point Gauss-Kronrod + { + 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[] // 21-point Gauss-Kronrod + { + 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[] // 31-point Gauss-Kronrod + { + 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[] // 41-point Gauss-Kronrod + { + 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[] // 51-point Gauss-Kronrod + { + 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[] // 61-point Gauss-Kronrod + { + 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, + } + }; + + /// + /// precomputed weight vector per order 15, 21, 31, 41, 51 and 61 + /// + static readonly double[][] PrecomputedKronrodWeights = + { + new[] // 15-point Gauss-Kronrod integration + { + 2.09482141084727828e-01, + 2.04432940075298892e-01, + 1.90350578064785410e-01, + 1.69004726639267903e-01, + 1.40653259715525919e-01, + 1.04790010322250184e-01, + 6.30920926299785533e-02, + 2.29353220105292250e-02, + }, + new[] // 21-point Gauss-Kronrod integration + { + 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, + }, + new[] // 31-point Gauss-Kronrod integration + { + 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, + }, + new[] // 41-point Gauss-Kronrod integration + { + 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, + }, + new[] // 51-point Gauss-Kronrod integration + { + 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, + }, + new[] // 61-point Gauss-Kronrod integration + { + 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, + }, + }; + + /// + /// precomputed Gauss weight vector per order 7, 10, 15, 20, 25 and 30 + /// + static readonly double[][] PrecomputedGaussWeights = + { + new [] // 7-point Gauss + { + 4.17959183673469388e-01, + 3.81830050505118945e-01, + 2.79705391489276668e-01, + 1.29484966168869693e-01, + }, + new[] // 10-point Gauss + { + 2.95524224714752870e-01, + 2.69266719309996355e-01, + 2.19086362515982044e-01, + 1.49451349150580593e-01, + 6.66713443086881376e-02, + }, + new[] // 15-point Gauss + { + 2.02578241925561273e-01, + 1.98431485327111576e-01, + 1.86161000015562211e-01, + 1.66269205816993934e-01, + 1.39570677926154314e-01, + 1.07159220467171935e-01, + 7.03660474881081247e-02, + 3.07532419961172684e-02, + }, + new[] // 20-point Gauss + { + 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, + }, + new[] // 25-point Gauss + { + 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, + }, + new[] // 30-point Gauss + { + 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, + } + }; + + #endregion Pre-computed Abscissa and weights + } +} From 9cb8bbea5f2d3d5b3bb43a8fce287a046041922e Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 11 Sep 2019 20:51:47 +0900 Subject: [PATCH 04/13] Add Gauss-Legendre to Integrate class --- .../IntegrationTests/IntegrationTest.cs | 41 ++++- src/Numerics/Integrate.cs | 156 +++++++++++++++++- src/Numerics/Integration/GaussLegendreRule.cs | 43 +++++ 3 files changed, 230 insertions(+), 10 deletions(-) diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index cc2e3b79..c8422144 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -430,13 +430,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests expected, Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b), 1e-10, - "DET Integral e^(-x^2 /2) from {0} to {1}", a, b); + "DET Integral of e^(-x^2 /2) from {0} to {1}", a, b); Assert.AreEqual( expected, Integrate.GaussKronrod((x) => Math.Exp(-x * x / 2), a, b), 1e-10, - "GK Integral e^(-x^2 /2) from {0} to {1}", a, b); + "GK Integral of e^(-x^2 /2) from {0} to {1}", a, b); + + Assert.AreEqual( + expected, + Integrate.GausLegendre((x) => Math.Exp(-x * x / 2), a, b, order: 128), + 1e-10, + "GL Integral of e^(-x^2 /2) from {0} to {1}", a, b); } // integral_(-oo)^(oo) sin(pi x) / (pi x) dx = 1 / pi integral_(-oo)^(oo) sin(x) / x dx @@ -453,13 +459,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests expected, factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), 1e-10, - "DET Integral sin(pi*x)/(pi*x) from -oo to oo"); + "DET Integral of sin(pi*x)/(pi*x) from -oo to oo"); Assert.AreEqual( expected, factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b), 1e-10, - "GK Integral sin(pi*x)/(pi*x) from -oo to oo"); + "GK Integral of sin(pi*x)/(pi*x) from -oo to oo"); + + Assert.AreEqual( + expected, + factor * Integrate.GausLegendre((x) => 1 / (1 + x * x), a, b, order: 128), + 1e-10, + "GL Integral of sin(pi*x)/(pi*x) from -oo to oo"); } // integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð @@ -473,30 +485,43 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests var expected = new Complex(r, i); var actualDET = ContourIntegrate.DoubleExponential((x) => 1 / new Complex(1, x * x), a, b); var actualGK = ContourIntegrate.GaussKronrod((x) => 1 / new Complex(1, x * x), a, b); + var actualGL = ContourIntegrate.GaussLegendre((x) => 1 / new Complex(1, x * x), a, b, order: 128); Assert.AreEqual( expected.Real, actualDET.Real, 1e-10, - "DET Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "DET Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Imaginary, actualDET.Imaginary, 1e-10, - "DET Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "DET Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Real, actualGK.Real, 1e-10, - "GK Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "GK Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); Assert.AreEqual( expected.Imaginary, actualGK.Imaginary, 1e-10, - "GK Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + "GK Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + + Assert.AreEqual( + expected.Real, + actualGL.Real, + 1e-10, + "GL Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); + + Assert.AreEqual( + expected.Imaginary, + actualGL.Imaginary, + 1e-10, + "GL Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b); } } } diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index cc7944b3..6057ed6a 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -156,6 +156,82 @@ namespace MathNet.Numerics } } + /// + /// Approximation of the definite integral of an analytic smooth function by Gauss-Legendre quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. + /// + /// The analytic smooth function to integrate. + /// Where the interval starts. + /// Where the interval stops. + /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. + /// Approximation of the finite integral in the given interval. + public static double GausLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) + { + // Reference: + // Formula used for variable subsitution from + // 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140. + // 2. quadgk.m, GNU Octave + + if (intervalBegin > intervalEnd) + { + return -GaussLegendreRule.Integrate(f, intervalEnd, intervalBegin, order); + } + + // (-oo, oo) => [-1, 1] + // + // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = t / (1 - t^2) + // g'(t) = (1 + t^2) / (1 - t^2)^2 + if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) + { + Func u = (t) => + { + return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); + }; + return GaussLegendreRule.Integrate(u, -1, 1, order); + } + // [a, oo) => [0, 1] + // + // integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt + // = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 - s) + // g'(s) = 1 / (1 - s)^2 + else if (double.IsInfinity(intervalEnd)) + { + Func u = (s) => + { + return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); + }; + return GaussLegendreRule.Integrate(u, 0, 1, order); + } + // (-oo, b] => [-1, 0] + // + // integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt + // = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 + s) + // g'(s) = 1 / (1 + s)^2 + else if (double.IsInfinity(intervalBegin)) + { + Func u = (s) => + { + return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); + }; + return GaussLegendreRule.Integrate(u, -1, 0, order); + } + // [a, b] => [-1, 1] + // + // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2 + // g'(t) = 3 / 4 * (b - a) * (1 - t^2) + else + { + Func u = (t) => + { + return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t); + }; + return GaussLegendreRule.Integrate(u, -1, 1, order); + } + } + /// /// Approximation of the definite integral of an analytic smooth function by Gauss-Kronrod quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. /// @@ -163,8 +239,8 @@ namespace MathNet.Numerics /// Where the interval starts. /// Where the interval stops. /// The expected relative accuracy of the approximation. - /// The maximum number of interval splittings permitted before stopping - /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points + /// The maximum number of interval splittings permitted before stopping. + /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points. /// Approximation of the finite integral in the given interval. public static double GaussKronrod(Func f, double intervalBegin, double intervalEnd, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15) { @@ -270,6 +346,82 @@ namespace MathNet.Numerics } } + /// + /// Approximation of the definite integral of an analytic smooth complex function by double-exponential quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts. + /// Where the interval stops. + /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. + /// Approximation of the finite integral in the given interval. + public static Complex GaussLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) + { + // Reference: + // Formula used for variable subsitution from + // 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140. + // 2. quadgk.m, GNU Octave + + if (intervalBegin > intervalEnd) + { + return -GaussLegendre(f, intervalEnd, intervalBegin, order); + } + + // (-oo, oo) => [-1, 1] + // + // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = t / (1 - t^2) + // g'(t) = (1 + t^2) / (1 - t^2)^2 + if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) + { + Func u = (t) => + { + return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t)); + }; + return GaussLegendreRule.ContourIntegrate(u, -1, 1, order); + } + // [a, oo) => [0, 1] + // + // integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt + // = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 - s) + // g'(s) = 1 / (1 - s)^2 + else if (double.IsInfinity(intervalEnd)) + { + Func u = (s) => + { + return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s)); + }; + return GaussLegendreRule.ContourIntegrate(u, 0, 1, order); + } + // (-oo, b] => [-1, 0] + // + // integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt + // = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds + // g(s) = s / (1 + s) + // g'(s) = 1 / (1 + s)^2 + else if (double.IsInfinity(intervalBegin)) + { + Func u = (s) => + { + return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s)); + }; + return GaussLegendreRule.ContourIntegrate(u, -1, 0, order); + } + // [a, b] => [-1, 1] + // + // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2 + // g'(t) = 3 / 4 * (b - a) * (1 - t^2) + else + { + Func u = (t) => + { + return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t); + }; + return GaussLegendreRule.ContourIntegrate(u, -1, 1, order); + } + } + /// /// Approximation of the definite integral of an analytic smooth function by Gauss-Kronrod quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity. /// diff --git a/src/Numerics/Integration/GaussLegendreRule.cs b/src/Numerics/Integration/GaussLegendreRule.cs index 60bed83f..fab03a4e 100644 --- a/src/Numerics/Integration/GaussLegendreRule.cs +++ b/src/Numerics/Integration/GaussLegendreRule.cs @@ -28,6 +28,7 @@ // using System; +using System.Numerics; using MathNet.Numerics.Integration.GaussRule; namespace MathNet.Numerics.Integration @@ -166,6 +167,48 @@ namespace MathNet.Numerics.Integration return a*sum; } + /// + /// Approximates a definite integral using an Nth order Gauss-Legendre rule. + /// + /// The analytic smooth complex function to integrate, defined on the real domain. + /// Where the interval starts, exclusive and finite. + /// Where the interval ends, exclusive and finite. + /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. + /// Approximation of the finite integral in the given interval. + public static Complex ContourIntegrate(Func f, double invervalBegin, double invervalEnd, int order) + { + GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.GetGaussPoint(order); + + Complex sum; + double ax; + int i; + int m = (order + 1) >> 1; + + double a = 0.5 * (invervalEnd - invervalBegin); + double b = 0.5 * (invervalEnd + invervalBegin); + + if (order.IsOdd()) + { + sum = gaussLegendrePoint.Weights[0] * f(b); + for (i = 1; i < m; i++) + { + ax = a * gaussLegendrePoint.Abscissas[i]; + sum += gaussLegendrePoint.Weights[i] * (f(b + ax) + f(b - ax)); + } + } + else + { + sum = 0.0; + for (i = 0; i < m; i++) + { + ax = a * gaussLegendrePoint.Abscissas[i]; + sum += gaussLegendrePoint.Weights[i] * (f(b + ax) + f(b - ax)); + } + } + + return a * sum; + } + /// /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. /// From ab1f6083458c50d969bd35367d13524b3297212f Mon Sep 17 00:00:00 2001 From: diluculo Date: Mon, 16 Sep 2019 11:44:00 +0900 Subject: [PATCH 05/13] Fix typos --- .../IntegrationTests/IntegrationTest.cs | 27 ++++++++++++++----- src/Numerics/Integrate.cs | 2 +- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index c8422144..dd988098 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -62,7 +62,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests } /// - /// Test Function: f(x,y) = 1 / (1 + x^2) + /// Test Function: f(x) = 1 / (1 + x^2) /// /// First input value. /// Function result. @@ -72,7 +72,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests } /// - /// Test Function: f(x,y) = log(x) + /// Test Function: f(x) = log(x) /// /// First input value. /// Function result. @@ -120,7 +120,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// Test Function Stop point. /// private const double StopD = 1; - + /// /// Target area square. /// @@ -140,7 +140,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// Target area. /// private const double TargetAreaD = -1; - + /// /// Test Integrate facade for simple use cases. /// @@ -183,6 +183,21 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests 1e-10, "DoubleExponential, Target 1e-10"); + // integrate_(0)^(1) log(x) dx = -1 + // Note that DoubleExponential returns -oo. + + Assert.AreEqual( + TargetAreaD, + Integrate.OnClosedInterval(TargetFunctionD, StartD, StopD), + 1e-10, + "Interval"); + + Assert.AreEqual( + TargetAreaD, + Integrate.GaussLegendre(TargetFunctionD, StartD, StopD, order: 1024), + 1e-10, + "GaussLegendre, order 128"); + Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 15), @@ -440,7 +455,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests Assert.AreEqual( expected, - Integrate.GausLegendre((x) => Math.Exp(-x * x / 2), a, b, order: 128), + 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); } @@ -469,7 +484,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests Assert.AreEqual( expected, - factor * Integrate.GausLegendre((x) => 1 / (1 + x * x), a, b, order: 128), + factor * Integrate.GaussLegendre((x) => 1 / (1 + x * x), a, b, order: 128), 1e-10, "GL Integral of sin(pi*x)/(pi*x) from -oo to oo"); } diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index 6057ed6a..9c7e0370 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -164,7 +164,7 @@ namespace MathNet.Numerics /// Where the interval stops. /// Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2-20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calculated on the fly. /// Approximation of the finite integral in the given interval. - public static double GausLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) + public static double GaussLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) { // Reference: // Formula used for variable subsitution from From e0cd83b64894cfcf1bd32dbce6d2cb1021f52ea6 Mon Sep 17 00:00:00 2001 From: diluculo Date: Mon, 16 Sep 2019 15:52:51 +0900 Subject: [PATCH 06/13] Remove substitution rule for [a, b] => [-1, 1] from DoubleExponential method. --- .../IntegrationTests/IntegrationTest.cs | 15 ++--- src/Numerics/Integrate.cs | 66 +++++++------------ src/Numerics/Integration/GaussKronrodRule.cs | 24 +++---- 3 files changed, 42 insertions(+), 63 deletions(-) diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index dd988098..bc2e967b 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -181,16 +181,13 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests TargetAreaC, Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10), 1e-10, - "DoubleExponential, Target 1e-10"); - - // integrate_(0)^(1) log(x) dx = -1 - // Note that DoubleExponential returns -oo. + "DoubleExponential"); Assert.AreEqual( TargetAreaD, - Integrate.OnClosedInterval(TargetFunctionD, StartD, StopD), + Integrate.DoubleExponential(TargetFunctionD, StartD, StopD), 1e-10, - "Interval"); + "DoubleExponential"); Assert.AreEqual( TargetAreaD, @@ -474,19 +471,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests expected, factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), 1e-10, - "DET Integral of sin(pi*x)/(pi*x) from -oo to oo"); + "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 -oo to oo"); + "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 -oo to oo"); + "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) ¥ð diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index 9c7e0370..5761f6c9 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -102,7 +102,7 @@ namespace MathNet.Numerics // (-oo, oo) => [-1, 1] // - // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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)) @@ -115,8 +115,8 @@ namespace MathNet.Numerics } // [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 + // 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)) @@ -129,8 +129,8 @@ namespace MathNet.Numerics } // (-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 + // 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)) @@ -141,18 +141,9 @@ namespace MathNet.Numerics }; return DoubleExponentialTransformation.Integrate(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.Integrate(u, -1, 1, targetAbsoluteError); + return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError); } } @@ -178,7 +169,7 @@ namespace MathNet.Numerics // (-oo, oo) => [-1, 1] // - // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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)) @@ -191,8 +182,8 @@ namespace MathNet.Numerics } // [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 + // 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)) @@ -205,8 +196,8 @@ namespace MathNet.Numerics } // (-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 + // 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)) @@ -219,7 +210,7 @@ namespace MathNet.Numerics } // [a, b] => [-1, 1] // - // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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 @@ -292,7 +283,7 @@ namespace MathNet.Numerics // (-oo, oo) => [-1, 1] // - // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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)) @@ -305,8 +296,8 @@ namespace MathNet.Numerics } // [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 + // 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)) @@ -319,8 +310,8 @@ namespace MathNet.Numerics } // (-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 + // 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)) @@ -331,18 +322,9 @@ namespace MathNet.Numerics }; 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); + return DoubleExponentialTransformation.ContourIntegrate(f, intervalBegin, intervalEnd, targetAbsoluteError); } } @@ -368,7 +350,7 @@ namespace MathNet.Numerics // (-oo, oo) => [-1, 1] // - // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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)) @@ -381,8 +363,8 @@ namespace MathNet.Numerics } // [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 + // 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)) @@ -395,8 +377,8 @@ namespace MathNet.Numerics } // (-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 + // 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)) @@ -409,7 +391,7 @@ namespace MathNet.Numerics } // [a, b] => [-1, 1] // - // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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 diff --git a/src/Numerics/Integration/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs index 82e1bf36..3d5224be 100644 --- a/src/Numerics/Integration/GaussKronrodRule.cs +++ b/src/Numerics/Integration/GaussKronrodRule.cs @@ -240,7 +240,7 @@ namespace MathNet.Numerics.Integration // (-oo, oo) => [-1, 1] // - // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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)) @@ -253,8 +253,8 @@ namespace MathNet.Numerics.Integration } // [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 + // 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) @@ -267,8 +267,8 @@ namespace MathNet.Numerics.Integration } // (-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 + // 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) @@ -281,7 +281,7 @@ namespace MathNet.Numerics.Integration } // [a, b] => [-1, 1] // - // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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 @@ -326,7 +326,7 @@ namespace MathNet.Numerics.Integration // (-oo, oo) => [-1, 1] // - // integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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)) @@ -339,8 +339,8 @@ namespace MathNet.Numerics.Integration } // [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 + // 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) @@ -353,8 +353,8 @@ namespace MathNet.Numerics.Integration } // (-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 + // 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) @@ -367,7 +367,7 @@ namespace MathNet.Numerics.Integration } // [a, b] => [-1, 1] // - // integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt + // 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 From c18ceade6444485199a31cff24f1b8c869b3fa95 Mon Sep 17 00:00:00 2001 From: diluculo Date: Mon, 16 Sep 2019 18:43:20 +0900 Subject: [PATCH 07/13] Add more integration tests --- .../IntegrationTests/IntegrationTest.cs | 182 ++++++++++++++++-- src/Numerics/Integrate.cs | 2 +- src/Numerics/Integration/GaussKronrodRule.cs | 8 +- 3 files changed, 171 insertions(+), 21 deletions(-) diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index bc2e967b..fc61f29c 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -51,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. @@ -80,7 +80,37 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests { 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. /// @@ -121,6 +151,36 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// 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. /// @@ -141,17 +201,35 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// 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, @@ -159,72 +237,81 @@ 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"); + "DoubleExponential, 1/(1 + x^2)"); Assert.AreEqual( TargetAreaC, Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10), 1e-10, - "DoubleExponential"); + "DoubleExponential, 1/(1 + x^2)"); + + // TargetFunctionD + // integral_(0)^(1) log(x) dx = -1 Assert.AreEqual( TargetAreaD, Integrate.DoubleExponential(TargetFunctionD, StartD, StopD), 1e-10, - "DoubleExponential"); - + "DoubleExponential, log(x)"); + Assert.AreEqual( TargetAreaD, Integrate.GaussLegendre(TargetFunctionD, StartD, StopD, order: 1024), 1e-10, - "GaussLegendre, order 128"); + "GaussLegendre, log(x), order 1024"); Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 15), 1e-10, - "GaussKronrod, Target 1e-10, order 15"); + "GaussKronrod, log(x), order 15"); Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 21), 1e-10, - "GaussKronrod, Target 1e-10, order 21"); + "GaussKronrod, log(x), order 21"); Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 31), 1e-10, - "GaussKronrod, Target 1e-10, order 31"); + "GaussKronrod, log(x), order 31"); Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 41), 1e-10, - "GaussKronrod, Target 1e-10, order 41"); + "GaussKronrod, log(x), order 41"); Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 51), 1e-10, - "GaussKronrod, Target 1e-10, order 51"); + "GaussKronrod, log(x), order 51"); Assert.AreEqual( TargetAreaD, Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, 1e-10, order: 61), 1e-10, - "GaussKronrod, Target 1e-10, order 61"); + "GaussKronrod, log(x), order 61"); double error, L1; var Q = Integrate.GaussKronrod(TargetFunctionD, StartD, StopD, out error, out L1, 1e-10, order: 15); @@ -233,6 +320,69 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests 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"); } /// diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index 5761f6c9..e0e503df 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -164,7 +164,7 @@ namespace MathNet.Numerics if (intervalBegin > intervalEnd) { - return -GaussLegendreRule.Integrate(f, intervalEnd, intervalBegin, order); + return -GaussLegendre(f, intervalEnd, intervalBegin, order); } // (-oo, oo) => [-1, 1] diff --git a/src/Numerics/Integration/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs index 3d5224be..9d798275 100644 --- a/src/Numerics/Integration/GaussKronrodRule.cs +++ b/src/Numerics/Integration/GaussKronrodRule.cs @@ -231,13 +231,13 @@ namespace MathNet.Numerics.Integration throw new ArgumentNullException(nameof(f)); } - Order = order; - if (intervalBegin > intervalEnd) { return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); } + Order = order; + // (-oo, oo) => [-1, 1] // // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt @@ -317,13 +317,13 @@ namespace MathNet.Numerics.Integration throw new ArgumentNullException(nameof(f)); } - Order = order; - if (intervalBegin > intervalEnd) { return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); } + Order = order; + // (-oo, oo) => [-1, 1] // // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt From 3cd5bd80e3c7995504009bcea548b2fb5ad8fa88 Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 18 Sep 2019 19:29:35 +0900 Subject: [PATCH 08/13] Add GaussKronrodPointFactory, but not yet support other than the precomputed order. --- src/Numerics/Integration/GaussKronrodRule.cs | 721 +++++------------- .../GaussRule/GaussKronrodPoint.cs | 382 ++++++++++ .../GaussRule/GaussKronrodPointFactory.cs | 36 + .../Integration/GaussRule/GaussPoints.cs | 35 + 4 files changed, 626 insertions(+), 548 deletions(-) create mode 100644 src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs create mode 100644 src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs create mode 100644 src/Numerics/Integration/GaussRule/GaussPoints.cs diff --git a/src/Numerics/Integration/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs index 9d798275..f3e95d6c 100644 --- a/src/Numerics/Integration/GaussKronrodRule.cs +++ b/src/Numerics/Integration/GaussKronrodRule.cs @@ -35,178 +35,63 @@ // 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 static class GaussKronrodRule + public class GaussKronrodRule { - const double epsilon = 2.2204460492503131e-016; + private readonly GaussPoints gaussKronrodPoint; /// - /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points. + /// Getter for the order. /// - static int Order = 15; - - static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1) + public int Order { - int gauss_start = 2; - int kronrod_start = 1; - int gauss_order = ((int)Order - 1) / 2; - - double kronrod_result = 0d; - double gauss_result = 0d; - double fp, fm; - - var KAbscissa = KronrodAbscissa(); - var KWeights = KronrodWeights(); - var GWeights = GaussWeights(); - - 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) + get { - 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]; + return gaussKronrodPoint.Order; } - 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 * epsilon * 2d)); - return kronrod_result; } - static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1) + /// + /// Getter that returns a clone of the array containing the Kronrod abscissas. + /// + public double[] KronrodAbscissas { - int gauss_start = 2; - int kronrod_start = 1; - int gauss_order = ((int)Order - 1) / 2; - - Complex kronrod_result = new Complex(); - Complex gauss_result = new Complex(); - Complex fp, fm; - - var KAbscissa = KronrodAbscissa(); - var KWeights = KronrodWeights(); - var GWeights = GaussWeights(); - - if (gauss_order.IsOdd()) + get { - 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; + return gaussKronrodPoint.Abscissas.Clone() as double[]; } - 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 * epsilon * 2d)); - return kronrod_result; } - 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) + /// + /// Getter that returns a clone of the array containing the Kronrod weights. + /// + public double[] KronrodWeights { - 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); - 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)) + get { - 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); - estimate += recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local); - error += error_local; - L1 += L1_local; - return estimate; + return gaussKronrodPoint.Weights.Clone() as double[]; } - L1 *= scale; - error = error_local; - return estimate; } - 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) + /// + /// Getter that returns a clone of the array containing the Gauss weights. + /// + public double[] GaussWeights { - 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); - var estimate = scale * r1; - - var tmp = estimate * rel_tol; - var abs_tol1 = Complex.Abs(tmp); - if (abs_tol == 0) + get { - abs_tol = abs_tol1; + return gaussKronrodPoint.SecondWeights.Clone() as double[]; } + } - 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); - estimate += contour_recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local); - error += error_local; - L1 += L1_local; - return estimate; - } - L1 *= scale; - error = error_local; - return estimate; + public GaussKronrodRule(int order) + { + gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); } /// @@ -236,7 +121,7 @@ namespace MathNet.Numerics.Integration return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); } - Order = order; + GaussPoints gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); // (-oo, oo) => [-1, 1] // @@ -249,7 +134,7 @@ namespace MathNet.Numerics.Integration { 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); + return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); } // [a, oo) => [0, 1] // @@ -263,7 +148,7 @@ namespace MathNet.Numerics.Integration { 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); + return recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); } // (-oo, b] => [-1, 0] // @@ -277,7 +162,7 @@ namespace MathNet.Numerics.Integration { 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); + return recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); } // [a, b] => [-1, 1] // @@ -290,7 +175,7 @@ namespace MathNet.Numerics.Integration { 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); + return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm, gaussKronrodPoint); } } @@ -322,7 +207,7 @@ namespace MathNet.Numerics.Integration return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); } - Order = order; + GaussPoints gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); // (-oo, oo) => [-1, 1] // @@ -335,7 +220,7 @@ namespace MathNet.Numerics.Integration { 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); + return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); } // [a, oo) => [0, 1] // @@ -349,7 +234,7 @@ namespace MathNet.Numerics.Integration { 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); + return contour_recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); } // (-oo, b] => [-1, 0] // @@ -363,7 +248,7 @@ namespace MathNet.Numerics.Integration { 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); + return contour_recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint); } // [a, b] => [-1, 1] // @@ -376,428 +261,168 @@ namespace MathNet.Numerics.Integration { 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); + return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm, gaussKronrodPoint); } } - #region Pre-computed Abscissa and weights - - static double[] KronrodAbscissa() + private static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPoints gaussKronrodPoint) { - switch (Order) + 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) { - default: - case 15: - return PrecomputedKronrodAbscissas[0]; - case 21: - return PrecomputedKronrodAbscissas[1]; - case 31: - return PrecomputedKronrodAbscissas[2]; - case 41: - return PrecomputedKronrodAbscissas[3]; - case 51: - return PrecomputedKronrodAbscissas[4]; - case 61: - return PrecomputedKronrodAbscissas[5]; + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_result += fp * GWeights[0]; } - } - - static double[] KronrodWeights() - { - switch (Order) + else { - default: - case 15: - return PrecomputedKronrodWeights[0]; - case 21: - return PrecomputedKronrodWeights[1]; - case 31: - return PrecomputedKronrodWeights[2]; - case 41: - return PrecomputedKronrodWeights[3]; - case 51: - return PrecomputedKronrodWeights[4]; - case 61: - return PrecomputedKronrodWeights[5]; + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_start = 1; + kronrod_start = 2; } - } + double L1 = Math.Abs(kronrod_result); - static double[] GaussWeights() - { - switch (Order) + 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) { - default: - case 15: - return PrecomputedGaussWeights[0]; - case 21: - return PrecomputedGaussWeights[1]; - case 31: - return PrecomputedGaussWeights[2]; - case 41: - return PrecomputedGaussWeights[3]; - case 51: - return PrecomputedGaussWeights[4]; - case 61: - return PrecomputedGaussWeights[5]; + 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; } - /// - /// precomputed abscissa vector per order 15, 21, 31, 41, 51 and 61 - /// - static readonly double[][] PrecomputedKronrodAbscissas = + private static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPoints gaussKronrodPoint) { - new[] // 15-point Gauss-Kronrod - { - 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[] // 21-point Gauss-Kronrod - { - 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[] // 31-point Gauss-Kronrod + 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()) { - 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[] // 41-point Gauss-Kronrod + fp = f(0); + kronrod_result = fp * KWeights[0]; + gauss_result += fp * GWeights[0]; + } + else { - 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[] // 51-point Gauss-Kronrod + 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) { - 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[] // 61-point Gauss-Kronrod + 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) { - 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, + 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; + } - /// - /// precomputed weight vector per order 15, 21, 31, 41, 51 and 61 - /// - static readonly double[][] PrecomputedKronrodWeights = + 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, GaussPoints gaussKronrodPoint) { - new[] // 15-point Gauss-Kronrod integration - { - 2.09482141084727828e-01, - 2.04432940075298892e-01, - 1.90350578064785410e-01, - 1.69004726639267903e-01, - 1.40653259715525919e-01, - 1.04790010322250184e-01, - 6.30920926299785533e-02, - 2.29353220105292250e-02, - }, - new[] // 21-point Gauss-Kronrod integration - { - 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, - }, - new[] // 31-point Gauss-Kronrod integration - { - 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, - }, - new[] // 41-point Gauss-Kronrod integration - { - 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, - }, - new[] // 51-point Gauss-Kronrod integration + 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) { - 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, - }, - new[] // 61-point Gauss-Kronrod integration + abs_tol = abs_tol1; + } + + if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local)) { - 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, - }, - }; + 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; + } - /// - /// precomputed Gauss weight vector per order 7, 10, 15, 20, 25 and 30 - /// - static readonly double[][] PrecomputedGaussWeights = + 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, GaussPoints gaussKronrodPoint) { - new [] // 7-point Gauss - { - 4.17959183673469388e-01, - 3.81830050505118945e-01, - 2.79705391489276668e-01, - 1.29484966168869693e-01, - }, - new[] // 10-point Gauss - { - 2.95524224714752870e-01, - 2.69266719309996355e-01, - 2.19086362515982044e-01, - 1.49451349150580593e-01, - 6.66713443086881376e-02, - }, - new[] // 15-point Gauss - { - 2.02578241925561273e-01, - 1.98431485327111576e-01, - 1.86161000015562211e-01, - 1.66269205816993934e-01, - 1.39570677926154314e-01, - 1.07159220467171935e-01, - 7.03660474881081247e-02, - 3.07532419961172684e-02, - }, - new[] // 20-point Gauss - { - 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, - }, - new[] // 25-point Gauss - { - 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, - }, - new[] // 30-point Gauss + 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) { - 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, + abs_tol = abs_tol1; } - }; - #endregion Pre-computed Abscissa and weights + 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/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs new file mode 100644 index 00000000..8659663a --- /dev/null +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs @@ -0,0 +1,382 @@ +using System; +using System.Collections.Generic; + +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 GaussPoints(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 GaussPoints(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 GaussPoints(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 GaussPoints(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 GaussPoints(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 GaussPoints(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 GaussPoints Generate(int order, double targetAbsoluteTolerance = 1E-10) + { + throw new NotSupportedException(string.Format("The order of Gauss Kronrod, {0} is not yet supported", order)); + } + } +} diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs new file mode 100644 index 00000000..ee294d37 --- /dev/null +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs @@ -0,0 +1,36 @@ +using System; + +namespace MathNet.Numerics.Integration.GaussRule +{ + /// + /// Creates a Gauss-Kronrod point. + /// + internal static class GaussKronrodPointFactory + { + [ThreadStatic] + private static GaussPoints 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. + /// Required precision to compute the abscissas/weights. 1e-10 is usually fine. + /// Object containing the non-negative abscissas/weights, and order. + public static GaussPoints GetGaussPoint(int order, double targetAbsoluteTolerance = 1E-10) + { + // 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)) + { + //Not yet supported! + gaussKronrodPoint = GaussKronrodPoint.Generate(order, targetAbsoluteTolerance); + } + } + + return gaussKronrodPoint; + } + } +} diff --git a/src/Numerics/Integration/GaussRule/GaussPoints.cs b/src/Numerics/Integration/GaussRule/GaussPoints.cs new file mode 100644 index 00000000..309a7247 --- /dev/null +++ b/src/Numerics/Integration/GaussRule/GaussPoints.cs @@ -0,0 +1,35 @@ +namespace MathNet.Numerics.Integration.GaussRule +{ + /// + /// Contains two set of abscissas, weights, and order. + /// + internal class GaussPoints + { + 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 GaussPoints(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondAbscissas, double[] secondWeights) + { + Order = order; + Abscissas = abscissas; + Weights = weights; + SecondOrder = secondOrder; + SecondAbscissas = secondAbscissas; + SecondWeights = secondWeights; + + } + + internal GaussPoints(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondWeights) + : this(order, abscissas, weights, secondOrder, null, secondWeights) + { } + } +} From 0f501660a632533560e2cf869a1a11e0959cc1ba Mon Sep 17 00:00:00 2001 From: diluculo Date: Fri, 20 Sep 2019 20:44:58 +0900 Subject: [PATCH 09/13] Add StieltjesPolynomial to support other than precomputed orders in Gauss-Kronrod rule. --- .../IntegrationTests/IntegrationTest.cs | 17 +- src/Numerics/Integration/GaussKronrodRule.cs | 14 +- .../GaussRule/GaussKronrodPoint.cs | 201 +++++++++++++++++- .../GaussRule/GaussKronrodPointFactory.cs | 8 +- .../Integration/GaussRule/GaussPointPair.cs | 40 ++++ .../Integration/GaussRule/GaussPoints.cs | 35 --- 6 files changed, 257 insertions(+), 58 deletions(-) create mode 100644 src/Numerics/Integration/GaussRule/GaussPointPair.cs delete mode 100644 src/Numerics/Integration/GaussRule/GaussPoints.cs diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index fc61f29c..a00d2af1 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -549,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]); } } @@ -576,6 +576,21 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests 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(7)] + [TestCase(8)] + [TestCase(9)] + [TestCase(10)] + 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) diff --git a/src/Numerics/Integration/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs index f3e95d6c..467dc076 100644 --- a/src/Numerics/Integration/GaussKronrodRule.cs +++ b/src/Numerics/Integration/GaussKronrodRule.cs @@ -43,7 +43,7 @@ namespace MathNet.Numerics.Integration { public class GaussKronrodRule { - private readonly GaussPoints gaussKronrodPoint; + private readonly GaussPointPair gaussKronrodPoint; /// /// Getter for the order. @@ -121,7 +121,7 @@ namespace MathNet.Numerics.Integration return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); } - GaussPoints gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); + GaussPointPair gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); // (-oo, oo) => [-1, 1] // @@ -207,7 +207,7 @@ namespace MathNet.Numerics.Integration return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); } - GaussPoints gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); + GaussPointPair gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); // (-oo, oo) => [-1, 1] // @@ -265,7 +265,7 @@ namespace MathNet.Numerics.Integration } } - private static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPoints 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; @@ -314,7 +314,7 @@ namespace MathNet.Numerics.Integration return kronrod_result; } - private static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPoints gaussKronrodPoint) + 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; @@ -363,7 +363,7 @@ namespace MathNet.Numerics.Integration 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, GaussPoints gaussKronrodPoint) + 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; @@ -394,7 +394,7 @@ namespace MathNet.Numerics.Integration 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, GaussPoints gaussKronrodPoint) + 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; diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs index 8659663a..fad3c88a 100644 --- a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs @@ -1,5 +1,6 @@ using System; using System.Collections.Generic; +using System.Linq; namespace MathNet.Numerics.Integration.GaussRule { @@ -11,9 +12,9 @@ namespace MathNet.Numerics.Integration.GaussRule /// /// Precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61. /// - internal static readonly Dictionary PreComputed = new Dictionary + internal static readonly Dictionary PreComputed = new Dictionary { - { 15, new GaussPoints(15, + { 15, new GaussPointPair(15, new[] // 15-point Gauss-Kronrod Abscissa { 0.00000000000000000e+00, @@ -44,7 +45,7 @@ namespace MathNet.Numerics.Integration.GaussRule 1.29484966168869693e-01, }) }, - { 21, new GaussPoints(21, + { 21, new GaussPointPair(21, new[] // 21-point Gauss-Kronrod Abscissa { 0.00000000000000000e+00, @@ -82,7 +83,7 @@ namespace MathNet.Numerics.Integration.GaussRule 6.66713443086881376e-02, }) }, - { 31, new GaussPoints(31, + { 31, new GaussPointPair(31, new[] // 31-point Gauss-Kronrod Abscissa { 0.00000000000000000e+00, @@ -133,7 +134,7 @@ namespace MathNet.Numerics.Integration.GaussRule 3.07532419961172684e-02, }) }, - { 41, new GaussPoints(41, + { 41, new GaussPointPair(41, new[] // 41-point Gauss-Kronrod Abscissa { 0.00000000000000000e+00, @@ -196,7 +197,7 @@ namespace MathNet.Numerics.Integration.GaussRule 1.76140071391521183e-02, }) }, - { 51, new GaussPoints(51, + { 51, new GaussPointPair(51, new[] // 51-point Gauss-Kronrod Abscissa { 0.00000000000000000e+00, @@ -272,7 +273,7 @@ namespace MathNet.Numerics.Integration.GaussRule 1.13937985010262879e-02, }) }, - { 61, new GaussPoints(61, + { 61, new GaussPointPair(61, new[] // 61-point Gauss-Kronrod Abscissa { 0.00000000000000000e+00, @@ -372,11 +373,191 @@ namespace MathNet.Numerics.Integration.GaussRule /// 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 GaussPoints Generate(int order, double targetAbsoluteTolerance = 1E-10) + internal static GaussPointPair Generate(int order) { - throw new NotSupportedException(string.Format("The order of Gauss Kronrod, {0} is not yet supported", order)); + 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; + + // Build polynomials + + var E = StieltjesPolynomial(gaussOrder + 1); // Stieltjes polynomials + var Eprime = E.Differentiate(); // derivative of E + var L = LegendrePolynomial(gaussOrder); // Legendre polynomials + var Lprime = L.Differentiate(); // derivative of L + + // Calculate Abscissa for Kronrod polynomial, E + + var roots = E.Roots(); + var kronrodAbscissas = roots + .Where(v => v.Imaginary == 0 && v.Real >= 0) + .Select(v => v.Real) + .Distinct() + .ToArray(); + + if (Math.Abs(kronrodAbscissas.Length - gaussAbscissas.Length) > 1) + throw new NotSupportedException("Fail to calculate Abscissas of Gauss-Kronrod rule."); + + // 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 p = Lprime.Evaluate(x); + var w2 = 2.0 / ((1.0 - x * x) * p * p); // Gauss weight + weights[i] = w2 + 2.0 / ((gaussOrder + 1.0) * p * E.Evaluate(x)); + } + for (int i = kronrodStart; i < abscissas.Length; i += 2) + { + var x = abscissas[i]; + weights[i] = 2.0 / ((gaussOrder + 1.0) * L.Evaluate(x) * Eprime.Evaluate(x)); + } + + return new GaussPointPair(order, abscissas, weights, gaussOrder, gaussWeights); + } + + /// + /// Returns the Stieltjes Polynomial of order. + /// + internal static Polynomial StieltjesPolynomial(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.org + // + // 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) + return LegendrePolynomial(1); + else if (order == 2) + return LegendrePolynomial(2) - 2 / 5 * LegendrePolynomial(0); + else if (order == 3) + return LegendrePolynomial(3) - 9 / 14 * LegendrePolynomial(1); + else if (order == 4) + return LegendrePolynomial(4) - 30 / 27 * LegendrePolynomial(2) + 14 / 891 * LegendrePolynomial(0); + else if (order == 5) + return LegendrePolynomial(5) - 35 / 44 * LegendrePolynomial(3) + 135 / 12584 * LegendrePolynomial(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 = 1d; + a[r - k] = 0d; + 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; + } + } + + // First few Legendre polynomials. + Polynomial[] legendrePolynomials = new Polynomial[] + { + new Polynomial(new[] { 1.0 }), + new Polynomial(new[] { 0.0, 1.0 }), + }; + var x = legendrePolynomials[1]; + + var P2 = n.IsOdd() ? legendrePolynomials[0] : legendrePolynomials[1]; + var P1 = n.IsOdd() ? legendrePolynomials[0] : legendrePolynomials[1]; + var P0 = legendrePolynomials[0]; + int degree = P2.Degree; + + Polynomial E = a[1] * P2; + for (int i = 2; i <= r; i++) + { + // Calculate Legendre Polynomial of degree (2 * i - 1 - q) + for (int k = 0; k < 2; k++) + { + degree++; + P2 = ((2 * degree - 1) * x * P1 - (degree - 1) * P0) / degree; + P0 = P1; + P1 = P2; + } + E += a[i] * P2; + } + + return E; + } + + /// + /// Returns the Legendre polynomial of order. + /// + internal static Polynomial LegendrePolynomial(int order) + { + // Calculate P[n, x] by recursion relations + // + // (n + 1) * P[n + 1, x] = (2 * n + 1) * x * P[n, x] - n * P[n - 1, x] + + Polynomial[] legendrePolynomials = new Polynomial[] + { + new Polynomial(new double[] { 1 }), + new Polynomial(new double[] { 0, 1 }), + }; + + if (order < legendrePolynomials.Length) + return legendrePolynomials[order]; + + var x = legendrePolynomials[1]; + + var P2 = legendrePolynomials[0]; + var P1 = legendrePolynomials[0]; + var P0 = legendrePolynomials[0]; + for (int i = 1; i <= order; i++) + { + var n = i - 1; + P2 = ((2 * n + 1) * x * P1 - n * P0) / (n + 1); + P0 = P1; + P1 = P2; + } + var L = P2; + + return L; } } } diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs index ee294d37..fdce236f 100644 --- a/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs @@ -8,15 +8,14 @@ namespace MathNet.Numerics.Integration.GaussRule internal static class GaussKronrodPointFactory { [ThreadStatic] - private static GaussPoints gaussKronrodPoint; + 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. - /// Required precision to compute the abscissas/weights. 1e-10 is usually fine. /// Object containing the non-negative abscissas/weights, and order. - public static GaussPoints GetGaussPoint(int order, double targetAbsoluteTolerance = 1E-10) + public static GaussPointPair GetGaussPoint(int order) { // Try to get the GaussKronrodPoint from the cached static field. bool gaussKronrodPointIsCached = gaussKronrodPoint != null && gaussKronrodPoint.Order == order; @@ -25,8 +24,7 @@ namespace MathNet.Numerics.Integration.GaussRule // Try to find the GaussKronrodPoint in the precomputed dictionary. if (!GaussKronrodPoint.PreComputed.TryGetValue(order, out gaussKronrodPoint)) { - //Not yet supported! - gaussKronrodPoint = GaussKronrodPoint.Generate(order, targetAbsoluteTolerance); + gaussKronrodPoint = GaussKronrodPoint.Generate(order); } } 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/GaussRule/GaussPoints.cs b/src/Numerics/Integration/GaussRule/GaussPoints.cs deleted file mode 100644 index 309a7247..00000000 --- a/src/Numerics/Integration/GaussRule/GaussPoints.cs +++ /dev/null @@ -1,35 +0,0 @@ -namespace MathNet.Numerics.Integration.GaussRule -{ - /// - /// Contains two set of abscissas, weights, and order. - /// - internal class GaussPoints - { - 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 GaussPoints(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondAbscissas, double[] secondWeights) - { - Order = order; - Abscissas = abscissas; - Weights = weights; - SecondOrder = secondOrder; - SecondAbscissas = secondAbscissas; - SecondWeights = secondWeights; - - } - - internal GaussPoints(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondWeights) - : this(order, abscissas, weights, secondOrder, null, secondWeights) - { } - } -} From eda7f568a6e7e4056473ded67edb2b3d57ea8f16 Mon Sep 17 00:00:00 2001 From: diluculo Date: Tue, 24 Sep 2019 19:13:13 +0900 Subject: [PATCH 10/13] Apply the Clenshaw algorithm to more accurately calculate the sum of Legendre series. --- .../IntegrationTests/IntegrationTest.cs | 10 +- .../GaussRule/GaussKronrodPoint.cs | 225 +++++++++++------- .../GaussRule/GaussKronrodPointFactory.cs | 2 +- 3 files changed, 152 insertions(+), 85 deletions(-) diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs index a00d2af1..546a46df 100644 --- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs +++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs @@ -580,10 +580,12 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests /// 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(7)] - [TestCase(8)] - [TestCase(9)] - [TestCase(10)] + [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); diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs index fad3c88a..49289e0a 100644 --- a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs @@ -1,6 +1,7 @@ using System; using System.Collections.Generic; using System.Linq; +using System.Numerics; namespace MathNet.Numerics.Integration.GaussRule { @@ -373,8 +374,9 @@ namespace MathNet.Numerics.Integration.GaussRule /// 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) + internal static GaussPointPair Generate(int order, double eps) { int gaussOrder = (order - 1) / 2; int gaussStart = gaussOrder.IsOdd() ? 0 : 1; @@ -384,24 +386,34 @@ namespace MathNet.Numerics.Integration.GaussRule var gaussAbscissas = gaussPoint.Abscissas; var gaussWeights = gaussPoint.Weights; - // Build polynomials + // Calculate Kronrod polynomial in terms of Legendre polynomials + // K(x) = c0*P(0, x) + c1*P(1, x) + ... - var E = StieltjesPolynomial(gaussOrder + 1); // Stieltjes polynomials - var Eprime = E.Differentiate(); // derivative of E - var L = LegendrePolynomial(gaussOrder); // Legendre polynomials - var Lprime = L.Differentiate(); // derivative of L + var c = StieltjesP(gaussOrder + 1); - // Calculate Abscissa for Kronrod polynomial, E + // Calculate Abscissas for Kronrod polynomial - var roots = E.Roots(); - var kronrodAbscissas = roots - .Where(v => v.Imaginary == 0 && v.Real >= 0) - .Select(v => v.Real) - .Distinct() - .ToArray(); + int r = gaussOrder.IsOdd() ? (gaussOrder - 1) / 2 + 1 : gaussOrder / 2 + 1; + var kronrodAbscissas = new double[r]; - if (Math.Abs(kronrodAbscissas.Length - gaussAbscissas.Length) > 1) - throw new NotSupportedException("Fail to calculate Abscissas of Gauss-Kronrod rule."); + 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); + + kronrodAbscissas[(k - 1) / 2] = x0; + } // Concatenate two abscissas @@ -416,34 +428,42 @@ namespace MathNet.Numerics.Integration.GaussRule for (int i = gaussStart; i < abscissas.Length; i += 2) { var x = abscissas[i]; - var p = Lprime.Evaluate(x); + + 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.Evaluate(x)); + weights[i] = w2 + 2.0 / ((gaussOrder + 1.0) * p * E.Item1); } for (int i = kronrodStart; i < abscissas.Length; i += 2) { var x = abscissas[i]; - weights[i] = 2.0 / ((gaussOrder + 1.0) * L.Evaluate(x) * Eprime.Evaluate(x)); + + 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 the Stieltjes Polynomial of order. + /// Returns coefficients of a Stieltjes polynomial in terms of Legendre polynomials. /// - internal static Polynomial StieltjesPolynomial(int order) + 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.org + // 3. Legendre-Stieltjes Polynomials, Boost.Math // - // Here, we are using Patterson algorithm, expanding the Stieltjes polynomial in terms of legendre polynomials. + // 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] + // 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, @@ -452,16 +472,16 @@ namespace MathNet.Numerics.Integration.GaussRule // // The added n + 1 Kronrod abscissae is the roots of the Kronrod polynomial. - if (order == 1) - return LegendrePolynomial(1); - else if (order == 2) - return LegendrePolynomial(2) - 2 / 5 * LegendrePolynomial(0); - else if (order == 3) - return LegendrePolynomial(3) - 9 / 14 * LegendrePolynomial(1); - else if (order == 4) - return LegendrePolynomial(4) - 30 / 27 * LegendrePolynomial(2) + 14 / 891 * LegendrePolynomial(0); - else if (order == 5) - return LegendrePolynomial(5) - 35 / 44 * LegendrePolynomial(3) + 135 / 12584 * LegendrePolynomial(1); + 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) - 30/27 * P(2, x) + P(4, x) + return new double[] { 0.0157126823793490460157126823793, 0, -1.11111111111111111111111111111, 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; @@ -485,79 +505,124 @@ namespace MathNet.Numerics.Integration.GaussRule for (int k = 1; k < r; k++) { double ratio = 1d; + BigInteger num = 1; + BigInteger den = 1; + double rat = 1d; a[r - k] = 0d; for (int i = r + 1 - k; i <= r; i++) { + num = num * (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); + den = den * (n - q + 2 * (i - k)) * (2 * (k + i - 1) - q - n) * (n + 1 + q + 2 * (k - i)) * (n - 1 - q + 2 * (i + k)); + + var gcd = Euclid.GreatestCommonDivisor(num, den); + num = num / gcd; + den = den / gcd; + + rat = (double)num / (double)den; 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; + + a[r - k] -= a[i] * rat; // ratio; } } - // First few Legendre polynomials. - Polynomial[] legendrePolynomials = new Polynomial[] - { - new Polynomial(new[] { 1.0 }), - new Polynomial(new[] { 0.0, 1.0 }), - }; - var x = legendrePolynomials[1]; - - var P2 = n.IsOdd() ? legendrePolynomials[0] : legendrePolynomials[1]; - var P1 = n.IsOdd() ? legendrePolynomials[0] : legendrePolynomials[1]; - var P0 = legendrePolynomials[0]; - int degree = P2.Degree; - - Polynomial E = a[1] * P2; - for (int i = 2; i <= r; i++) + // K = sum c[k] P[k, x] + + double[] c = new double[2 * r - q]; + for (int i = 1; i < a.Length; i++) { - // Calculate Legendre Polynomial of degree (2 * i - 1 - q) - for (int k = 0; k < 2; k++) - { - degree++; - P2 = ((2 * degree - 1) * x * P1 - (degree - 1) * P0) / degree; - P0 = P1; - P1 = P2; - } - E += a[i] * P2; + c[2 * i - 1 - q] = a[i]; } - return E; + return c; } /// - /// Returns the Legendre polynomial of order. + /// Return value and derivative of a Legendre series at given points. /// - internal static Polynomial LegendrePolynomial(int order) + internal static Tuple LegendreSeries(double[] a, double x) { - // Calculate P[n, x] by recursion relations + // 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 // - // (n + 1) * P[n + 1, x] = (2 * n + 1) * x * P[n, x] - n * P[n - 1, x] + // 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; - Polynomial[] legendrePolynomials = new Polynomial[] + for (int k = a.Length - 1; k >= 1; k--) { - new Polynomial(new double[] { 1 }), - new Polynomial(new double[] { 0, 1 }), - }; + 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; - if (order < legendrePolynomials.Length) - return legendrePolynomials[order]; + 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; - var x = legendrePolynomials[1]; + return new Tuple( value, derivative ); + } - var P2 = legendrePolynomials[0]; - var P1 = legendrePolynomials[0]; - var P0 = legendrePolynomials[0]; - for (int i = 1; i <= order; i++) + /// + /// 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++) { - var n = i - 1; - P2 = ((2 * n + 1) * x * P1 - n * P0) / (n + 1); - P0 = P1; - P1 = P2; + 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 L = P2; - return L; + 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 index fdce236f..372fb577 100644 --- a/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs @@ -24,7 +24,7 @@ namespace MathNet.Numerics.Integration.GaussRule // Try to find the GaussKronrodPoint in the precomputed dictionary. if (!GaussKronrodPoint.PreComputed.TryGetValue(order, out gaussKronrodPoint)) { - gaussKronrodPoint = GaussKronrodPoint.Generate(order); + gaussKronrodPoint = GaussKronrodPoint.Generate(order, 1E-10); } } From 8b449c76d84780357df67c1d41f644ddd6b206f1 Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 25 Sep 2019 10:55:44 +0900 Subject: [PATCH 11/13] Clean up codes. --- .paket/Paket.Restore.targets | 2 +- .../GaussRule/GaussKronrodPoint.cs | 20 ++++--------------- 2 files changed, 5 insertions(+), 17 deletions(-) diff --git a/.paket/Paket.Restore.targets b/.paket/Paket.Restore.targets index 952ad42f..ca9842d8 100644 --- a/.paket/Paket.Restore.targets +++ b/.paket/Paket.Restore.targets @@ -203,7 +203,7 @@ $([System.String]::Copy('%(PaketReferencesFileLines.Identity)').Split(',')[0]) $([System.String]::Copy('%(PaketReferencesFileLines.Identity)').Split(',')[1]) $([System.String]::Copy('%(PaketReferencesFileLines.Identity)').Split(',')[4]) - $([System.String]::Copy('%(PaketReferencesFileLines.Identity)').Split(',')[5]) + $([System.String]::Copy('%(PaketReferencesFileLines.Identity)').Split(',')[5]) %(PaketReferencesFileLinesInfo.PackageVersion) diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs index 49289e0a..8b682b86 100644 --- a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs @@ -414,7 +414,7 @@ namespace MathNet.Numerics.Integration.GaussRule kronrodAbscissas[(k - 1) / 2] = x0; } - + // Concatenate two abscissas var abscissas = new double[gaussAbscissas.Length + kronrodAbscissas.Length]; @@ -504,26 +504,14 @@ namespace MathNet.Numerics.Integration.GaussRule a[r] = 1.0; for (int k = 1; k < r; k++) { - double ratio = 1d; - BigInteger num = 1; - BigInteger den = 1; - double rat = 1d; - a[r - k] = 0d; + double ratio = 1.0; + a[r - k] = 0.0; for (int i = r + 1 - k; i <= r; i++) { - num = num * (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); - den = den * (n - q + 2 * (i - k)) * (2 * (k + i - 1) - q - n) * (n + 1 + q + 2 * (k - i)) * (n - 1 - q + 2 * (i + k)); - - var gcd = Euclid.GreatestCommonDivisor(num, den); - num = num / gcd; - den = den / gcd; - - rat = (double)num / (double)den; 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] * rat; // ratio; + a[r - k] -= a[i] * ratio; } } From c3b180fd99a2f0f7e7ce4d6126dd527a9c1d8b08 Mon Sep 17 00:00:00 2001 From: diluculo Date: Wed, 25 Sep 2019 12:00:30 +0900 Subject: [PATCH 12/13] Fix typos. --- .../Integration/GaussRule/GaussKronrodPoint.cs | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs index 8b682b86..e1f4e93a 100644 --- a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs +++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs @@ -412,9 +412,11 @@ namespace MathNet.Numerics.Integration.GaussRule } 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]; @@ -475,13 +477,13 @@ namespace MathNet.Numerics.Integration.GaussRule 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 }; + 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) - 30/27 * P(2, x) + P(4, x) - return new double[] { 0.0157126823793490460157126823793, 0, -1.11111111111111111111111111111, 0, 1 }; + 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 }; + return new double[] { 0, 0.0107279084551811824539097266370, 0, -0.795454545454545454545454545455, 0, 1 }; int n = order - 1; int q = n.IsOdd() ? 1 : 0; From a7f18f291c262c9886c4aeb686c663787d1372ea Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 13 Oct 2019 14:49:07 +0200 Subject: [PATCH 13/13] Revert breaking change --- src/Numerics/Integrate.cs | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs index e0e503df..b7218248 100644 --- a/src/Numerics/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -46,11 +46,23 @@ namespace MathNet.Numerics /// 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 double OnClosedInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) + public static double OnClosedInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError) { return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError); } + /// + /// Approximation of the definite integral of an analytic smooth function on a closed interval. + /// + /// The analytic smooth function to integrate. + /// 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 double OnClosedInterval(Func f, double intervalBegin, double intervalEnd) + { + return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, 1e-8); + } + /// /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d]. /// @@ -91,7 +103,7 @@ namespace MathNet.Numerics public static double DoubleExponential(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) { // Reference: - // Formula used for variable subsitution from + // 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 @@ -158,7 +170,7 @@ namespace MathNet.Numerics public static double GaussLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) { // Reference: - // Formula used for variable subsitution from + // 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 @@ -272,7 +284,7 @@ namespace MathNet.Numerics public static Complex DoubleExponential(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8) { // Reference: - // Formula used for variable subsitution from + // 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 @@ -339,7 +351,7 @@ namespace MathNet.Numerics public static Complex GaussLegendre(Func f, double intervalBegin, double intervalEnd, int order = 128) { // Reference: - // Formula used for variable subsitution from + // 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