Parabolas, numerically

Recap, and background

Two posts ago we showed how, given four points in the plane in general position, but with a few restrictions, it was possible to find two parabolas through those points. We used computer algebra.

The steps were:

1. Create four equations

   \[ (Ax_i+By_i)^2+Cx_i+Dy_i+E=0 \]

   for each of the four \((x_i,y_i)\) coordinates.

2. Solve the first three equations for \(C\), \(D\) and \(E\): the results will be expressions in \(A\) and \(B\).

3. Substitute the values from the previous step into the last equation and solve for \(A\) and \(B\) - there will in general be two solutions.

4. Substitute those \(A\) and \(B\) values into the expressions for \(C\), \(D\) and \(E\) to obtain the parabola equations.

An example

As an example, with the four points:

\[ (2,3),\quad (2,1),\quad (-4,1),\quad (1,0). \]

we have the equations (which we note are linear in \(C\), \(D\) and \(E\)):

\[\begin{gathered} (2A+3B)^2+2C+3D+E=0\\ (2A+B)^2+2C+D+E=0\\ (-4A+B)^2-4C+D+E=0\\ (A)^2+C+E=0\end{gathered}\]

for which we find from the first three that

\[\begin{align*} C&=2A^2-2AB\\ D&=-4AB-4B^2\\ E&=-8A^2+4AB+3B^2 \end{align*}\]

or alternatively, that

\[\begin{bmatrix}C\\ D\\ E\end{bmatrix}=\begin{bmatrix*}[r]2&-1&0\\ 0&-2&-4\\ -8&2&3\end{bmatrix*}\begin{bmatrix}A^2\\ 2AB\\ B^2\end{bmatrix}\]

Substituting these expressions into the last equation produces (after a bit of algebraic simplification):

\[ -5A^2+2AB+3B^2=0 \]

which has the two solutions

\[ A = s, B = -\frac{5}{4}s\qquad\mathrm{and}\qquad A=t,\;B=t \]

for arbitrary \(s\) and \(t\). Substituting these back into the expressions for \(C\), \(D\) and \(E\):

\[\begin{align*} C&=\frac{16}{3}s^2&D&=-\frac{40}{9}s^2&E&=-\frac{19}{3}s^2\\[2mm] C&=0&D&=-8t^2&E&=-t^2 \end{align*}\]

These expressions can now be used for the parabola equations. In each equation, the squared parameter \(s^2\) or \(t^2\) can be factored out.

Numerical approach

Here we show how to do the same thing numerically.

Rather than explain in full algebraic detail, we'll simply work through an example, from which the general method will be obvious.

We start with the four points \((x_i,y_i)\):

\[ (2,3),\quad (2,1),\quad (-4,1),\quad (1,0). \]

We first create a \(3\times 3\) matrix \(M\) for which the first two columns are the first three \(x\) and \(y\) values, and the last column is all ones:

\[ M=\begin{bmatrix*}[r]2&3&1\\ 2&1&1\\ -4&1&1\end{bmatrix*} \]

Next is a \(4\times 3\) matrix \(N\) whose rows consist of the values \(x^2,\;2xy,\; y^2\) for each \((x,y)\) coordinates:

\[ N = \begin{bmatrix*}[r]4&12&9\\ 4&4&1\\ 16&-8&1\\ 1&0&0\end{bmatrix*} \]

Next multiply the inverse of \(M\) by the top three rows of \(N\):

\[ P=\begin{bmatrix*}[r]2&3&1\\ 2&1&1\\ -4&1&1\end{bmatrix*}^{-1}\begin{bmatrix*}[r]4&12&9\\ 4&4&1\\16&-8&1\end{bmatrix*}=\begin{bmatrix*}[r]2&-2&0\\ 0&-4&-4\\ -8&4&3\end{bmatrix*} \]

What this matrix contains, of course, are the coefficients of \(A^2\), \(AB\) and \(B^2\) in the initial expressions for \(C\), \(D\) and \(E\).

Then compute

