Disk Triangle Picking

Author

Eric W. Weisstein
February 6, 2005

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

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

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

Code

RandomTriangleInDisk[]:=Table[Sqrt[#1]{Cos[2Pi#2],Sin[2Pi#2]}&@@Table[Random[],{2}],{3}]
RandomTriangles[n_]:=Table[{{Hue[Random[]],Polygon[#]},Line[Append[#,First[#]]]}&@RandomTriInDisk[],{n}]

Random Plots

l=RandomTriangles[50];
movie=Table[Graphics[
{Circle[{0,0},1],Take[l,n]},
AspectRatio->Automatic],{n,0,Length[l]}];

Export[gifs<>"distripk.gif",Show[#,PlotRange->{{-1,1},{-1,1}}]&/@movie]//Timing

Show[Graphics[{Circle[{0,0},1],l},AspectRatio->Automatic]]

Exact

N[(35/(48(π^2))),10]

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

areas=Area/@Triangle/@Table[RandomTriangleInDisk[],{10^6}]/Pi;//Timing

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

<<Graphics`Graphics`

BarChart[Length/@Split[Ceiling[700Sort[areas]]],Ticks->None,BarEdges->None,BarSpacing->-.2,PlotRange->{{0,250},All},AxesLabel->TraditionalForm/@{A,P[A]}]

2-D

Monte Carlo

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

Length[tris=Partition[pts,3]]

TriangleArea2D[l_]:=Abs[Det[Transpose[Append[Transpose[l],{1,1,1}]]]/2]

Mean[data=TriangleArea2D/@tris]/Pi//Timing

Histogram[data,{.02},"Probability",AxesLabel->{A,P[A]}]

Exact

N[35/(48Pi^2)]

CAD

Area[Triangle[{{x1,y1},{x2,y2},{x3,y3}}]]

Simplify`PiecewiseIntegrate;
PiecewiseDump`$ShowIntegrals=True;

Length[i=Integrate[Abs[-x2 y1+x3 y1+x1 y2-x3 y2-x1 y3+x2 y3]/2Boole[x1^2+y1^2<1&&x2^2+y2^2<1&&x3^2+y3^2<1],{x1,-∞,∞},{y1,-∞,∞},{x2,-∞,∞},{y2,-∞,∞},{x3,-∞,∞},{y3,-∞,∞}]]//Timing

Exact

i[[1]]

ReleaseHold[i[[1]]]//Timing

...crash...

i[[1]]/.HoldForm[Integrate[x__]]:>IntegrateWithAssumptions[x]//Timing

Numerical All at Once

(i/.HoldForm[Integrate[x__]]:>NIntegrate[x,Method->QuasiMonteCarlo,MaxPoints->10^6]/Pi^4)//Timing

Numeric One at a Time

l={};
Do[
Print[{n,int=i[[n]]/.HoldForm[Integrate[x__,opts___?OptionQ]]:>NIntegrate[x,Method->QuasiMonteCarlo,MaxPoints->10^6]//Timing}];
AppendTo[l,int],
{n,Length[i]}
];

Plus@@l

int=Last[Plus@@l]/Pi^4

N[(35/(48(π^2))),10]

Break up the Integral

r[u_,v_]:=Sqrt[u]{Cos[v],Sin[v]}

{r[u1,0],r[uu,vv]}

vol=Area[Evaluate[Triangle[{r[u1,0],r[u2,v2],r[u3,v3]}]]]

vol=(1/2) Abs[(Sqrt[u1]) (Sqrt[u2]) Sin[v2]-(Sqrt[u2]) (Sqrt[u3]) Cos[v3] Sin[v2]-(Sqrt[u1]) (Sqrt[u3]) Sin[v3]+(Sqrt[u2]) (Sqrt[u3]) Cos[v2] Sin[v3]];

(1/(2(π^3)))NIntegrate[vol,{u1,0,1},{u2,0,1},{v2,0,π},{u3,0,1},{v3,0,2π},Method->QuasiMonteCarlo]

(1/(2(π^3)))Integrate[vol,{u1,0,1},{u2,0,1},{v2,0,π},{u3,0,1},{v3,0,2π}]//Timing

vol1=(Sqrt[u1 u2(1-(w2^2))]) -(Sqrt[u2 u3(1-(w2^2))])  w3-(Sqrt[u1 u3(1-(w3^2))])+(Sqrt[u2 u3(1-(w3^2))])  w2;
vol2=(Sqrt[u1 u2(1-(w2^2))]) -(Sqrt[u2 u3(1-(w2^2))])  w3+(Sqrt[u1 u3(1-(w3^2))])-(Sqrt[u2 u3(1-(w3^2))])  w2;

(1/(4(π^3)))NIntegrate[Evaluate[((Abs[vol1]+Abs[vol2])/(Sqrt[(1-(w2^2))(1-(w3^2))]))],{u1,0,1},{u2,0,1},{w2,-1,1},{u3,0,1},{w3,-1,1},Method->QuasiMonteCarlo]

(1/(4(π^3)))Integrate[Evaluate[((Abs[vol1])/(Sqrt[(1-(w2^2))(1-(w3^2))]))],{u1,0,1},{u2,0,1},{w2,-1,1},{u3,0,1},{w3,-1,1}]//Timing

(1/(4(π^3)))Integrate[Evaluate[((Abs[vol2])/(Sqrt[(1-(w2^2))(1-(w3^2))]))],{u1,0,1},{u2,0,1},{w2,-1,1},{u3,0,1},{w3,-1,1}]//Timing

(1/(4(π^3)))NIntegrate[((Abs[#])/(Sqrt[(1-(w2^2))(1-(w3^2))])),{u1,0,1},{u2,0,1},{w2,-1,1},{u3,0,1},{w3,-1,1},Method->QuasiMonteCarlo,MaxPoints->10^6]&/@{vol1,vol2}

Plus@@%

Reduce[vol1>0&&0<u1<1&&0<u2<1&&0<u3<1&&-1<w2<1&&-1<w3<1,{u1,u2,u3,w1,w2}]

Expanding 2nd part

Expand[((Sqrt[u1 u2(1-(w2^2))]) -(Sqrt[u2 u3(1-(w2^2))])  w3+(Sqrt[u1 u3(1-(w3^2))])-(Sqrt[u2 u3(1-(w3^2))])  w2)^2]

i=u1 u2+u1 u3-u1 u2 (w2^2)+u2 u3 (w2^2)-2 u2(1-(w2^2))(Sqrt[u1 u3 ])  w3-u1 u3 (w3^2)+u2 u3 (w3^2)-2 u2 u3 (w2^2) (w3^2)+2 u1(Sqrt[u2 (1-(w2^2)) u3 (1-(w3^2))])-2u3 (Sqrt[u2 (1-(w2^2))u1 (1-(w3^2))]) w3 -2 w2  u2(Sqrt[u1 (1-(w2^2)) u3 (1-(w3^2))])+2 w2 u2 u3(Sqrt[(1-(w2^2))(1-(w3^2))]) w3 -2 w2 (1-(w3^2))u3(Sqrt[u1 u2])

i2=(1/(2(π^3)))NIntegrate[vol2,{u1,0,1},{u2,0,1},{w2,-1,1},{u3,0,1},{w3,-1,1},Method->QuasiMonteCarlo]

(1/(2(π^3)))NIntegrate[Evaluate[((Sqrt[i])/(2 (Sqrt[(1-(w2^2)) (1-(w3^2))])))],{u1,0,1},{u2,0,1},{w2,-1,1},{u3,0,1},{w3,-1,1},Method->QuasiMonteCarlo]

Integrate[((Sqrt[i])/(2 (Sqrt[(1-(w2^2)) (1-(w3^2))]))),{u1,0,1}]

Integrate[((Sqrt[i])/(2 (Sqrt[(1-(w2^2)) (1-(w3^2))]))),{w2,-1,1}]

Acute

N[(4/(π^2))-(1/8),10]

With[{n=10^5},
Count[Triangle/@Table[RandomTriInDisk[],{n}],_?AcuteQ]/n
]//N//Timing