2

I have written a function to get the canonical form of a 3D quadric surface. This involves finding a unitary transformation of the coordinates which eliminates the cross-terms such as $xy$ or $xz$. For example, let $q(x,y,z)=2xy-x+y+z-1$. Then by a rotation of axes, this equation is converted to $y^2- x^2+ \sqrt{2}x+z-1$.

The code is as follows

canonicalForm[polynomial_] := Module[{a, b, c, svd, p, v, transform},
  v = Normal@ CoefficientArrays[polynomial];
  a = Last[v];
  b = v[[2]];
  c = First[v]; (** form of x\[Transpose].a.x + b.x + c = 0 **)
  a = (a + Transpose[a])/2;           (** make it symmetric **)
  svd = SingularValueDecomposition[a];
  If[Equal @@ Drop[svd, {2}],
    p = First[svd],
    p = Normalize/@ Last@ Eigensystem[a]; (* THIS IS LINE 10 *)
    If[Expand[p.Transpose@ p] != IdentityMatrix[3],
      p = Chop@ First@ SchurDecomposition@ N@ a]
  ];
  transform = p.{{x}, {y}, {z}};
  a = Simplify@Flatten[Transpose[transform].a.transform];
  b = Together[b.transform];
  Chop@ First[a + b + c]]

Three equations are used to test the code: $$\begin{align} q_1&=3 x^2-2 x y+2 x z+5 y^2-2 y z+3 z^2+2 z\\ q_2&=2 x y-\sqrt{2} x-\sqrt{2} y+z-1\\ q_3&=y z-x y+z^2+z-2 \end{align}$$ The code works perfectly for $q_1$ and $q_2$, but for $q_3$ it takes a long time and then gives a very lengthy output in terms of the roots of some polynomials.

canonicalForm[3x^2 -2x y +2x z +5y^2 -2y z +3z^2 +2z]

$6 x^2+3 y^2+2 z^2+\frac{1}{3} \left(\sqrt{6} x+2 \sqrt{3} y+3 \sqrt{2} z\right)$

canonicalForm[2x y -Sqrt[2]x -Sqrt[2]y +z -1]

$-x^2+y^2-2 y+z-1$

canonicalForm[y z -x y +z^2 +z -2]

a vary lengthy output which is also mathematically wrong

For $q_1$ singular value decomposition results in a unitary transform. For $q_2$ the normalized eigenvectors give the desired unitary transform. But for $q_3$ none of these works and the code switches to SchurDecomposition, which only accepts approximate numerical values.

Now, if I skip calculating the normalized eigenvectors and replace the 10th, 11th and 12th lines with this:

p = Chop@ First@ SchurDecomposition@ N@ a

The code will produce the desired correct results, but in approximate numerical form.

canonicalForm[y z -x y +z^2 +z -2]

$-0.585043 x^2+0.233192 x+0.344446 y^2+0.397113 y+1.2406 z^2+0.88765 z-2$

So it seems calculating the eigenvectors before using SchurDecomposition somehow affects its performance. I wasn't able to suppress this error by defining temporary variables and any other method that I could think of. The only way is deleting those three lines and directly using SchurDecomposition.

Questions

  1. Is there any way other than circumventing eigenvectors, that I could make this work? Is this a ?
  2. Is there any other method that results in a unitary transform and gives the exact values, rather than approximate numerical ones given by SchurDecomposition?
