Diagonalización

En este documento ilustraremos cómo usar Python para resolver los problemas tipo propuestos por L. Merino y E. Santos en página de resolución de ejercicios tipo correspondientes al bloque “Diagonalización”.

Ejemplo

Dada la matriz

\[ A=\left( \begin{array}{ccc} 5 & 1 & -3 \\ 1 & 5 & -3 \\ 2 & 2 & -2 \end{array}\right). \]

Determinar si es o no diagonalizable y, en caso de que lo sea, determinar su forma diagonal y una matriz de paso.

from sympy import Matrix, eye

A=Matrix([(5,1,-3),(1,5,-3),(2,2,-2)])
A

\(\displaystyle \left[\begin{matrix}5 & 1 & -3\\1 & 5 & -3\\2 & 2 & -2\end{matrix}\right]\)

Podemos calcular los valores propios con eigenvals.

A.eigenvals()

{4: 2, 0: 1}

O ver toda la información (incluyendo multiplicidades y vectores propios) con eigenvects.

A.eigenvects()

[(0,
  1,
  [Matrix([
   [1/2],
   [1/2],
   [  1]])]),
 (4,
  2,
  [Matrix([
   [-1],
   [ 1],
   [ 0]]),
   Matrix([
   [3],
   [0],
   [1]])])]

De la salida que obtenemos, observamos que tenemos dos valores propios: 0 con multiplicidad algebraica y geométrica 1, y 4 con multiplicidad algebraica y geométrica 2. También podemos construir la matriz de paso con los datos que hemos obtenido, pero es más fácil usar diagonalize.

A.diagonalize()

(Matrix([
 [1, -1, 3],
 [1,  1, 0],
 [2,  0, 1]]),
 Matrix([
 [0, 0, 0],
 [0, 4, 0],
 [0, 0, 4]]))

P, D = A.diagonalize()

P

\(\displaystyle \left[\begin{matrix}1 & -1 & 3\\1 & 1 & 0\\2 & 0 & 1\end{matrix}\right]\)

D

\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 4 & 0\\0 & 0 & 4\end{matrix}\right]\)

Comprobemos que efectivamente \(P^{-1}AP\) es diagonal.

P.inv()*A*P==D

True

Hagamos este proceso paso a paso usando calculando el polinomio característico y los subespacios propios.

Empezamos calculando el polinomio característico de A y sus raíces (valores propios).

from sympy.abc import x

A.charpoly()

\(\displaystyle \operatorname{PurePoly}{\left( \lambda^{3} - 8 \lambda^{2} + 16 \lambda, \lambda, domain=\mathbb{Z} \right)}\)

A.charpoly(28)

\(\displaystyle \operatorname{PurePoly}{\left( 28^{3} - 8 28^{2} + 16 28, 28, domain=\mathbb{Z} \right)}\)

p=A.charpoly(x)
p

\(\displaystyle \operatorname{PurePoly}{\left( x^{4} - 8 x^{3} + 13 x^{2} - 6 x, x, domain=\mathbb{Z} \right)}\)

p.all_roots()

[0, 1, 1, 6]

Por lo que 0 tiene multiplicidad 1, y 4 multiplicidad 2.

Calculemos ahora los subespacios propios asociados a cada uno de esos valores propios.

V0=A.nullspace()
V0

[Matrix([
 [-1],
 [ 0],
 [ 2],
 [ 1]])]

V4=(A-4*eye(3)).nullspace()
V4

---------------------------------------------------------------------------
ShapeError                                Traceback (most recent call last)
Cell In[40], line 1
----> 1 V4=(A-4*eye(3)).nullspace()
      2 V4

File /usr/local/lib/python3.9/site-packages/sympy/core/decorators.py:106, in call_highest_priority.<locals>.priority_decorator.<locals>.binary_op_wrapper(self, other)
    104         if f is not None:
    105             return f(self)
--> 106 return func(self, other)

