Merge pull request #371 from larzw/feature/Gauss-Legendre

Gauss-Legendre Rule Feature

src/Numerics/Integrate.cs

@@ -62,5 +62,33 @@ namespace MathNet.Numerics
         {
             return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, 1e-8);
         }
+
+        /// <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 GaussLegendre(Func<double, double> f, double invervalBegin, double invervalEnd, int order)
+        {
+            return GaussLegendreRule.Integrate(f, invervalBegin, invervalEnd, order);
+        }
+
+        /// <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 GaussLegendre(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
+        {
+            return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order);
+        }
     }
 }

src/Numerics/Integration/GaussLegendreRule.cs

@@ -0,0 +1,252 @@
+// <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;
+       
+        /// <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>
+        public int Order
+        {
+            get
+            {
+                return gaussLegendrePoint.Order;
+            }
+        }
+
+        /// <summary>
+        /// Getter for the InvervalBegin. 
+        /// </summary>
+        public double IntervalBegin
+        {
+            get
+            {
+                return gaussLegendrePoint.IntervalBegin;
+            }
+        }
+
+        /// <summary>
+        /// Getter for the InvervalEnd. 
+        /// </summary>
+        public double IntervalEnd
+        {
+            get
+            {
+                return gaussLegendrePoint.IntervalEnd;
+            }
+        }
+
+        /// <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)
+        {
+            GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.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 = 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;
+        }
+
+        /// <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)
+        {
+            GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.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 = gaussLegendrePoint.Weights[0]*gaussLegendrePoint.Weights[0]*f(b, d);
+
+                double t;
+                for (j = 1, t = 0.0; j < m; j++)
+                {
+                    cy = c*gaussLegendrePoint.Abscissas[j];
+                    t += gaussLegendrePoint.Weights[j]*(f(b, d + cy) + f(b, d - cy));
+                }
+
+                sum += gaussLegendrePoint.Weights[0]*t;
+
+                for (i = 1, t = 0.0; i < m; i++)
+                {
+                    ax = a*gaussLegendrePoint.Abscissas[i];
+                    t += gaussLegendrePoint.Weights[i]*(f(b + ax, d) + f(b - ax, d));
+                }
+
+                sum += gaussLegendrePoint.Weights[0]*t;
+
+                for (i = 1; i < m; i++)
+                {
+                    ax = a*gaussLegendrePoint.Abscissas[i];
+                    for (j = 1; j < m; j++)
+                    {
+                        cy = c*gaussLegendrePoint.Abscissas[j];
+                        sum += gaussLegendrePoint.Weights[i]*gaussLegendrePoint.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*gaussLegendrePoint.Abscissas[i];
+                    for (j = 0; j < m; j++)
+                    {
+                        cy = c*gaussLegendrePoint.Abscissas[j];
+                        sum += gaussLegendrePoint.Weights[i]*gaussLegendrePoint.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 TargetFunctionB(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,109 @@ 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 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 TestIntegrateGaussLegendre(int order)
+        {
+            double appoximateArea = Integrate.GaussLegendre(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(TargetFunctionB, StartA, StopA, StartB, StopB, order);
+            double relativeError = Math.Abs(TargetAreaB - appoximateArea) / TargetAreaB;
+            Assert.Less(relativeError, 1e-15);
+        }
+
+        /// <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 TestIntegrateGaussLegendre2D(int order)
+        {
+            double appoximateArea = Integrate.GaussLegendre(TargetFunctionB, 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 TestGaussLegendreRuleGetAbscissasGetWeightsOrderViaIntegration(int order)
+        {
+            GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
+
+            double appoximateArea = 0;
+            for (int i = 0; i < gaussLegendre.Order; 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.IntervalBegin, StartA);
+        }
+
+        /// <summary>
+        /// Gauss-Legendre rule supports obtaining IntervalEnd.
+        /// </summary>
+        [TestCase]
+        public void TestGaussLegendreRuleIntervalEnd()
+        {
+            const int order = 19;
+            GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
+            Assert.AreEqual(gaussLegendre.IntervalEnd, StopA);
+        }
     }
-}
+}