polfosol
  • 952
  • 5
  • 20
  • If the matrix is symmetric (after all, you are dealing with quadratic forms), then you shouldn't need SchurDecomposition[]. – J. M.'s missing motivation Apr 21 '20 at 13:16
  • @J.M. that was another ironic issue. Sometimes it seems neither SVD nor eigenvectors give the desired result for a symmetric matrix. – polfosol Apr 21 '20 at 13:18
  • I just "plugged" the polynomial 1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2 into canonicalForm[], and lo and behold!, I got 1. What does this mean? (I'm somewhat of a novice in this area.) See my preceding questions https://math.stackexchange.com/questions/3660652/to-which-of-the-seventeen-standard-quadrics-do-these-two-equations-reduce and https://mathoverflow.net/questions/359459/interpret-certain-expressions-in-terms-of-classical-quadratic-surfaces. I used Q_1=x, Q_2= y, Q_3 =z. I'll now try the code of J. M. in his answer. – Paul B. Slater May 06 '20 at 00:00

1 Answers1

3

Writing a robust routine for dealing with quadric surfaces (and conic sections for that matter) is a lot of work, so I'll give a skeleton from which I hope the OP (or other more motivated users) can build on further.

I'll take the first quadric as an example:

quad = 3 x^2 - 2 x y + 2 x z + 5 y^2 - 2 y z + 3 z^2 + 2 z;

The first important step is to ensure that CoefficientArrays[] returns a symmetric matrix:

coefs = Normal[CoefficientArrays[quad, Variables[quad], "Symmetric" -> True]]
   {0, {0, 0, 2}, {{3, -1, 1}, {-1, 5, -1}, {1, -1, 3}}}

Then, one should also recall that Eigensystem does not yield orthonormal eigenvectors by default, so one has to do some extra work:

{vals, vecs} = MapAt[Orthogonalize, Eigensystem[coefs[[3]]], {2}]
   {{6, 3, 2}, {{1/Sqrt[6], -Sqrt[(2/3)], 1/Sqrt[6]},
                {1/Sqrt[3], 1/Sqrt[3], 1/Sqrt[3]},
                {-(1/Sqrt[2]), 0, 1/Sqrt[2]}}}

Now that you have an orthogonal matrix, you can do this:

quad /. Thread[{x, y, z} -> Transpose[vecs].{xt, yt, zt}] // Simplify
   Sqrt[2/3] xt + 6 xt^2 + (2 yt)/Sqrt[3] + 3 yt^2 + zt (Sqrt[2] + 2 zt)

which should now yield a result that is amenable to completing the square.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
  • Your method is indeed simpler and more efficient than mine. But still, the issue with SchurDecomposition is unaddressed, and furthermore, it still doesn't work on the third equation. – polfosol Apr 21 '20 at 14:33
  • I just tried it on the third equation; you need to do MapAt[FullSimplify @* Orthogonalize, Eigensystem[coefs[[3]]], {2}] instead, and quad /. Thread[{x, y, z} -> Transpose[vecs].{xt, yt, zt}] ought to give a result that has Root[] objects in the coefficients, but is otherwise free of cross-terms. My point anyway was that you don't need to use SchurDecomposition[] here, unless your coefficients are already inexact numbers. – J. M.'s missing motivation Apr 21 '20 at 14:39
  • Per my just posted comment to the question itself, I took quad =1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2, and fed it through the three indicated steps and got 1 - ((61 + 13 Sqrt[61]) Sqrt[3538 + 272 Sqrt[61]] xt)/1098 + 3/2 (14 + Sqrt[61]) xt^2 + ( Sqrt[3538 - 272 Sqrt[61]] (61 - 13 Sqrt[61]) yt)/1098 - 3/2 (-14 + Sqrt[61]) yt^2 and Sqrt[2/3] xt + 6 xt^2 + (2 yt)/Sqrt[3] + 3 yt^2 + zt (Sqrt[2] + 2 zt) Clear? Signifigance? – Paul B. Slater May 06 '20 at 00:08
  • This is why I said writing a robust routine is hard work, @Paul. A good routine ought to be able to detect degenerate cases like yours (cf. the criteria listed here and many other references). – J. M.'s missing motivation May 06 '20 at 00:41
  • Thanks--and thanks for the quadric surfaces link also, J. M. But maybe you could further "spoon-feed" me as to the nature of the degeneracy? What conclusions can one make as to the nature of quad =1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2? What transformation can be employed to reduce it to a presumably simpler form? Actually, as the previous links of mine indicate, I'm primarily interested in the inequality (!) constraint $ x^2+3 x y +\left(3 y+z\right){}^2<3 y+2 x z$. . So, what can I conclude as to the volume that the quadric surface encloses. – Paul B. Slater May 06 '20 at 01:44
  • @Paul, if I proceed with the analysis in the link I gave, your "quadric" is actually a straight line, which of course does not enclose any volume. – J. M.'s missing motivation May 06 '20 at 02:32
  • OK, thanks! So, what's the equation of the line? And can we say something about the surface, quad <1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2? – Paul B. Slater May 06 '20 at 04:02
  • @Paul, since I already told you that it's a line, you can use FindInstance[] to generate two points on it (say, $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$), and then consider the line with the parametric equations $((1-t)x_1+t x_2,(1-t)y_1+t y_2,(1-t)z_1+t z_2)$. – J. M.'s missing motivation May 06 '20 at 04:08