File /usr/local/lib/python3.9/site-packages/sympy/matrices/common.py:2937, in MatrixArithmetic.__sub__(self, a)
   2935 @call_highest_priority('__rsub__')
   2936 def __sub__(self, a):
-> 2937     return self + (-a)

File /usr/local/lib/python3.9/site-packages/sympy/core/decorators.py:106, in call_highest_priority.<locals>.priority_decorator.<locals>.binary_op_wrapper(self, other)
    104         if f is not None:
    105             return f(self)
--> 106 return func(self, other)

File /usr/local/lib/python3.9/site-packages/sympy/matrices/common.py:2642, in MatrixArithmetic.__add__(self, other)
   2640 if hasattr(other, 'shape'):
   2641     if self.shape != other.shape:
-> 2642         raise ShapeError("Matrix size mismatch: %s + %s" % (
   2643             self.shape, other.shape))
   2645 # honest SymPy matrices defer to their class's routine
   2646 if getattr(other, 'is_Matrix', False):
   2647     # call the highest-priority class's _eval_add

ShapeError: Matrix size mismatch: (4, 4) + (3, 3)

Podemos usar hstack para combinar los vectores que hemos encontrado.

P=Matrix.hstack(V0[0],V4[0],V4[1])
P

\(\displaystyle \left[\begin{matrix}\frac{1}{2} & -1 & 3\\\frac{1}{2} & 1 & 0\\1 & 0 & 1\end{matrix}\right]\)

O de una forma más compacta.

Matrix.hstack(*(V0+V4))

\(\displaystyle \left[\begin{matrix}\frac{1}{2} & -1 & 3\\\frac{1}{2} & 1 & 0\\1 & 0 & 1\end{matrix}\right]\)

P.inv()*A*P

\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 4 & 0\\0 & 0 & 4\end{matrix}\right]\)

Ejercicio

Estudiar si la siguiente matriz es o no diagonalizable y, en tal caso, determinar su forma diagonal y una matriz de paso. \[ A=\left( \begin{array}{ccc} 2 & 1 & -2 \\ 0 & 2 & -1 \\ 0 & 0 & 1 \end{array}\right) \]

A=Matrix([(2,1,-2),(0,2,-1),(0,0,1)])
A

\(\displaystyle \left[\begin{matrix}2 & 1 & -2\\0 & 2 & -1\\0 & 0 & 1\end{matrix}\right]\)

Al ser triangular, ya sabemos cuáles van a ser los valores propios (2 y 1).

p=A.charpoly(x)
p.all_roots()

[1, 2, 2]

El valor propio 2 tiene multiplicidad algebraica 2, pero sin embargo, su subespacio propio asociado tiene dimensión uno, por lo que \(A\) no es diagonalizable.

V2=(A-2*eye(3)).nullspace()
V2

[Matrix([
 [1],
 [0],
 [0]])]

try:
    A.diagonalize()
except:
    print("No es diagonalizable")

No es diagonalizable

Ejercicio

Razonar que la siguiente matriz \(A\) es diagonalizable y determinar su forma diagonal y una matriz de paso. Determinar \(A^k\) en función de \(k\). \[ A=\left( \begin{array}{cccc} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 \\ 1 & 0 & -2 & 5 \end{array}\right) \]

A=Matrix([(1,0,0,1),(0,1,0,0),(0,0,1,-2),(1,0,-2,5)])
A

\(\displaystyle \left[\begin{matrix}1 & 0 & 0 & 1\\0 & 1 & 0 & 0\\0 & 0 & 1 & -2\\1 & 0 & -2 & 5\end{matrix}\right]\)

A.eigenvects()

[(0,
  1,
  [Matrix([
   [-1],
   [ 0],
   [ 2],
   [ 1]])]),
 (1,
  2,
  [Matrix([
   [0],
   [1],
   [0],
   [0]]),
   Matrix([
   [2],
   [0],
   [1],
   [0]])]),
 (6,
  1,
  [Matrix([
   [ 1/5],
   [   0],
   [-2/5],
   [   1]])])]

