El polinomio caracteristico de un matriz es
\begin{equation*} p(\lambda) = c_{0} \lambda^{n} + c_{1} \lambda^{n-1} + c_{2} \lambda^{n-2} + \dots + c_{n-1} \lambda + c_{n} \end{equation*}function Leverrier_Faddeev(A)
    I = 1 // matriz identidad
    c[0] = 1
    B[0] = I 
    for k=1 to n-1 do
        c[k] = - 1/k * traza(A * B[k-1])
        B[k] = A * B[k-1] + c[k] * I
    end for 
    c[n] = - 1/n * traza(A * B[n-1])
    determinante = (-1)^n * c[n]
    if c[0]!=0
        inversa = -1/c[n] * B[n-1]
    else
        inversa = 'No tiene inversa'
    end if
    mostrar c, determinante, inversa
end function
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*} c_{0} &= 1 \\ B_{0} &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{align*}Bucle 1
\begin{align*} c_{1} &= -\frac{1}{1} \ \mathrm{tr}(A B_{0}) = -\frac{1}{1} (10) = -10 \\ B_{1} &= AB_{0} + c_{1} 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} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} -7 & 1 & 5 \\ 3 & -7 & 1 \\ 4 & 6 & -6 \end{bmatrix} \end{align*}Bucle 2
\begin{align*} c_{2} &= -\frac{1}{2} \ \mathrm{tr}(A B_{1}) = -\frac{1}{2} (-8) = 4 \\ B_{2} &= AB_{1} + c_{2} 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} \end{align*}Calculando $c_{3}$
\begin{align*} c_{3} &= -\frac{1}{3} \ \mathrm{tr}(A B_{2}) = -\frac{1}{3} (120) = -40 \end{align*}El polinomio caracteristico es
\begin{equation*} p(\lambda) = c_{0} \lambda^{3} + c_{1} \lambda^{2} + c_{2} \lambda + c_{3} = \lambda^{3} - 10 \lambda^{2} + 4 \lambda - 40 \end{equation*}El determinante es
\begin{equation*} \det(A) = (-1)^{3} c_{n} = (-1)^{3} (-40) = 40 \end{equation*}Su inversa es
\begin{equation*} A^{-1} = -\frac{1}{c_{n}} B_{2} = -\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*}
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[0] = 1.0
    B[0] = I
    
    for k in range(1,n):
        c[k] = - 1/k * np.trace(np.dot(A,B[k-1]))
        B[k] = np.dot(A,B[k-1]) + c[k]*I
        
    c[n] = - 1/n * np.trace(np.dot(A,B[n-1]))
    determinante = (-1)**n * c[n]
    
    if c[0]!=0:
        inversa = -1/c[n] * B[n-1]
    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)
    
    
In [3]:
    
#revisando
print('p(x) =', np.poly(A))
print('determinante =', np.linalg.det(A))
print('inversa =')
print(np.linalg.inv(A))
    
    
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)
    
    
In [5]:
    
#revisando
print('p(x) =', np.poly(B))
print('determinante =', np.linalg.det(B))
print('inversa =')
print(np.linalg.inv(B))
    
    
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)
    
    
In [7]:
    
#revisando
print('p(x) =', np.poly(C))
print('determinante =', np.linalg.det(C))
print('inversa =')
print(np.linalg.inv(C))