Gauss-Legendre Rule Feature

- Support for obtaining the Gauss-Legendre abscissas/weights
- Integration in 1D and 2D.
- Added Unit tests

src/Numerics/Integration/GaussLegendreRule.cs

@@ -0,0 +1,263 @@
+// <copyright file="GaussLegendreRule.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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.
+// </copyright>
+
+using System;
+using MathNet.Numerics.Integration.GaussRule;
+
+namespace MathNet.Numerics.Integration
+{
+    /// <summary>
+    /// Approximates a definite integral using an Nth order Gauss-Legendre rule. Precomputed Gauss-Legendre abscissas/weights for orders 2,. . ., 20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calulcated on the fly.
+    /// </summary>
+    public class GaussLegendreRule
+    {
+        private readonly GaussPoint gaussLegendrePoint;
+        private static GaussPoint nonNegativeGaussLegendrePoint;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="GaussLegendreRule"/> class.
+        /// </summary>
+        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
+        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
+        /// <param name="order">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 calulcated on the fly.</param>
+        public GaussLegendreRule(double intervalBegin, double intervalEnd, int order)
+        {
+            gaussLegendrePoint = Map(GaussLegendrePointFactory.GetGaussPoint(order), intervalBegin, intervalEnd);
+        }
+
+        /// <summary>
+        /// Gettter for the ith abscissa.
+        /// </summary>
+        /// <param name="index">Index of the ith abscissa.</param>
+        /// <returns>The ith abscissa.</returns>
+        public double GetAbscissa(int index)
+        {
+            return gaussLegendrePoint.Abscissas[index];
+        }
+
+        /// <summary>
+        /// Getter for the ith weight.
+        /// </summary>
+        /// <param name="index">Index of the ith weight.</param>
+        /// <returns>The ith weight.</returns>
+        public double GetWeight(int index)
+        {
+            return gaussLegendrePoint.Weights[index];
+        }
+
+        /// <summary>
+        /// Getter for the order.
+        /// </summary>
+        /// <returns>The order, which defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</returns>
+        public int GetOrder()
+        {
+            return gaussLegendrePoint.Order;
+        }
+
+        /// <summary>
+        /// Getter for the InvervalBegin. 
+        /// </summary>
+        /// <returns>Where the interval starts, inclusive and finite.</returns>
+        public double GetIntervalBegin()
+        {
+            return gaussLegendrePoint.IntervalBegin;
+        }
+
+        /// <summary>
+        /// Getter for the InvervalEnd. 
+        /// </summary>
+        /// <returns>Where the interval stops, inclusive and finite.</returns>
+        public double GetIntervalEnd()
+        {
+            return gaussLegendrePoint.IntervalEnd;
+        }
+
+        /// <summary>
+        /// Getter for the GaussPoint.
+        /// </summary>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
+        /// <returns>Object containing the non-negative abscissas/weights, order, and intervalBegin/intervalEnd. The non-negative abscissas/weights are generated over the interval [-1,1] for the given order.</returns>
+        private static GaussPoint GetGaussPoint(int order)
+        {
+            // Try to get the point from the static field first. If it's not there, or it's the wrong order, then generate it on the fly and store it in the static field.
+            if (nonNegativeGaussLegendrePoint == null || nonNegativeGaussLegendrePoint.Order != order)
+            {
+                return GaussLegendrePointFactory.GetGaussPoint(order);
+            }
+
+            return nonNegativeGaussLegendrePoint;
+        }
+
+        /// <summary>
+        /// Maps the non-negative abscissas/weights from the interval [-1, 1] to the interval [intervalBegin, intervalEnd].
+        /// </summary>
+        /// <param name="gaussPoint">Object containing the non-negative abscissas/weights, order, and intervalBegin/intervalEnd. The non-negative abscissas/weights are generated over the interval [-1,1] for the given order.</param>
+        /// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
+        /// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
+        /// <returns>Object containing the abscissas/weights, order, and intervalBegin/intervalEnd.</returns>
+        private static GaussPoint Map(GaussPoint gaussPoint, double intervalBegin, double intervalEnd)
+        {
+            double[] abscissas = new double[gaussPoint.Order];
+            double[] weights = new double[gaussPoint.Order];
+
+            double a = 0.5*(intervalEnd - intervalBegin);
+            double b = 0.5*(intervalEnd + intervalBegin);
+
+            int m = (gaussPoint.Order + 1) >> 1;
+
+            for (int i = 1; i <= m; i++)
+            {
+                int index1 = gaussPoint.Order - i;
+                int index2 = i - 1;
+                int index3 = m - i;
+
+                abscissas[index1] = gaussPoint.Abscissas[index3]*a + b;
+                abscissas[index2] = -gaussPoint.Abscissas[index3]*a + b;
+
+                weights[index1] = gaussPoint.Weights[index3]*a;
+                weights[index2] = gaussPoint.Weights[index3]*a;
+            }
+
+            return new GaussPoint(intervalBegin, intervalEnd, gaussPoint.Order, abscissas, weights);
+        }
+
+        /// <summary>
+        /// Approximates a definite integral using an Nth order Gauss-Legendre rule.
+        /// </summary>
+        /// <param name="f">The analytic smooth function to integrate.</param>
+        /// <param name="invervalBegin">Where the interval starts, exclusive and finite.</param>
+        /// <param name="invervalEnd">Where the interval ends, exclusive and finite.</param>
+        /// <param name="order">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 calulcated on the fly.</param>
+        /// <returns>Approximation of the finite integral in the given interval.</returns>
+        public static double Integrate(Func<double, double> f, double invervalBegin, double invervalEnd, int order)
+        {
+            nonNegativeGaussLegendrePoint = GetGaussPoint(order);
+
+            double sum, ax;
+            int i;
+            int m = (order + 1) >> 1;
+
+            double a = 0.5*(invervalEnd - invervalBegin);
+            double b = 0.5*(invervalEnd + invervalBegin);
+
+            if (order.IsOdd())
+            {
+                sum = nonNegativeGaussLegendrePoint.Weights[0]*f(b);
+                for (i = 1; i < m; i++)
+                {
+                    ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
+                    sum += nonNegativeGaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
+                }
+            }
+            else
+            {
+                sum = 0.0;
+                for (i = 0; i < m; i++)
+                {
+                    ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
+                    sum += nonNegativeGaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
+                }
+            }
+
+            return a*sum;
+        }
+
+        /// <summary>
+        /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d].
+        /// </summary>
+        /// <param name="f">The 2-dimensional analytic smooth function to integrate.</param>
+        /// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
+        /// <param name="invervalEndA">Where the interval ends for the first (inside) integral, exclusive and finite.</param>
+        /// <param name="invervalBeginB">Where the interval starts for the second (outside) integral, exclusive and finite.</param>
+        /// /// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
+        /// <param name="order">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 calulcated on the fly.</param>
+        /// <returns>Approximation of the finite integral in the given interval.</returns>
+        public static double Integrate(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
+        {
+            nonNegativeGaussLegendrePoint = GetGaussPoint(order);
+
+            double ax, cy, sum;
+            int i, j;
+            int m = (order + 1) >> 1;
+
+            double a = 0.5*(invervalEndA - invervalBeginA);
+            double b = 0.5*(invervalEndA + invervalBeginA);
+            double c = 0.5*(invervalEndB - invervalBeginB);
+            double d = 0.5*(invervalEndB + invervalBeginB);
+
+            if (order.IsOdd())
+            {
+                sum = nonNegativeGaussLegendrePoint.Weights[0]*nonNegativeGaussLegendrePoint.Weights[0]*f(b, d);
+
+                double t;
+                for (j = 1, t = 0.0; j < m; j++)
+                {
+                    cy = c*nonNegativeGaussLegendrePoint.Abscissas[j];
+                    t += nonNegativeGaussLegendrePoint.Weights[j]*(f(b, d + cy) + f(b, d - cy));
+                }
+
+                sum += nonNegativeGaussLegendrePoint.Weights[0]*t;
+
+                for (i = 1, t = 0.0; i < m; i++)
+                {
+                    ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
+                    t += nonNegativeGaussLegendrePoint.Weights[i]*(f(b + ax, d) + f(b - ax, d));
+                }
+
+                sum += nonNegativeGaussLegendrePoint.Weights[0]*t;
+
+                for (i = 1; i < m; i++)
+                {
+                    ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
+                    for (j = 1; j < m; j++)
+                    {
+                        cy = c*nonNegativeGaussLegendrePoint.Abscissas[j];
+                        sum += nonNegativeGaussLegendrePoint.Weights[i]*nonNegativeGaussLegendrePoint.Weights[j]*(f(b + ax, d + cy) + f(ax + b, d - cy) + f(b - ax, d + cy) + f(b - ax, d - cy));
+                    }
+                }
+            }
+            else
+            {
+                sum = 0.0;
+                for (i = 0; i < m; i++)
+                {
+                    ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
+                    for (j = 0; j < m; j++)
+                    {
+                        cy = c*nonNegativeGaussLegendrePoint.Abscissas[j];
+                        sum += nonNegativeGaussLegendrePoint.Weights[i]*nonNegativeGaussLegendrePoint.Weights[j]*(f(b + ax, d + cy) + f(ax + b, d - cy) + f(b - ax, d + cy) + f(b - ax, d - cy));
+                    }
+                }
+            }
+
+            return c*a*sum;
+        }
+    }
+}

