Homework 10

INSERT YOUR NAME HERE

Directions: Add work to this notebook to solve the problems below.

Be careful to name the objects you create as described in the problem. Some problems may include some tests you can run to test your code.

This assignment will not be collected.

Some problems are derived from problems in these books:

  • 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.

Standard imports

In [ ]:
from scipy import linalg
import numpy as np
import math as m
import matplotlib.pyplot as plt

1. The dominant eigenpair

Consider the matrix defined below:

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

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 [ ]:

Tests for your code:

This vector should be very small:

In [ ]:
A @ eigenvector_0 - eigenvalue_0 * eigenvector_0

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. (This is an example of a tridiagonal matrix.)

In [2]:
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 [ ]: