Método de Leverrier-Faddeev

El polinomio caracteristico de un matriz es

\begin{equation*} p(\lambda) = c_{0} + c_{1} \lambda + \dots + c_{n-2} \lambda^{n-2} + c_{n-1} \lambda^{n-1} + c_{n} \lambda^{n} \end{equation*}

Algoritmo

\begin{align*} B_{0} &= 0 \\ c_{n} &= 1 \\ k &= 1,\ldots ,n \\ &\quad B_{k} = AB_{k-1} + c_{n-k+1} I \\ &\quad c_{n-k}= -\frac{1}{k} \ \mathrm{tr}(AB_{k}) \\ \det A &= (-1)^n c_{0} \\ A^{-1} &= - \frac 1 {c_0} B_n \end{align*}

Seudocódigo

function Leverrier_Faddeev(A)
    I = 1 // matriz identidad
    B[0] = 0 // matriz nula
    c[n] = 1 

    for k=1 to n do
        B[k] = A * B[k-1] + c[n-k+1] * I
        c[n-k] = - 1/k * traza(A * B[k])
    end for  

    determinante = (-1)^n * c[0]

    if c[0]!=0
        inversa = -1/c[0] * B[n]
    else
        inversa = 'No tiene inversa'
    end if

    mostrar c, determinante, inversa
end function

Ejemplo 1

Determinar el polinomio caracteristico de $A$

\begin{equation*} A = \begin{bmatrix} 3 & 1 & 5 \\ 3 & 3 & 1 \\ 4 & 6 & 4 \end{bmatrix} \end{equation*}

Iniciando valores

\begin{align*} B_{0} &= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \\ c_3 &= 1 \end{align*}

Bucle 1

\begin{align*} B_{1} &= AB_{0} + c_{3} I = \begin{bmatrix} 3 & 1 & 5 \\ 3 & 3 & 1 \\ 4 & 6 & 4 \end{bmatrix} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \\ c_{2} &= -\frac{1}{1} \ \mathrm{tr}(A B_{1}) = -\frac{1}{1} (10) = -10 \end{align*}

Bucle 2

\begin{align*} B_{2} &= AB_{1} + c_{2} I = \begin{bmatrix} 3 & 1 & 5 \\ 3 & 3 & 1 \\ 4 & 6 & 4 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} + \begin{bmatrix} -10 & 0 & 0 \\ 0 & -10 & 0 \\ 0 & 0 & -10 \end{bmatrix} = \begin{bmatrix} -7 & 1 & 5 \\ 3 & -7 & 1 \\ 4 & 6 & -6 \end{bmatrix} \\ c_{1} &= -\frac{1}{2} \ \mathrm{tr}(A B_{2}) = -\frac{1}{2} (-8) = 4 \end{align*}

Bucle 3

