# Computing surface integral on the unit sphere

I have a function $$f$$ which takes three unit vectors of $$\mathbb R^3$$ and returns a number, so $$f \colon \mathbb R^3 \times \mathbb R^3 \times \mathbb R^3 \to \mathbb R$$ which is defined as $$f(\mathbf x,\mathbf y, \mathbf z) := |\det(\mathbf x, \mathbf y, \mathbf z)|$$ i.e. determinant of the 3 by 3 matrix whose columns are the three coordinates (wrt canonical base) of the three vectors.

I would like to compute the surface integral $$\int_{S^2 \times S^2 \times S^2} f(\mathbf x, \mathbf y, \mathbf z) \, d\sigma(x) d\sigma(y) d\sigma(z)$$ where $$S^2 := \{(a,b,c) \in \mathbb R^3: a^2 + b^2 + c^2 = 1\}$$ is the unit sphere and $$\sigma$$ is the surface measure over it.

I have used

region = ImplicitRegion[x^2 + y^2 + z^2 == 1, {x,y,z}];
Integrate[Abs[Det[{{a, d, g}, {b,e,h}, {c, f, i}}]], {a,b,c} ∈ region, {d,e,f} ∈ region,{g,h,i} ∈ region]


where, of course, "function" is replaced by the above definition ($$a,b,c$$ are the coordinates of $$x$$ and so on). Mathematica works for quite a lot and returns... 0.

This is clearly meaningless, as the function is non-negative (but not identically vanishing). Maybe in this way I am computing a volume integral (on a region of measure 0)? What is the correct way to write this? Thanks.

EDIT: I have tried also spherical coordinates: I have assumed (rotational invariance) that $$\mathbf = (0,0,1)$$. Hence I gave to compute

Integrate[Abs[Det[( {
{Sin[a] Cos[b], Sin[c] Cos[d], 0},
{Sin[a] Sin[b], Sin[c] Sin[d], 0},
{Cos[a], Cos[c], 1}
} )]]*Sin[a]*Sin[c], {a, 0, 2*Pi}, {c, 0, 2*Pi}, {b, 0, Pi}, {d,
0, Pi}, Assumptions -> a ∈ Reals && b ∈ Reals && c
∈ Reals && d ∈ Reals]


returning 0 after a while. I have also tried the following:

Block[{t, a},
F = {t, a} \[Function] {Sin[t] Cos[a], Sin[t] Sin[a], Cos[t]};
DF = {t, a} \[Function] Evaluate[D[F[t, a], {{t, a}, 1}]];
jacobidet = {t, a} \[Function]
Evaluate[Simplify[Sqrt[Det[Transpose[DF[t, a]].DF[t, a]]]]];];
Integrate[Abs[Det[( {
{Sin[t] Cos[a], d, g},
{Sin[t] Sin[a], e, h},
{Cos[t], f, i}
} )]]* jacobidet[t, a], {a, 0, 2*Pi}, {t, 0, Pi}]


hoping that I could then integrate again the resultign expression in the other variables. But also this did not work - as it could be expected actually.

I am stuck...

• Just to be clear, are you trying to compute the portion of the surface of the unit sphere between the 3 points that you're feeding into the function? (i.e. your answer is proportional to N/4Pi) – Matt Stein Mar 1 '19 at 15:39
• @MattStein Not really. I am trying to compute the (average of the) volume (without sign) of the tetrahedron spanned by the 3 points feeding into the function. In other words, the function is indeed the abs of the volume of the tetrahedron. Then I integrate this over the sphere w.r.t. uniform distribution. Do you see what I mean? Thanks for your interest. – Romeo Mar 1 '19 at 17:41
• The Det computes the parallelpiped volume. I guess the tetrahedron volume would be dividing it by 3! = 6, correct? – MikeY Mar 1 '19 at 18:35
• Also, you should be able to fix one of the vectors in your Cartesian formulation to just {1,0,0}, correct? If you are looking for the expected volume of the tetrahedron defined by 3 random vectors. – MikeY Mar 1 '19 at 18:47
• @MikeY Yep, you are right, I forgot the factor 6, but I think it is not relevant for the computation of the integral. I have also fixed one vector to be (0,0,1) (see last attempt). However even this reduced model does not work. Any ideas of what is going wrong? – Romeo Mar 1 '19 at 19:00

