# Solving a system of equations with conditions related to the number of solutions

The equation below describes a conic with oblique axis:

$$9 + 22 x + 9 x^2 + 46 y + 24 x y + 16 y^2=0$$

It is a parabola, as the coefficients in $x^2$, $y^2$, $xy$ form a perfect square. To find the coordinates of its vertex, it's either the hard way using algebra only on paper (I summarise that at the end of my entry) or possibly an easy way with Mathematica.
Let me describe my Mathematica approach:

Looking at the $x^2$ and $y^2$ coefficients I know that the slope of the tangent to the parabola at its vertex is $4/3$ so this tangent is of the form $y = 4x/3 + b$ and it intercepts the parabola at its vertex $(x,y)$ exactly one time. When I translate this in Mathematica terms, I have the following system of equations to solve:

 {x, y, b} /.
Solve[{9 + 22 x + 9 x^2 + 46 y + 24 x y + 16 y^2 == 0,
y - (4/3)x - b == 0}, {x, y, b}]


which is not enough of course to find unique values for $(x,y,b)$.
Is it possible to add conditions to the above Solve expressions such as only one solution for $(x,y,b)$ is to be returned (or more exactly as there are squares in the expression) two solutions but identical?.

On paper, I transformed with factors the equation of the parabola. End result put into Mathematica:

eq1 = (3 x + 4 y + 5)^2  - 2 (4 x - 3 y + 8)
eqc = 9 + 22 x + 9 x^2 + 46 y + 24 x y + 16 y^2
eqc == eq1 // FullSimplify
tangentvertex = Reduce[4 x - 3 y + 8 == 0, y]
vertex = {x, y}/.Solve[{eqc == 0, y - 8/3 - (4 x)/3 == 0}, {x, y}] // FullSimplify


Vertex coordinates at last:

{-47/25, 4/25}

-

Here's a higher calculus way: maximize the curvature.

eqn = 9 + 22 x + 9 x^2 + 46 y + 24 x y + 16 y^2;
(* 5/(26 + 9 x^2 + 40 y + 16 y^2 + 6 x (5 + 4 y))^(3/2) *)

Maximize[{curvature, eqn == 0}, {x, y}]
(* {5, {x -> -(47/25), y -> 4/25}} *)


Another solution:

