Tetrahedron Tetrahedron Picking

Author

Eric W. Weisstein
February 6, 2005

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

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

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

Initialization

Exact

c=13/720-Pi^2/15015

N[c,10]

Inequalities

ineqs[x_,y_,z_]:=And@@{x+y+z≤2 ,x-y-z≤0,-x+y-z≤0,-x-y+z≤0};

<<Graphics`InequalityGraphics`

g=InequalityPlot3D[ineqs[x,y,z],{x},{y},{z}]

Numeric Monte Carlo

StandardDeviationOfMean[l_]:=Sqrt[Variance[l]/Length[l]]
StandardDeviationOfVariance[l_]:=Module[{n=Length[l]},
Sqrt[(n-1)^2/n^3(Plus@@((l-Mean[l])^4)/n)-(n-1)(n-3)Variance[l]^2/n^3]
]

Length[pts=Partition[Select[Table[(3^(1/3))Random[],{10^7},{3}],ineqs@@#&],4]]//Timing

Show[{g/.Polygon[x_]:>Line[Append[x,First[x]]],Graphics3D[{Red,Tetrahedron3D/@Take[pts,15]}]}]

dist=3Abs/@TetrahedronVolume/@pts;//Timing

Partition[Through[{Mean,StandardDeviationOfMean,Variance,StandardDeviationOfVariance}[dist]],2]

<<Graphics`Graphics`

BarChart[Length/@Split[Ceiling[600Sort[dist]]],Ticks->None,BarEdges->None,BarSpacing->-.2,PlotRange->{{0,70},All},AxesLabel->TraditionalForm/@{V,P[V]}]

Numerical Integration

Setup

a=Array[x,{4,3}];

{#,0,1}&/@Flatten[Array[x,{4,3}]]

TetrahedronVolume[{{0,0,0},{1,0,0},{0,1,0},{0,0,1}}]

And@@PolyhedronInequalities[Orient[Tetrahedron3D[{{0,0,0},{1,0,0},{0,1,0},{0,0,1}}],{1,1,1}/3],{x,y,z}]//FullSimplify

ineq[{x_,y_,z_}]:=x>0&&y>0&&z>0&&x+y+z<1

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

ineqs1=And@@ineq/@vars

ineqs1=x1>0&&y1>0&&z1>0&&x1+y1+z1<1&&x2>0&&y2>0&&z2>0&&x2+y2+z2<1&&x3>0&&y3>0&&z3>0&&x3+y3+z3<1&&x4>0&&y4>0&&z4>0&&x4+y4+z4<1;

vol=TetrahedronVolume[vars]

vol=(1/6) (-x3 y2 z1+x4 y2 z1+x2 y3 z1-x4 y3 z1-x2 y4 z1+x3 y4 z1+x3 y1 z2-x4 y1 z2-x1 y3 z2+x4 y3 z2+x1 y4 z2-x3 y4 z2-x2 y1 z3+x4 y1 z3+x1 y2 z3-x4 y2 z3-x1 y4 z3+x2 y4 z3+x2 y1 z4-x3 y1 z4-x1 y2 z4+x3 y2 z4+x1 y3 z4-x2 y3 z4);

Normalization

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

v5

<<Calculus`Integration`

Integrate[Boole[ineqs],Evaluate[Sequence@@List/@Flatten[vars]]]//Timing

v6

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

QuasiMonteCarlo

v5.0

<<Calculus`Integration`

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

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

6 1296NIntegrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@List/@Flatten[vars]],Method->QuasiMonteCarlo,MaxPoints->10^7,PrecisionGoal->3]//Timing

v5.1

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

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

Integral

v5.0

<<Calculus`Integration`

i=Calculus`Integration`Private`DebugIntegrate[Evaluate[Abs[TetrahedronVolume[vars]]]Boole[ineqs],Evaluate[Sequence@@List/@Flatten[vars]]];//Timing

...Uses more than 1.6 GB...

v5.1

(i=6 1296 Integrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])]])//Timing

...crashed?...

Simplify`PiecewiseIntegrate;
PiecewiseDump`$ShowIntegrals=True;
(i=Integrate[Abs[vol]Boole[ineqs1],Sequence@@({#,-∞,∞}&/@Flatten[vars])])//Length//Timing

v6

(i=6 1296 Integrate[Abs[vol]Boole[ineqs],Evaluate[Sequence@@({#,-Infinity,Infinity}&/@Flatten[vars])]])//Timing

Simplify`PiecewiseIntegrate;
PiecewiseDump`$ShowIntegrals=True;
i=Integrate[Abs[vol]Boole[ineqs],Sequence@@({#,-∞,∞}&/@Flatten[vars])];//Timing

...Used more than 1.6 GB and ran out of memory...