I'll use the spherical-coordinates approach, and I'll assume for now that by "surface measure" you mean the Haar measure of $$4\pi$$ total area covering the unit sphere uniformly. If instead you are looking for an average, see further below.

The full integral would be

Integrate[Abs[Det[{{Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]},
{Sin[θ2] Cos[φ2], Sin[θ2] Sin[φ2], Cos[θ2]},
{Sin[θ3] Cos[φ3], Sin[θ3] Sin[φ3], Cos[θ3]}}]]*
Sin[θ1] Sin[θ2] Sin[θ3],
{θ1, 0, π}, {φ1, 0, 2 π}, {θ2, 0, π}, {φ2, 0, 2 π}, {θ3, 0, π}, {φ3, 0, 2 π}]


but doesn't seem to evaluate.

Rotational invariance means that we can set $$\theta_3=\phi_3=0$$ (but keep a factor of $$4\pi$$ for the corresponding surface integral measure), as well as $$\phi_2=0$$ (but keep a factor of $$2\pi$$ for the corresponding line integral measure):

4π*2π*Integrate[Abs[Det[{{Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]},
{Sin[θ2], 0, Cos[θ2]},
{0, 0, 1}}]]*
Sin[θ1] Sin[θ2],
{θ1, 0, π}, {φ1, 0, 2 π}, {θ2, 0, π}]


8 π^4

The numerical integral agrees with this result of $$8\pi^4\approx779.273$$:

NIntegrate[Abs[Det[{{Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]},
{Sin[θ2] Cos[φ2], Sin[θ2] Sin[φ2], Cos[θ2]},
{Sin[θ3] Cos[φ3], Sin[θ3] Sin[φ3], Cos[θ3]}}]]*
Sin[θ1] Sin[θ2] Sin[θ3],
{θ1, 0, π}, {φ1, 0, 2 π}, {θ2, 0, π}, {φ2, 0, 2 π}, {θ3, 0, π}, {φ3, 0, 2 π}]


765 ± 17

Playing with the Method option of NIntegrate will give better precision.

If you're looking for the average instead, you need to divide by the used surface measure, $$(4\pi)^3$$ (three times a spherical surface of $$4\pi$$):

8 π^4/(4 π)^3


π/8

and if you are looking further for the average tetrahedral volume instead of the parallelepiped, then as @MikeY says you further divide by 6:

8 π^4/(4 π)^3/6


π/48

This result agrees with @MikeY's numerical answer.

• How are you getting from a surface area to a volume? That is not just a matter of dividing by 6. – Matt F. Mar 2 '19 at 19:21
• @MattF. But that is actually what you have to integrate: the integrand is the volume (up to a factor 6). And it has to be integrated wrt a surface measure. The area formula then allows to represent the surface integral as a volume integral (roughly $d\sigma = (something) du dv$ where u,v are coordinates on the surface and "something" is the Jacobi determinant). Do you see what I mean? – Romeo Mar 2 '19 at 22:06
• @MattF. As @Romeo says the integrand f=Abs[Det[...]] is the volume of the parallelepiped, and integrating it over the three unit spheres with $d\theta_1\,d\phi_1\,d\theta_2\,d\phi_2\,d\theta_3\,d\phi_3$ gives a result in "units" of Volume*Area^3. Dividing this by $(4\pi)^3$ gives "units" of average parallelepiped Volume. Dividing this by 6 (unitless) gives "units" of average tetrahedral Volume. – Roman Mar 2 '19 at 23:25

Not a symbolic answer, but gives a target to shoot for.

