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

`> `