Demo # 10
Determinants, Roots, and Eigensystems in Mathematica

Some of the material contained in this demo has already been touched on in Demo # 7, so I suggest that you review that notebook before you proceed.

Root Finding

There are two principle tools for finding roots of polynomial equations in Mathematica:
		
Solve -> generates analytical expressions for the n roots
                      of an nth order equation
NSolve -> generates numerical expressions for the n roots

For algebraic equations that are not polynomial in form, e.g. transcendental equations, numerical solutions can be
found with:
		
FindRoot -> Finds a numerical approximation to the closest
                              root of some equation

Let's start with the quadratic equation

Solve[ 
      a λ^2  +  b λ  +  c  ==  
        0,  λ]

This is a very well-known result to you! Next, the cubic equation:

Solve[ 
      a λ^3  +  b λ^2  + c λ  +  d  ==  
        0, λ]

Although this is the general formula, it is clearly a mess, so the Cardan-Tartaglia method described in class is a better way to proceed, at least on paper!

Here is a simple problem:

Solve[ x^3  - x^2  - 4 x  -  6  ==  0,  x]

We get two complex conjugate roots and one real root. Notice that the complex roots always appear as complex conjugate pairs. If we change the equation:

Solve[ 5 x^3  -  2 x^2  + 3 x  - 7  == 0,  x]

Now things are really a mess. NSolve will get us out of this:

NSolve[  5 x^3  -  2 x^2  + 3 x  - 7  == 0,  x]

Again, one real and two complex roots. Let's finish with something you really don't want to try by hand!!

NSolve[ 
      57  x^7  -  105 x^6 + 9 x^5  - 2 x^4  + 98 x^3  + 3 x^2  - 
          5 x  - 89  == 0, x]

Here we have six complex roots, appearing as three complex conjugate pairs and one real root.

Let's finish with an example of solving a transcendental equation. Suppose we wanted to solve:

Tan[x] - x ==0

First, let's plot the function to get an idea of where the roots lie:

Plot[ Tan[x] - x, {x, 0, 10}]

There are a whole set of zero crossings. To get the one at the origin:

Clear[x];
FindRoot[ Tan[x]  -  x  == 0, {x,  .1}]

Actually, I think the solution is at 0, not 0.013 as Mathematica seems to imply.  Let's try improving the the precision:

\(FindRoot[\ Tan[x]\  - \ x\  == 0, \ {x,  .1}, 
	
      WorkingPrecision\  -> \ $MachinePrecision\  + 3]\)

Oops!   The answer changed when I changed the precision.  To make sure that the solution is at 0 and not some small but finite number, let's expand the plot about zero

Plot[ Tan[x] - x, {x, 0, 0.0010}]

Lets expand the function in a Taylor series

Series[Tan[x] - x, {x, 0, 9}]

It looks like the only solution near 0 is 0.

Now let's try to find the root near x=4:

FindRoot[ Tan[x]  -  x  == 0, {x, 4. }]

Clearly, Mathematica's Plot and FindRoot functions make root-finding a very easy process!  But, you must be careful not to be misled.  Plotting the function near the root in question is usually a good idea.

Determinants

Let us start by entering a 4 by 4 matrix and declaring that all matrix output should be in MatrixForm:

$Post := If[MatrixQ[#1],  MatrixForm[#1],  #1] &;
matrix1  =  {{1, 4, 5, 8}, {4, -2, 3,  8}, {4, -7, 9, 
          11}, {6, -3, 5, 7}}

The Determinant of this matrix is obtained simply by:

Det[matrix1]

Notice how easy it is to test the various theorems about determinants that were discussed in your lectures. For example, suppose that we interchange rows 1 and 2:

matrix2  =  {{4, -2, 3,  8}, {1, 4, 5, 8}, {4, -7, 9, 
          11}, {6, -3, 5, 7}}

Det[matrix2]

Thus, the determinant changes sign if we interchange any two rows! What about interchanging two columns? Here is matrix1 with columns 1 and 2 reversed:

matrix3  =  {{4, 1, 5, 8}, {-2, 4, 3, 8}, {-7, 4, 9, 
          11}, {-3, 6, 5, 7}}

Det[matrix3]

Again, the sign is inverted by a single column exchange!
		
We can prove many other things by explicit example. Here is a general formula for the 4 by 4 determinant:

matrix4  =  Array[ a, {4, 4}]

Det[matrix4]

I am sure that you are glad that you did not have to work this one out by hand!
		
Let's now zero out the lower half of the matrix:

matrix5  =  
      matrix4  /.  {a[2, 1] -> 0, a[3, 1] -> 0, a[3, 2] -> 0, a[4, 1] -> 0, 
          a[4, 2] -> 0, a[4, 3] -> 0}

Det[matrix5]

Thus, as stated in class, the determinant of a "triangular" matrix with zeros filling either the upper or lower half is equal to the product of the diagonal elements.

Eigensystems

Next, let's work through some exercises on matrix eigenvalue problems.

This is a simple 3 by 3 matrix:

a  =  {{2, 1, 1}, {1, 3, 1}, {1, 0, 1}}

We get the eigenvalues simply from:

Eigenvalues[a]

These are a mess, because Mathematica is insisting on analytical solutions of the cubic characteristic equation. We can get a simpler result by:

Eigenvalues[N[a]]

or, by:

Eigenvalues[a] // N

The eigenvectors, expressed as rows of a matrix, are also easily obtained:

evecmat  = Eigenvectors[N[a]]

This is a matrix in which the three eigenvectors appear as rows of the matrix.  Thus, the eigenvector corresponding to eigenvalue 3.87939 is (.53699, .822716, .186495).

The individual eigenvectors can be extracted for later use from the rows of the eigenvector matrix:

evec1 = evecmat[[1]]
evec2 = evecmat[[2]]
evec3 = evecmat[[3]]

The first (and subsequent) column(s) of a matix can also be easily extracted as follows:

evecmat[[All, 1]]

Recall that in class, it was shown that if we define a matrix "b" whose columns are the eigenvectors, i.e.

b  =  Transpose[Eigenvectors[N[a]]]

then the transformed matrix: Inverse[b] a b, has the eigenvalues as its diagonal elements. Let's check this:

newa  =  Inverse[b]  .  a  .  b

Indeed, to within numerical precision, this is a matrix with all zeros, except the eigenvalues on the diagonal.

Actually, we can suppress the very small numbers that result from roundoff errors characteristic of numerical work using Mathematic'a Chop function which sets to 0 all numbers smaller than 1×(10^(-10)).

newa  =  Chop[Inverse[b]  .  a  .  b]

Finally, we note that we could have found the eigenvalues of our matrix "a" in the same way that we would on paper:

newnewa  =   a  -  λ IdentityMatrix[3]

chareqn  =  Det[newnewa]

NSolve[ chareqn  ==  0 ,  λ]

Again, we find the same three eigenvalues!

In summary, we see that Mathematica has very powerful capabilities for evaluating determinants, finding roots, and obtaining eigenvalues and eigenvectors of matrices.