Skip to content

Commit bf5c6f2

Browse files
diluculodiluculo
committed
Integration: Add triple integral
1 parent 306fb06 commit bf5c6f2

File tree

3 files changed

+164
-3
lines changed

3 files changed

+164
-3
lines changed

‎src/Numerics.Tests/IntegrationTests/IntegrationTest.cs‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -111,6 +111,18 @@ private static double TargetFunctionG(double x)
111111
return Math.Sqrt(x) / Math.Sqrt(1 - x * x);
112112
}
113113

114+
/// <summary>
115+
/// Test Function: f(x, y, z) = y sin(x) + z cos(x)
116+
/// </summary>
117+
/// <param name="x">First input value.</param>
118+
/// <param name="y">Second input value.</param>
119+
/// <param name="z">Third input value.</param>
120+
/// <returns>Function result.</returns>
121+
private static double TargetFunctionH(double x, double y, double z)
122+
{
123+
return y * Math.Sin(x) + z * Math.Cos(x);
124+
}
125+
114126
/// <summary>
115127
/// Test Function Start point.
116128
/// </summary>
@@ -216,6 +228,41 @@ private static double TargetFunctionG(double x)
216228
/// </summary>
217229
private const double TargetAreaG = 1.1981402347355922074;
218230

231+
/// <summary>
232+
/// Target volume.
233+
/// </summary>
234+
private const double TargetVolumeH = 2.0;
235+
236+
/// <summary>
237+
/// Test Function Start point.
238+
/// </summary>
239+
private const double Xmin_H = 0;
240+
241+
/// <summary>
242+
/// Test Function end point.
243+
/// </summary>
244+
private const double Xmax_H = Math.PI;
245+
246+
/// <summary>
247+
/// Test Function Start point.
248+
/// </summary>
249+
private const double Ymin_H = 0;
250+
251+
/// <summary>
252+
/// Test Function end point.
253+
/// </summary>
254+
private const double Ymax_H = 1;
255+
256+
/// <summary>
257+
/// Test Function Start point.
258+
/// </summary>
259+
private const double Zmin_H = -1;
260+
261+
/// <summary>
262+
/// Test Function end point.
263+
/// </summary>
264+
private const double Zmax_H = 1;
265+
219266
/// <summary>
220267
/// Test Integrate facade for simple use cases.
221268
/// </summary>
@@ -383,6 +430,21 @@ public void TestIntegrateFacade()
383430
Integrate.GaussKronrod(TargetFunctionG, StartG, StopG, 1e-10, order: 15),
384431
1e-10,
385432
"GaussKronrod, sqrt(x)/sqrt(1 - x^2), order 15");
433+
434+
// TargetFunctionH
435+
// integral_(0)^(¥ð) integral_(0)^(1) integral_(-1)^(1) y sin(x) + z cos(x) dz dy dx = 2.0
436+
437+
Assert.AreEqual(
438+
TargetVolumeH,
439+
Integrate.OnCuboid(TargetFunctionH, Xmin_H, Xmax_H, Ymin_H, Ymax_H, Zmin_H, Zmax_H),
440+
1e-15,
441+
"Cuboid, order 32");
442+
443+
Assert.AreEqual(
444+
TargetVolumeH,
445+
Integrate.OnCuboid(TargetFunctionH, Xmin_H, Xmax_H, Ymin_H, Ymax_H, Zmin_H, Zmax_H, order: 22),
446+
1e-15,
447+
"Cuboid, Order 22");
386448
}
387449

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

579+
/// <summary>
580+
/// Gauss-Legendre rule supports 3-dimensional integration over the cuboid.
581+
/// </summary>
582+
/// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
583+
[TestCase(19)]
584+
[TestCase(20)]
585+
[TestCase(21)]
586+
[TestCase(22)]
587+
public void TestGaussLegendreRuleIntegrate3D(int order)
588+
{
589+
double appoximateVolume = GaussLegendreRule.Integrate(TargetFunctionH, Xmin_H, Xmax_H, Ymin_H, Ymax_H, Zmin_H, Zmax_H, order);
590+
double relativeError = Math.Abs(TargetVolumeH - appoximateVolume) / TargetVolumeH;
591+
Assert.Less(relativeError, 1e-15);
592+
}
593+
594+
[TestCase(19)]
595+
[TestCase(20)]
596+
[TestCase(21)]
597+
[TestCase(22)]
598+
public void TestIntegrateOverTetrahedron(int order)
599+
{
600+
// Variable change from (u, v, w) to (x, y, z):
601+
// x = u;
602+
// y = (1 - u)*v;
603+
// z = (1 - u)*(1 - v)*w;
604+
// Jacobian determinant of the transform
605+
// |J| = (1 - u)^2 (1 - v)
606+
// integrate[f, {x, 0, 1}, {y, 0, 1 - x}, {z, 0, 1 - x - y}]
607+
// = integrate[f |J|, {u, 0, 1}, {v, 0, 1}, {w, 0, 1}]
608+
609+
double J(double u, double v, double w) => (1.0 - u) * (1.0 - u) * (1.0 - v);
610+
611+
double f1(double x, double y, double z) => Math.Sqrt(x + y + z);
612+
double expected1 = 0.1428571428571428571428571; // 1/7
613+
Assert.AreEqual(
614+
expected1,
615+
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),
616+
1e-10,
617+
"Integral 3D of sqrt(x + y + z) on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");
618+
619+
double f2(double x, double y, double z) => Math.Pow(1 + x + y + z, -4);
620+
double expected2 = 0.02083333333333333333333333; // 1/48
621+
Assert.AreEqual(
622+
expected2,
623+
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),
624+
1e-15,
625+
"Integral 3D of (1 + x + y + z)^(-4) on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");
626+
627+
double f3(double x, double y, double z) => Math.Sin(x + 2 * y + 4 * z);
628+
double expected3 = 0.1319023268901816727730723; // 2/3 sin^4(1/2) (2 + 4 cos(1) + cos(2))
629+
Assert.AreEqual(
630+
expected3,
631+
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),
632+
1e-15,
633+
"Integral 3D of sin(x + 2 y + 4 z) on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");
634+
635+
double f4(double x, double y, double z) => x * x + y;
636+
double expected4 = 0.05833333333333333333333333; // 7/120
637+
Assert.AreEqual(
638+
expected4,
639+
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),
640+
1e-15,
641+
"Integral 3D of x^2 + y on [0, 1] x [0, 1 - x] x [0, 1 - x - y]");
642+
}
643+
517644
/// <summary>
518645
/// Gauss-Legendre rule supports obtaining the ith abscissa/weight. In this case, they're used for integration.
519646
/// </summary>

