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.
from sympy import *
init_printing(use_latex='mathjax')
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
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:
31/100 == S(31)/100
That means they're not the same. On the other hand:
Rational(31,100) == S(31)/100
OK, that's enough of that; let's move on to matrix computations. The problem tells us to
Let's get to that. First, the eigenvalues:
A.eigenvals()
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.
A.eigenvects()
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:
evects_ll=[Matrix.orthogonalize(*x[2],normalize=True) for x in A.eigenvects()]
evects_ll
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):
cols=cols=[y for x in evects_ll for y in x]
cols
Finally, I can make my matrix $P$ with these guys as columns by using numpy's hstack
function:
import numpy as np
P=Matrix(np.hstack(cols))
P
And finally, let's recover our diagonal matrix $D$ that will satisfy $A=PDP^{-1}$ as $D=P^{-1}AP$:
D=P.inv()*A*P
D
Indeed diagonal with the eigenvalues of $A$ (namely $-\frac 54$, $0$ and $\frac 34$ twice) as its diagonal entries, as expected.