Two dimensional quadratic approximation

Failure of a linear approximation to a surface will occur whenever the triangulation of the coordinate grid is insufficiently fine, which is normally compensated by reducing the grid size. Near critical points, however, a smaller grid only makes matters worse. The remedy is to use a quadratic approximation, unless the critical point is of still higher order. The most general second degree polynomial in two variables has the form

Ax²+2Bxy+Cy²+2Dx+2Ey+F,

which implies a basis composed of six monomials and therefore a $7\times7$ basic determinant

$$\left\vert \begin{array}{ccccccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 ... ...z_2 & z_3 & z_4 & z_5 & z_6 & z \end{array} \right\vert = 0.$$

It would not seem to be advisable to try to express such large Vandermonde determinants in symbolic form, so any particular instance should be inverted numerically to obtain the row of coefficients Y^TM^-1. An inhomogeneous second order polynomial can be written as a quadratic form,

z = U^TQU,

where

$$U^T = \left[ \begin{array}{ccc} 1 & x & y \end{array} \right],$$

and

$$Q=\left[ \begin{array}{ccc} F&E&D\\ E&C&B\\ D&B&A \end{array} \right].$$

The contour lines for z can be represented by y as a function of x if the terms in y are gathered together and expressed as a quadratic form. Write the previous equation in the form

Cy²+2(Bx+E)y+(Ax²+2Dx+F-z)=0,

using the coefficients gotten from inverting the Vandermonde matrix, and solve it for y using the quadratic formula,

$$y=\frac{-(Bx+E)+\sqrt{(Bx+E)^2-C(Ax^2+2Dx+F-z)}}{C}.$$

The other half of the curve could be traced out with the negative square root.

Nevertheless there is a better way to procede. Traditionally the cross term is eliminated by rotating x and y,

$$x' \cos\ \theta + y' \sin\ \theta,$$ x = (42)

$$-x' \sin\ \theta + y' \cos\ \theta.$$ y = (43)

Substituting these values into 2Bxy then leaves, after some trigonometric manipulation, a coefficient for x'y',

$$(A-C)\sin\ 2\theta+2B\cos\ 2\theta,$$

which will vanish if $\theta$ is chosen to be

$$2\theta = {\rm tan}^{-1}\frac{2B}{C-A}.$$

The expression becomes

$$(A\cos^2\ \theta-2B\sin\ \theta\cos\ \theta+C\cos^2\ \theta){x'}^2+ 2(D\cos\ \theta-E\sin\ \theta)x'+$$

$$(A\cos^2\ \theta+2B\sin\ \theta\cos\ \theta+C\cos^2\ \theta){y'}^2+ 2(D\sin\ \theta+E\cos\ \theta)y'+F.$$

With further trigonometric manipulation,

$$\sin\ 2\theta=\frac{2B}{\sqrt{A^2+4B^2+C^2-2AC}},$$ (44)

$$\cos\ 2\theta=\frac{C-A}{\sqrt{A^2+4B^2+C^2-2AC}},$$ (45)

from which the half-angle formulas lead to

$$\sin\ \theta=\sqrt{\frac{1-\cos\ 2\theta}{2}}$$ (46)

$$\cos\ \theta=\sqrt{\frac{1+\cos\ 2\theta}{2}}.$$ (47)

The second step is to complete the squares, using the formula

$$\alpha x^2+\beta x=\alpha(x+\frac{\beta}{2\alpha})^2- (\frac{\beta}{2\alpha})^2,$$

to get

A'x'² + C'y'² + F' = 0.

Then, we might have

$$(ZU)^T=\left[\begin{array}{ccc}\xi&\eta&\zeta\end{array}\right]$$

with

$$Z^T Q Z = \Lambda.$$

Q has not been diagonalized because diagonalization is a transformation of the type Z^-1QZ rather than Z^TQZ; but the translation which completes the square is not unitary, as would be necessary for the two transformations to coincide.

The two remaining diagonal elements of $\Lambda$ may have the same sign, or opposite signs, according to whether the critical point is extremal or a saddle point. Accordingly coefficients of their eigenvectors could be chosen according to one of the following two forms,

$$(ZU)^T=\left[\begin{array}{ccc}0&\cos\ t&\sin\ t\end{array}\right],$$

$$(ZU)^T=\left[\begin{array}{ccc}0&\cosh\ t&\sinh\ t\end{array}\right],$$

depending upon whether the eigenvalues agree or disagree in sign. Returning to the original coordinate system gives contour lines uniformly parameterized by t.

It still remains to ascertain how much of the conic section lies within the simplex that will be contoured; one way to avoid the problem would be to graph it through some parameter range, with a visibility function to determine whether the segments to be graphed lie within the chosen simplex or not.

Another important point concerns the location of the six data points needed for quadratic interpolation; one way to obtain them is to choose two consecutive squares from the coordinate grid, and another would be to calculate supplemental points within the square iteslf. Yet another solution would be to work with a triangular or hexagonal grid. Such alternatives tend to be incompatible with normal data acquisition techniques, but might be feasible whenever data could be generated on demand at specified coordinates.