You want the average volume of a tetrahedron with edges defined by points randomly picked from the sphere. So first the routine to randomly pick from the sphere, define a multivariate normal with mean {0,0,0} and variance = IdentityMatrix[3]. Then sample randomly from this distribution, which is evenly dispersed in all directions, and normalize so they are of length = 1.

mnd = MultinormalDistribution[{0, 0, 0}, IdentityMatrix[3]];
spherePoint := (sp = RandomVariate[mnd])/Norm[sp];


Now just compute the expected value of your function.

n = 100000;
res = 1/n Sum[Abs@Det[{spherePoint, spherePoint, spherePoint}]/3!, {n}]


0.065389

By comparison, with n=100 the estimate is .0643, so it converges rapidly.

No obvious connection to $$\pi$$.

• That's $\pi/48\approx0.0654498$, see my answer for the derivation. – Roman Mar 2 '19 at 5:40
• Ahhh, I should have picked that up. – MikeY Mar 2 '19 at 15:04
• @MikeY A stupid question, sorry: is the normal distribution you are considering the uniform distribution on the sphere? I do not think so, am I wrong? But we want to pick the points random wrt uniform distribution, not normal one. Btw I did not know Mathematica can also do probability, good point :-) – Romeo Mar 3 '19 at 13:55
• It’s the multivariate normal in 3 dimensions. The beauty of it is when you have covariance matrix $\Sigma$ equal the identity matrix, the direction that points are located from the mean is uniformly distributed over a sphere. The distance from the mean looks like the bell curve you are thinking of. That’s why we normalize it to 1. Works in arbitrary n-dimensional problems too. Try it in 2D to confirm. – MikeY Mar 3 '19 at 14:04
• @MikeY You can also use RandomPoint[Sphere[],n] instead. – Henrik Schumacher Mar 7 '19 at 6:58

For you spherical coordinates approach, I think you just used the wrong measure. It should be Sin[b] Sin[d] not Sin[a] Sin[c]. So:

Integrate[
Abs[Det[{
{Sin[a] Cos[b], Sin[c] Cos[d], 0},
{Sin[a] Sin[b], Sin[c] Sin[d], 0},
{Cos[a], Cos[c], 1}
}]] Sin[b] Sin[d],
{a, 0, 2 Pi}, {c, 0, 2 Pi}, {b, 0, Pi}, {d, 0, Pi}
]


12 π

• Looking here it seems I got the right jacobian, didn't I? Btw it is curious because the volume of the unit sphere is $4/3 \pi$ and you get the average volume of "something inside the sphere" is $12 \pi \cdot1/6 = 2\pi > 4/3 \pi$. So that cannot be, for sure. Thanks for your time and interest (btw I am not the downvoter). – Romeo Mar 2 '19 at 9:46

Picking up on the interpretation of the meaning of the integral discussed in comments:

trying to compute the average of the volume of a tetrahedron spanned by the origin and 3 points on the surface of a sphere, chosen randomly from a uniform distribution on the surface

and taking a brute-force approach:

Mean@
Table[
Volume@ Tetrahedron[{{0, 0, 0}} ~ Join ~ RandomPoint[Sphere[], 3]],
{100000}
]

(* Out: 0.0655448 *)


The result is in accordance with the symbolic one ($$\pi /48$$) derived in Roman's answer above.

• As an OBTW, while fiddling with our solutions, I found that the fastest method is a mix of our approaches,res = 1/n Sum[Abs@Det[RandomPoint[Sphere[], 3]]/3!, {n}], which is $10\times$ faster than yours and over twice as fast as mine. The RandomPoint[ ] is faster, but the Volume calculation is a big time sink. Surprisingly, a Join is too. – MikeY Mar 4 '19 at 14:36
• @Mike and Marco, you can just generate all the random points for your Monte Carlo in one go before computation: With[{n = 1*^6}, Total[Abs[Det /@ RandomPoint[Sphere[], {n, 3}]]/(6 n), Method -> "CompensatedSummation"]]. Of course, this consumes more memory. – J. M.'s technical difficulties Mar 8 '19 at 7:19