Sistemas lineales, factorización y condicionamiento

Los sistemas lineales aparecen dentro de casi todos los algoritmos del curso. Newton resuelve un jacobiano por un paso; un método implícito resuelve una ecuación en cada instante; una regresión ajusta parámetros; la estabilidad local depende de una matriz. Escribir $x=A^{-1}b$ es correcto en álgebra, pero rara vez describe el cálculo que conviene hacer. Numéricamente se factoriza $A$ y se resuelven sistemas triangulares.

En dinámica, estos sistemas aparecen al clasificar equilibrios y linealizar campos. Su implementación se contrasta con la documentación oficial de NumPy y se verifica mediante residuo, condicionamiento y experimentos de perturbación.

Eliminación gaussiana

Para resolver $Ax=b$, la eliminación usa operaciones de fila para convertir $A$ en una matriz triangular superior $U$. En forma factorizada,

$$PA=LU,$$

donde $P$ registra intercambios, $L$ los multiplicadores y $U$ la matriz triangular. Se resuelve primero

$$Ly=Pb$$

por sustitución hacia adelante y después $Ux=y$ hacia atrás. La factorización cuesta del orden de $n^3$ operaciones; cada nuevo vector $b$ cuesta del orden de $n^2$. Por eso conviene reutilizarla cuando la matriz no cambia.

El pivoteo parcial elige en cada columna un elemento de magnitud grande como pivote. Evita divisiones por números pequeños y limita crecimiento en muchos problemas prácticos. No arregla una matriz singular y no elimina el condicionamiento inherente; mejora la estabilidad del algoritmo.

Considere

$$A=\begin{pmatrix}10^{-12}&1\\1&1\end{pmatrix}, \qquad b=\begin{pmatrix}1\\2\end{pmatrix}.$$

Eliminar sin intercambiar usa el multiplicador $10^{12}$ y forma diferencias enormes antes de recuperar una respuesta moderada. Intercambiar filas coloca el pivote $1$ arriba y evita ese crecimiento. La solución matemática no cambió; sí cambió la trayectoria aritmética.

numpy.linalg.solve resuelve sistemas cuadrados de rango completo mediante rutinas LAPACK. Requiere revisar forma, tipo y excepciones. Construir np.linalg.inv(A) @ b suele costar más, almacenar una matriz innecesaria y ocultar que el objetivo era resolver.

Criterios de residuo y condicionamiento

Después de calcular $\widehat x$, forme

$$r=b-A\widehat x.$$

Usa los datos originales, no una matriz ya sobrescrita por la factorización. El residuo relativo

$$\eta=\frac{\|r\|}{\|A\|\,\|\widehat x\|+\|b\|}$$

es una medida de error hacia atrás: indica cuánto habría que perturbar el problema para que $\widehat x$ fuera exacta. Una rutina estable suele producir $\eta$ cercano a la precisión de máquina, incluso cuando el error hacia adelante es grande por mala condición.

El número de condición en una norma es

$$\kappa(A)=\|A\|\,\|A^{-1}\|.$$

En norma dos, es la razón entre el mayor y el menor valor singular. Si $\kappa_2(A)=10^8$, perturbaciones relativas del orden de $10^{-10}$ pueden amplificarse hasta el orden de $10^{-2}$ en el peor caso. La cota no predice que toda perturbación alcance ese máximo, pero señala cuántas cifras pueden estar en riesgo.

numpy.linalg.cond estima la condición con varias normas. cond(A) grande no dice que solve esté mal implementado; dice que el problema transmite con fuerza incertidumbre de datos o redondeo hacia la solución.

Sensibilidad ante perturbaciones

Sea

$$A_\varepsilon=\begin{pmatrix}1&1\\1&1+\varepsilon\end{pmatrix}, \qquad b=\begin{pmatrix}2\\2+\varepsilon\end{pmatrix}.$$

La solución es $(1,1)^T$. Reste la primera ecuación de la segunda: $\varepsilon x_2=\varepsilon$. Cuando $\varepsilon$ se aproxima a la precisión de los datos, esa resta extrae información de una diferencia diminuta.

Perturba sólo $b_2$ a $2+\varepsilon+\delta$. Entonces

$$x_2=1+\frac{\delta}{\varepsilon}, \qquad x_1=1-\frac{\delta}{\varepsilon}.$$

