Disk Covering Problem

Author

Eric W. Weisstein
September 7, 2007

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

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

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

Code

FractionCovered[disks_,dx_:.01]:=Module[
{
x,y,
r=disks[[1,2]],
c=Compile[{x,y},Evaluate[Boole[Or@@(Norm[{x,y}-#[[1]]]^2<#[[2]]^2&/@disks)]]]
},
Count[Flatten[Table[c[x,y],{x,-1,1,dx},{y,-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)]),dx}]],1]dx^2/Pi
]

DrawDisks[disks_]:=Graphics[{
{Blue,Disk[]},
Opacity[.5],White,EdgeForm[Black],disks
},
Axes->True,Ticks->False,
PlotLabel->"covered fraction = "<>ToString[Round[100FractionCovered[disks]]]<>"%"]

f[x_,dx_:.01]:=FractionCovered[x,dx]/;NumericQ[x[[1,1,1]]]

Numeric Solution

r(1)

disks=Disk[#,1]&/@(xs=Array[x,{1,2}])

(soln=NMaximize[f[disks],Flatten[xs],Method->"NelderMead"])//Timing

DrawDisks[disks/.Last[soln]]

r(2)

disks=Disk[#,1]&/@(xs=Array[x,{2,2}])

(soln=NMaximize[f[disks],Flatten[xs],Method->"NelderMead"])//Timing

DrawDisks[disks/.Last[soln]]

r(3)

disks=Disk[#,(Sqrt[3])/2]&/@(xs=Array[x,{3,2}])

(soln=NMaximize[f[disks],Flatten[xs],Method->"NelderMead"])//Timing

DrawDisks[disks/.Last[soln]]

r(4)

disks=Disk[#,(Sqrt[2])/2]&/@(xs=Array[x,{4,2}])

(soln=NMaximize[f[disks],Flatten[xs],Method->"NelderMead"])//Timing

DrawDisks[disks/.Last[soln]]

r(5)

disks=Disk[#,.609382864]&/@(xs=Array[x,{5,2}]);

(soln=NMaximize[f[disks],Flatten[xs],Method->"NelderMead"])//Timing

DrawDisks[disks/.Last[soln]]

(soln=NMaximize[f[disks,.005],Flatten[xs],Method->"NelderMead"])//Timing

DrawDisks[disks/.Last[soln]]

Exact Solution

System of Equations

eqns={
2Sin[t]-Sin[t+phi/2+psi]-Sin[psi-t-phi/2]==0,
2Sin[phi]-Sin[t+phi/2+chi]-Sin[chi-t-phi/2]==0,
2Sin[t]+Sin[chi+t]-Sin[chi-t]-Sin[psi+phi]-Sin[psi-phi]-2Sin[psi-2t]==0,
Cos[2psi-chi+phi]-Cos[2psi+chi-phi]-2Cos[chi]+Cos[2psi+chi-2t]+Cos[2psi-chi-2t]==0
};

Approximate Solution

soln=FindRoot[
eqns,
{t,36 Degree},{phi, 36 Degree},{psi,72 Degree},{chi,72 Degree},
WorkingPrecision->20
]

Sin/@{t,phi/2}/.soln

Cos[t+phi/2]/.soln

N[Root[625-800 #1+22 ((#1)^2)-1776 ((#1)^3)+1665 ((#1)^4)+1360 ((#1)^5)-3480 ((#1)^6)+2112 ((#1)^7)+1296 ((#1)^8)&,3],20]

Exact Solution

Attempt to solve directly

Solve[eqns,{t,phi,psi,chi}]

eqn2=Join[TrigExpand//@(eqns/.phi->2halfphi)/.{Cos->cos,Sin->sin},
sin[#]^2+cos[#]^2==1&/@{halfphi,t,psi,chi}]

Solve[eqn2,{sin[t],sin[halfphi]},{cos[t],cos[halfphi],sin[psi],cos[psi],sin[chi],cos[chi]}]

Took 2.05 GB and entered swap death

Use Gröbner bases

solns=Factor/@GroebnerBasis[eqn2,{sin[t],sin[halfphi]},{cos[t],cos[halfphi],sin[psi],cos[psi],sin[chi],cos[chi]},MonomialOrder->EliminationOrder]

solns/.halfphi->phi/2/.soln/.sin:>Sin

GroebnerBasis[{(38484 sin[halfphi]-199557 ((sin[halfphi])^3)+168228 ((sin[halfphi])^5)-47244 sin[t]+66171 ((sin[halfphi])^2) sin[t]-35628 ((sin[halfphi])^4) sin[t]+24573 sin[halfphi] ((sin[t])^2)+14036 ((sin[halfphi])^3) ((sin[t])^2)+85037 ((sin[t])^3)-113862 ((sin[halfphi])^2) ((sin[t])^3)-11260 sin[halfphi] ((sin[t])^4)+11022 ((sin[t])^5))==0,(329664-1059045 ((sin[halfphi])^2)+680160 ((sin[halfphi])^4)+192816 ((sin[halfphi])^6)-592926 sin[halfphi] sin[t]+1713220 ((sin[halfphi])^3) sin[t]-1346928 ((sin[halfphi])^5) sin[t]-617309 ((sin[t])^2)+1764300 ((sin[halfphi])^2) ((sin[t])^2)-1030416 ((sin[halfphi])^4) ((sin[t])^2)+52084 sin[halfphi] ((sin[t])^3)-85620 ((sin[t])^4))==0},sin[t],sin[halfphi]]

Factor[%]

RootReduce/@Last/@Reduce[(-675+4482 ((sin[t])^2)-12267 ((sin[t])^4)+22580 ((sin[t])^6)-32756 ((sin[t])^8)+31744 ((sin[t])^10)-16720 ((sin[t])^12)+3136 ((sin[t])^14)+576 ((sin[t])^16))==0&&sin[t]>0,sin[t]]

N[%]

r1=Root[-675+4482 ((#1)^2)-12267 ((#1)^4)+22580 ((#1)^6)-32756 ((#1)^8)+31744 ((#1)^10)-16720 ((#1)^12)+3136 ((#1)^14)+576 ((#1)^16)&,4];

t=ArcSin[r1]

-675+4482 ((#1)^2)-12267 ((#1)^4)+22580 ((#1)^6)-32756 ((#1)^8)+31744 ((#1)^10)-16720 ((#1)^12)+3136 ((#1)^14)+576 ((#1)^16)&[x]//TraditionalForm

N[%,20]

Factor[GroebnerBasis[{(38484 sin[halfphi]-199557 ((sin[halfphi])^3)+168228 ((sin[halfphi])^5)-47244 sin[t]+66171 ((sin[halfphi])^2) sin[t]-35628 ((sin[halfphi])^4) sin[t]+24573 sin[halfphi] ((sin[t])^2)+14036 ((sin[halfphi])^3) ((sin[t])^2)+85037 ((sin[t])^3)-113862 ((sin[halfphi])^2) ((sin[t])^3)-11260 sin[halfphi] ((sin[t])^4)+11022 ((sin[t])^5))==0,(329664-1059045 ((sin[halfphi])^2)+680160 ((sin[halfphi])^4)+192816 ((sin[halfphi])^6)-592926 sin[halfphi] sin[t]+1713220 ((sin[halfphi])^3) sin[t]-1346928 ((sin[halfphi])^5) sin[t]-617309 ((sin[t])^2)+1764300 ((sin[halfphi])^2) ((sin[t])^2)-1030416 ((sin[halfphi])^4) ((sin[t])^2)+52084 sin[halfphi] ((sin[t])^3)-85620 ((sin[t])^4))==0},sin[halfphi],sin[t]]]

r2=Select[sin[halfphi]/.{ToRules[RootReduce/@Roots[(-75-22 ((sin[halfphi])^2)+36149 ((sin[halfphi])^4)-435956 ((sin[halfphi])^6)+1931140 ((sin[halfphi])^8)-4241024 ((sin[halfphi])^10)+4995456 ((sin[halfphi])^12)-3032064 ((sin[halfphi])^14)+746496 ((sin[halfphi])^16))==0,sin[halfphi]]]},Positive][[1]]

-75-22 ((#1)^2)+36149 ((#1)^4)-435956 ((#1)^6)+1931140 ((#1)^8)-4241024 ((#1)^10)+4995456 ((#1)^12)-3032064 ((#1)^14)+746496 ((#1)^16)&[x]//TraditionalForm

phi=2ArcSin[r2]

N[Cos[t+phi/2],20]

RootReduce[TrigExpand[Cos[t+phi/2]]]

N[%,20]

625-800 #1+22 ((#1)^2)-1776 ((#1)^3)+1665 ((#1)^4)+1360 ((#1)^5)-3480 ((#1)^6)+2112 ((#1)^7)+1296 ((#1)^8)&[x]//TraditionalForm

Minimize

a[y_]:=80y^2+64y
b[y_]:=416y^3+384y^2+64y
c[y_]:=848y^4+928y^3+352y^2+32y
d[y_]:=768y^5+992y^4+736y^3+288y^2+96y
e[y_]:=256y^6+384y^5+592y^4+480y^3+336y^2+96y+16
f[y_]:=128y^5+192y^4+256y^3+160y^2+96y+32
g[y_]:=64y^2+64y+16

List@@(1||2||3)

NMaximize[{Max[Last/@(List@@Reduce[a[y]x^6-b[y]x^5+c[y]x^4-d[y]x^3+e[y]x^2-f[y]x+g[y]==0&&2y+1>x>(Sqrt[2]),x])],-1<y<1},y]

Limiting Area Ratio

N[r=((3Sqrt[3])/(2π)),20]