src/Numerics/Integration/GaussRule/GaussPoint.cs

@@ -0,0 +1,61 @@
+// <copyright file="GaussLegendreRule.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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.
+// </copyright>
+
+namespace MathNet.Numerics.Integration.GaussRule
+{
+    /// <summary>
+    /// Contains the abscissas/weights, order, and intervalBegin/intervalEnd.
+    /// </summary>
+    class GaussPoint
+    {
+        public double[] Abscissas { get; private set; }
+
+        public double[] Weights { get; private set; }
+
+        public double IntervalBegin { get; private set; }
+
+        public double IntervalEnd { get; private set; }
+
+        public int Order { get; private set; }
+
+        public GaussPoint(double intervalBegin, double intervalEnd, int order, double[] abscissas, double[] weights)
+        {
+            Abscissas = abscissas;
+            Weights = weights;
+            IntervalBegin = intervalBegin;
+            IntervalEnd = intervalEnd;
+            Order = order;
+        }
+
+        public GaussPoint(int order, double[] abscissas, double[] weights) : this(-1, 1, order, abscissas, weights)
+        {
+        }
+    }
+}

src/Numerics/Numerics.csproj

@@ -96,6 +96,9 @@
     <Compile Include="GoodnessOfFit.cs" />
     <Compile Include="IntegralTransforms\Fourier.cs" />
     <Compile Include="IntegralTransforms\Hartley.cs" />