Vemos que las multiplicidades algebraicas y geométricas coinciden. Por tanto, A es diagonalizable.

P,D = A.diagonalize()
P

\(\displaystyle \left[\begin{matrix}-1 & 0 & 2 & 1\\0 & 1 & 0 & 0\\2 & 0 & 1 & -2\\1 & 0 & 0 & 5\end{matrix}\right]\)

P.inv()*A*P

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 6\end{matrix}\right]\)

Calculemos P de nuevo “a mano”.

A.charpoly().all_roots()

[0, 1, 1, 6]

V0=A.nullspace()
V0

[Matrix([
 [-1],
 [ 0],
 [ 2],
 [ 1]])]

V1=(A-eye(4)).nullspace()
V1

[Matrix([
 [0],
 [1],
 [0],
 [0]]),
 Matrix([
 [2],
 [0],
 [1],
 [0]])]

V6=(A-6*eye(4)).nullspace()
V6

[Matrix([
 [ 1/5],
 [   0],
 [-2/5],
 [   1]])]

P=Matrix.hstack(*(V0+V1+V6))
P

\(\displaystyle \left[\begin{matrix}-1 & 0 & 2 & \frac{1}{5}\\0 & 1 & 0 & 0\\2 & 0 & 1 & - \frac{2}{5}\\1 & 0 & 0 & 1\end{matrix}\right]\)

P.inv()*A*P

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 6\end{matrix}\right]\)

Calculemos ahora las potencias de \(A\). Para todo \(k\) entero positivo, \(A^k=(P^{-1}DP)^k=P^{-1}D^kP\).

from sympy import var, powsimp

k=var("k", integer=True, nonzero=True)

D=Matrix.diag([0,1,1,6])
D

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 6\end{matrix}\right]\)

D**k

\(\displaystyle \left[\begin{matrix}0^{k} & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 6^{k}\end{matrix}\right]\)

Ak=P.inv()*Matrix.diag([0,1,1,6**k])*P
Ak

\(\displaystyle \left[\begin{matrix}\frac{6^{k}}{6} + \frac{2}{3} & 0 & \frac{1}{3} & \frac{6^{k}}{6} - \frac{2}{15}\\0 & 1 & 0 & 0\\\frac{2}{5} & 0 & \frac{1}{5} & - \frac{2}{25}\\\frac{5 \cdot 6^{k}}{6} - \frac{2}{3} & 0 & - \frac{1}{3} & \frac{5 \cdot 6^{k}}{6} + \frac{2}{15}\end{matrix}\right]\)

powsimp(Ak)

\(\displaystyle \left[\begin{matrix}\frac{6^{k}}{6} + \frac{2}{3} & 0 & \frac{1}{3} & \frac{6^{k}}{6} - \frac{2}{15}\\0 & 1 & 0 & 0\\\frac{2}{5} & 0 & \frac{1}{5} & - \frac{2}{25}\\\frac{5 \cdot 6^{k}}{6} - \frac{2}{3} & 0 & - \frac{1}{3} & \frac{5 \cdot 6^{k}}{6} + \frac{2}{15}\end{matrix}\right]\)

Ejercicio

Estudiar para qué valores de los parámetros \(a\) y \(b\) es diagonalizable la matriz \[ A=\left( \begin{array}{ccc} 1 & a & 0 \\ 0 & 0 & b \\ 0 & b & 0 \end{array}\right). \]

from sympy.abc import a,b

A=Matrix([[1,a,0],[0,0,b],[0,b,0]])
A

\(\displaystyle \left[\begin{matrix}1 & a & 0\\0 & 0 & b\\0 & b & 0\end{matrix}\right]\)

p=A.charpoly(x)
p

\(\displaystyle \operatorname{PurePoly}{\left( x^{3} - x^{2} - b^{2} x + b^{2}, x, domain=\mathbb{Z}\left[b\right] \right)}\)

