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

Add complex version of double-exponential integral

pull/655/head
diluculo 7 years ago
parent
92ae8e358a
commit
067ec6fed3
4 changed files with 392 additions and 34 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. 70
      src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
  2. 143
      src/Numerics/Integrate.cs
  3. 20
      src/Numerics/Integration/DoubleExponentialTransformation.cs
  4. 193
      src/Numerics/Integration/NewtonCotesTrapeziumRule.cs

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

@ -27,9 +27,10 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
using MathNet.Numerics.Integration;
using NUnit.Framework;
using System;
using System.Numerics;
namespace MathNet.Numerics.UnitTests.IntegrationTests
{
@ -59,7 +60,17 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
{
return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y));
}
/// <summary>
/// Test Function: f(x,y) = 1 / (1 + x^2)
/// </summary>
/// <param name="x">First input value.</param>
/// <returns>Function result.</returns>
private static double TargetFunctionC(double x)
{
return 1 / (1 + x * x);
}
/// <summary>
/// Test Function Start point.
/// </summary>
@ -80,6 +91,16 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
/// </summary>
private const double StopB = 1;
/// <summary>
/// Test Function Start point.
/// </summary>
private const double StartC = double.NegativeInfinity;
/// <summary>
/// Test Function Stop point.
/// </summary>
private const double StopC = double.PositiveInfinity;
/// <summary>
/// Target area square.
/// </summary>
@ -90,6 +111,11 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
/// </summary>
private const double TargetAreaB = 11.7078776759298776163;
/// <summary>
/// Target area.
/// </summary>
private const double TargetAreaC = Constants.Pi;
/// <summary>
/// Test Integrate facade for simple use cases.
/// </summary>
@ -119,6 +145,18 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
TargetAreaB,
1e-10,
"Rectangle, Gauss-Legendre Order 22");
Assert.AreEqual(
TargetAreaC,
Integrate.DoubleExponential(TargetFunctionC, StartC, StopC),
1e-5,
"Integral by substitution");
Assert.AreEqual(
TargetAreaC,
Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10),
1e-10,
"Integral by substitution, Target 1e-10");
}
/// <summary>
@ -326,7 +364,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
{
Assert.AreEqual(
expected,
Integrate.OnOpenInterval((x) => Math.Exp(-x * x / 2), a, b),
Integrate.DoubleExponential((x) => Math.Exp(-x * x / 2), a, b),
1e-10,
"Integral e^(-x^2 /2) from {0} to {1}", a, b);
}
@ -343,9 +381,33 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
{
Assert.AreEqual(
expected,
factor * Integrate.OnOpenInterval((x) => 1 / (1 + x * x), a, b),
factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b),
1e-10,
"Integral sin(pi*x)/(pi*x) from -oo to oo");
}
// integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð
[TestCase(double.NegativeInfinity, double.PositiveInfinity, 2.2214414690791831235, -2.2214414690791831235)]
// integral_(0)^(oo) 1/(1 + j x^2) dx = -1/2 (-1)^(3/4) ¥ð
[TestCase(0, double.PositiveInfinity, 1.1107207345395915618, -1.1107207345395915618)]
// integral_(-oo)^(0) 1/(1 + j x^2) dx = -1/2 (-1)^(3/4) ¥ð
[TestCase(double.NegativeInfinity, 0, 1.1107207345395915618, -1.1107207345395915618)]
public void TestContourIntegralBySubstitution(double a, double b, double r, double i)
{
var expected = new Complex(r, i);
var actual = ContourIntegrate.DoubleExponential((x) => 1 / new Complex(1, x * x), a, b);
Assert.AreEqual(
expected.Real,
actual.Real,
1e-10,
"Integral e^(-x^2 /2) / (1 + j e^x) from {0} to {1}", a, b);
Assert.AreEqual(
expected.Imaginary,
actual.Imaginary,
1e-10,
"Integral e^(-x^2 /2) / (1 + j e^x) from {0} to {1}", a, b);
}
}
}

