Resolver un sistema grande no siempre significa terminar una eliminación completa. Podemos empezar con una aproximación, medir lo que falta y corregirla. Haremos las primeras iteraciones a mano, explicaremos cuándo se acercan a la solución y veremos cómo cambiar la escala del problema ayuda. Cada método necesita sus propias hipótesis; obtener un número en pantalla no significa que haya convergido.
Si necesitas repasar:
1. Una aproximación se comprueba sustituyéndola en todas las ecuaciones
Partimos de Ax=b, con A cuadrada e invertible, y llamamos x* a su solución. Un método iterativo construye x⁽⁰⁾, x⁽¹⁾,… usando cálculos más baratos que una factorización completa. El coste total depende del trabajo por paso y del número de pasos necesarios.
El residual se calcula con datos conocidos; el error normalmente no.
No basta con observar que las componentes parecen estabilizarse.
Una primera aproximación no tiene que satisfacer todavía el sistema
4x+y=6, 2x+5y=12; probar (x, y)=(1,5;2,4).
Primera fila
4·1,5+2,4=8,4; el residual de esa fila es 6−8,4=−2,4.Segunda
2·1,5+5·2,4=15; residual 12−15=−3.Referencia
La solución (1,2) sí da seis y doce. El error de la aproximación es (0,5;0,4), distinto del residual.
Residual y error tienen números y significados diferentes.
Comprobar solo la primera ecuación no certifica una solución del sistema.
Una iteración barata puede ahorrar memoria, pero no garantiza éxito
A es dispersa y una factorización puede generar mucho relleno.
Operación base
Un método iterativo puede usar principalmente productos Av y unas pocas soluciones con matrices sencillas.
Coste total
Si un paso cuesta c y se hacen k, el trabajo incluye aproximadamente kc más preparación, comunicación y comprobaciones.
Elegir
Un método que necesite demasiados pasos puede ser peor que una factorización. Se consideran estructura, tolerancia y reutilización de varios segundos miembros.
La palabra iterativo no es por sí sola una garantía de rapidez.
También se decide cómo tratar falta de convergencia y límite de recursos.
2. Una descomposición convierte el sistema en una regla de actualización
En el programa resolvemos con M; no hace falta almacenar una inversa explícita.
Jacobi toma M igual a la diagonal y usa todas las componentes antiguas. Gauss–Seidel utiliza las nuevas según se calculan, lo que corresponde a otra descomposición. Cambiar ese orden cambia el método y su matriz de error.
Despejar cada componente para Jacobi
A=[[4,1],[2,5]], b=(6,12), inicio (0,0).
Regla
x nuevo=(6−y antiguo)/4; y nuevo=(12−2x antiguo)/5.M=diag(4,5) y N tiene filas (0;−1) y (−2;0).Primer paso
Usando ambos valores antiguos cero: (x, y)=(3/2;12/5)=(1,5;2,4).Segundo
Con ese vector antiguo completo: x=(6−2,4)/4=0,9; y=(12−3)/5=1,8.
Los dos componentes del segundo paso son (0,9;1,8).
La matriz B tiene filas (0;−1/4) y (−2/5;0), con signos negativos.
Gauss–Seidel ya usa el primer componente nuevo
El mismo sistema y el mismo inicio (0,0).
Primero x
x nuevo=(6−0)/4=1,5.Después y
Se usa ese x nuevo: y nuevo=(12−2·1,5)/5=1,8.Siguiente paso
x=(6−1,8)/4=1,05 y y=(12−2·1,05)/5=1,98.No son los iterados de Jacobi.
Usar datos nuevos no es un detalle intercambiable.
Aquí B de Gauss–Seidel tiene filas (0;−1/4) y (0;1/10), con radio espectral 1/10.
3. La ecuación del error distingue convergencia general y crecimiento transitorio
En dimensión finita, el método estacionario converge desde todo inicio exactamente cuando el radio espectral ρ(B), máximo módulo de sus autovalores, es menor que uno.
Una norma inducida con ‖B‖<1 es suficiente y da una cota geométrica. No es necesario que esa desigualdad se cumpla en una norma concreta elegida. Una matriz no normal puede amplificar algunos errores durante un tiempo antes de que terminen decayendo.
Obtener una cota con las sumas absolutas de fila
Para Jacobi del ejemplo, B tiene filas (0;−1/4) y (−2/5;0).
Norma infinito
La mayor suma absoluta de fila es max(1/4,2/5)=2/5.Cota
‖e⁽ᵏ⁾‖∞≤(2/5)ᵏ‖e⁽⁰⁾‖∞.Como el factor es menor que uno, hay convergencia desde todo inicio.
Radio
B²=(1/10)I, así sus autovalores son ±1/√10.El radio es aproximadamente 0,3162, menor que la cota de norma 0,4.
El radio y la norma no tienen por qué ser el mismo número.
El punto fijo del ejemplo sigue siendo (1,2), no el vector cero salvo para el error.
Un método convergente puede amplificar los primeros errores
B tiene filas (1/2;10) y (0;1/2), con e⁽⁰⁾=(0,1).
Autovalores
Los dos son 1/2, así Bᵏ→0.Sin embargo ‖B‖∞=10,5.Primeros errores
e⁽¹⁾=(10;1/2), e⁽²⁾=(10;1/4), e⁽³⁾=(7,5;1/8).La norma crece al principio.
Forma general
e⁽ᵏ⁾=(10k(1/2)ᵏ⁻¹;(1/2)ᵏ) para k≥1; ambas componentes terminan tendiendo a cero.
ρ(B)<1 no obliga a una caída monótona en cualquier norma.
No se usa 10,5 como factor de contracción; esa cota no serviría para demostrar decaimiento.
Un inicio especial puede converger sin que exista garantía general
B=diag(1/2,2).
Primer inicio
Con error (1,0), Bᵏe=(2⁻ᵏ,0) tiende a cero.Otro inicio
Con error (0,1), Bᵏe=(0,2ᵏ) crece.Criterio
ρ(B)=2.Se refuta convergencia desde todo inicio, aunque una trayectoria particular se acerque al punto fijo.
Una ejecución afortunada no demuestra una propiedad universal.
El criterio espectral se refiere a todos los errores iniciales.
4. Jacobi necesita diagonal no nula y una garantía de convergencia propia
Si |aᵢᵢ|>Σⱼ≠ᵢ|aᵢⱼ| en cada fila, Jacobi tiene ‖B‖∞<1: la dominancia diagonal estricta por filas es suficiente. No es necesaria. Si alguna diagonal es cero, su despeje no está definido aunque el sistema completo sea invertible.
No se trasladan automáticamente garantías entre métodos. En particular, ser simétrica definida positiva permite el marco clásico de gradiente conjugado, pero no asegura convergencia de Jacobi sin más condiciones.
Fallarle a una condición suficiente no demuestra divergencia
A tiene filas (1;2) y (1/10;1).
Dominancia
La primera fila no es estrictamente dominante: uno no supera dos.
Matriz de Jacobi
B tiene filas (0;−2) y (−1/10;0), luego B²=(1/5)I.Radio
ρ(B)=1/√5<1, así Jacobi converge desde todo inicio. Su norma infinito dos no certifica esa convergencia, pero tampoco la niega.
Una condición suficiente no es una condición necesaria.
A es invertible porque su determinante es 4/5, distinto de cero.
Una matriz positiva puede hacer divergir Jacobi
A es 3×3 con diagonal uno y todas las entradas fuera de diagonal iguales a 0,9.
Positividad
A=(1/10)I+(9/10)11ᵀ tiene autovalores 14/5 y 1/10, este último doble: todos positivos.
Jacobi
Como D=I, B=I−A tiene autovalores −9/5 y 9/10, este último doble.
Diagnóstico
ρ(B)=9/5>1.Por ejemplo, un error en dirección (1,1,1) aumenta por factor de módulo 1,8 en cada paso.
SPD no basta para garantizar este método estacionario.
El mismo problema sí cumple la hipótesis estructural de gradiente conjugado clásico, que es otro algoritmo.
5. Los métodos de Krylov aprovechan subespacios y no una matriz B fija
Desde un residual inicial r₀, un espacio de Krylov reúne combinaciones de r₀, Ar₀, A²r₀,… Los métodos eligen aproximaciones en x₀ más ese espacio. No se analiza cada uno como si repitiera necesariamente una misma B estacionaria.
Gradiente conjugado clásico usa A real simétrica definida positiva y minimiza el error en norma A en el subespacio correspondiente. GMRES minimiza la norma euclídea del residual y admite matrices no simétricas. El almacenamiento de vectores y los reinicios importan. Terminación en dimensión finita se refiere al algoritmo apropiado en aritmética exacta; el redondeo cambia ese ideal.
Dos pasos de gradiente conjugado en un sistema diagonal positivo
A=diag(1,4), b=(1,1), x₀=(0,0).
Primera dirección
r₀=p₀=(1,1). α₀=(r₀ᵀr₀)/(p₀ᵀAp₀)=2/5, luego x₁=(2/5;2/5), r₁=(3/5;−3/5).Nueva dirección
β₀=(r₁ᵀr₁)/(r₀ᵀr₀)=9/25 y p₁=r₁+β₀p₀=(24/25;−6/25).Segundo paso
α₁=5/8.Entonces x₂=x₁+α₁p₁=(1;1/4), solución exacta.Las dos direcciones son conjugadas: p₀ᵀAp₁=0.
El método combina direcciones de acuerdo con la geometría de A.
A tiene autovalores uno y cuatro positivos; esa hipótesis no se omite al trasladar el algoritmo.
GMRES escoge el residual mínimo del subespacio disponible
A=[[1,1],[0,1]], b=(0,1), x₀=0.
Primer espacio
Es span{b}; probar x=(0, α) da residual (−α;1−α).Su norma al cuadrado es α²+(1−α)².
Minimizar
La derivada 4α−2 se anula en α=1/2.El primer iterado es (0;1/2) y el residual (−1/2;1/2).
Segundo espacio
b y Ab=(1,1) son independientes.La solución (−1,1)=2b−Ab pertenece a ese espacio completo y el residual mínimo es cero en aritmética exacta, sin reiniciar.
El objetivo se expresa sobre el residual, no sobre una regla Jacobi.
Reiniciar elimina parte de la información del espacio; no se conserva sin más la garantía de esos dos pasos.
6. Precondicionar transforma el cálculo conservando la solución
En ambos se aplica una solución con M; no es necesario formar su inversa. M debe ser útil por coste y por propiedades del problema transformado.
En gradiente conjugado precondicionado se usa M simétrica definida positiva con una interpretación simétrica: M⁻¹ᐟ²AM⁻¹ᐟ². M⁻¹A puede no ser simétrica en el producto euclídeo ordinario. No se ignora esa distinción al aplicar un método que requiere simetría.
Izquierda y derecha conservan x, pero usan incógnitas y residuos diferentes
A=[[4,1],[2,5]], b=(6,12), M=diag(4,5).
Izquierda
M⁻¹A tiene filas (1;1/4) y (2/5;1); M⁻¹b=(3/2;12/5).La solución x=(1,2) permanece.Derecha
AM⁻¹ tiene filas (1;1/5) y (1/2;1).
Se resuelve para z=Mx=(4,10), y después se recupera x.Residual
En la izquierda, el residual transformado es M⁻¹r; en la derecha, b−AM⁻¹z coincide con el residual original al recuperar x. Hay que decir cuál se está midiendo.
No se compara un residual precondicionado con uno original como si fueran la misma cifra.
Usar la diagonal es barato aquí; la transformación no garantiza por sí sola la calidad de cualquier método.
Una escala diagonal equilibra un problema con dos magnitudes muy distintas
A=diag(1,100), M=diag(1,100).
Antes
El número de condición euclídeo de A es cien y la forma xᵀAx=x₁²+100x₂² tiene niveles elípticos.
Cambio simétrico
u=M¹ᐟ²x=(x₁;10x₂).El operador M⁻¹ᐟ²AM⁻¹ᐟ² es I; su condición es uno.
Alcance
Aquí M=A es fácil porque es diagonal.En una matriz general elegir M=A trasladaría a cada aplicación la dificultad de resolver el sistema original.
El ejemplo enseña el efecto de la escala, no una receta universal gratuita.
La solución física se recupera en las coordenadas x, no se interpreta u como si fuera el mismo vector.
7. Un precondicionador se valora por el ahorro total y el método que lo utiliza
Un precondicionador necesita preparación, memoria y trabajo por aplicación. Una reducción de iteraciones puede no compensar esos gastos. Si se reutiliza para varios segundos miembros, su preparación se amortiza. Se comparan tiempos o costes medidos bajo el mismo objetivo de precisión.
En problemas SPD, cotas de gradiente conjugado dependen de autovalores y condición del operador simétrico apropiado. La distribución espectral también importa. En problemas no normales, mirar solo autovalores o una única condición no predice toda la trayectoria; ya vimos un crecimiento transitorio con radio menor que uno.
Menos iteraciones no siempre basta, pero reutilizar preparación cambia la cuenta
Costes ilustrativos: sin precondicionar, 200 pasos de coste uno; con M, preparación 120 y 20 pasos de coste dos.
Una resolución
Sin M: 200.
Con M: 120+20·2=160.Cinco resoluciones
Si la misma preparación sirve para las cinco, el total con M es 120+5·40=320, frente a 1000 sin M.
Otro problema
Si sin M bastaran diez pasos y con M cinco, una resolución costaría diez frente a 130. Esos recuentos deben medirse; no son una garantía por usar el nombre precondicionador.
La comparación incluye preparación y coste por paso.
Si los segundos miembros o cambios de matriz impiden reutilizar M, se vuelve a contar su preparación.
Una condición grande no obliga a agotar una cota pesimista
A diagonal SPD con solo dos valores distintos: uno y cien, repetidos en cualquier multiplicidad.
Condición
κ₂(A)=100.Polinomio
p(λ)=(1−λ)(1−λ/100) tiene p(0)=1 y anula ambos autovalores.Por tanto p(A)=0.Krylov
En aritmética exacta, gradiente conjugado puede resolver desde cualquier inicio en como máximo dos pasos en este caso, por la estructura espectral. Eso no convierte en innecesarias las cotas de condición para otros espectros.
La distribución de autovalores puede dar información que una sola razón no muestra.
Con redondeo, perder ortogonalidad modifica la terminación exacta ideal.
8. La parada informa qué se ha medido y qué garantía de error existe
Para b y x* no nulos, una cota relativa es ‖e‖/‖x*‖≤κ(A)‖r‖/‖b‖ en normas compatibles. Un residual pequeño exige interpretar escala y condición.
Se conservan tolerancias, máximo de pasos, motivo de salida y residual del sistema original. Un incremento pequeño puede ser estancamiento. Con una contracción conocida, en cambio, el último incremento sí produce una cota: error del punto nuevo ≤q/(1−q) veces ese incremento.
Una dirección mal escalada puede esconder un error grande
A=diag(1,10⁻⁶), b=(1;10⁻⁶), aproximación x=(1,0).
Solución
x*=(1,1); el error es (0,−1), de norma infinito uno.Residual
r=(0;10⁻⁶), muy pequeño en norma absoluta.A⁻¹ tiene diagonal uno y un millón, explicando la amplificación.
Interpretar
κ∞=10⁶ y ‖r‖∞/‖b‖∞=10⁻⁶.Su producto da uno, coherente con el error relativo infinito uno.
El residual con precondicionador M=A sería (0,1).
La cifra pequeña del residual original no certifica seis decimales de solución.
La solución satisface bien una ecuación mal escalada en sentido residual, pero no está cerca en esa segunda componente.
La contracción escalar permite traducir una parada en una cota
x siguiente=x/2+1, desde x₀=0.
Fijo y error
x*=2; x⁽ᵏ⁾=2(1−2⁻ᵏ), con error 2¹⁻ᵏ.Incremento
|x⁽ᵏ⁾−x⁽ᵏ⁻¹⁾|=2¹⁻ᵏ para k≥1.Con q=1/2, la cota q/(1−q) coincide aquí con el incremento.
Tolerancia
Para error no mayor que 10⁻³, k=11 da 2⁻¹⁰=1/1024≈0,0009765625.La garantía usa la contracción, no solo dos decimales parecidos.
El informe puede declarar una precisión certificada por la hipótesis demostrada.
Si se alcanza el máximo de iteraciones sin satisfacer el criterio, se informa límite alcanzado en lugar de convergencia.
Ahora lo resolvemos de principio a fin.
EJEMPLO RESUELTO 1
Jacobi una vez
A=[[4,1],[1,3]], b=(1,2), x₀=(0,0). Calcula un paso Jacobi.
Datos
A=[[4,1],[1,3]], b=(1,2), x₀=(0,0).Se solicita un paso Jacobi.
Despejes
La primera componente nueva es (1−segunda antigua)/4. La segunda nueva es (2−primera antigua)/3.
Primera
Con segunda antigua cero: x₁, primera=1/4.Segunda
Se usa primera antigua cero, no el cuarto recién calculado: x₁, segunda=2/3.
Resultado
El vector tras un paso es (1/4;2/3). Aún no tiene por qué resolver exactamente el sistema.
Control
Gauss–Seidel habría usado el nuevo cuarto y obtenido segunda componente 7/12. Esa diferencia comprueba que aquí sí se aplicó Jacobi.
x₁=(1/4,2/3).
En el primer paso ambas ecuaciones usan solo los valores antiguos cero.
EJEMPLO RESUELTO 2
Una contracción escalar
xₖ₊₁=0,5xₖ+1. Halla su punto fijo.
Punto fijo
En el límite buscado, x siguiente=x; se resuelve x=0,5x+1.Agrupar
x−0,5x=1, es decir, 0,5x=1.Resolver
El único punto fijo es dos.
Restar el fijo
Con eₖ=xₖ−2, eₖ₊₁=0,5xₖ+1−2=0,5eₖ.Iterar el error
eₖ=(0,5)ᵏe₀, que tiende a cero desde cualquier inicio real finito.
Conclusión
La secuencia converge a dos. La existencia del punto fijo se obtuvo con la ecuación; la convergencia se justificó aparte con el factor de error de módulo menor que uno.
Converge a x=2 desde cualquier x₀.
eₖ₊₁=0,5eₖ.
Los símbolos que hemos usado
- x* / xₖ
- Solución exacta / aproximación en el paso k.
- eₖ=xₖ−x*
- Error firmado de la aproximación.
- rₖ=b−Axₖ
- Residual: incumplimiento de la ecuación.
- B
- Matriz que propaga errores en una iteración estacionaria.
- ρ(B)
- Radio espectral: máximo módulo de los autovalores de B.
- M
- Precondicionador o parte fácilmente resoluble de una descomposición, según el contexto.
Ahora te toca a ti.
Tienes cinco ejercicios. Inténtalos; las pistas te ayudan si te atascas.
Practicar este tema →