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 := MATRIX([[2, 4], [4, 1]])

A general vector in R^2.

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

v := VECTOR([x, y])

Verifying that A corresponds to this quadratic expression.

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

expn1 := 2*x^2+8*x*y+y^2

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));

VECTOR([-2.531128876, 5.531128873])

The columns here are the corresponding eigenvectors.

> print(A_evects);

MATRIX([[.6618025631, -.7496781760], [-.7496781760,...

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



(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)));

MATRIX([[0, 0], [0, 0]])

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);

B := MATRIX([[-2.531128875, .3e-8], [.2e-8, 5.53112...

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

v_new := VECTOR([x_new, y_new])

The same expression in the new coordinates.

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

expn2 := -2.531128875*x_new^2+.5e-8*x_new*y_new+5.5...

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

[Maple Plot]

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

[Maple Plot]