+    <Compile Include="Integration\GaussLegendreRule.cs" />
+    <Compile Include="Integration\GaussRule\GaussLegendrePointFactory.cs" />
+    <Compile Include="Integration\GaussRule\GaussPoint.cs" />
     <Compile Include="Interpolation\Barycentric.cs" />
     <Compile Include="Interpolation\CubicSpline.cs" />
     <Compile Include="Interpolation\LogLinear.cs" />

src/UnitTests/IntegrationTests/IntegrationTest.cs

@@ -46,6 +46,17 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
             return Math.Exp(-x / 5) * (2 + Math.Sin(2 * x));
         }
 
+        /// <summary>
+        /// Test Function: f(x,y) = exp(-x/5) (2 + sin(x * y))
+        /// </summary>
+        /// <param name="x">First input value.</param>
+        /// <param name="y">Second input value.</param>
+        /// <returns>Function result.</returns>
+        private static double TargeFunctionB(double x, double y)
+        {
+            return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y));
+        }
+
         /// <summary>
         /// Test Function Start point.
         /// </summary>
@@ -56,11 +67,26 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
         /// </summary>
         private const double StopA = 10;
 
+        /// <summary>
+        /// Test Function Start point.
+        /// </summary>
+        private const double StartB = 0;
+
+        /// <summary>
+        /// Test Function Stop point.
+        /// </summary>
+        private const double StopB = 1;
+
         /// <summary>
         /// Target area square.
         /// </summary>
         private const double TargetAreaA = 9.1082396073229965070;
 