\begin{align*} B_{3} &= AB_{2} + c_{1} I = \begin{bmatrix} 3 & 1 & 5 \\ 3 & 3 & 1 \\ 4 & 6 & 4 \end{bmatrix} \begin{bmatrix} -7 & 1 & 5 \\ 3 & -7 & 1 \\ 4 & 6 & -6 \end{bmatrix} + \begin{bmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{bmatrix} = \begin{bmatrix} 6 & 26 & -14 \\ -8 & -8 & 12 \\ 6 & -14 & 6 \end{bmatrix} \\ c_{0} &= -\frac{1}{3} \ \mathrm{tr}(A B_{3}) = -\frac{1}{3} (120) = -40 \end{align*}

El polinomio caracteristico es

\begin{equation*} p(\lambda) = c_{0} + c_{1} \lambda + c_{2} \lambda^{2} + c_{3} \lambda^{3} = - 40 + 4 \lambda - 10 \lambda^{2} + \lambda^{3} \end{equation*}

El determinante es

\begin{equation*} \det(A) = (-1)^{3} c_0 = (-1)^{3} (-40) = 40 \end{equation*}

Su inversa es

\begin{equation*} A^{-1} = -\frac{1}{c_{0}} B_{3} = -\frac{1}{-40} \begin{bmatrix} 6 & 26 & -14 \\ -8 & -8 & 12 \\ 6 & -14 & 6 \end{bmatrix} = \begin{bmatrix} 0.15 & 0.65 & -0.35 \\ -0.20 & -0.20 & 0.30 \\ 0.15 & -0.35 & 0.15 \end{bmatrix} \end{equation*}

Implementación


In [1]:
import numpy as np

def Leverrier_Faddeev(A):
    m, n = A.shape 
    I = np.eye(n)
    B = np.zeros((n+1,n,n))
    c = np.zeros(n+1)
    
    c[n] = 1.0
    
    for k in range(1,n+1):
        B[k] = np.dot(A,B[k-1]) + c[n-k+1]*I
        c[n-k] = - 1/k * np.trace(np.dot(A,B[k]))
        
    determinante = (-1)**n * c[0]
    
    if c[0]!=0:
        inversa = -1/c[0] * B[n]
    else:
        inversa = np.zeros((n,n)) * np.nan

    print('p(x) =', c)
    print('determinante =', determinante)
    print('inversa =')
    print(inversa)

In [2]:
A = np.array([[3,1,5],
              [3,3,1],
              [4,6,4]], float)
Leverrier_Faddeev(A)


p(x) = [-40.   4. -10.   1.]
determinante = 40.0
inversa =
[[ 0.15  0.65 -0.35]
 [-0.2  -0.2   0.3 ]
 [ 0.15 -0.35  0.15]]

In [3]:
#revisando
print('p(x) =', np.poly(A))
print('determinante =', np.linalg.det(A))
print('inversa =')
print(np.linalg.inv(A))


p(x) = [  1. -10.   4. -40.]
determinante = 40.0
inversa =
[[ 0.15  0.65 -0.35]
 [-0.2  -0.2   0.3 ]
 [ 0.15 -0.35  0.15]]

In [4]:
B = np.array([[1,-4,-1,-4],
              [2,0,5,-4],
              [-1,1,-2,3],
              [-1,4,-1,6]],float)
Leverrier_Faddeev(B)


p(x) = [ 2. -7.  9. -5.  1.]
determinante = 2.0
inversa =
[[  1.   -1.   -4.    2. ]
 [ -0.5  -4.  -11.    2.5]
 [ -0.    3.    8.   -2. ]
 [  0.5   3.    8.   -1.5]]

In [5]:
#revisando
print('p(x) =', np.poly(B))
print('determinante =', np.linalg.det(B))
print('inversa =')
print(np.linalg.inv(B))


p(x) = [ 1. -5.  9. -7.  2.]
determinante = 2.0
inversa =
[[  1.   -1.   -4.    2. ]
 [ -0.5  -4.  -11.    2.5]
 [ -0.    3.    8.   -2. ]
 [  0.5   3.    8.   -1.5]]

In [6]:
C = np.array([[0.01,0,0,0,0],
              [0,0.01,0,0,0],
              [0,0,0.99,0,0],
              [0,0,0,100,0],
              [0,0,0,0,10000]],float)
Leverrier_Faddeev(C)


p(x) = [  1.05845976e+01   1.99010181e+04  -1.01020099e+06   1.01020102e+06
  -1.01010100e+04   1.00000000e+00]
determinante = -10.5845976024
inversa =
[[ -9.35323058e+02  -0.00000000e+00  -0.00000000e+00  -0.00000000e+00
   -0.00000000e+00]
 [ -0.00000000e+00  -9.35323058e+02  -0.00000000e+00  -0.00000000e+00
   -0.00000000e+00]
 [ -0.00000000e+00  -0.00000000e+00  -9.44941026e+00  -0.00000000e+00
   -0.00000000e+00]
 [ -0.00000000e+00  -0.00000000e+00  -0.00000000e+00  -9.52868576e-02
   -0.00000000e+00]
 [ -0.00000000e+00  -0.00000000e+00  -0.00000000e+00  -0.00000000e+00
    4.25900631e-03]]

In [7]:
#revisando
print('p(x) =', np.poly(C))
print('determinante =', np.linalg.det(C))
print('inversa =')
print(np.linalg.inv(C))


p(x) = [  1.00000000e+00  -1.01010100e+04   1.01020102e+06  -1.01020099e+06
   1.99009999e+04  -9.90000000e+01]
determinante = 99.0
inversa =
[[  1.00000000e+02   0.00000000e+00   0.00000000e+00   0.00000000e+00
    0.00000000e+00]
 [  0.00000000e+00   1.00000000e+02   0.00000000e+00   0.00000000e+00
    0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   1.01010101e+00   0.00000000e+00
    0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.00000000e-02
    0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
    1.00000000e-04]]