# How can I reduce my solve (the number of the planes tangent to the three spheres)?

Let's be three spheres $(S_1)$ has centre $A(1,2,1)$, radius $R=2$, $(S_2)$ has centre $B(3, -1, 1)$, radius $R=1$ and $(S_3)$ has centre $C(-1, -1, 1)$, radius $R=1$. How many planes tangent to the three spheres?

I suppose the equation of the plane has the form $a x + by + c z + d=0$, where $a^2 + b^2 + c^2 = 1.$

The distance from the point $A$ to the plane is $|a + 2b + c+d|$. The plane tangent to the three sphere $(S_1)$ if and only if $|a + 2b + c+d|=2$. Similarly to the spheres $(S_2)$ and $(S_3)$. My code

f[x_, y_, z_] := Abs[a x + b y + c z + d];
Solve[{f[1, 2, 1] == 2, f[3, -1, 1] == 1, f[-1, -1, 1] == 1,
a^2 + b^2 + c^2 == 1}, {a, b, c, d}, Reals]


I got

{{a -> -(1/2), b -> -(2/3), c -> -(Sqrt[11]/6), d -> 1/6 (-1 + Sqrt[11])}, {a -> -(1/2), b -> -(2/3), c -> Sqrt[11]/6, d -> 1/6 (-1 - Sqrt[11])}, {a -> -(1/2), b -> 2/3, c -> -(Sqrt[11]/6), d -> 1/6 (7 + Sqrt[11])}, {a -> -(1/2), b -> 2/3, c -> Sqrt[11]/6, d -> 1/6 (7 - Sqrt[11])}, {a -> 0, b -> -1, c -> 0, d -> 0}, {a -> 0, b -> -(1/3), c -> -((2 Sqrt[2])/3), d -> 2/3 (-2 + Sqrt[2])}, {a -> 0, b -> -(1/3), c -> (2 Sqrt[2])/3, d -> -(2/3) (2 + Sqrt[2])}, {a -> 0, b -> 1/3, c -> -((2 Sqrt[2])/3), d -> 2/3 (2 + Sqrt[2])}, {a -> 0, b -> 1/3, c -> (2 Sqrt[2])/3, d -> -(2/3) (-2 + Sqrt[2])}, {a -> 0, b -> 1, c -> 0, d -> 0}, {a -> 1/2, b -> -(2/3), c -> -(Sqrt[11]/6), d -> 1/6 (-7 + Sqrt[11])}, {a -> 1/2, b -> -(2/3), c -> Sqrt[11]/6, d -> 1/6 (-7 - Sqrt[11])}, {a -> 1/2, b -> 2/3, c -> -(Sqrt[11]/6), d -> 1/6 (1 + Sqrt[11])}, {a -> 1/2, b -> 2/3, c -> Sqrt[11]/6, d -> 1/6 (1 - Sqrt[11])}

Now I count by hand and DeleteCases.

• If you just want to know the number of planes, can't you work on topological grounds? Assuming the spheres are far enough apart, we have: 1 above all 3 spheres; 3 where the plane is above 2 of the spheres and below the 3rd; 3 where the plane is above 1 of the spheres and below the other 2; and finally, 1 below all 3 spheres. Giving us 8 spheres. Commented Jan 27, 2018 at 0:54
• @CarlWoll Only 7 spheres. Commented Jan 27, 2018 at 1:22
• Oops, I meant 8 planes. Commented Jan 27, 2018 at 1:37
• Sorry. 7 planes. Commented Jan 27, 2018 at 1:54
• In general there are 8 planes. In your example 2 of them are degenerate. Commented Jan 27, 2018 at 2:02

I think the equations are not quite what is wanted.

For each sphere_j there is a point {x[j],y[j],z[j]} that must satisfy the equations that (i) it is on that particular sphere (ii) it is on the plane (iii) the gradient of the sphere, at that point, is parallel to the plane direction vector {a,b,c}. We also use a normalization equation for the plane parameters. The code below will set up this system.

centers = {{1, 2, 1}, {3, -1, 1}, {-1, -1, 1}};

plane[x_, y_, z_] := {a, b, c, d}.{x, y, z, -1}
spheres =
MapIndexed[
With[{coords = {x[#2[[1]]], y[#2[[1]]], z[#2[[1]]]} - #[[1]]},
coords.coords - #[[2]]^2] &, sphereVals];
spheres];
planeHitsSpherePolys = Map[plane[x[#], y[#], z[#]] &, Range[3]];
planeTangSpherePolys = Flatten[MapIndexed[#1 - k[#2[[1]]]*{a, b, c} &, grads]];


So here are the equations.

allPolys =
Join[spheres, planeHitsSpherePolys,
planeTangSpherePolys, {a^2 + b^2 + c^2 + d^2 - 1}]

(* {-4 + (-1 + x[1])^2 + (-2 + y[1])^2 + (-1 +
z[1])^2, -1 + (-3 + x[2])^2 + (1 + y[2])^2 + (-1 +
z[2])^2, -1 + (1 + x[3])^2 + (1 + y[3])^2 + (-1 + z[3])^2, -d +
a x[1] + b y[1] + c z[1], -d + a x[2] + b y[2] + c z[2], -d +
a x[3] + b y[3] + c z[3], -a k[1] + 2 (-1 + x[1]), -b k[1] +
2 (-2 + y[1]), -c k[1] + 2 (-1 + z[1]), -a k[2] +
2 (-3 + x[2]), -b k[2] + 2 (1 + y[2]), -c k[2] +
2 (-1 + z[2]), -a k[3] + 2 (1 + x[3]), -b k[3] +
2 (1 + y[3]), -c k[3] + 2 (-1 + z[3]), -1 + a^2 + b^2 + c^2 + d^2} *)


We don't really care about the ratios between direction vector and gradients, and from the wording of the question I will guess we also do not require the three intersection points.

solns = {a, b, c, d} /. Solve[allPolys == 0];


Numerically:

N[solns]

(* {{0., -1., 0., 0.}, {0., 1., 0., 0.}, {0.,
0.134077, -0.379228, -0.915538}, {0., -0.134077, 0.379228,
0.915538}, {0., -0.310496, -0.878216, 0.363769}, {0., 0.310496,
0.878216, -0.363769}, {-0.251372,
0.335162, -0.277902, -0.864436}, {0.251372, -0.335162, 0.277902,
0.864436}, {0.405875,
0.541167, -0.448712, -0.584004}, {-0.405875, -0.541167, 0.448712,
0.584004}, {0.426112, -0.568149, -0.471085, 0.523177}, {-0.426112,
0.568149,
0.471085, -0.523177}, {-0.46644, -0.62192, -0.515669, -0.360189}, \
{0.46644, 0.62192, 0.515669, 0.360189}} *)


These come in pairs since the negative of a solution is also a solution (giving the same plane of course). So taking that into account brings it to seven such planes, as was claimed in the original post. Also we can renormalize and simplify. The rationale for renormalizing is that maybe the original norm caused an excess of radicals, and dividing each solution by the same nonzero component could give a better form.

possolns = Select[solns, #[[2]] >= 0 &];
simpler = Simplify[Map[#/Max[Abs[#]] &, possolns]]

(* Out[79]= {{0, 1, 0, 0}, {0, 1/(4 + 2 Sqrt[2]), 1 - Sqrt[2], -1}, {0,
1/(2 Sqrt[2]), 1, 1 - Sqrt[2]}, {-(3/(7 + Sqrt[11])), 4/(
7 + Sqrt[11]), 1/38 (11 - 7 Sqrt[11]), -1}, {3/(1 + Sqrt[11]), 4/(
1 + Sqrt[11]), 1/10 (-11 + Sqrt[11]), -1}, {-(3/4), 1, Sqrt[11]/4,
1/4 (-7 + Sqrt[11])}, {3/4, 1, Sqrt[11]/4, 1/4 (-1 + Sqrt[11])}} *)


Someone with better graphics chops than myself might be able to show all these planes along with the spheres.

• Lichblau With my way, the equation y=0 is an answer. Your way is not. But It is easy to see that y = 0 is correct. Commented Jan 27, 2018 at 0:02
• There are 7 planes in my problem. Not eight. Commented Jan 27, 2018 at 3:42
• Thanks for pointing that out. There was an error in the gradient equations, now fixed. Commented Jan 27, 2018 at 15:41