Apollonius Spheres

Consider three points in the plane, each being the midpoint of a circle. Let the circles around these three points have the property that they touch each other pairwise. Find the coordinates of the midpoint and the radius of a fourth circle which touches all of the three original circles. This is the famous Apollonius circle problem.[Soddy 1936, Coxeter 1968]. For a solution of this problem using Mathematica see [Trott 1997].

Here we will treat a generalization of the problem in three dimensions. Lets us start with four points in (ℝ^3) arranged in such a manner that spheres around these points touch each other. (Between four points there are six different distances, so the described constellation is not generically possible.) We will calculate the radius and the coordinates of a fifth point such that a sphere around this fifth point touches the original four spheres (in case such a constellation exists).

Calculating the Radius

This is the set of equations under consideration:

\(\(ℰ = {
        \((x - x\_1)\)\^2 + \((y - y\_1)\)\^2 + \((z - z\_1)\)\^2 - 
          \((r\_1 + r)\)\^2, 
		\ \ \ 
        \((x - x\_2)\)\^2 + \((y - y\_2)\)\^2 + \((z - z\_2)\)\^2 - 
          \((r\_2 + r)\)\^2, 
	\ \ \ \ \ 
        \((x - x\_3)\)\^2 + \((y - y\_3)\)\^2 + \((z - z\_3)\)\^2 - 
          \((r\_3 + r)\)\^2, \ \ \ 
			\ 
        \((x - x\_4)\)\^2 + \((y - y\_4)\)\^2 + \((z - z\_4)\)\^2 - 
          \((r\_4 + r)\)\^2, \ 
			\ 
        \((x\_1 - x\_2)\)\^2 + \((y\_1 - y\_2)\)\^2 + \((z\_1 - z\_2)\)\^2 - 
          \((r\_1 + r\_2)\)\^2, 
		\ \ \ 
        \((x\_1 - x\_3)\)\^2 + \((y\_1 - y\_3)\)\^2 + \((z\_1 - z\_3)\)\^2 - 
          \((r\_1 + r\_3)\)\^2, 
		\ \ \ 
        \((x\_1 - x\_4)\)\^2 + \((y\_1 - y\_4)\)\^2 + \((z\_1 - z\_4)\)\^2 - 
          \((r\_1 + r\_4)\)\^2, 
		\ \ \ 
        \((x\_2 - x\_3)\)\^2 + \((y\_2 - y\_3)\)\^2 + \((z\_2 - z\_3)\)\^2 - 
          \((r\_2 + r\_3)\)\^2, 
		\ \ \ 
        \((x\_2 - x\_4)\)\^2 + \((y\_2 - y\_4)\)\^2 + \((z\_2 - z\_4)\)\^2 - 
          \((r\_2 + r\_4)\)\^2, 
		\ \ \ 
        \((x\_3 - x\_4)\)\^2 + \((y\_3 - y\_4)\)\^2 + \((z\_3 - z\_4)\)\^2 - 
          \((r\_3 + r\_4)\)\^2}; \)\)

{(Subscript[x, i]), (Subscript[y, i]), (Subscript[z, i])} is the midpoint if the i-th sphere with radius (Subscript[r, i]). {x, y, z} is the radius of the fifth sphere with radius r.The first four equations express the fact that the fifth sphere exactly touches the first four. The remaining six equations describe the fact that the original four spheres all touch each other pairwise.

Without loss of generality we assume the midpoint of the first sphere to be at the origin, the midpoint of the second sphere on the x-axis and the midpoint of the third sphere in the x-y\(\*"" plane"\).

\(x\_1 = 0; y\_1 = 0; z\_1 = 0; 
y\_2 = 0; z\_2 = 0; 
z\_3 = 0; \)

We will calculate the radius of the fifth sphere as a function of the radii of the four original spheres. Because our defining set ℰ of equations is purely polynomial, we make use of Mathematica's enhanced GroebnerBasis command.

If the radius of the fifth sphere is known, the calculation of its midpoint  is more or less straightforward, which is why we derive an equation for this radius r of the fifth sphere. As a first step, we eliminate the midpoint coordinates of the fifth sphere from the first four of the above equations.

ℊ𝒷 = 
      GroebnerBasis[Take[ℰ, 4], 
         {(Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 3]), r, (Subscript[x, 2]), (Subscript[x, 3]), (Subscript[y, 4]), (Subscript[x, 4]), (Subscript[y, 4]), (Subscript[z, 4])}, {x, y, 
          z}, MonomialOrder -> EliminationOrder];

The resulting polynomial contains the known data from the original four spheres and the unknown radius of the fifth equation.

Short[ℊ𝒷,6]

Using the last six of the above equations ℰ allows us to eliminate the midpoint coordinates of the first four equations in favour of their radii, again using the GroebnerBasis command In the following calculation most of the nontrivial work is done. We also use PolynomialReduce (which requires for its second arguments the polynomials in the form of a Gröbner basis).

