Ball Line Picking

Author

Eric W. Weisstein
March 14, 2017

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

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

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

Analytic

Circle

Subtract@@(Sqrt[#1]{Cos[#2],Sin[#2]}&@@@{{(Subscript[r, 1]),(Subscript[θ, 1])},{(Subscript[r, 2]),(Subscript[θ, 2])}})

Simplify[Sqrt[#.#&[Subtract@@(Sqrt[#1]{Cos[#2],Sin[#2]}&@@@{{(Subscript[r, 1]),(Subscript[θ, 1])},{(Subscript[r, 2]),(Subscript[θ, 2])}})]]]//TraditionalForm

Timing[i=(1/(2π))Integrate[(Sqrt[r1+r2-2 Cos[θ] (Sqrt[r1 r2])]),{θ,0,2π},{r1,0,1},{r2,0,1}]]

Integrate[s P[2,s,1],{s,0,2}]

Monte Carlo

Length[lines=Norm/@Subtract@@@Partition[Select[Table[Random[Real,{-1,1}],{10^6},{3}],Norm[#]<1&],2]]//Timing

bins=Length/@Split[Sort[Round[lines 40]]];

<<Graphics`Graphics`

BarChart[bins,Ticks->None,BarEdges->None,BarSpacing->-.2]

N[36/35]

RandomPoint

Generate

lines=RandomPoint[Ball[],{(10^7),2}];//Timing

Length[lens=EuclideanDistance@@@lines]//Timing

Various ways of subtracting

Length[lens=EuclideanDistance@@@lines]//Timing

Length[lens=Norm/@Subtract@@@lines]//Timing

Length[lens=Norm[Subtract@@#]&/@lines]//Timing

Distribution

Length[lens]

With[{bins=100},Show[{Histogram[lens,bins],Plot[Evaluate[(2/bins)Length[lens]P[3,s,1]],{s,0,2},PlotStyle->Opacity[.5]]}]]

Mean[lens]//InputForm

Mean[Take[lens,(10^6)]]//InputForm

Integrate[s P[3,s,1],{s,0,2}]

N[%]

Distribution

P[n_,s_,R_]:=n s^(n-1)/R^n BetaRegularized[1-s^2/4/R^2,(n+1)/2,1/2]

<<Graphics`Colors`

Plot[Evaluate[P[#,s]&/@Range[8]],{s,0,2},AxesLabel->TraditionalForm/@{(Subscript[P, n])[s/R],s/R},
PlotStyle->{Red,Orange,Yellow,Green,Cyan,Blue,Violet,Purple}]

Inability to Integarte BetaRegularized

$Version

Integrate[BetaRegularized[x,a,b],x]

FunctionExpand[%]

Normalization

Table[Integrate[FunctionExpand[P[n,s]],{s,0,2R},Assumptions->0<s<2R],{n,5}]//Timing

Table[Integrate[FunctionExpand[P[n,s]],{s,0,2R},Assumptions->0<s<2R],{n,10}]//Timing

Table[Integrate[FullSimplify[FunctionExpand[P[n,s]]],{s,0,2R}],{n,10}]

Small n

{#,Together//@FullSimplify[FunctionExpand[P[#,s]],s>0&&R>0]}&/@Range[5]//ColumnForm

{#,FullSimplify[FunctionExpand[P[#,s]],R>0&&s>0]/.{(Sqrt[4 (R^2) (s^4)-(s^6)]):>(s^2)(Sqrt[4(R^2)-(s^2)])}}&/@Range[5]//ColumnForm

Means

Exact

Integrate[s P[n,s,R],{s,0,2R}]//Timing

FullSimplify[s FunctionExpand[P[n,s,R]],0<s<2R]

v4.2

Timing[m=With[{R=1},Table[Integrate[FullSimplify[s FunctionExpand[P[n,s,R]]],{s,0,2R}],{n,20}]]]

Integrate[FunctionExpand[s^20*Beta[1 - s^2/4, 21/2, 1/2]],s]

With[{R=1,n=21},Table[Integrate[FullSimplify[s FunctionExpand[P[n,s,R]]],{s,0,2R}]]]//Timing

v5

Integrate[((n ((((s/R)))^n) Beta[1-((s^2)/(4 (R^2))),((1+n)/2),(1/2)] Gamma[1+(n/2)])/((Sqrt[π]) Gamma[((1+n)/2)])),{s,0,2R}]//Timing

v6

With[{R=1},Integrate[((n ((((s/R)))^n) Beta[1-((s^2)/(4 (R^2))),((1+n)/2),(1/2)] Gamma[1+(n/2)])/((Sqrt[π]) Gamma[((1+n)/2)])),{s,0,2R}]]//Timing

Timing[m=With[{R=1},Table[Integrate[FullSimplify[s FunctionExpand[P[n,s,R]]],{s,0,2R}],{n,20}]]]

eo=Transpose[Partition[Most[m],2]]

Through[{Numerator,Denominator}[eo[[1]]]]

Through[{Numerator,Denominator}[Pi eo[[2]]]]

Numeric

N[%/R]

NIntegrate[s FunctionExpand[P[#,s,1]],{s,0,2}]&/@Range[10]//Chop

Plot

<<Utilities`Plot`

means={#,NIntegrate[s P[#,s,1],{s,0,2}]}&/@Range[10]//Chop;//Timing

ListPlot[means,PlotRange->All,AxesLabel->TraditionalForm/@{n,Subscript[OverBar[s],n]},
PlotStyle->{PointSize[.02],Red},AxesOrigin->{1,Automatic}]

2-Ball

PDF

FullSimplify[P[2,s,a]//FunctionExpand,s>0&&a>0]

Integrate[((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]

Mean

Numerical

With[{a=1},NIntegrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]]

(128/(45 π))//N

V5.2 (wrong)

$Version

Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]//Timing

FullSimplify[Last[%],a>0]

...wrong...

V6 (wrong)

$Version

Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]//Timing

...wrong...

With[{a=1},Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]]/.s->s/a

...correct..

V7 (wrong)

$Version

Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]//Timing

((5-5 (Sqrt[1-16 (a^4)])-40 (a^4) (Sqrt[1-16 (a^4)])+96 (a^6) (-2+5 ArcCos[4 (a^2)]))/(45 (a^5) π))/.{a->1}//N

...wrong...

Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a},Assumptions->a>0]//Timing

N[{((5-5 I (Sqrt[-1+16 (a^4)])-40 I (a^4) (Sqrt[-1+16 (a^4)])+96 (a^6) (-2+5 I Log[4 (a^2)+(Sqrt[-1+16 (a^4)])]))/(45 (a^5) π)),((5-192 (a^6))/(45 (a^5) π))}/.a->1]

V8 (wrong)

$Version

Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]//Timing

...wrong...

Trick to get correct result

With[{a=1},Integrate[s((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,2a}]]//Timing

FullSimplify[Last[%1],0<a<1/2]

CDF

$Version

Integrate[((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,s}]//Timing

...hangs...

With[{a=1},Integrate[((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),{s,0,s}]]/.s->s/a//Timing

FullSimplify[Last[%],0<s<2a&&a>0]

%/.s->2a

Plots

Plot[#/.a->1,{s,0,2}]&/@{((s (-s (Sqrt[4 (a^2)-(s^2)])+4 (a^2) ArcSec[((2 a)/s)]))/((a^4) π)),((-((s (2 (a^2)+(s^2)) (Sqrt[4-((s^2)/(a^2))]))/(a^3))+8 ArcCsc[((2 a)/s)]+((8 (s^2) ArcSec[((2 a)/s)])/(a^2)))/(4 π))}

3-Ball

$Version

PDF

FullSimplify[P[3,s,a]//FunctionExpand,s>0&&a>0]//Apart

Integrate[((3 (s^2))/(a^3))-((9 (s^3))/(4 (a^4)))+((3 (s^5))/(16 (a^6))),{s,0,2a}]

Mean

Integrate[s(((3 (s^2))/(a^3))-((9 (s^3))/(4 (a^4)))+((3 (s^5))/(16 (a^6)))),{s,0,2a}]

CDF

Integrate[((3 (s^2))/(a^3))-((9 (s^3))/(4 (a^4)))+((3 (s^5))/(16 (a^6))),{s,0,s}]

((s^3)/(a^3))-((9 (s^4))/(16 (a^4)))+((s^6)/(32 (a^6)))/.s->2a

Plots

Plot[#/.a->1,{s,0,2}]&/@{((3 (s^2))/(a^3))-((9 (s^3))/(4 (a^4)))+((3 (s^5))/(16 (a^6))),((s^3)/(a^3))-((9 (s^4))/(16 (a^4)))+((s^6)/(32 (a^6)))}