Factorizamos el polinomio característico.

p.factor_list()

(1,
 [(PurePoly(x - 1, x, domain='ZZ[b]'), 1),
  (PurePoly(x - b, x, domain='ZZ[b]'), 1),
  (PurePoly(x + b, x, domain='ZZ[b]'), 1)])

O bien como expresión.

p.as_expr().factor()

\(\displaystyle \left(- b + x\right) \left(b + x\right) \left(x - 1\right)\)

O bien con solve.

from sympy import solve

solve(p.as_expr(),x)

[1, -b, b]

Los valores propios son \(\{1,b,-b\}\). Así si \(b\not\in\{0,1,-1\}\), los tres valores propios son diferentes y la matriz es diagonalizable.

Estudiemos ahora qué ocurre para \(b=0\), en cuyo caso los valores propios son \(1\) y \(0\) con multiplicidad dos. Nos interesa por tanto ver si el subespacio propio asociado a \(0\) tiene dimensión dos.

Ab0=A.subs({b:0})
Ab0

\(\displaystyle \left[\begin{matrix}1 & a & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\)

Que tiene rango 1, por lo que que el valor propio 0 tiene multiplicidad geométrica 2.

Ab0.nullspace()

[Matrix([
 [-a],
 [ 1],
 [ 0]]),
 Matrix([
 [0],
 [0],
 [1]])]

Por lo que para \(b=0\), la multiplicidad algebraica y geométrica de cada valor propio coincide, siendo así Ab0 diagonalizable.

P0, D= Ab0.diagonalize()
P0

\(\displaystyle \left[\begin{matrix}- a & 0 & 1\\1 & 0 & 0\\0 & 1 & 0\end{matrix}\right]\)

P0.inv()*Ab0*P0

\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 1\end{matrix}\right]\)

Veamos ahora el caso \(b=1\). Los valores propios son \(1\), con multiplicidad dos, y \(-1\).

Ab1=A.subs({b:1})
Ab1

\(\displaystyle \left[\begin{matrix}1 & a & 0\\0 & 0 & 1\\0 & 1 & 0\end{matrix}\right]\)

Ab1-eye(3)

\(\displaystyle \left[\begin{matrix}0 & a & 0\\0 & -1 & 1\\0 & 1 & -1\end{matrix}\right]\)

En este caso el rango de \(A-I\) depende de \(a\). - Si \(a\neq 0\), el rango es dos, y la dimension del subespacio propio asociado a 1 sería 1, por lo que la matriz no sería diagonalizable. - Si \(a=0\), entonces la matriz es diagonalizable, pues el rango de \(A-I\) es 1, y por tanto su núcleo tiene dimensión 2.

Hay que tener cuidado con usar rank o rref directamente, pues sympy considera a como un símbolo.

(Ab1-eye(3)).rank()

2

(Ab1-eye(3)).rref()

(Matrix([
 [0, 1, 0],
 [0, 0, 1],
 [0, 0, 0]]),
 (1, 2))

Podemos considerar el menor \(2\times 2\) superior derecho, que nos proporciona la distinción de los valores de \(a\) que hemos hecho arriba.

M=(Ab1-eye(3))[:2,1:]
M

\(\displaystyle \left[\begin{matrix}a & 0\\-1 & 1\end{matrix}\right]\)

M.det()

\(\displaystyle a\)

Para \(b=-1\), tenemos como valores propios a 1 (multiplicidad 2) y -1.

Abm1=A.subs({b:1})
Abm1

\(\displaystyle \left[\begin{matrix}1 & a & 0\\0 & 0 & 1\\0 & 1 & 0\end{matrix}\right]\)

(Abm1-eye(3))

\(\displaystyle \left[\begin{matrix}0 & a & 0\\0 & -1 & 1\\0 & 1 & -1\end{matrix}\right]\)

Repitiéndose la distinción que hicimos para \(b=1\).