ℯ = 
      Collect[Last[
          PolynomialReduce[First[ℊ𝒷], 
            GroebnerBasis[
              Take[ℰ, {5, 10}], {(Subscript[x, 2]), (Subscript[x, 3]), (Subscript[y, 3]), (Subscript[x, 4]), 
                (Subscript[y, 4]), (Subscript[z, 4]), (Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 3]), (Subscript[r, 4]), r}], {(Subscript[x, 2]), (Subscript[x, 3]), (Subscript[y, 3]), 
              (Subscript[x, 4]), (Subscript[y, 4]), (Subscript[z, 4]), (Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 3]), (Subscript[r, 4]), r}]], r, Simplify]

This quadratic equation in r is the equation we looked for and a generalization of Söddy's formula of the corresponding two-dimensional problem.

An Example

Let as demonstrate the formula using an explicit example.
Here are the midpoints of the given four spheres.

\(p\_1 = {0, 0, 0}; 
p\_2 = {1, 0, 0}; 
p\_3 = {0, 1, 0}; 

    p\_4 = {0, 0, 1}; \)

The radii of these four spheres at (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]), and (Subscript[p, 4])are:

\(ρ = 
      Simplify[Solve[

          			{r\_1 + r\_2 == 
              \@\(\((p\_1 - p\_2)\).\((p\_1 - p\_2)\)\), 
            r\_1 + r\_3 == \@\(\((p\_1 - p\_3)\).\((p\_1 - p\_3)\)\), 

            			\ r\_1 + r\_4 == 
              \@\(\((p\_1 - p\_4)\).\((p\_1 - p\_4)\)\), 
            r\_2 + r\_3 == \@\(\((p\_2 - p\_3)\).\((p\_2 - p\_3)\)\), 

            				r\_2 + r\_4 == 
              \@\(\((p\_2 - p\_4)\).\((p\_2 - p\_4)\)\), 
            r\_3 + r\_4 == \@\(\((p\_3 - p\_4)\).\((p\_3 - p\_4)\)\)}, 
          
			{r\_1, r\_2, r\_3, r\_4}]]\)

Using the above equation ℯ  we get for the radius r of the fifth sphere:

𝓇 = 
      First[Select[Solve[ℯ == 0 /. ρ, r], 
          (r /. #) > 0&]]

Using this result we can calculate the coordinates of the midpoint of the fifth sphere.

(Subscript[p, 5]) = {x, y, z};

\(\(\ xyz = 
      Solve[\({\((p\_5 - p\_1)\).\((p\_5 - p\_1)\) - \((r\_1 + r)\)\^2 
                  == 0, 
					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                \((p\_5 - p\_2)\).\((p\_5 - p\_2)\) - \((r\_2 + r)\)\^2 
                  == 0, \ \ \ \ \ \ \ \ \ \ \ 

                							\ \ \ \ \ \ \ \ \ \ \ 
                \((p\_5 - p\_3)\).\((p\_5 - p\_3)\) - \((r\_3 + r)\)\^2 
                  == 0} /. ρ\) /. 𝓇, 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {x, y, z}] // Simplify\)\)

The second of these solutions is the one contained in the tetrahedron formed by the midpoints of the original four spheres.

Finally let us construct a picture of the five spheres. The four original spheres are shown in red and the fifth sphere in green.

Needs["Graphics`Shapes`"]

\(\(Show[
      Graphics3D[{{SurfaceColor[RGBColor[1, 0, 0], RGBColor[1, 0, 0]], 
            TranslateShape[
              Sphere[r\_1 /. 
                  ρ[[1]], 40, 30], 
              p\_1], TranslateShape[
              Sphere[r\_2 /. 
                  ρ[[1]], 40, 30], 
              p\_2], TranslateShape[
              Sphere[r\_3 /. 
                  ρ[[1]], 40, 30], 
              p\_3], TranslateShape[
              Sphere[r\_4 /. 
                  ρ[[1]], 40, 30], 
              p\_4]}, {SurfaceColor[RGBColor[0, 1, 0], RGBColor[0, 1, 0]], 

            				TranslateShape[
              Sphere[r /. 𝓇, 20, 15], {x, y, z} /. 
                xyz[[2]]]}}], 
      ViewPoint -> {0, \(-2\), 0.35}, \ Boxed\  -> \ False]; \)\)

Making holes in the polygons makes it easier to see how the fifth sphere fits (we also change the viewpoint in comparision to the last picture):

\(holePolygons[𝓅_Polygon, f_]\  := 

      \(Function[m, 
          \(Function[n, 
              
{EdgeForm[], \ Thickness[0.0001], 

                MapThread[Polygon[Join[#1, Reverse[#2]]]&, \ 

                  					{
                    Partition[\(Append[#, First[#]]&\)[First[𝓅]], 2, 
                      1], 
					\ 
                    Partition[\(Append[#, First[#]]&\)[n], 2, 1]}], 

                \(Line[Append[#, First[#]]]&\)[n]}]\)[\ \ 
            \(\((m + f \((# - m)\))\)&\)/@\ First[𝓅]]]\)[
        Plus\ @@\ First[𝓅]/Length[First[𝓅]]]\)

Show[%% /. 
        𝓅_Polygon :> holePolygons[𝓅, 0.8], 
      ViewPoint -> {2, 2.3, 1.7}, Boxed -> False];