Skip to content

Commit

Permalink
Merge pull request mathnet#1013 from diluculo/integral3
Browse files Browse the repository at this point in the history 
Integration: Add triple integral
  • Loading branch information
@cdrnet
cdrnet authored Aug 1, 2023
2 parents 78fd5ad + bf5c6f2 commit ad8ef49
Show file tree
Hide file tree
Showing 3 changed files with 164 additions and 3 deletions.
127 changes: 127 additions & 0 deletions src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,18 @@ private static double TargetFunctionG(double x)
return Math.Sqrt(x) / Math.Sqrt(1 - x * x);
}

/// <summary>
/// Test Function: f(x, y, z) = y sin(x) + z cos(x)
/// </summary>
/// <param name="x">First input value.</param>
/// <param name="y">Second input value.</param>
/// <param name="z">Third input value.</param>
/// <returns>Function result.</returns>
private static double TargetFunctionH(double x, double y, double z)
{
return y * Math.Sin(x) + z * Math.Cos(x);
}

/// <summary>
/// Test Function Start point.
/// </summary>
Expand Down Expand Up @@ -216,6 +228,41 @@ private static double TargetFunctionG(double x)
/// </summary>
private const double TargetAreaG = 1.1981402347355922074;

/// <summary>
/// Target volume.
/// </summary>
private const double TargetVolumeH = 2.0;

/// <summary>
/// Test Function Start point.
/// </summary>
private const double Xmin_H = 0;

/// <summary>
/// Test Function end point.
/// </summary>
private const double Xmax_H = Math.PI;

/// <summary>
/// Test Function Start point.
/// </summary>
private const double Ymin_H = 0;

/// <summary>
/// Test Function end point.
/// </summary>
private const double Ymax_H = 1;

/// <summary>
/// Test Function Start point.
/// </summary>
private const double Zmin_H = -1;

/// <summary>
/// Test Function end point.
/// </summary>
private const double Zmax_H = 1;

/// <summary>
/// Test Integrate facade for simple use cases.
/// </summary>
Expand Down Expand Up @@ -383,6 +430,21 @@ public void TestIntegrateFacade()
Integrate.GaussKronrod(TargetFunctionG, StartG, StopG, 1e-10, order: 15),
1e-10,
"GaussKronrod, sqrt(x)/sqrt(1 - x^2), order 15");

// TargetFunctionH
// integral_(0)^(¥ð) integral_(0)^(1) integral_(-1)^(1) y sin(x) + z cos(x) dz dy dx = 2.0

Assert.AreEqual(
TargetVolumeH,
Integrate.OnCuboid(TargetFunctionH, Xmin_H, Xmax_H, Ymin_H, Ymax_H, Zmin_H, Zmax_H),
1e-15,
"Cuboid, order 32");

Assert.AreEqual(
TargetVolumeH,
Integrate.OnCuboid(TargetFunctionH, Xmin_H, Xmax_H, Ymin_H, Ymax_H, Zmin_H, Zmax_H, order: 22),
1e-15,
"Cuboid, Order 22");
}

/// <summary>
Expand Down Expand Up @@ -514,6 +576,71 @@ public void TestGaussLegendreRuleIntegrate2D(int order)
Assert.Less(relativeError, 1e-15);
}

/// <summary>
/// Gauss-Legendre rule supports 3-dimensional integration over the cuboid.
/// </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 TestGaussLegendreRuleIntegrate3D(int order)
{
double appoximateVolume = GaussLegendreRule.Integrate(TargetFunctionH, Xmin_H, Xmax_H, Ymin_H, Ymax_H, Zmin_H, Zmax_H, order);
double relativeError = Math.Abs(TargetVolumeH - appoximateVolume) / TargetVolumeH;
Assert.Less(relativeError, 1e-15);
}

