Homework 8

Student's Name: PLEASE INSERT YOUR NAME HERE

Directions: Add work to this notebook to solve the problems below. Some other cells have been included which are important for testing your work and should give you some feeback.

Check your work. I give some tests that should be passed if you get it correct.

Your notebook should run without printing errors and without user input. If this is not the case, points will be deducted from your grade.

You may not use eigenvalue/eigenvector methods built into python (except possibly to check your work). You must use the methods discussed in class.

Problem Sources:

  • LL: Programming for Computations - Python by Svein Linge and Hans Petter Langtangen, 2nd edition.
  • L: A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
  • TAK: Applied Scientific Computing With Python by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh.
In [ ]:
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv
from scipy import linalg

1. The dominant eigenpair

Consider the matrix defined below:

In [ ]:
A = np.array([[max(5-row,col+1) for col in range(5)] for row in range(5)])
A

This matrix has five distinct real eigenvalues.

Use the power method to find a dominant eigenvalue and eigenvector pair. Be sure that the eigenvalue is accurate to 5 decimal places. Store the eigenvalue and eigenvector in the variables eigenvalue_0 and eigenvector_0, respectively. Check that your answer is likely correct by verifying that

A @ eigenvector_0 - eigenvalue_0 * eigenvector_0

is small.

(This exercise is based on TAK §4.6 #1.)

In [ ]:

2. The eigenvalue of smallest absolute value

Find the eigenvalue of smallest absolute value. Make sure the result is accurate to 5 decimal places and store the result in eigenvalue_1. Store the corresponding eigenvector in eigenvector_1. Check that your answer corresponds to an eigenvalue.

You may use code from the notebook Eigenvalues and Eigenvectors used in class or code from § 4.6 of TAK.

(This is problem # 2 from TAK § 4.6.)

In [ ]:

3. Some other eigenpairs

Find the eigenvalue which is closest to $3$. Store the result in eigenvalue_2, and be sure the result is accurate to 5 decimal places. Store the corresponding eigenvector in eigenvector_2. Similarly store the eigenvalue nearest to $-1$ in the variables eigenvalue_3 and store the corresponding eigenvector in eigenvector_3.

Use inverse iteration with an origin shift as discussed in TAK § 4.6 and in the Eigenvalues and Eigenvectors notebook from class.

In [ ]:

4. The last eigenvector and eigenvalue pair

The trace of a matrix is the sum of the diagonal entries. The trace of a matrix also is the sum of the eigenvalues (counting algebraic multiplicity).

Compute the trace of the matrix $A$ and store the result in the variable trace.

Use trace to approximate the fifth eigenvalue of $A$ (using the four other eigenvalues you have found). Since you know it is near the approximate eigenvalue, you can use inverse iteration with an origin shift to get a better approximate. Compute the fifth eigenvalue accurate to 5 decimal places and store the result in eigenvalue_4. Store the corresponding eigenvector in eigenvector_4. Check your work.

In [ ]:

5. Diagonalization

Use the work from the problems above to find an invertible matrix $P$ and a diagonal matrix $D$ such that $$A = P D P^{-1}.$$ Store $P$ in the variable P and $D$ in the variable D. Compute the inverse of $P$ and store it in the variable P_inv. Check that the equation is satisfied (up to error of order $10^{-4}$.

To compute the inverse you can use np.linalg.inv in Numpy's linear algebra module.

In [ ]:

6. The shape of a dominant eigenvector

The function below returns the matrix $B_n$, which is the $n \times n$ matrix, which has ones on the diagonal and in the entries adjacent to the diagonal but zeros everywhere else.

In [ ]:
def B(n):
    B_n = np.zeros((n,n))
    for j in range(n):
        B_n[j,j] = 1
    for j in range(n-1):
        B_n[j,j+1] = 1
        B_n[j+1,j] = 1
    return B_n

Write a function dominant_eigenvector(n, epsilon) which returns the an approximate dominant eigenvector of $B_n$. This eigenvector must be computed using the power method. Compute until the assoicated approximate eigenvalue changes by less than epsilon (which must be positive).

Note that there is not a single dominant eigenvector: scalar multiples of a dominant eigenvector are also dominant eigenvectors.

Guess a formula for a dominant eigenvector and verify it agrees with what your function computes. Write a function eigenvector_guess(n) which returns your guess for the dominant eigenvector. Verify that it is correct by comparing with what your dominant_eigenvector(n, epsilon) and by working with the matrix $B_n$.

Hints: First concentrate on the case that $n$ is odd. Normalize the dominant eigenvector so that the middle entry is one. Try plotting the eigenvector using plt.plot(eigenvector). What does it look like? (It relates to waves...)

In [ ]: