Sphere Tetrahedron Picking

Author

Eric W. Weisstein
March 14, 2017

This notebook downloaded from http://mathworld.wolfram.com/notebooks/ComputationalGeometry/SphereTetrahedronPicking.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/SphereTetrahedronPicking.html.

©2017 Wolfram Research, Inc. except for portions noted otherwise

Area of a Random Triangle

TetrahedronVolume[{x1_,x2_,x3_,x4_}]:=Module[{i={1,1,1,1}},Det[Join[Transpose[{x1,x2,x3,x4}],{i}]]/3!]

Exact

(4/105)π//N

mu

Attempt to Guess

v=RootApproximant[(mu/π),1]

(4/105)π

N[v-(4/105)]

v=Rationalize[(mu/π),.001]

N[v-(4/105)]

{4/105.,1/26.}

WolframAlpha["0.11964",{{"PossibleClosedForm",3},"FormattedData"}]

RandomPoint

tet=RandomPoint[Sphere[],{(10^6),4}];//Timing

Length[vols=Abs[TetrahedronVolume/@tet]]//Timing

TimeConstrained[Length[vols=Volume[Tetrahedron[#]]&/@tet],60]//Timing

TimeConstrained[Length[vols=RegionMeasure[Tetrahedron[#]]&/@tet],60]//Timing

bug(330535): Area, Volume, and RegionMeasure should be faster on shape primitives

Through[{Mean,Variance}[vols]]

Histogram[vols,100]

Monte Carlo

SphericalDistribution[num_]:=Table[{Cos[#1]Sqrt[1-#2^2],Sin[#1]Sqrt[1-#2^2],#2}&[Random[Real,{0,2Pi}],Random[Real,{-1,1}]],{num}]

Length[pts=Abs/@TetrahedronVolume/@Partition[SphericalDistribution[10^7],4]]//Timing

{mu,var}=Through[{Mean,Variance}[pts]]

mu^2-(2/81)^2

Histogram[pts]

Distances using (u, v)

r[u_,v_]:=With[{θ=2π u,ϕ=ArcCos[2v-1]},
{Cos[θ]Sin[ϕ],Sin[θ]Sin[ϕ],Cos[ϕ]}
]

{r[0,1],r[0,v1]}//Simplify

voluv=FullSimplify[TetrahedronVolume[{r[0,1],r[0,v2],r[u3,v3],r[u4,v4]}],And@@(0<#<1&/@{v2,u3,v3,u4,v4})]

voluv=-(4/3) ((Sqrt[(-1+v2) v2 (-1+v3) v3]) (-1+v4) Sin[2 π u3]-(-1+v2) (Sqrt[(-1+v3) v3 (-1+v4) v4]) Sin[2 π (u3-u4)]-(-1+v3) (Sqrt[(-1+v2) v2 (-1+v4) v4]) Sin[2 π u4]);

<<Utilities`Simplify`

Flatten[t/.Table[Solve[Recognize[a,1,t,10^n]==0,t],{n,5,15}]]

N[%,10]

Numerical

v5.2

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->QuasiMonteCarlo]//Timing

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->QuasiMonteCarlo,PrecisionGoal->3,MaxPoints->10^6]//Timing

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->QuasiMonteCarlo,PrecisionGoal->4,MaxPoints->10^7]//Timing

v6

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1}]//Timing

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->{SymbolicPreprocessing,Method->QuasiMonteCarlo}]//Timing

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->{SymbolicPreprocessing,Method->QuasiMonteCarlo},PrecisionGoal->3,MaxPoints->10^6]//Timing

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->{SymbolicPreprocessing,Method->AdaptiveQuasiMonteCarlo},PrecisionGoal->3,MaxPoints->10^6]//Timing

NIntegrate[Abs[voluv],{v2,0,1},{u3,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->{SymbolicPreprocessing,Method->AdaptiveQuasiMonteCarlo},PrecisionGoal->4,MaxPoints->10^7]//Timing

Analytic

Timing[i1=Integrate[Sqrt[v^2],{v3,0,1}]]

Timing[i1=Integrate[Sqrt[v^2],{v3,0,1}]]

(4/3)Integrate[Sqrt[v^2],{v2,0,1}]

Timing[i1=Integrate[Sqrt[v^2],{u3,0,1}]]

NIntegrate[i1,{v2,0,1},{v3,0,1},{u4,0,1},{v4,0,1},Method->QuasiMonteCarlo]

...wrong..

Timing[i1=Integrate[Sqrt[v^2],{v3,0,1}]]

Timing[i1=Integrate[Sqrt[v^2],{v3,0,1}]]

(4/3)Integrate[Sqrt[v^2],{v2,0,1}]

Distances using (u, θ)

r[θ_,u_]:={Sqrt[1-u^2]Cos[θ],Sqrt[1-u^2]Sin[θ],u}

{r[0,1],r[0,u]}

v=TetrahedronVolume[{r[0,1],r[0,u2],r[θ3,u3],r[θ4,u4]}]

v//InputForm

v=(1/6) (-(Sqrt[1-(u3^2)]) (-(Sqrt[1-(u2^2)]) (1-u4)+(1-u2) (Sqrt[1-(u4^2)]) Cos[θ4]) Sin[θ3]+(Sqrt[1-(u4^2)]) (-(Sqrt[1-(u2^2)]) (1-u3)+(1-u2) (Sqrt[1-(u3^2)]) Cos[θ3]) Sin[θ4]);

vol=(1/(6(Sqrt[1-v1^2]))) (-(Sqrt[1-(u3^2)]) (-(Sqrt[1-(u2^2)]) (1-u4)+(1-u2) (Sqrt[1-(u4^2)]) Cos[θ4]) (Sqrt[1-v1^2])+(Sqrt[1-(u4^2)]) (-(Sqrt[1-(u2^2)]) (1-u3)+(1-u2) (Sqrt[1-(u3^2)])v1)Sin[θ4]);

Numerical

v5.x

(1/(16(π^2)))NIntegrate[Abs[vol],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6]//Timing

(1/(16(π^2)))NIntegrate[Abs[vol],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->3]//Timing

v6

(1/(16(π^2)))NIntegrate[Abs[vol],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{θ4,0,2π}]//Timing

(1/(16(π^2)))NIntegrate[Abs[vol],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{θ4,0,2π},Method->{SymbolicPreprocessing,Method->QuasiMonteCarlo}]//Timing

(1/(16(π^2)))NIntegrate[Abs[vol],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{θ4,0,2π},Method->{SymbolicPreprocessing,Method->AdaptiveQuasiMonteCarlo},PrecisionGoal->3]//Timing

Numerical splitting into Two parts

vol1=(1/(6(Sqrt[1-v2^2])(Sqrt[1-v1^2]))) (-(Sqrt[1-(u3^2)]) (-(Sqrt[1-(u2^2)]) (1-u4)+(1-u2) (Sqrt[1-(u4^2)]) v2) (Sqrt[1-v1^2])+(Sqrt[1-(u4^2)]) (-(Sqrt[1-(u2^2)]) (1-u3)+(1-u2) (Sqrt[1-(u3^2)])v1)(Sqrt[1-v2^2]));

vol2=(1/(6(Sqrt[1-v2^2])(Sqrt[1-v1^2]))) (-(Sqrt[1-(u3^2)]) (-(Sqrt[1-(u2^2)]) (1-u4)+(1-u2) (Sqrt[1-(u4^2)]) v2) (Sqrt[1-v1^2])+(Sqrt[1-(u4^2)]) (-(Sqrt[1-(u2^2)]) (1-u3)+(1-u2) (Sqrt[1-(u3^2)])v1)(-1)(Sqrt[1-v2^2]));

v5.x

i1=(1/(16(π^2)))NIntegrate[Abs[vol1],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{v2,-1,1},Method->QuasiMonteCarlo,MaxPoints->10^6]

i2=(1/(16(π^2)))NIntegrate[Abs[vol2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{v2,-1,1},Method->QuasiMonteCarlo,MaxPoints->10^6]

i1+i2

v6

i1=(1/(16(π^2)))NIntegrate[Abs[vol1],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{v2,-1,1}]

i1=(1/(16(π^2)))NIntegrate[Abs[vol1],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{v2,-1,1},Method->{SymbolicPreprocessing,Method->QuasiMonteCarlo}]

i2=(1/(16(π^2)))NIntegrate[Abs[vol2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{v2,-1,1}]

i2=(1/(16(π^2)))NIntegrate[Abs[vol2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{v1,-1,1},{v2,-1,1},Method->{SymbolicPreprocessing,Method->QuasiMonteCarlo}]

Normalization

Integrate[1,{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,2Pi},{θ4,0,2Pi}]

(Subscript[μ, 1]^′)

(1/(32(π^2)))NIntegrate[Evaluate[Abs[v]],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,2Pi},{θ4,0,2Pi},Method->QuasiMonteCarlo,PrecisionGoal->3,MaxPoints->10^6]//Timing

(1/(16(π^2)))NIntegrate[Evaluate[Abs[v]],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi},Method->QuasiMonteCarlo,PrecisionGoal->3,MaxPoints->10^6]//Timing

(1/(16(π^2)))NIntegrate[Evaluate[Sqrt[v^2]],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi},Method->QuasiMonteCarlo,PrecisionGoal->3,MaxPoints->10^6]//Timing

(1/(16(π^2)))NIntegrate[Evaluate[Abs[v]],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi},Method->QuasiMonteCarlo,PrecisionGoal->4,MaxPoints->10^7]//Timing

(1/(16(π^2)))NIntegrate[Evaluate[Abs[v]],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Timing

(1/(16(π^2)))NIntegrate[Evaluate[Abs[v]],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi},Method->MultiDimensional]//Timing

v4.2

(1/(16(π^2)))Integrate[Sqrt[v^2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Timing

FullSimplify[(1/(192 (π^2)))(I ((π^2) Log[-I (Sqrt[2])]-(π^2) Log[I (Sqrt[2])]+3 (π^2) Log[-2 I (Sqrt[2])]-3 (π^2) Log[2 I (Sqrt[2])]))]

N[(π/48)]

...wrong...

v5

(1/(16(π^2)))Integrate[Sqrt[v^2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Timing

N[Pi/192]

...wrong...

v6

(1/(16(π^2)))Integrate[Sqrt[v^2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Timing

N[Pi/96]

...wrong...

(1/(16(π^2)))NIntegrate[Sqrt[v^2],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,π},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6]

By Pieces: (Subscript[u, 2])

Timing[i1=(1/(16(π^2)))Integrate[Sqrt[v^2],{u2,-1,1},Assumptions->-1<u3<1&&-1<u4<1&&0<θ3<π&&0<θ4<2π]]

NIntegrate[Evaluate[i1],{u3,-1,1},{u4,-1,1},{θ3,0,π},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6]//Timing

By Pieces: (Subscript[u, 3])(v5)

Timing[i1=Simplify[(1/(16(π^2)))Integrate[Sqrt[v^2],{u3,-1,1},Assumptions->-1<u2<1&&-1<u4<1&&0<θ3<π&&0<θ4<2π]]]

NIntegrate[Evaluate[i1],{u2,-1,1},{u4,-1,1},{θ3,0,π},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6]//Timing

By Pieces: (Subscript[u, 4])(v5)

Timing[i1=Simplify[(1/(16(π^2)))Integrate[Sqrt[v^2],{u4,-1,1},Assumptions->-1<u2<1&&-1<u3<1&&0<θ3<π&&0<θ4<2π]]]

NIntegrate[Evaluate[i1],{u2,-1,1},{u3,-1,1},{θ3,0,π},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6]//Timing

By Pieces: (Subscript[θ, 3])

Timing[i1=Simplify[(1/(16(π^2)))Integrate[Sqrt[v^2],{θ3,0,Pi},Assumptions->-1<u2<1&&-1<u3<1&&-1<u4<1&&0<θ4<2π]]]

NIntegrate[Evaluate[i1],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ4,0,2π},Method->QuasiMonteCarlo,MaxPoints->10^6]//Timing

By Pieces: (Subscript[θ, 4])

i1=(1/(16(π^2)))Integrate[Sqrt[v^2],{θ4,0,2Pi},Assumptions->-1<u2<1&&-1<u3<1&&-1<u4<1&&0<θ3<π]

NIntegrate[Evaluate[i1],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,π},Method->QuasiMonteCarlo,MaxPoints->10^6]

..wrong..

Attempt to CAD

Simplify`PiecewiseIntegrate;
PiecewiseDump`$ShowIntegrals=True;
i=Integrate[Abs[v],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,2Pi},{θ4,0,2Pi}];//Timing

...ran for 10 hours using 246 MB with no result...

Simplify`PiecewiseIntegrate;
PiecewiseDump`$ShowIntegrals=True;
i=Integrate[Abs[v],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}];//Timing

Timing[
i=Integrate[Abs[v],{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,2Pi},{θ4,0,2Pi}]
]

...crash...

(Subscript[μ, 2]^′)

(1/(16(π^2)))Integrate[v^2,{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Simplify//Timing

(Subscript[μ, 4]^′)

(1/(16(π^2)))Integrate[v^4,{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Simplify//Timing

(Subscript[μ, 6]^′)

(1/(16(π^2)))Integrate[v^6,{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Simplify//Timing

(Subscript[μ, 8]^′)

(1/(16(π^2)))Integrate[v^8,{u2,-1,1},{u3,-1,1},{u4,-1,1},{θ3,0,Pi},{θ4,0,2Pi}]//Simplify//Timing