# Is there a way to solve the Apollonius Circle problem in Mathematica?

Assuming x, a, and b are complex numbers, is there a way to reduce the equation Abs[x - a] == k Abs[x - b] to something like Abs[x - ...] = ...?

Or if it can't be done by complex number, can it be solved by specified the coordinate explicitely? I've try something like this:

d[x_, y_, xx_, yy_] := (x - xx)^2  + (y - yy)^2
d[x, y, x1, y1] - k d[x, y, x2, y2] == 0 // Expand


and got:

x^2 - k x^2 - 2 x x1 + x1^2 + 2 k x x2 - k x2^2 + y^2 - k y^2 - 2 y y1 + y1^2 + 2 k y y2 - k y2^2 == 0


but how can I transform the above equation to some form like (x-...)^2 + (y-...)^2 = ... by mathematica?

-
Solve[Abs[x - a] - k Abs[x - b] == 0, x] ? – Dr. belisarius Sep 24 '12 at 11:55
@belisarius No, I want something like Abs[x - ...] = .... But the Solve treat x as Real not Complex. – chyx Sep 24 '12 at 11:58
Try Solve[Abs[x - (ar + I ai)] - k Abs[x - (br + I bi)] == 0, x] ... why do you say that Solve treat them as Reals? – Dr. belisarius Sep 24 '12 at 12:01
It gives result:{{x -> (I (ai - I ar + bi k - I br k))/(1 + k)}, {x -> ( I (-ai + I ar + bi k - I br k))/(-1 + k)}}. But I want to get the centre and radius of $x$. – chyx Sep 24 '12 at 12:03
@Artes: I understood the OP as asking "given $k$, $a$, and $b$ in $|z-a|=k|z-b|$, how can I use Mathematica to find the center and radius of the Apollonian circle?" – J. M. Sep 24 '12 at 13:29

Thanks to @Artes and after glancing http://stackoverflow.com/questions/8462244/controlling-measure-zero-sets-of-solutions-with-manipulate-a-case-study (the math seems too hard for me...). I do it with more details to mathematica and finally got the result.

