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...

  • 1
    $\begingroup$ 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) $\endgroup$ – Matt Stein Mar 1 '19 at 15:39
  • 1
    $\begingroup$ @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. $\endgroup$ – Romeo Mar 1 '19 at 17:41
  • 1
    $\begingroup$ The Det computes the parallelpiped volume. I guess the tetrahedron volume would be dividing it by 3! = 6, correct? $\endgroup$ – MikeY Mar 1 '19 at 18:35
  • 1
    $\begingroup$ 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. $\endgroup$ – MikeY Mar 1 '19 at 18:47
  • 1
    $\begingroup$ @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? $\endgroup$ – 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


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


This result agrees with @MikeY's numerical answer.

  • 1
    $\begingroup$ How are you getting from a surface area to a volume? That is not just a matter of dividing by 6. $\endgroup$ – Matt F. Mar 2 '19 at 19:21
  • $\begingroup$ @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? $\endgroup$ – Romeo Mar 2 '19 at 22:06
  • $\begingroup$ @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. $\endgroup$ – 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}]


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

No obvious connection to $\pi$.

  • $\begingroup$ That's $\pi/48\approx0.0654498$, see my answer for the derivation. $\endgroup$ – Roman Mar 2 '19 at 5:40
  • 1
    $\begingroup$ Ahhh, I should have picked that up. $\endgroup$ – MikeY Mar 2 '19 at 15:04
  • $\begingroup$ @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 :-) $\endgroup$ – Romeo Mar 3 '19 at 13:55
  • 2
    $\begingroup$ 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. $\endgroup$ – MikeY Mar 3 '19 at 14:04
  • 1
    $\begingroup$ @MikeY You can also use RandomPoint[Sphere[],n] instead. $\endgroup$ – 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:

        {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 π

  • $\begingroup$ 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). $\endgroup$ – 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:

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

(* Out: 0.0655448 *)

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

  • $\begingroup$ 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. $\endgroup$ – MikeY Mar 4 '19 at 14:36
  • $\begingroup$ @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. $\endgroup$ – J. M.'s ennui Mar 8 '19 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.