[TestCase(19)]
[TestCase(20)]
[TestCase(21)]
[TestCase(22)]
public void TestIntegrateOverTetrahedron(int order)
{
// Variable change from (u, v, w) to (x, y, z):
// x = u;
// y = (1 - u)*v;
// z = (1 - u)*(1 - v)*w;
// Jacobian determinant of the transform
// |J| = (1 - u)^2 (1 - v)
// integrate[f, {x, 0, 1}, {y, 0, 1 - x}, {z, 0, 1 - x - y}]
// = integrate[f |J|, {u, 0, 1}, {v, 0, 1}, {w, 0, 1}]

double J(double u, double v, double w) => (1.0 - u) * (1.0 - u) * (1.0 - v);

double f1(double x, double y, double z) => Math.Sqrt(x + y + z);
double expected1 = 0.1428571428571428571428571; // 1/7
Assert.AreEqual(
expected1,
Integrate.OnCuboid((u, v, w) => f1(u, (1 - u) * v, (1 - u) * (1 - v) * w) * J(u, v, w), 0, 1, 0, 1, 0, 1, order: order),
1e-10,
"Integral 3D of sqrt(x + y + z) on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");

double f2(double x, double y, double z) => Math.Pow(1 + x + y + z, -4);
double expected2 = 0.02083333333333333333333333; // 1/48
Assert.AreEqual(
expected2,
Integrate.OnCuboid((u, v, w) => f2(u, (1 - u) * v, (1 - u) * (1 - v) * w) * J(u, v, w), 0, 1, 0, 1, 0, 1, order: order),
1e-15,
"Integral 3D of (1 + x + y + z)^(-4) on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");

double f3(double x, double y, double z) => Math.Sin(x + 2 * y + 4 * z);
double expected3 = 0.1319023268901816727730723; // 2/3 sin^4(1/2) (2 + 4 cos(1) + cos(2))
Assert.AreEqual(
expected3,
Integrate.OnCuboid((u, v, w) => f3(u, (1 - u) * v, (1 - u) * (1 - v) * w) * J(u, v, w), 0, 1, 0, 1, 0, 1, order: order),
1e-15,
"Integral 3D of sin(x + 2 y + 4 z) on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");

double f4(double x, double y, double z) => x * x + y;
double expected4 = 0.05833333333333333333333333; // 7/120
Assert.AreEqual(
expected4,
Integrate.OnCuboid((u, v, w) => f4(u, (1 - u) * v, (1 - u) * (1 - v) * w) * J(u, v, w), 0, 1, 0, 1, 0, 1, order: order),
1e-15,
"Integral 3D of x^2 + y on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");
}

/// <summary>
/// Gauss-Legendre rule supports obtaining the ith abscissa/weight. In this case, they're used for integration.
/// </summary>
Expand Down
21 changes: 19 additions & 2 deletions src/Numerics/Integrate.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ public static double OnClosedInterval(Func<double, double> f, double intervalBeg
/// <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="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)
Expand All @@ -85,13 +85,30 @@ public static double OnRectangle(Func<double, double, double> f, double inverval
/// <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="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>
/// Approximates a 3-dimensional definite integral using an Nth order Gauss-Legendre rule over the cuboid or rectangular prism [a1,a2] x [b1,b2] x [c1,c2].
/// </summary>
/// <param name="f">The 3-dimensional analytic smooth function to integrate.</param>
/// <param name="invervalBeginA">Where the interval starts for the first integral, exclusive and finite.</param>
/// <param name="invervalEndA">Where the interval ends for the first integral, exclusive and finite.</param>
/// <param name="invervalBeginB">Where the interval starts for the second integral, exclusive and finite.</param>
/// <param name="invervalEndB">Where the interval ends for the second integral, exclusive and finite.</param>
/// <param name="invervalBeginC">Where the interval starts for the third integral, exclusive and finite.</param>
/// <param name="invervalEndC">Where the interval ends for the third 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 OnCuboid(Func<double, double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, double invervalBeginC, double invervalEndC, int order = 32)
{
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, invervalBeginC, invervalEndC, order);
}

/// <summary>
/// Approximation of the definite integral of an analytic smooth function by double-exponential quadrature. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
/// </summary>
Expand Down
19 changes: 18 additions & 1 deletion src/Numerics/Integration/GaussLegendreRule.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,7 +192,7 @@ public static Complex ContourIntegrate(Func<double, Complex> f, double invervalB
/// <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="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 Integrate(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
Expand Down Expand Up @@ -258,5 +258,22 @@ public static double Integrate(Func<double, double, double> f, double invervalBe

return c*a*sum;
}

/// <summary>
/// Approximates a 3-dimensional definite integral using an Nth order Gauss-Legendre rule over the cuboid [a1,a2] x [b1,b2] x [c1,c2].
/// </summary>
/// <param name="f">The 3-dimensional analytic smooth function to integrate.</param>
/// <param name="invervalBeginA">Where the interval starts for the first integral, exclusive and finite.</param>
/// <param name="invervalEndA">Where the interval ends for the first integral, exclusive and finite.</param>
/// <param name="invervalBeginB">Where the interval starts for the second integral, exclusive and finite.</param>
/// <param name="invervalEndB">Where the interval ends for the second integral, exclusive and finite.</param>
/// <param name="invervalBeginC">Where the interval starts for the third integral, exclusive and finite.</param>
/// <param name="invervalEndC">Where the interval ends for the third 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 intervals.</returns>
public static double Integrate(Func<double, double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, double invervalBeginC, double invervalEndC, int order)
{
return Integrate((z) => Integrate((x, y) => f(x, y, z), invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order), invervalBeginC, invervalEndC, order);
}
}
}

0 comments on commit ad8ef49

Please sign in to comment.