1. ¿Por qué los autovalores son críticos en sistemas dinámicos?

En la teoría de sistemas lineales, la ecuación de estado suele escribirse como:

\dot{x}(t) = A\,x(t)

donde A es una matriz cuadrada que describe la dinámica del sistema. Los autovalores (eigenvalues) de A determinan:

• Estabilidad: si la parte real de todos los autovalores es negativa, el sistema es asintóticamente estable.
• Frecuencias naturales: la parte imaginaria indica oscilaciones sinusoidales.
• Respuesta transitoria y velocidad de convergencia.

Por ello, el cálculo preciso y rápido de los autovalores es una tarea esencial en control, análisis modal y simulación de procesos.

2. Algoritmos clásicos para obtener autovalores

Existen varios enfoques, cada uno con ventajas y limitaciones. A continuación, una tabla comparativa (Bootstrap table) que ayuda a decidir cuál usar según el caso de uso.

En la práctica, numpy.linalg.eig y scipy.linalg.eig emplean variantes del algoritmo QR optimizado, mientras que scipy.sparse.linalg.eigs usa Arnoldi para casos dispersos.

3. Implementación paso a paso en Python

3.1. Preparación del entorno

# Instalación de dependencias (Python 3.10+ recomendado)
python -m pip install numpy scipy matplotlib

3.2. Ejemplo completo: análisis de estabilidad de un sistema de segundo orden

Consideremos la siguiente matriz de estado:

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
# Matriz del sistema (ejemplo clásico de masa‑resorte‑amortiguador)
A = np.array([[0, 1],
              [-2, -3]])

3.2.1. Cálculo de autovalores con numpy.linalg.eig

# Autovalores y autovectores
w, v = np.linalg.eig(A)
print("Autovalores:", w)
print("Autovectores (columnas):\n", v)

Salida esperada:

Autovalores: [-1.5+1.32287566j -1.5-1.32287566j]

Como la parte real es -1.5, el sistema es asintóticamente estable (todos los modos decaen).

3.2.2. Visualización del espectro en el plano complejo

plt.figure(figsize=(6,4))
plt.axhline(0, color='gray', linewidth=0.5)
plt.axvline(0, color='gray', linewidth=0.5)
plt.scatter(w.real, w.imag, color='red', s=80, marker='x')
plt.title('Espectro de autovalores de A')
plt.xlabel('Parte real')
plt.ylabel('Parte imaginaria')
plt.grid(True, which='both', linestyle='--', alpha=0.7)
plt.show()

3.2.3. Power Iteration (autovalor dominante)

def power_iteration(M, num_iter=1000, tol=1e-10):
    n = M.shape[0]
    b_k = np.random.rand(n)
    b_k = b_k / np.linalg.norm(b_k)
    for _ in range(num_iter):
        b_k1 = M @ b_k
        b_k1_norm = np.linalg.norm(b_k1)
        b_k_next = b_k1 / b_k1_norm
        if np.linalg.norm(b_k_next - b_k) < tol:
            break
        b_k = b_k_next
    eigenvalue = b_k.T @ M @ b_k
    return eigenvalue, b_k
lambda_max, vec_max = power_iteration(A)
print('Autovalor dominante (Power Iteration):', lambda_max)

En matrices no simétricas, el método converge al autovalor con mayor módulo, que puede no coincidir con el de mayor parte real.

3.2.4. Cálculo rápido de algunos autovalores en matrices dispersas

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs
# Convertimos A a formato disperso (solo a modo de ejemplo)
A_sparse = csr_matrix(A)
# Calculamos los 2 autovalores con mayor módulo
w_sparse, _ = eigs(A_sparse, k=2, which='LM')
print('Autovalores (sparse, Arnoldi):', w_sparse)

4. Buenas prácticas, rendimiento y estabilidad numérica

• Escalado de la matriz: antes de aplicar QR o Jacobi, escalar la matriz (norma 1 o infinito) reduce el riesgo de overflow/underflow.
• Shift estratégico en QR: el uso del Wilkinson shift acelera la convergencia y mejora la precisión de los autovalores reales.
• Precisión de punto flotante: para problemas sensibles, prefiera np.float64 o np.longdouble si la plataforma lo permite.
• Validación de resultados: compare la suma de los autovalores con el rastro (trace(A)) y el producto con el determinante (det(A)) para detectar errores de cálculo.
• Uso de bibliotecas probadas: siempre que sea posible, delegue a numpy.linalg o scipy.linalg, que están optimizadas con LAPACK/BLAS.
• Memoria en matrices grandes: para n > 10⁴, prefiera algoritmos basados en subespacios Krylov (Arnoldi/Lanczos) y estructuras dispersas.

5. Troubleshooting común

5.1. Convergencia lenta en Power Iteration

Posibles causas y soluciones:

1. El cociente entre el mayor y el segundo mayor módulo de autovalor está cercano a 1. Solución: use Rayleigh Quotient Iteration o deflación.
2. Matriz mal condicionada. Solución: aplicar pre‑condicionamiento (por ejemplo, LU o QR).

5.2. Valores complejos inesperados

Si la matriz debería ser simétrica pero aparecen componentes imaginarias, verifique:

• Redondeo numérico: use np.real_if_close para eliminar pequeñas partes imaginarias.
• Errores de entrada: una pequeña asimetría rompe la hermiticidad.

5.3. Fallos de memoria en QR para matrices muy grandes

Utilice versiones “out‑of‑core” o algoritmos de subespacio (Arnoldi) y asegúrese de trabajar con tipos float32 cuando la precisión lo permite.