Browse Source

Merge pull request #423 from larzw/gauss-legendre

Gauss-Legendre
netstandard
Christoph Ruegg 10 years ago
committed by GitHub
parent
commit
656af9870f
  1. 126
      docs/content/Integration.md
  2. 2
      docs/tools/templates/template.cshtml
  3. 2
      src/Numerics/Integrate.cs
  4. 65
      src/Numerics/Integration/GaussLegendreRule.cs
  5. 187
      src/Numerics/Integration/GaussRule/GaussLegendrePoint.cs
  6. 174
      src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs
  7. 16
      src/Numerics/Integration/GaussRule/GaussPoint.cs
  8. 1
      src/Numerics/Numerics.csproj
  9. 38
      src/UnitTests/IntegrationTests/IntegrationTest.cs

126
docs/content/Integration.md

@ -9,14 +9,140 @@
Numerical Integration
=====================
The following double precision numerical integration or quadrature rules are supported in Math.NET Numerics under the `MathNet.Numerics.Integration` namespace. Unless stated otherwise, the examples below evaluate the integral $\int_0^{10} x^2 \, dx = \frac{1000}{3} \approx 333.\overline{3}$.
Simpson's Rule
--------------
[lang=csharp]
// Composite approximation with 4 partitions
double composite = SimpsonRule.IntegrateComposite(x => x * x, 0.0, 10.0, 4);
// Approximate value using IntegrateComposite with 4 partitions is: 333.33333333333337
Console.WriteLine("Approximate value using IntegrateComposite with 4 partitions is: " + composite);
// Three point approximation
double threePoint = SimpsonRule.IntegrateThreePoint(x => x * x, 0.0, 10.0);
// Approximate value using IntegrateThreePoint is: 333.333333333333
Console.WriteLine("Approximate value using IntegrateThreePoint is: " + threePoint);
Newton Cotes Trapezium Rule
---------------------------
[lang=csharp]
// Adaptive approximation with a relative error of 1e-5
double adaptive = NewtonCotesTrapeziumRule.IntegrateAdaptive(x => x * x, 0.0, 10.0, 1e-5);
// Approximate value of the integral using IntegrateAdaptive with a relative error of 1e-5 is: 333.333969116211
Console.WriteLine("Approximate value using IntegrateAdaptive with a relative error of 1e-5: " + adaptive);
// Composite approximation with 15 partitions
double composite = NewtonCotesTrapeziumRule.IntegrateComposite(x => x * x, 0.0, 10.0, 15);
//Approximate value of the integral using IntegrateComposite with 15 partitions is: 334.074074074074
Console.WriteLine("Approximate value using IntegrateComposite with 15 partitions is: " + composite);
// Two point approximation
double twoPoint = NewtonCotesTrapeziumRule.IntegrateTwoPoint(x => x * x, 0.0, 10.0);
//Approximate value using IntegrateTwoPoint is: 500
Console.WriteLine("Approximate value using IntegrateTwoPoint is: " + twoPoint);
Double-Exponential Transformation
---------------------------------
The Double-Exponential Transformation is suited for integration of smooth functions with no discontinuities, derivative discontinuities, and poles inside the interval.
[lang=csharp]
// Approximate using a relative error of 1e-5.
double integrate = DoubleExponentialTransformation.Integrate(x => x * x, 0.0, 10.0, 1e-5);
// Approximate value using a relative error of 1e-5 is: 333.333333333332
Console.WriteLine("Approximate value using a relative error of 1e-5 is: " + integrate);
Gauss-Legendre Rule
-------------------
A fixed-order Gauss-Legendre integration routine is provided for fast integration of smooth functions with known polynomial order. The N-point Gauss-Legendre rule is exact for polynomials of order $2N-1$ or less. For example, these rules are useful when integrating basis functions to form mass matrices for the Galerkin method [[GSL]](https://www.gnu.org/software/gsl/).
The basic idea of Gauss-Legendre integration is to approximate the integral of a function $f(x)$ using $N$ Weights $w_i$ and abscissas (or nodes) $x_i$.
$$$
\int_a^b f(x) \, dx \approx \sum_{i = 0}^{N - 1} w_i f(x_i)
This algorithm calculates the abscissas and weights for a given order and integration interval. For efficiency, pre-computed abscissas and weights for the orders $ N = 2 - 20, \, 32, \, 64, \, 96, 100, \, 128, \, 256, \, 512, \, 1024$ are used. Otherwise, they are calculated on the fly using Newton's method. For more information on the algorithm see [[Holoborodko, Pavel] ](http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/).
### Abscissas and Weights
We'll first use the abscissas and weights to approximate an integral using a 5-point Gauss-Legendre rule
[lang=csharp]
// Create a 5-point Gauss-Legendre rule over the integration interval [0, 10]
GaussLegendreRule rule = new GaussLegendreRule(0.0, 10.0, 5);
double sum = 0; // Will hold the approximate value of the integral
for (int i = 0; i < rule.Order; i++) // rule.Order = 5
{
// Access the ith abscissa and weight
sum += rule.GetWeight(i) * rule.GetAbscissa(i) * rule.GetAbscissa(i);
}
// Approximate value is: 333.333333333333
Console.WriteLine("Approximate value is: " + sum);
If you prefer direct access to the abscissas and weights, as opposed to using the methods
- ```double GetAbscissa(int i)```
- ```double GetWeight(int i)```
then use the properties `Abscissas` and `Weights`
[lang=csharp]
// Create a 5-point Gauss-Legendre rule over the integration interval [0, 10]
GaussLegendreRule rule = new GaussLegendreRule(0.0, 10.0, 5);
double[] x = rule.Abscissas; // Creates a clone and returns array of abscissas
double[] w = rule.Weights; // Creates a clone and returns array of weights
double sum = 0; // Will hold the approximate value of the integral
for (int i = 0; i < rule.Order; i++) // rule.Order = 5
{
// Access the ith abscissa and weight
sum += w[i] * x[i] * x[i];
}
// Approximate value is: 333.333333333333
Console.WriteLine("Approximate value is: " + sum);;
In addition to obtaining the abscissas and weights, the order and integration interval can be obtained
[lang=csharp]
// Create a 5-point Gauss-Legendre rule over the integration interval [0, 10]
GaussLegendreRule rule = new GaussLegendreRule(0.0, 10.0, 5);
// The order of the rule is: 5
Console.WriteLine("The order of the rule is: " + rule.Order);
// The lower integral bound is 0
Console.WriteLine("The lower integral bound is: " + rule.IntervalBegin);
// The upper integral bound is 10
Console.WriteLine("The upper integral bound is: " + rule.IntervalEnd);
### Integrate Method
For convenience, we provide an overloaded static method `double Integrate(...)` which preforms 1D and 2D integration of a function. The first parameter to the method is a delegate of type `Func<double, double>` or `Func<double, double, double>` for 1D and 2D integration respectively. So for example
[lang=csharp]
// 1D integration using a 5-point Gauss-Legendre rule over the integration interval [0, 10]
double integrate1D = GaussLegendreRule.Integrate(x => x * x, 0.0, 10.0, 5);
// Approximate value of the 1D integral is: 333.333333333333
Console.WriteLine("Approximate value of the 1D integral is: " + integrate1D);
// 2D integration using a 5-point Gauss-Legendre rule over the integration interval [0, 10] X [1, 2]
double integrate2D = GaussLegendreRule.Integrate((x, y) => (x * x) * (y * y), 0.0, 10.0, 1.0, 2.0, 5);
// Approximate value of the 2D integral is: 777.777777777778
Console.WriteLine("Approximate value of the 2D integral is: " + integrate2D);
where we used $\int_0^{10}\int_1^2 x^2 y^2 \,dydx = \frac{7000}{9} \approx 777.\overline{7}$ for the 2D integral example.

2
docs/tools/templates/template.cshtml

@ -92,7 +92,7 @@
<li class="nav-header">Evaluation</li>
<li><a href="@Root/Functions.html">Special Functions</a></li>
<li>Differentiation</li>
<li>Integration</li>
<li><a href="@Root/Integration.html">Integration</a></li>
<li class="nav-header">Statistics/Probability</li>
<li><a href="@Root/DescriptiveStatistics.html">Descriptive Statistics</a></li>

2
src/Numerics/Integrate.cs

@ -70,7 +70,7 @@ namespace MathNet.Numerics
/// <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>
/// <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 OnRectangle(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
{

65
src/Numerics/Integration/GaussLegendreRule.cs

@ -33,7 +33,7 @@ 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.
/// 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
{
@ -44,10 +44,10 @@ namespace MathNet.Numerics.Integration
/// </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>
/// <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);
_gaussLegendrePoint = GaussLegendrePointFactory.GetGaussPoint(intervalBegin, intervalEnd, order);
}
/// <summary>
@ -60,6 +60,17 @@ namespace MathNet.Numerics.Integration
return _gaussLegendrePoint.Abscissas[index];
}
/// <summary>
/// Getter that returns a clone of the array containing the abscissas.
/// </summary>
public double[] Abscissas
{
get
{
return _gaussLegendrePoint.Abscissas.Clone() as double[];
}
}
/// <summary>
/// Getter for the ith weight.
/// </summary>
@ -70,6 +81,17 @@ namespace MathNet.Numerics.Integration
return _gaussLegendrePoint.Weights[index];
}
/// <summary>
/// Getter that returns a clone of the array containing the weights.
/// </summary>
public double[] Weights
{
get
{
return _gaussLegendrePoint.Weights.Clone() as double[];
}
}
/// <summary>
/// Getter for the order.
/// </summary>
@ -103,46 +125,13 @@ namespace MathNet.Numerics.Integration
}
}
/// <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>
/// <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)
{
@ -185,7 +174,7 @@ namespace MathNet.Numerics.Integration
/// <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>
/// <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)
{

187
src/Numerics/Integration/GaussRule/GaussLegendrePoint.cs

File diff suppressed because one or more lines are too long

174
src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs

File diff suppressed because one or more lines are too long

16
src/Numerics/Integration/GaussRule/GaussPoint.cs

@ -32,19 +32,19 @@ namespace MathNet.Numerics.Integration.GaussRule
/// <summary>
/// Contains the abscissas/weights, order, and intervalBegin/intervalEnd.
/// </summary>
class GaussPoint
internal class GaussPoint
{
public double[] Abscissas { get; private set; }
internal double[] Abscissas { get; private set; }
public double[] Weights { get; private set; }
internal double[] Weights { get; private set; }
public double IntervalBegin { get; private set; }
internal double IntervalBegin { get; private set; }
public double IntervalEnd { get; private set; }
internal double IntervalEnd { get; private set; }
public int Order { get; private set; }
internal int Order { get; private set; }
public GaussPoint(double intervalBegin, double intervalEnd, int order, double[] abscissas, double[] weights)
internal GaussPoint(double intervalBegin, double intervalEnd, int order, double[] abscissas, double[] weights)
{
Abscissas = abscissas;
Weights = weights;
@ -53,7 +53,7 @@ namespace MathNet.Numerics.Integration.GaussRule
Order = order;
}
public GaussPoint(int order, double[] abscissas, double[] weights) : this(-1, 1, order, abscissas, weights)
internal GaussPoint(int order, double[] abscissas, double[] weights) : this(-1, 1, order, abscissas, weights)
{
}
}

1
src/Numerics/Numerics.csproj

@ -97,6 +97,7 @@
<Compile Include="IntegralTransforms\Fourier.cs" />
<Compile Include="IntegralTransforms\Hartley.cs" />
<Compile Include="Integration\GaussLegendreRule.cs" />
<Compile Include="Integration\GaussRule\GaussLegendrePoint.cs" />
<Compile Include="Integration\GaussRule\GaussLegendrePointFactory.cs" />
<Compile Include="Integration\GaussRule\GaussPoint.cs" />
<Compile Include="Interpolation\Barycentric.cs" />

38
src/UnitTests/IntegrationTests/IntegrationTest.cs

@ -251,26 +251,14 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
}
/// <summary>
/// Gauss-Legendre rule supports 2-dimensional integration over the rectangle.
/// Gauss-Legendre rule supports obtaining the ith abscissa/weight. 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 TestIntegrateGaussLegendre2D(int order)
{
}
/// <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)
public void TestGaussLegendreRuleGetAbscissaGetWeightOrderViaIntegration(int order)
{
GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
@ -284,10 +272,28 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
Assert.Less(relativeError, 5e-16);
}
/// <summary>
/// Gauss-Legendre rule supports obtaining array of abscissas/weights.
/// </summary>
[Test]
public void TestGaussLegendreRuleAbscissasWeightsViaIntegration()
{
const int order = 19;
GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
double[] abscissa = gaussLegendre.Abscissas;
double[] weight = gaussLegendre.Weights;
for (int i = 0; i < gaussLegendre.Order; i++)
{
Assert.AreEqual(gaussLegendre.GetAbscissa(i),abscissa[i]);
Assert.AreEqual(gaussLegendre.GetWeight(i), weight[i]);
}
}
/// <summary>
/// Gauss-Legendre rule supports obtaining IntervalBegin.
/// </summary>
[TestCase]
[Test]
public void TestGetGaussLegendreRuleIntervalBegin()
{
const int order = 19;
@ -298,7 +304,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
/// <summary>
/// Gauss-Legendre rule supports obtaining IntervalEnd.
/// </summary>
[TestCase]
[Test]
public void TestGaussLegendreRuleIntervalEnd()
{
const int order = 19;

Loading…
Cancel
Save