143
src/Numerics/Integrate.cs
View File

@ -28,6 +28,7 @@
// </copyright>
using System;
using System.Numerics;
using MathNet.Numerics.Integration;
namespace MathNet.Numerics
@ -50,15 +51,44 @@ namespace MathNet.Numerics
return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError);
}
/// <summary>
/// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d].
/// </summary>
/// <param name="f">The 2-dimensional analytic smooth function to integrate.</param>
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
/// <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 calculated 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)
{
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order);
}
/// <summary>
/// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d].
/// </summary>
/// <param name="f">The 2-dimensional analytic smooth function to integrate.</param>
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
/// <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>
/// <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)
{
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, 32);
}
/// <summary>
/// Approximation of the definite integral of an analytic smooth function by substitution. 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, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="intervalBegin">Where the interval starts.</param>
/// <param name="intervalEnd">Where the interval stops.</param>
/// <param name="targetAbsoluteError">The expected relative accuracy of the approximation.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static double OnOpenInterval(Func<double, double> f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
public static double DoubleExponential(Func<double, double> f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
{
// Reference:
// Formula used for variable subsitution from
@ -67,7 +97,7 @@ namespace MathNet.Numerics
if (intervalBegin > intervalEnd)
{
return -OnOpenInterval(f, intervalEnd, intervalBegin, targetAbsoluteError);
return -DoubleExponential(f, intervalEnd, intervalBegin, targetAbsoluteError);
}
// (-oo, oo) => [-1, 1]
@ -81,7 +111,7 @@ namespace MathNet.Numerics
{
return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
};
return OnClosedInterval(u, -1d, 1d, targetAbsoluteError);
return DoubleExponentialTransformation.Integrate(u, -1, 1, targetAbsoluteError);
}
// [a, oo) => [0, 1]
//
@ -95,7 +125,7 @@ namespace MathNet.Numerics
{
return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
};
return OnClosedInterval(u, 0d, 1d, targetAbsoluteError);
return DoubleExponentialTransformation.Integrate(u, 0, 1, targetAbsoluteError);
}
// (-oo, b] => [-1, 0]
//
@ -109,7 +139,7 @@ namespace MathNet.Numerics
{
return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
};
return OnClosedInterval(u, -1d, 0d, targetAbsoluteError);
return DoubleExponentialTransformation.Integrate(u, -1, 0, targetAbsoluteError);
}
// [a, b] => [-1, 1]
//
@ -122,37 +152,90 @@ namespace MathNet.Numerics
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
return OnClosedInterval(u, -1d, 1d, targetAbsoluteError);
return DoubleExponentialTransformation.Integrate(u, -1, 1, targetAbsoluteError);
}
}
}
/// <summary>
/// Numerical Contour Integration over a real variable, of a complex-valued function.
/// </summary>
public static class ContourIntegrate
{
/// <summary>
/// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d].
/// Approximation of the definite integral of an analytic smooth complex function by substitution. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>
/// <param name="f">The 2-dimensional analytic smooth function to integrate.</param>
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
/// <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 calculated on the fly.</param>
/// <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="targetAbsoluteError">The expected relative accuracy of the approximation.</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)
public static Complex DoubleExponential(Func<double, Complex> f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
{
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order);
}
// 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
/// <summary>
/// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d].
/// </summary>
/// <param name="f">The 2-dimensional analytic smooth function to integrate.</param>
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
/// <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>
/// <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)
{
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, 32);
if (intervalBegin > intervalEnd)
{
return -DoubleExponential(f, intervalEnd, intervalBegin, targetAbsoluteError);
}
// (-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 DoubleExponentialTransformation.ContourIntegrate(u, -1, 1, targetAbsoluteError);
}
// [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 DoubleExponentialTransformation.ContourIntegrate(u, 0, 1, targetAbsoluteError);
}
// (-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 DoubleExponentialTransformation.ContourIntegrate(u, -1, 0, targetAbsoluteError);
}
// [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 DoubleExponentialTransformation.ContourIntegrate(u, -1, 1, targetAbsoluteError);
}
}
}
}