\[ Q=\begin{bmatrix}1&0&1\end{bmatrix}\begin{bmatrix*}[r]2&-2&0\\ 0&-4&-4\\ -8&4&3\end{bmatrix*}+\begin{bmatrix}1&0&0\end{bmatrix}=\begin{bmatrix*}[r]-5&2&3\end{bmatrix*} \]

Here, the first matrix on the right consists of the last (so far unused) coordinate values and a 1; the central matrix is \(P\) (just computed) and the last matrix is the bottom row of \(N\). Thus:

\[ Q=\begin{bmatrix}x_4&y_4&1\end{bmatrix}P+\begin{bmatrix}n_{41}& n_{42}& n_{43}\end{bmatrix} \]

The values of \(Q\) are the coefficients \(a,b,c\) for the final equation \(aA^2+2bAB+cB^2=0\) for \(A\) and \(B\).

Now, it is easy to show that the solutions of the equation

\[ aA^2+2bAB+cB^2=0 \]

are

\[ A=s,\quad B = \frac{-b+\sqrt{b^2-ac}}{c}s\qquad\mathrm{and}\qquad A=t,\quad B=\frac{-b-\sqrt{b^2-ac}}{c}t \]

for arbitrary values of \(s\) and \(t\). Since in our case we have \(a=-5\), \(b=1\) and \(c=3\):

\[ A=s,\quad B = \frac{-1+\sqrt{1+15}}{3}s=s\qquad\mathrm{and}\qquad A=t,\quad B=\frac{-1-\sqrt{1+15}}{3}t=-\frac{5}{3}t \]

We can set \(s=t=1\), or we can set \(s=1\) and \(t=3\) (which eliminates fractions) and substitute those values for \(A\) and \(B\) into the equations for \(C\), \(D\) and \(E\). So putting \(s=1\):

\[\begin{bmatrix}C\\D\\ E\end{bmatrix}=\begin{bmatrix*}[r]2&-2&0\\ 0&-4&-4\\ -8&4&3\end{bmatrix*}\begin{bmatrix*}[r]1\\ 1\\ 1\end{bmatrix*}=\begin{bmatrix*}[r]0\\ -8\\ -1\end{bmatrix*}\]

If we put \(t=1\), in the second solution, then \(A = 1\) and \(B = -5/3\), from which \(C\), \(D\) and \(E\) can be easily obtained.

Finishing off

From above, we have the two solutions, for

\[\begin{align*} A,B,C,D,E&=1,1,0,-8,1\\ & = 1,\dfrac{5}{3},\dfrac{16}{3},-\dfrac{40}{9},-\dfrac{19}{3} \end{align*}\]

This produces the two parabolas:

\[\begin{align*} &(x+y)^2-8y+1=0\\ &(3x+5y)^2+48x-40y-57=0 \end{align*}\]

where in the last equation every coefficient has been multiplied by 9 to clear all the fractions.

Using Python

With NumPy:

[]  import numpy as np
[]  import numpy.linalg as la

[]  xs = np.array([2,2,-4,1])
[]  ys = np.array([3,1,1,0])

[]  M = np.matrix([xs[:-1],ys[:-1],[1,1,1]]).T

[]  N = np.matrix([xs**2,2*xs*ys,ys**2]).T

[]  P = -la.inv(M1)*N[:-1,:]

[]  Q = np.matrix([[xs[-1],ys[-1],1]])*P+N[-1,:]

[]  [a,b,c] = Q.tolist()[0]

[]  b/=2
[]  b0 = (-b+np.sqrt(b*b-a*c))/c
[]  b1 = (-b-np.sqrt(b*b-a*c))/c

[]  AB = np.matrix([[1,b0,b0*b0],[1,b1,b1*b1]]).T

[]  R = P*AB

[]  coeffs = np.vstack([AB[:-1,:],R]).T

[]  display(coeffs)

  matrix([[ 1.00000000e+00,  1.00000000e+00,  2.77555756e-17, -8.00000000e+00, -1.00000000e+00],
          [ 1.00000000e+00, -1.66666667e+00,  5.33333333e+00, -4.44444444e+00, -6.33333333e+00]])

The parabolas can then be plotted using the parameterization given in the previous post.