Problem 39, Section 7.1

I'll be working through this one, essentially following the instructions given there. I'll be making comments in between the computations.

First things first: let's load up the relevant Python packages and then define the matrix we're working with.

In [1]:
from sympy import *
init_printing(use_latex='mathjax')
In [2]:
A=Matrix([[S(31)/100,S(58)/100,S(8)/100,S(44)/100],
          [S(58)/100,-S(56)/100,S(44)/100,-S(58)/100],
          [S(8)/100,S(44)/100,S(19)/100,-S(8)/100],
          [S(44)/100,-S(58)/100,-S(8)/100,S(31)/100]])
A
Out[2]:
$\displaystyle \left[\begin{matrix}\frac{31}{100} & \frac{29}{50} & \frac{2}{25} & \frac{11}{25}\\\frac{29}{50} & - \frac{14}{25} & \frac{11}{25} & - \frac{29}{50}\\\frac{2}{25} & \frac{11}{25} & \frac{19}{100} & - \frac{2}{25}\\\frac{11}{25} & - \frac{29}{50} & - \frac{2}{25} & \frac{31}{100}\end{matrix}\right]$

That's the matrix problem 39 works with. Note what I did to ensure that the numbers truly are fractions rather than floating-point numbers subject to rounding errors: instead of, say, 0.31 or 31/100, I entered S(31)/100. That tells Sympy that I want the actual rational number, with precisely that numerator and denominator, to be used as such in computations (rather than rounded, etc.).

In the past I've used expression like Rational(31,100); that's equivalent to what I did here, but not to 31/100. Check this:

In [3]:
31/100 == S(31)/100
Out[3]:
False

That means they're not the same. On the other hand:

In [4]:
Rational(31,100) == S(31)/100
Out[4]:
True

OK, that's enough of that; let's move on to matrix computations. The problem tells us to

  • find the eigenvalues of $A$;
  • for each eigenvalue $\lambda$ find an orthonormal basis for the corresponding eigenspace (which is the nullspace of $A-\lambda I_4$);
  • then, as in examples 2 and 3 in the same section, write $A=PDP^{-1}$ where $P$ is the matrix whose columns are the orthonormal eigenvectors we constructed and $D$ is a diagonal matrix.

Let's get to that. First, the eigenvalues:

In [5]:
A.eigenvals()
Out[5]:
$\displaystyle \left\{ - \frac{5}{4} : 1, \ 0 : 1, \ \frac{3}{4} : 2\right\}$

OK, eigenvalues $-\frac 54$ and $0$ each with multiplicity $1$ and $\frac 34$ with multiplicity $2$ (so that guy's eigenspace will have a two-vector orthonormal basis). Let's have a look at eigenvectors too.

In [6]:
A.eigenvects()
Out[6]:
$\displaystyle \left[ \left( - \frac{5}{4}, \ 1, \ \left[ \left[\begin{matrix}-1\\2\\- \frac{1}{2}\\1\end{matrix}\right]\right]\right), \ \left( 0, \ 1, \ \left[ \left[\begin{matrix}-1\\- \frac{1}{2}\\2\\1\end{matrix}\right]\right]\right), \ \left( \frac{3}{4}, \ 2, \ \left[ \left[\begin{matrix}\frac{3}{2}\\1\\1\\0\end{matrix}\right], \ \left[\begin{matrix}1\\0\\0\\1\end{matrix}\right]\right]\right)\right]$

As before, they (the eigenvectors) are nested inside a list of lists: for each eigenvalue we have a list of the form (eigenvalue,multiplicity,[list of eigenvectors]). Now, I want to take that [list of eigenvectors] and orthonormalize it. SymPy makes that dead-simple:

In [7]:
evects_ll=[Matrix.orthogonalize(*x[2],normalize=True) for x in A.eigenvects()]
evects_ll
Out[7]:
$\displaystyle \left[ \left[ \left[\begin{matrix}- \frac{2}{5}\\\frac{4}{5}\\- \frac{1}{5}\\\frac{2}{5}\end{matrix}\right]\right], \ \left[ \left[\begin{matrix}- \frac{2}{5}\\- \frac{1}{5}\\\frac{4}{5}\\\frac{2}{5}\end{matrix}\right]\right], \ \left[ \left[\begin{matrix}\frac{3 \sqrt{17}}{17}\\\frac{2 \sqrt{17}}{17}\\\frac{2 \sqrt{17}}{17}\\0\end{matrix}\right], \ \left[\begin{matrix}\frac{8 \sqrt{17}}{85}\\- \frac{6 \sqrt{17}}{85}\\- \frac{6 \sqrt{17}}{85}\\\frac{\sqrt{17}}{5}\end{matrix}\right]\right]\right]$

This is a list of lists (which is why I called it evects_ll)! But each of the innermost lists is orthonormal, at least. The innermost elements are our vectors. First, I want to "flatten" the list (I've done this before somewhere), so as to get just a plain list of vectors that I'll then make a matrix $P$ out of (as columns):

In [8]:
cols=cols=[y for x in evects_ll for y in x]
cols
Out[8]:
$\displaystyle \left[ \left[\begin{matrix}- \frac{2}{5}\\\frac{4}{5}\\- \frac{1}{5}\\\frac{2}{5}\end{matrix}\right], \ \left[\begin{matrix}- \frac{2}{5}\\- \frac{1}{5}\\\frac{4}{5}\\\frac{2}{5}\end{matrix}\right], \ \left[\begin{matrix}\frac{3 \sqrt{17}}{17}\\\frac{2 \sqrt{17}}{17}\\\frac{2 \sqrt{17}}{17}\\0\end{matrix}\right], \ \left[\begin{matrix}\frac{8 \sqrt{17}}{85}\\- \frac{6 \sqrt{17}}{85}\\- \frac{6 \sqrt{17}}{85}\\\frac{\sqrt{17}}{5}\end{matrix}\right]\right]$

Finally, I can make my matrix $P$ with these guys as columns by using numpy's hstack function:

In [9]:
import numpy as np
P=Matrix(np.hstack(cols))
P
Out[9]:
$\displaystyle \left[\begin{matrix}- \frac{2}{5} & - \frac{2}{5} & \frac{3 \sqrt{17}}{17} & \frac{8 \sqrt{17}}{85}\\\frac{4}{5} & - \frac{1}{5} & \frac{2 \sqrt{17}}{17} & - \frac{6 \sqrt{17}}{85}\\- \frac{1}{5} & \frac{4}{5} & \frac{2 \sqrt{17}}{17} & - \frac{6 \sqrt{17}}{85}\\\frac{2}{5} & \frac{2}{5} & 0 & \frac{\sqrt{17}}{5}\end{matrix}\right]$

And finally, let's recover our diagonal matrix $D$ that will satisfy $A=PDP^{-1}$ as $D=P^{-1}AP$:

In [10]:
D=P.inv()*A*P
D
Out[10]:
$\displaystyle \left[\begin{matrix}- \frac{5}{4} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & \frac{3}{4} & 0\\0 & 0 & 0 & \frac{3}{4}\end{matrix}\right]$

Indeed diagonal with the eigenvalues of $A$ (namely $-\frac 54$, $0$ and $\frac 34$ twice) as its diagonal entries, as expected.