Ball Tetrahedron Picking

Author

Eric W. Weisstein
March 14, 2017

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

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

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

Random Picking

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

Exact Mean

N[((12π)/715),10]

Monte Carlo

$Version

Length[pts=Select[Table[Random[Real,{-1,1}],{10^6},{3}],Norm[#]<1&]]//Timing

tet=Partition[pts,4];

Through[{Mean,Variance}[Abs/@TetrahedronVolume/@tet]]

RandomPoint

tet=RandomPoint[Ball[],{(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]

Symbolic Integration

Setup

ineq[{x_,y_,z_}]:=x^2+y^2+z^2<1

vars={{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4}};

ineqs=And@@ineq/@vars

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

vol=TetrahedronVolume[vars]

Normalization

v5.0

<<Calculus`Integration`

Integrate[Boole[ineqs],Evaluate[Sequence@@({#,-∞,∞}&/@Flatten[vars])]]//Timing

v5.2

Integrate[Boole[ineqs],Evaluate[Sequence@@({#,-∞,∞}&/@Flatten[vars])]]//Timing

...negative of correct result...

Integrate[Boole[ineqs],Evaluate[Sequence@@({#,-∞,∞}&/@Flatten[vars])],GenerateConditions->True]//Timing

v6

Integrate[Boole[ineqs],Evaluate[Sequence@@({#,-∞,∞}&/@Flatten[vars])]]//Timing

...0 times the correct result ;)...

Integrate[Boole[ineqs],Evaluate[Sequence@@({#,-∞,∞}&/@Flatten[vars])],GenerateConditions->True]//Timing

Numerical

v5.0

<<Calculus`Integration`

((((4/3)π))^(-4))NIntegrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@List/@Flatten[vars]],Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->2]//Timing

v5.2

((((4/3)π))^(-4))NIntegrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])],Method->QuasiMonteCarlo,MaxPoints->10^6,PrecisionGoal->2]//Timing

...hangs...

((((4/3)π))^(-4))NIntegrate[Evaluate[Abs[vol]]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])],Method->{SymbolicPiecewiseSubdivision,Method->QuasiMonteCarlo,ExpandSpecialPiecewise->{Boole}},MaxPoints->10^6,PrecisionGoal->2]//Timing

v6

(((((4/3)π))^(-4))NIntegrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])],Method->{SymbolicPiecewiseSubdivision,Method->QuasiMonteCarlo,ExpandSpecialPiecewise->{Boole}},MaxPoints->10^7,PrecisionGoal->#]//Timing)&/@{2,3}

(((((4/3)π))^(-4))NIntegrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])],Method->{SymbolicPiecewiseSubdivision,Method->AdaptiveQuasiMonteCarlo,ExpandSpecialPiecewise->{Boole}},MaxPoints->10^7,PrecisionGoal->#]//Timing)&/@{2,3}

((((4/3)π))^(-4))NIntegrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])],Method->{SymbolicPiecewiseSubdivision,Method->QuasiMonteCarlo,ExpandSpecialPiecewise->{Boole}},MaxPoints->10^6,PrecisionGoal->3]//Timing