20
src/Numerics/Integration/DoubleExponentialTransformation.cs
View File

@ -29,6 +29,7 @@
using System;
using System.Linq;
using System.Numerics;
namespace MathNet.Numerics.Integration
{
@ -63,6 +64,25 @@ namespace MathNet.Numerics.Integration
targetRelativeError);
}
/// <summary>
/// Approximate the integral by the double exponential transformation
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, defined on the real domain.</param>
/// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="targetRelativeError">The expected relative accuracy of the approximation.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex ContourIntegrate(Func<double, Complex> f, double intervalBegin, double intervalEnd, double targetRelativeError)
{
return NewtonCotesTrapeziumRule.ContourIntegrateAdaptiveTransformedOdd(
f,
intervalBegin, intervalEnd,
Enumerable.Range(0, NumberOfMaximumLevels).Select(EvaluateAbcissas),
Enumerable.Range(0, NumberOfMaximumLevels).Select(EvaluateWeights),
1.0,
targetRelativeError);
}
/// <summary>
/// Compute the abscissa vector for a single level.
/// </summary>

193
src/Numerics/Integration/NewtonCotesTrapeziumRule.cs
View File

@ -29,6 +29,7 @@
using System;
using System.Collections.Generic;
using System.Numerics;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Integration
@ -58,6 +59,23 @@ namespace MathNet.Numerics.Integration
return (intervalEnd - intervalBegin)/2*(f(intervalBegin) + f(intervalEnd));
}
/// <summary>
/// Direct 2-point approximation of the definite integral in the provided interval by the trapezium rule.
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, defined on real domain.</param>
/// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex ContourIntegrateTwoPoint(Func<double, Complex> f, double intervalBegin, double intervalEnd)
{
if (f == null)
{
throw new ArgumentNullException(nameof(f));
}
return (intervalEnd - intervalBegin) / 2 * (f(intervalBegin) + f(intervalEnd));
}
/// <summary>
/// Composite N-point approximation of the definite integral in the provided interval by the trapezium rule.
/// </summary>
@ -92,6 +110,40 @@ namespace MathNet.Numerics.Integration
return step*sum;
}
/// <summary>
/// Composite N-point approximation of the definite integral in the provided interval by the trapezium rule.
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, defined on real domain.</param>
/// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="numberOfPartitions">Number of composite subdivision partitions.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex ContourIntegrateComposite(Func<double, Complex> f, double intervalBegin, double intervalEnd, int numberOfPartitions)
{
if (f == null)
{
throw new ArgumentNullException(nameof(f));
}
if (numberOfPartitions <= 0)
{
throw new ArgumentOutOfRangeException(nameof(numberOfPartitions), Resources.ArgumentPositive);
}
double step = (intervalEnd - intervalBegin) / numberOfPartitions;
double offset = step;
Complex sum = 0.5 * (f(intervalBegin) + f(intervalEnd));
for (int i = 0; i < numberOfPartitions - 1; i++)
{
// NOTE (ruegg, 2009-01-07): Do not combine intervalBegin and offset (numerical stability!)
sum += f(intervalBegin + offset);
offset += step;
}
return step * sum;
}
/// <summary>
/// Adaptive approximation of the definite integral in the provided interval by the trapezium rule.
/// </summary>
@ -132,6 +184,47 @@ namespace MathNet.Numerics.Integration
return sum;
}
/// <summary>
/// Adaptive approximation of the definite integral in the provided interval by the trapezium rule.
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, define don real domain.</param>
/// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="targetError">The expected accuracy of the approximation.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex ContourIntegrateAdaptive(Func<double, Complex> f, double intervalBegin, double intervalEnd, double targetError)
{
if (f == null)
{
throw new ArgumentNullException(nameof(f));
}
int numberOfPartitions = 1;
double step = intervalEnd - intervalBegin;
Complex sum = 0.5 * step * (f(intervalBegin) + f(intervalEnd));
for (int k = 0; k < 20; k++)
{
Complex midpointsum = 0;
for (int i = 0; i < numberOfPartitions; i++)
{
midpointsum += f(intervalBegin + ((i + 0.5) * step));
}
midpointsum *= step;
sum = 0.5 * (sum + midpointsum);
step *= 0.5;
numberOfPartitions *= 2;
if (sum.AlmostEqualRelative(midpointsum, targetError))
{
break;
}
}
return sum;
}
/// <summary>
/// Adaptive approximation of the definite integral by the trapezium rule.
/// </summary>
@ -230,5 +323,105 @@ namespace MathNet.Numerics.Integration
return sum*linearSlope;
}
}
/// <summary>
/// Adaptive approximation of the definite integral by the trapezium rule.
/// </summary>
/// <param name="f">The analytic smooth complex function to integrate, defined on the real domain.</param>
/// <param name="intervalBegin">Where the interval starts, inclusive and finite.</param>
/// <param name="intervalEnd">Where the interval stops, inclusive and finite.</param>
/// <param name="levelAbscissas">Abscissa vector per level provider.</param>
/// <param name="levelWeights">Weight vector per level provider.</param>
/// <param name="levelOneStep">First Level Step</param>
/// <param name="targetRelativeError">The expected relative accuracy of the approximation.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static Complex ContourIntegrateAdaptiveTransformedOdd(
Func<double, Complex> f,
double intervalBegin, double intervalEnd,
IEnumerable<double[]> levelAbscissas, IEnumerable<double[]> levelWeights,
double levelOneStep, double targetRelativeError)
{
if (f == null)
{
throw new ArgumentNullException(nameof(f));
}
if (levelAbscissas == null)
{
throw new ArgumentNullException(nameof(levelAbscissas));
}
if (levelWeights == null)
{
throw new ArgumentNullException(nameof(levelWeights));
}
double linearSlope = 0.5 * (intervalEnd - intervalBegin);
double linearOffset = 0.5 * (intervalEnd + intervalBegin);
targetRelativeError /= 5 * linearSlope;
using (var abcissasIterator = levelAbscissas.GetEnumerator())
using (var weightsIterator = levelWeights.GetEnumerator())
{
double step = levelOneStep;
// First Level
abcissasIterator.MoveNext();
weightsIterator.MoveNext();
double[] abcissasL1 = abcissasIterator.Current;
double[] weightsL1 = weightsIterator.Current;
Complex sum = f(linearOffset) * weightsL1[0];
for (int i = 1; i < abcissasL1.Length; i++)
{
sum += weightsL1[i] * (f((linearSlope * abcissasL1[i]) + linearOffset) + f(-(linearSlope * abcissasL1[i]) + linearOffset));
}
sum *= step;
// Additional Levels
double previousDelta = double.MaxValue;
for (int level = 1; abcissasIterator.MoveNext() && weightsIterator.MoveNext(); level++)
{
double[] abcissas = abcissasIterator.Current;
double[] weights = weightsIterator.Current;
Complex midpointsum = 0;
for (int i = 0; i < abcissas.Length; i++)
{
midpointsum += weights[i] * (f((linearSlope * abcissas[i]) + linearOffset) + f(-(linearSlope * abcissas[i]) + linearOffset));
}
midpointsum *= step;
sum = 0.5 * (sum + midpointsum);
step *= 0.5;
double delta = Complex.Abs(sum - midpointsum);
if (level == 1)
{
previousDelta = delta;
continue;
}
double r = Math.Log(delta) / Math.Log(previousDelta);
previousDelta = delta;
if (r > 1.9 && r < 2.1)
{
// convergence region
delta = Math.Sqrt(delta);
}
if (sum.Real.AlmostEqualNormRelative(midpointsum.Real, delta, targetRelativeError)
&& sum.Imaginary.AlmostEqualNormRelative(midpointsum.Imaginary, delta, targetRelativeError))
{
break;
}
}
return sum * linearSlope;
}
}
}
}

Write Preview
Loading…
Cancel
Save