‎src/Numerics/Integrate.cs‎

Lines changed: 19 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ public static double OnClosedInterval(Func<double, double> f, double intervalBeg
7070
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
7171
/// <param name="invervalEndA">Where the interval ends for the first (inside) integral, exclusive and finite.</param>
7272
/// <param name="invervalBeginB">Where the interval starts for the second (outside) integral, exclusive and finite.</param>
73-
/// /// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
73+
/// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
7474
/// <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>
7575
/// <returns>Approximation of the finite integral in the given interval.</returns>
7676
public static double OnRectangle(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
@@ -85,13 +85,30 @@ public static double OnRectangle(Func<double, double, double> f, double inverval
8585
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
8686
/// <param name="invervalEndA">Where the interval ends for the first (inside) integral, exclusive and finite.</param>
8787
/// <param name="invervalBeginB">Where the interval starts for the second (outside) integral, exclusive and finite.</param>
88-
/// /// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
88+
/// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
8989
/// <returns>Approximation of the finite integral in the given interval.</returns>
9090
public static double OnRectangle(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB)
9191
{
9292
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, 32);
9393
}
9494

95+
/// <summary>
96+
/// 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].
97+
/// </summary>
98+
/// <param name="f">The 3-dimensional analytic smooth function to integrate.</param>
99+
/// <param name="invervalBeginA">Where the interval starts for the first integral, exclusive and finite.</param>
100+
/// <param name="invervalEndA">Where the interval ends for the first integral, exclusive and finite.</param>
101+
/// <param name="invervalBeginB">Where the interval starts for the second integral, exclusive and finite.</param>
102+
/// <param name="invervalEndB">Where the interval ends for the second integral, exclusive and finite.</param>
103+
/// <param name="invervalBeginC">Where the interval starts for the third integral, exclusive and finite.</param>
104+
/// <param name="invervalEndC">Where the interval ends for the third integral, exclusive and finite.</param>
105+
/// <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>
106+
/// <returns>Approximation of the finite integral in the given interval.</returns>
107+
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)
108+
{
109+
return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, invervalBeginC, invervalEndC, order);
110+
}
111+
95112
/// <summary>
96113
/// 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.
97114
/// </summary>

‎src/Numerics/Integration/GaussLegendreRule.cs‎

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -192,7 +192,7 @@ public static Complex ContourIntegrate(Func<double, Complex> f, double invervalB
192192
/// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
193193
/// <param name="invervalEndA">Where the interval ends for the first (inside) integral, exclusive and finite.</param>
194194
/// <param name="invervalBeginB">Where the interval starts for the second (outside) integral, exclusive and finite.</param>
195-
/// /// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
195+
/// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
196196
/// <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>
197197
/// <returns>Approximation of the finite integral in the given interval.</returns>
198198
public static double Integrate(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
@@ -258,5 +258,22 @@ public static double Integrate(Func<double, double, double> f, double invervalBe
258258

259259
return c*a*sum;
260260
}
261+
262+
/// <summary>
263+
/// Approximates a 3-dimensional definite integral using an Nth order Gauss-Legendre rule over the cuboid [a1,a2] x [b1,b2] x [c1,c2].
264+
/// </summary>
265+
/// <param name="f">The 3-dimensional analytic smooth function to integrate.</param>
266+
/// <param name="invervalBeginA">Where the interval starts for the first integral, exclusive and finite.</param>
267+
/// <param name="invervalEndA">Where the interval ends for the first integral, exclusive and finite.</param>
268+
/// <param name="invervalBeginB">Where the interval starts for the second integral, exclusive and finite.</param>
269+
/// <param name="invervalEndB">Where the interval ends for the second integral, exclusive and finite.</param>
270+
/// <param name="invervalBeginC">Where the interval starts for the third integral, exclusive and finite.</param>
271+
/// <param name="invervalEndC">Where the interval ends for the third integral, exclusive and finite.</param>
272+
/// <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>
273+
/// <returns>Approximation of the finite integral in the given intervals.</returns>
274+
public static double Integrate(Func<double, double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, double invervalBeginC, double invervalEndC, int order)
275+
{
276+
return Integrate((z) => Integrate((x, y) => f(x, y, z), invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order), invervalBeginC, invervalEndC, order);
277+
}
261278
}
262279
}

0 commit comments

Comments
 (0)