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

Moved the caching to the factory, thereby making the rule stateless

pull/371/head
Larz White 10 years ago
parent
4f64aa52ab
commit
e60cb1af1b
2 changed files with 31 additions and 42 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. 59
      src/Numerics/Integration/GaussLegendreRule.cs
  2. 14
      src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs

59
src/Numerics/Integration/GaussLegendreRule.cs
View File

@ -39,8 +39,7 @@ namespace MathNet.Numerics.Integration
public class GaussLegendreRule
{
private readonly GaussPoint gaussLegendrePoint;
private static GaussPoint nonNegativeGaussLegendrePoint;
/// <summary>
/// Initializes a new instance of the <see cref="GaussLegendreRule"/> class.
/// </summary>
@ -99,22 +98,6 @@ namespace MathNet.Numerics.Integration
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>
@ -158,7 +141,7 @@ namespace MathNet.Numerics.Integration
/// <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);
GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.GetGaussPoint(order);
double sum, ax;
int i;
@ -169,11 +152,11 @@ namespace MathNet.Numerics.Integration
if (order.IsOdd())
{
sum = nonNegativeGaussLegendrePoint.Weights[0]*f(b);
sum = gaussLegendrePoint.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));
ax = a*gaussLegendrePoint.Abscissas[i];
sum += gaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
}
}
else
@ -181,8 +164,8 @@ namespace MathNet.Numerics.Integration
sum = 0.0;
for (i = 0; i < m; i++)
{
ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
sum += nonNegativeGaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
ax = a*gaussLegendrePoint.Abscissas[i];
sum += gaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
}
}
@ -201,7 +184,7 @@ namespace MathNet.Numerics.Integration
/// <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);
GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.GetGaussPoint(order);
double ax, cy, sum;
int i, j;
@ -214,32 +197,32 @@ namespace MathNet.Numerics.Integration
if (order.IsOdd())
{
sum = nonNegativeGaussLegendrePoint.Weights[0]*nonNegativeGaussLegendrePoint.Weights[0]*f(b, d);
sum = gaussLegendrePoint.Weights[0]*gaussLegendrePoint.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));
cy = c*gaussLegendrePoint.Abscissas[j];
t += gaussLegendrePoint.Weights[j]*(f(b, d + cy) + f(b, d - cy));
}
sum += nonNegativeGaussLegendrePoint.Weights[0]*t;
sum += gaussLegendrePoint.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));
ax = a*gaussLegendrePoint.Abscissas[i];
t += gaussLegendrePoint.Weights[i]*(f(b + ax, d) + f(b - ax, d));
}
sum += nonNegativeGaussLegendrePoint.Weights[0]*t;
sum += gaussLegendrePoint.Weights[0]*t;
for (i = 1; i < m; i++)
{
ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
ax = a*gaussLegendrePoint.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));
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));
}
}
}
@ -248,11 +231,11 @@ namespace MathNet.Numerics.Integration
sum = 0.0;
for (i = 0; i < m; i++)
{
ax = a*nonNegativeGaussLegendrePoint.Abscissas[i];
ax = a*gaussLegendrePoint.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));
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));
}
}
}

14
src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs
View File

@ -38,19 +38,25 @@ namespace MathNet.Numerics.Integration.GaussRule
/// </summary>
static class GaussLegendrePointFactory
{
private static GaussPoint gaussLegendrePoint;
/// <summary>
/// Getter for the GaussPoint.
/// </summary>
/// <param name="order">Defines 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.</param>
/// <param name="order">Defines 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. Furthermore, caching in a private static field is used to speed up the algorithm.</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>
public static GaussPoint GetGaussPoint(int order)
{
GaussPoint gaussLegendrePoint;
// Try to get the GaussPoint from the cached static field.
if (gaussLegendrePoint != null && gaussLegendrePoint.Order == order)
{
return gaussLegendrePoint;
}
// Try to find it in the precomputed dictionary first.
// Try to find the GaussPoint in the precomputed dictionary.
if (!GaussLegendrePoints.TryGetValue(order, out gaussLegendrePoint))
{
gaussLegendrePoint = GenerateGaussLegendrePoint(order, 1e-10); // Generate the abscissas/weights on the fly.
gaussLegendrePoint = GenerateGaussLegendrePoint(order, 1e-10); // Generate the GaussPoint on the fly.
}
return gaussLegendrePoint;

Write Preview
Loading…
Cancel
Save