Quadratic Forms and Eigenvalues

This worksheet explores the relationship between symmetric matrices and
quadratic forms. It shows how the eigenvalues of such a matrix relate to
the geometric character of the graph of the quadratic form. It also discusses
in the context of an example how the eigenvectors of the symmetric matrix
determine a rotation of coordinates making the quadratic form diagonal.

> with(plots): with(plottools):

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

The symmetric matrix corresponding to 2 x^2 + 8 x y + y^2

> A := matrix(2,2,[2,4,4,1]);

A general vector in R^2.

> v := vector(2,[x,y]);

Verifying that A corresponds to this quadratic expression.

> expn1 := expand(dotprod(v,evalm(A &* v)));

Compute the eigenvalues (returned as a list below) and eigenvectors
(returned in the columns of A_evects) of the matrix A.
Recall a nonzero v is eigenvector with eigenvalue lambda if A v = lambda v.

> evalf(Eigenvals(A,A_evects));

The columns here are the corresponding eigenvectors.

> print(A_evects);

Check that the columns of A_evects have been returned as unit vectors.
If not, we'd normalize them now.

> for j from 1 to 2 do
print(sqrt(sum((A_evects[i,j])^2,i=1..2)));
od;

(printf(`Length of column %d is %g\n`,j,sqrt(sum((A_evects[i,j])^2,i=1..2))) would be nicer.)

The eigenvectors of a symmetric matrix can always be chosen to be an orthonormal
basis for R^n, and this reflects itself in the identity below.
This identity also means A_evects is the matrix of a rotation.

> print(evalm(inverse(A_evects) - transpose(A_evects)));

In general, if an n x n matrix A has n independent eigenvectors in the columns of
A, then inverse (A_evects) &* A &* A_evects will be a diagonal matrix.

If we transform coordinates by [x_new y_new] = A_evects &* [x y], then B will
be the symmetric matrix corresponding to expn1 in coordinates x_new,y_new.

> B := evalm(transpose(A_evects) &* A &* A_evects);

> v_new := vector(2,[x_new,y_new]);

The same expression in the new coordinates.

> expn2 := expand(dotprod(v_new,evalm(B &* v_new)));

> plot3d(expn1,x=-5..5,y=-5..5,axes=boxed,shading=zhue,view=-10..10);

> plot3d(expn2,x_new=-5..5,y_new=-5..5,axes=boxed,shading=zhue,view=-10..10);

>