radius[x_] := Sqrt[x . x]
poly = (radius[{x, y} - {x1, y1}]^2 - k^2 radius[{x, y} - {x2, y2}]^2 )/ (1 - k^2)
{x0, y0} = -Coefficient[D[poly, #]/2, #, 0] & /@ {x, y} // FullSimplify
r = -(poly - (D[poly, x]/2)^2 - (D[poly, y]/2)^2) // FullSimplify


And the result reveals to be $|z - \frac{a - k^2 b}{1 - k^2}| = \frac{|a - b||k|}{|1 - k^2|}$

-
This approach seems to assume (much) more than was given in the question. A truly satisfactory solution would not presuppose the solution is a circle, for instance (or even a conic section): that ought to emerge from the answer. – whuber Sep 24 '12 at 16:48
Yes, I also want that kind of answer. An automatic way to transform the quadratic polynomial to conic section would be cool. – chyx Sep 25 '12 at 2:41

Well, I don't really have the time to explain this right now and I did it so long ago, that I'm not sure that I could anyway. The ideas are explained here: http://mathworld.wolfram.com/BowlofIntegers.html

d[A_, B_] := Sqrt[(A - B).(A - B)];
tangentCircle[{Circle[{x1_, y1_}, r1_],
Circle[{x2_, y2_}, r2_], Circle[{x3_, y3_}, r3_]}] :=

Module[{sols, x, y, r},
sols = NSolve[{
(x - x1)^2 + (y - y1)^2 == (r - r1)^2,
(x - x2)^2 + (y - y2)^2 == (r + r2)^2,
(x - x3)^2 + (y - y3)^2 == (r + r3)^2}, {x, y, r}];
Circle[{x, y}, r] /. sols] /;
r1 == r2 + r3  &&  r1 > d[{x1, y1}, {x2, y2}];
tangentCircle[{Circle[{x1_, y1_}, r1_],
Circle[{x2_, y2_}, r2_], Circle[{x3_, y3_}, r3_]}] :=

Module[{a2, b2, c2, d2, a3, b3, c3, d3,
x, y, r, sign},
If[r1 > d[{x1, y1}, {x2, y2}], sign = -1, sign = 1];
a2 = 2 (x1 - x2); b2 = 2 (y1 - y2); c2 = 2 (sign r1 - r2);
d2 = (x1^2 + y1^2 - r1^2) - (x2^2 + y2^2 - r2^2);
a3 = 2 (x1 - x3); b3 = 2 (y1 - y3); c3 = 2 (sign r1 - r3);
d3 = (x1^2 + y1^2 - r1^2) - (x3^2 + y3^2 - r3^2);
x = (b3 d2 - b2 d3 - b3 c2 r + b2 c3 r)/(a2 b3 - b2 a3) // N;
y = (-a3 d2 + a2 d3 + a3 c2 r - a2 c3 r)/(a2 b3 - a3 b2) // N;
r = Min[
Abs[r /. Solve[(x - x1)^2 + (y - y1)^2 == (r + sign r1)^2, r]]];
Circle[{x, y}, r]];
tangentCircles[{Circle[{x1_, y1_}, r1_],
Circle[{x2_, y2_}, r2_], Circle[{x3_, y3_}, r3_]}] :=

Module[{a2, b2, c2, d2, a3, b3, c3, d3,
x, y, r, sign},
If[r1 > d[{x1, y1}, {x2, y2}], sign = -1, sign = 1];
a2 = 2 (x1 - x2); b2 = 2 (y1 - y2); c2 = 2 (sign r1 - r2);
d2 = (x1^2 + y1^2 - r1^2) - (x2^2 + y2^2 - r2^2);
a3 = 2 (x1 - x3); b3 = 2 (y1 - y3); c3 = 2 (sign r1 - r3);
d3 = (x1^2 + y1^2 - r1^2) - (x3^2 + y3^2 - r3^2);
x = (b3 d2 - b2 d3 - b3 c2 r + b2 c3 r)/(a2 b3 - b2 a3) // N;
y = (-a3 d2 + a2 d3 + a3 c2 r - a2 c3 r)/(a2 b3 - a3 b2) // N;
Circle[{x, y}, Abs[r]] /.
Solve[(x - x1)^2 + (y - y1)^2 == (r + sign r1)^2, r]
];
triplets[{a_, b_, c_}, d_] :=
{{a, b, d}, {a, c, d}, {b, c, d}};
triplets[{a_, b_, c_Circle}, l_List] :=
triplets[{a, b, c}, #] & /@ l;
apollonianStep[{a_Circle, b_, c_}] :=
triplets[{a, b, c}, tangentCircle[{a, b, c}]];
apollonianStep[l_List] := apollonianStep /@ l;

{circ1, circ2, circ3, circ4} = {
Circle[{0, 0}, 1/6],
Circle[{1/6 - 1/11, 0}, 1/11],
Circle[{-8/105, -2/35}, 1/14],
Circle[{-3/50, 2/25}, 1/15]};
init = {{circ1, circ2, circ3}, {circ1, circ2, circ4},
{circ1, circ3, circ4}, {circ2, circ3, circ4}};

circs = TimeConstrained[
Union[Flatten[init //. {a_, b_, Circle[p_, r_]} :>
apollonianStep[{a, b, Circle[p, r]}] /; r > .01]], 10];
moreCircs = Union[Flatten[init //. {a_, b_, Circle[p_, r_]} :>
apollonianStep[{a, b, Circle[p, r]}] /; r > .0005]];

number[Circle[c_, r_]] :=
Text[Style[Round[1/r], FontSize -> 1400 r], c];
Graphics[{circs, moreCircs,
Map[number, DeleteCases[circs, Circle[{0, 0}, _]]]},
ImageSize -> 700]


-
...these are different Apollonian circles, I think. – J. M. Sep 24 '12 at 16:21
@J.M. Right, next time I'll read the question. – Mark McClure Sep 24 '12 at 18:39