+        /// <summary>
+        /// Target area.
+        /// </summary>
+        private const double TargetAreaB = 11.7078776759298776163;
+
         /// <summary>
         /// Test integrate portal.
         /// </summary>
@@ -178,5 +204,79 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
                 "Composite {0} Partitions",
                 partitions);
         }
+
+        /// <summary>
+        /// Gauss-Legendre rule supports integration.
+        /// </summary>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
+        [TestCase(19)]
+        [TestCase(20)]
+        [TestCase(21)]
+        [TestCase(22)]
+        public void TestGaussLegendreRuleIntegration(int order)
+        {
+            double appoximateArea = GaussLegendreRule.Integrate(TargetFunctionA, StartA, StopA, order);
+            double relativeError = Math.Abs(TargetAreaA - appoximateArea) / TargetAreaA;
+            Assert.Less(relativeError, 5e-16);
+        }
+
+        /// <summary>
+        /// Gauss-Legendre rule supports 2-dimensional integration over the rectangle.
+        /// </summary>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
+        [TestCase(19)]
+        [TestCase(20)]
+        [TestCase(21)]
+        [TestCase(22)]
+        public void TestGaussLegendreRuleIntegrate2D(int order)
+        {
+            double appoximateArea = GaussLegendreRule.Integrate(TargeFunctionB, StartA, StopA, StartB, StopB, order);
+            double relativeError = Math.Abs(TargetAreaB - appoximateArea) / TargetAreaB;
+            Assert.Less(relativeError, 1e-15);
+        }
+
+        /// <summary>
+        /// Gauss-Legendre rule supports obtaining the abscissas/weights. In this case, they're used for integration. 
+        /// </summary>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
+        [TestCase(19)]
+        [TestCase(20)]
+        [TestCase(21)]
+        [TestCase(22)]
+        public void TestGaussLegendreRuleGetAbscissasGetWeightsGetOrderViaIntegration(int order)
+        {
+            GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
+
+            double appoximateArea = 0;
+            for (int i = 0; i < gaussLegendre.GetOrder(); i++)
+            {
+                appoximateArea += gaussLegendre.GetWeight(i) * TargetFunctionA(gaussLegendre.GetAbscissa(i));
+            }
+
+            double relativeError = Math.Abs(TargetAreaA - appoximateArea) / TargetAreaA;
+            Assert.Less(relativeError, 5e-16);
+        }
+
+        /// <summary>
+        /// Gauss-Legendre rule supports obtaining IntervalBegin.
+        /// </summary>
+        [TestCase]
+        public void TestGetGaussLegendreRuleIntervalBegin()
+        {
+            const int order = 19;
+            GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
+            Assert.AreEqual(gaussLegendre.GetIntervalBegin(), StartA);
+        }
+
+        /// <summary>
+        /// Gauss-Legendre rule supports obtaining IntervalEnd.
+        /// </summary>
+        [TestCase]
+        public void TestGaussLegendreRuleGetIntervalEnd()
+        {
+            const int order = 19;
+            GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
+            Assert.AreEqual(gaussLegendre.GetIntervalEnd(), StopA);
+        }
     }
-}
+}

src/UnitTests/UnitTests.csproj

@@ -65,6 +65,10 @@
     <LangVersion>5</LangVersion>
   </PropertyGroup>
   <ItemGroup>
+    <Reference Include="nunit.framework, Version=3.0.5813.39035, Culture=neutral, PublicKeyToken=2638cd05610744eb, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\packages\NUnit\lib\nunit.framework.dll</HintPath>
+    </Reference>
     <Reference Include="System" />
     <Reference Include="System.Core" />
     <Reference Include="System.Numerics" />
@@ -416,11 +420,4 @@
     <None Include="paket.references" />
   </ItemGroup>
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
-  <ItemGroup>
-    <Reference Include="nunit.framework">
-      <HintPath>..\..\packages\test\NUnit\lib\nunit.framework.dll</HintPath>
-      <Private>True</Private>
-      <Paket>True</Paket>
-    </Reference>
-  </ItemGroup>
 </Project>