Browse Source

Cleaned up Gauss Legendre code.

pull/423/head
Larz White 10 years ago
parent
commit
d475314597
  1. 65
      src/Numerics/Integration/GaussLegendreRule.cs
  2. 187
      src/Numerics/Integration/GaussRule/GaussLegendrePoint.cs
  3. 174
      src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs
  4. 16
      src/Numerics/Integration/GaussRule/GaussPoint.cs
  5. 1
      src/Numerics/Numerics.csproj
  6. 38
      src/UnitTests/IntegrationTests/IntegrationTest.cs

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