Solve[{eqn == 0,
Pick[eqn, Exponent[#, x] + Exponent[#, y] & /@ List @@ eqn, 2] == 0 /.
Thread[{x, y} -> D[eqn, {{x, y}}]]}, {x, y}]
(* {{x -> -(47/25), y -> 4/25}} *)


where Pick yields the quadratic part:

Pick[eqn, Exponent[#, x] + Exponent[#, y] & /@ List @@ eqn, 2]
(* 9 x^2 + 24 x y + 16 y^2 *)


This works, given that the equation is a parabola, because the quadratic terms factor into $(a\,x+b\,y)^2$, where the vector $(a,b)$ is orthogonal to the tangent line at the vertex.

A third way:

generalForm = n (x - a) - m (y - b) - k (m (x - a) + n (y - b))^2;
Solve[CoefficientList[#, {x, y}] & /@ (eqn == generalForm), {a, b, k, m, n}]
(* {{a -> -(47/25), b -> 4/25, k -> -(1/4), m -> -6, n -> -8}} *)

-
Thanks, I knew about this way, see my note below that I gave to the answer given by J.M. I implemented in MMA the formula such as it is given in the quoted link where the numerator in the curvature expression is a bit different from yours . Your answer is shorter and very elegant. –  Sigis K Apr 21 '13 at 20:50
@SigismondKmiecik There must be a dozen different possible solutions, maybe two. –  Michael E2 Apr 21 '13 at 21:51

Here's a more-or-less mechanical approach:

conic = 9 + 22 x + 9 x^2 + 46 y + 24 x y + 16 y^2;

{co, li, qu} = Normal[CoefficientArrays[conic]];
{9, {22, 46}, {{9, 24}, {0, 16}}}

(* principal axis computation *)
b = qu[[1, 2]]; h = Subtract @@ Diagonal[qu];
{c, s} = Normalize[{1, b/(h + Sign[b] Norm[{h, b}])}]

FullSimplify[conic /. Thread[{x, y} -> {{c, -s}, {s, c}}.{x, y}]]
9 + 25 x (2 + x) + 10 y

(* calculus approach *)
{{c, -s}, {s, c}}.{x, y} /. First[Solve[{% == 0, D[%, x] == 0}, {x, y}]]
{-47/25, 4/25}


Another method:

(* perpendicular to principal axis *)
{cp, sp} = Normalize[{1, -1/(b/(h + Sign[b] Norm[{h, b}]))}]
{4/5, -3/5}

sols = Solve[{conic == 0, cp x + sp y + w == 0}, {x, y}] // FullSimplify;
Union[sols /. First[Solve[Apply[Equal, x /. sols], w]]]
{{x -> -47/25, y -> 4/25}}

-
As a short note, the formula I used for computing the slope of the principal axis of a conic is intimately related to the Jacobi method for diagonalizing a symmetric matrix. –  Guess who it is. Apr 18 '13 at 8:11
Thanks to both of you for these instructive answers. –  Sigis K Apr 20 '13 at 18:00
If you liked my solution, then you'll want to see this book. –  Guess who it is. Apr 20 '13 at 18:17
There is another way to find the apex of this conic, minimizing the ratio of curvature, for an implicit f(x,y) = 0 ) function thus a more general approach. I never found any textbook where this subject is broached in details and it was only on an internet that I found the relevant link - ann.jussieu.fr/~frey/papers/meshing/… - see 3.1 formula based on gradient and hessian matrices. It was quite mechanical to convert it into MMA but without it a pain in the neck to get this way the apex coordinates.SK –  Sigis K Apr 20 '13 at 18:21
Yes, a gradient/Hessian approach should work too, but I think it's a bit messier to set up. –  Guess who it is. Apr 20 '13 at 18:23

I will walk another road, so bear with me; this is something I remembered from analytic geometry. The purpose is to turn your conic section

eq = 9 + 22 x + 9 x^2 + 46 y + 24 x y + 16 y^2


to its normal form. For this, I need to make a change of variables. First of all, I take the second order terms and create a symmetric matrix (I suppose there is a better way to do that):

m = CoefficientArrays[eq, {x, y}][[3]]//Normal


{{9, 24}, {0, 16}}

temp = (m - DiagonalMatrix[Diagonal[m]])/2


{{0, 12}, {0, 0}}

symm = m - temp + Transpose[temp]


{{9, 12}, {12, 16}}

Having constructed this symmetric matrix, I find its eigenvalues and eigenvectors (the eigenvalues are always real because the matrix is symmetric):

eig = Eigensystem[symm]


{{25, 0}, {{3, 4}, {-4, 3}}}

Now we proceed with obtaining the needed rotation in order to straighten the axis:

change = Thread[{x, y} -> Transpose[eig[[2]]].{X, Y} // FullSimplify]


{x -> 3 X - 4 Y, y -> 4 X + 3 Y}

Through them the equation becomes

neweq = eq /. change // FullSimplify


9 + 125 X (2 + 5 X) + 50 Y

Now we look for the translations in order to put the vertex on $(0,0)$:

Solve[Coefficient[neweq /. X -> X + A, X] == 0, A][[1, 1]]


A -> -(1/5)

neweq2 = neweq /. X -> X - 1/5 // Simplify


-16 + 625 X^2 + 50 Y

For $Y$, it is even more obvious:

neweq2 = neweq2 /. Y -> Y + 16/50 // FullSimplify


625 X^2 + 50 Y

and this is the normal form of your conic. To find the vertex, we transform back the point $(0,0)$ to your initial coordinates:

change /. {X -> X - 1/5, Y -> Y + 16/50} /. {X -> 0, Y -> 0}


{x -> -(47/25), y -> 4/25}

-