6
$\begingroup$

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

Improve this question
$\endgroup$
12
  • 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$ Commented Mar 1, 2019 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$ Commented Mar 1, 2019 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$ Commented Mar 1, 2019 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$ Commented Mar 1, 2019 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$ Commented Mar 1, 2019 at 19:00

4 Answers 4

Reset to default
7
+50
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ How are you getting from a surface area to a volume? That is not just a matter of dividing by 6. $\endgroup$ Commented Mar 2, 2019 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$ Commented Mar 2, 2019 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$ Commented Mar 2, 2019 at 23:25
5
$\begingroup$

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

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

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 π

Improve this answer
$\endgroup$
1
  • $\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$ Commented Mar 2, 2019 at 9:46
1
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\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$ Commented Mar 4, 2019 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$ Commented Mar 8, 2019 at 7:19

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.