# 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. – Carl Woll Jan 27 '18 at 0:54
• @CarlWoll Only 7 spheres. – minhthien_2016 Jan 27 '18 at 1:22
• Oops, I meant 8 planes. – Carl Woll Jan 27 '18 at 1:37
• Sorry. 7 planes. – minhthien_2016 Jan 27 '18 at 1:54
• In general there are 8 planes. In your example 2 of them are degenerate. – Carl Woll Jan 27 '18 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. – minhthien_2016 Jan 27 '18 at 0:02
• There are 7 planes in my problem. Not eight. – minhthien_2016 Jan 27 '18 at 3:42
• Thanks for pointing that out. There was an error in the gradient equations, now fixed. – Daniel Lichtblau Jan 27 '18 at 15:41