tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Add StieltjesPolynomial to support other than precomputed orders in Gauss-Kronrod rule.

pull/655/head
diluculo 7 years ago
parent
3cd5bd80e3
commit
0f501660a6
6 changed files with 257 additions and 58 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 17
      src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
  2. 14
      src/Numerics/Integration/GaussKronrodRule.cs
  3. 201
      src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs
  4. 8
      src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs
  5. 40
      src/Numerics/Integration/GaussRule/GaussPointPair.cs
  6. 35
      src/Numerics/Integration/GaussRule/GaussPoints.cs

17
src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
View File

@ -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);
}
/// <summary>
/// Gauss-Kronrod rule supports integration.
/// </summary>
/// <param name="order">Defines an Nth order Gauss-Kronrod rule. The order also defines the number of abscissas and weights for the rule.</param>
[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)

14
src/Numerics/Integration/GaussKronrodRule.cs
View File

@ -43,7 +43,7 @@ namespace MathNet.Numerics.Integration
{
public class GaussKronrodRule
{
private readonly GaussPoints gaussKronrodPoint;
private readonly GaussPointPair gaussKronrodPoint;
/// <summary>
/// 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<double, double> f, out double error, out double pL1, GaussPoints gaussKronrodPoint)
private static double integrate_non_adaptive_m1_1(Func<double, double> 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<double, Complex> f, out double error, out double pL1, GaussPoints gaussKronrodPoint)
private static Complex contour_integrate_non_adaptive_m1_1(Func<double, Complex> 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<double, double> 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<double, double> 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<double, Complex> 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<double, Complex> 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;

201
src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs
View File

@ -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
/// <summary>
/// Precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61.
/// </summary>
internal static readonly Dictionary<int, GaussPoints> PreComputed = new Dictionary<int, GaussPoints>
internal static readonly Dictionary<int, GaussPointPair> PreComputed = new Dictionary<int, GaussPointPair>
{
{ 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.
/// </summary>
/// <param name="order">Defines an Nth order Gauss-Kronrod rule. The order also defines the number of abscissas and weights for the rule.</param>
/// <param name="targetAbsoluteTolerance">Required precision to compute the abscissas/weights.</param>
/// <returns>Object containing the non-negative abscissas/weights, order.</returns>
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);
}
/// <summary>
/// Returns the Stieltjes Polynomial of order.
/// </summary>
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;
}
/// <summary>
/// Returns the Legendre polynomial of order.
/// </summary>
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;
}
}
}

8
src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs
View File

@ -8,15 +8,14 @@ namespace MathNet.Numerics.Integration.GaussRule
internal static class GaussKronrodPointFactory
{
[ThreadStatic]
private static GaussPoints gaussKronrodPoint;
private static GaussPointPair gaussKronrodPoint;
/// <summary>
/// Getter for the GaussKronrodPoint.
/// </summary>
/// <param name="order">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.</param>
/// <param name="targetAbsoluteTolerance">Required precision to compute the abscissas/weights. 1e-10 is usually fine.</param>
/// <returns>Object containing the non-negative abscissas/weights, and order.</returns>
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);
}
}

40
src/Numerics/Integration/GaussRule/GaussPointPair.cs
View File

@ -0,0 +1,40 @@
namespace MathNet.Numerics.Integration.GaussRule
{
/// <summary>
/// Contains two GaussPoint.
/// </summary>
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)
{ }
}
}

35
src/Numerics/Integration/GaussRule/GaussPoints.cs
View File

@ -1,35 +0,0 @@
namespace MathNet.Numerics.Integration.GaussRule
{
/// <summary>
/// Contains two set of abscissas, weights, and order.
/// </summary>
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)
{ }
}
}
Write Preview
Loading…
Cancel
Save