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

Add Gauss-Legendre to Integrate class

pull/655/head
diluculo 7 years ago
parent
cca30a3a0e
commit
9cb8bbea5f
3 changed files with 230 additions and 10 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. 41
      src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
  2. 156
      src/Numerics/Integrate.cs
  3. 43
      src/Numerics/Integration/GaussLegendreRule.cs

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

@ -430,13 +430,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
expected,
Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b),
1e-10,
"DET Integral e^(-x^2 /2) from {0} to {1}", a, b);
"DET Integral of e^(-x^2 /2) from {0} to {1}", a, b);
Assert.AreEqual(
expected,
Integrate.GaussKronrod((x) => Math.Exp(-x * x / 2), a, b),
1e-10,
"GK Integral e^(-x^2 /2) from {0} to {1}", a, b);
"GK Integral of e^(-x^2 /2) from {0} to {1}", a, b);
Assert.AreEqual(
expected,
Integrate.GausLegendre((x) => Math.Exp(-x * x / 2), a, b, order: 128),
1e-10,
"GL Integral of e^(-x^2 /2) from {0} to {1}", a, b);
}
// integral_(-oo)^(oo) sin(pi x) / (pi x) dx = 1 / pi integral_(-oo)^(oo) sin(x) / x dx
@ -453,13 +459,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
expected,
factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b),
1e-10,
"DET Integral sin(pi*x)/(pi*x) from -oo to oo");
"DET Integral of sin(pi*x)/(pi*x) from -oo to oo");
Assert.AreEqual(
expected,
factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b),
1e-10,
"GK Integral sin(pi*x)/(pi*x) from -oo to oo");
"GK Integral of sin(pi*x)/(pi*x) from -oo to oo");
Assert.AreEqual(
expected,
factor * Integrate.GausLegendre((x) => 1 / (1 + x * x), a, b, order: 128),
1e-10,
"GL Integral of sin(pi*x)/(pi*x) from -oo to oo");
}
// integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð
@ -473,30 +485,43 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
var expected = new Complex(r, i);
var actualDET = ContourIntegrate.DoubleExponential((x) => 1 / new Complex(1, x * x), a, b);
var actualGK = ContourIntegrate.GaussKronrod((x) => 1 / new Complex(1, x * x), a, b);
var actualGL = ContourIntegrate.GaussLegendre((x) => 1 / new Complex(1, x * x), a, b, order: 128);
Assert.AreEqual(
expected.Real,
actualDET.Real,
1e-10,
"DET Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
"DET Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
Assert.AreEqual(
expected.Imaginary,
actualDET.Imaginary,
1e-10,
"DET Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
"DET Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
Assert.AreEqual(
expected.Real,
actualGK.Real,
1e-10,
"GK Integral Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
"GK Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
Assert.AreEqual(
expected.Imaginary,
actualGK.Imaginary,
1e-10,
"GK Integral Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
"GK Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
Assert.AreEqual(
expected.Real,
actualGL.Real,
1e-10,
"GL Integral of Re[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
Assert.AreEqual(
expected.Imaginary,
actualGL.Imaginary,
1e-10,
"GL Integral of Im[e^(-x^2 /2) / (1 + j e^x)] from {0} to {1}", a, b);
}
}
}

156
src/Numerics/Integrate.cs
View File

@ -156,6 +156,82 @@ namespace MathNet.Numerics
}
}
/// <summary>
/// Approximation of the definite integral of an analytic smooth function by Gauss-Legendre quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>
/// <param name="f">The analytic smooth function to integrate.</param>
/// <param name="intervalBegin">Where the interval starts.</param>
/// <param name="intervalEnd">Where the interval stops.</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 calculated on the fly.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static double GausLegendre(Func<double, double> f, double intervalBegin, double intervalEnd, int order = 128)
{
// Reference:
// Formula used for variable subsitution from
// 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140.
// 2. quadgk.m, GNU Octave
if (intervalBegin > intervalEnd)
{
return -GaussLegendreRule.Integrate(f, intervalEnd, intervalBegin, order);
}
// (-oo, oo) => [-1, 1]
//
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
{
Func<double, double> u = (t) =>
{
return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
};
return GaussLegendreRule.Integrate(u, -1, 1, order);
}
// [a, oo) => [0, 1]
//
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt
// = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd))
{
Func<double, double> u = (s) =>
{
return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
};
return GaussLegendreRule.Integrate(u, 0, 1, order);
}
// (-oo, b] => [-1, 0]
//
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt
// = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin))
{
Func<double, double> u = (s) =>
{
return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
};
return GaussLegendreRule.Integrate(u, -1, 0, order);
}
// [a, b] => [-1, 1]
//
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt
// g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else
{
Func<double, double> u = (t) =>
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
return GaussLegendreRule.Integrate(u, -1, 1, order);
}
}
/// <summary>
/// Approximation of the definite integral of an analytic smooth function by Gauss-Kronrod quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>
@ -163,8 +239,8 @@ namespace MathNet.Numerics
/// <param name="intervalBegin">Where the interval starts.</param>
/// <param name="intervalEnd">Where the interval stops.</param>
/// <param name="targetRelativeError">The expected relative accuracy of the approximation.</param>
/// <param name="maximumDepth">The maximum number of interval splittings permitted before stopping</param>
/// <param name="order">The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points</param>
/// <param name="maximumDepth">The maximum number of interval splittings permitted before stopping.</param>
/// <param name="order">The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static double GaussKronrod(Func<double, double> f, double intervalBegin, double intervalEnd, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15)
{
@ -270,6 +346,82 @@ namespace MathNet.Numerics
}
}
/// <summary>
/// Approximation of the definite integral of an analytic smooth complex function by double-exponential quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, defined on the real domain.</param>
/// <param name="intervalBegin">Where the interval starts.</param>
/// <param name="intervalEnd">Where the interval stops.</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 calculated on the fly.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex GaussLegendre(Func<double, Complex> f, double intervalBegin, double intervalEnd, int order = 128)
{
// Reference:
// Formula used for variable subsitution from
// 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140.
// 2. quadgk.m, GNU Octave
if (intervalBegin > intervalEnd)
{
return -GaussLegendre(f, intervalEnd, intervalBegin, order);
}
// (-oo, oo) => [-1, 1]
//
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
{
Func<double, Complex> u = (t) =>
{
return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
};
return GaussLegendreRule.ContourIntegrate(u, -1, 1, order);
}
// [a, oo) => [0, 1]
//
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt
// = integral_{0}^{1} f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd))
{
Func<double, Complex> u = (s) =>
{
return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
};
return GaussLegendreRule.ContourIntegrate(u, 0, 1, order);
}
// (-oo, b] => [-1, 0]
//
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt
// = -integral_{-1}^{0} f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin))
{
Func<double, Complex> u = (s) =>
{
return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
};
return GaussLegendreRule.ContourIntegrate(u, -1, 0, order);
}
// [a, b] => [-1, 1]
//
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt
// g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else
{
Func<double, Complex> u = (t) =>
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
return GaussLegendreRule.ContourIntegrate(u, -1, 1, order);
}
}
/// <summary>
/// Approximation of the definite integral of an analytic smooth function by Gauss-Kronrod quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>

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

@ -28,6 +28,7 @@
// </copyright>
using System;
using System.Numerics;
using MathNet.Numerics.Integration.GaussRule;
namespace MathNet.Numerics.Integration
@ -166,6 +167,48 @@ namespace MathNet.Numerics.Integration
return a*sum;
}
/// <summary>
/// Approximates a definite integral using an Nth order Gauss-Legendre rule.
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, defined on the real domain.</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 calculated on the fly.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex ContourIntegrate(Func<double, Complex> f, double invervalBegin, double invervalEnd, int order)
{
GaussPoint gaussLegendrePoint = GaussLegendrePointFactory.GetGaussPoint(order);
Complex sum;
double 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>

Write Preview
Loading…
Cancel
Save