Una perturbación $\delta$ pequeña se amplifica por $1/\varepsilon$. El residuo respecto del sistema perturbado puede ser diminuto y, aun así, la solución diferir mucho de la original. Éste es condicionamiento del problema, no inestabilidad necesaria del algoritmo numérico.

Escalar filas o columnas puede reducir disparidades de unidades y mejorar el comportamiento del cálculo, pero la transformación debe registrarse. Si $D_rAD_c\widetilde x=D_rb$ y $x=D_c\widetilde x$, la solución física se recupera al final. El escalamiento no inventa información que una matriz casi singular no contiene.

Matriz jacobiana y dinámica linealizada

Hirsch, Smale y Devaney desarrollan álgebra lineal y sistemas de dimensión $n$ antes de volver al caso no lineal. La matriz exponencial $e^{At}$ resuelve $x'=Ax$; autovalores y subespacios organizan contracción y expansión. Hirsch, 2013 Morris W. Hirsch, Stephen Smale y Robert L. Devaney (2013) Differential Equations, Dynamical Systems, and an Introduction to Chaos 3.ª ed. · Academic Press Ubicación consultada: caps. 5–6, pp. 73–138

En Newton, el sistema

$$J_F(x_k)s_k=-F(x_k)$$

puede volverse mal condicionado cerca de una bifurcación donde el jacobiano pierde rango. Un paso grande no debe recortarse sin diagnóstico: revisa valores singulares, condición, residuo y la geometría del problema. En un método implícito, matrices de la forma $I-hJ$ cambian con $h$; su factorización puede reutilizarse sólo si la implementación conoce cuándo la matriz sigue siendo válida.

Fuchs repasa sistemas lineales, autovalores y autovectores en el bloque matemático que sostiene sus modelos. Es una referencia para la interpretación, no para atribuirle los detalles internos de LAPACK. Fuchs, 2013 Armin Fuchs (2013) Nonlinear Dynamics in Complex Systems: Theory and Applications for the Life-, Neuro- and Natural Sciences Springer Ubicación consultada: cap. 10, secciones 10.2–10.3, pp. 185–193

Estructura matricial y múltiples términos independientes

Si $A$ es simétrica definida positiva, la factorización de Cholesky $A=LL^T$ usa aproximadamente la mitad del trabajo de LU y preserva estructura. Si hay múltiples términos independientes $B=[b_1,\ldots,b_m]$, se resuelve $AX=B$ en bloque. Para matrices grandes y esparsas, almacenar todos los ceros y usar factorización densa desperdicia memoria; se eligen formatos y solucionadores esparsos.

No toda matriz simétrica es definida positiva. Intentar Cholesky y recibir un fallo puede revelar que la hipótesis era falsa o que la matriz perdió definitud por redondeo. Antes de añadir arbitrariamente un múltiplo de la identidad, hay que entender qué propiedad del modelo se está modificando.

En problemas de mínimos cuadrados con más ecuaciones que incógnitas, las ecuaciones normales $A^TAx=A^Tb$ elevan aproximadamente al cuadrado el número de condición. QR o SVD suelen ser preferibles. numpy.linalg.solve no es la rutina para una matriz rectangular; la forma del problema decide el algoritmo.

Pruebas del solucionador y del problema

1. Construye una solución exacta $x_\star$ y define $b=Ax_\star$. Compara $\widehat x$ con $x_\star$ y calcula el residuo.
2. Perturba $b$ en direcciones aleatorias y en la dirección singular más sensible.
3. Repite en float32 y float64.
4. Escala filas y columnas y recupera las variables originales.
5. Compara solve con formar la inversa sólo como experimento de costo y error, no como alternativa recomendada.

Slotine y Li usan sistemas lineales y funciones cuadráticas dentro del análisis no lineal de estabilidad. Ese contexto muestra por qué una factorización fiable es parte de un argumento mayor: el cálculo debe sostener la interpretación del equilibrio, no reemplazarla. Slotine, 1991 Jean-Jacques E. Slotine y Weiping Li (1991) Applied Nonlinear Control Prentice Hall Ubicación consultada: caps. 2–3, pp. 17–94

Ejercicios de condicionamiento y validación

Resolver un sistema lineal no consiste en obtener un vector sin error de ejecución. La respuesta necesita residuo, condición, escala, estructura de la matriz y procedencia de los datos; juntos indican cuántas cifras y qué interpretación pueden defenderse.