Solución de sistemas de ecuaciones lineales

Solucionar un sistema de ecuaciones lineales es encontrar el valor de las variables que satisfacen todas las ecuaciones del sistema. En este capítulo se estudiarán los métodos de solución de sistemas de ecuaciones lineales.

Sea \(A\) una matriz de \(n \times n\) y \(b\) un vector de \(n \times 1\). El sistema de ecuaciones lineales \(Ax = b\) se puede escribir como:

\[\begin{bmatrix}a_{1,1}&a_{1,2}\dots&a_{1,n}\\ a_{2,1}&a_{2,2}\dots&a_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n,1}&a_{n,2}\dots&a_{n,n}\\ \end{bmatrix}\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix}\]

Definiremos \(A_{\cdot,j}\) como la columna \(j\) de la matriz \(A\) y \(A_{i,\cdot}\) como la fila \(i\) de la matriz \(A\).

Sistemas de ecuaciones lineales con matrices triangulares

Sea \(A\) una matriz triangular superior de \(n \times n\) y \(b\) un vector de \(n \times 1\). El sistema de ecuaciones lineales \(Ax = b\) se puede escribir como:

\[\begin{bmatrix}a_{1,1}&a_{1,2}\dots&a_{1,n}\\ 0&a_{2,2}\dots&a_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0\dots&a_{n,n}\\ \end{bmatrix}\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix}\]

El sistema de ecuaciones lineales se puede resolver por sustitución hacia atrás. El algoritmo para resolver el sistema de ecuaciones lineales es:

Inicializar \(x_n = \frac{b_n}{a_{n,n}}\).
Para \(i = n-1, n-2, \dots, 1\):
- Calcular \(x_i = \frac{b_i - \sum_{j=i+1}^n a_{i,j}x_j}{a_{i,i}}\).

Sea \(A\) una matriz triangular inferior de \(n \times n\) y \(b\) un vector de \(n \times 1\). El sistema de ecuaciones lineales \(Ax = b\) se puede escribir como:

\[\begin{bmatrix}a_{1,1}&0\dots&0\\ a_{2,1}&a_{2,2}\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ a_{n,1}&a_{n,2}\dots&a_{n,n}\\ \end{bmatrix}\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix}\]

El sistema de ecuaciones lineales se puede resolver por sustitución hacia adelante. El algoritmo para resolver el sistema de ecuaciones lineales es:

Inicializar \(x_1 = \frac{b_1}{a_{1,1}}\).
Para \(i = 2, 3, \dots, n\):
- Calcular \(x_i = \frac{b_i - \sum_{j=1}^{i-1} a_{i,j}x_j}{a_{i,i}}\).

Factorización LU

La factorización LU es una forma de factorizar una matriz \(A\) en el producto de dos matrices \(L\) y \(U\), donde la matriz \(L\) es una matriz triangular inferior y la matriz \(U\) es una matriz triangular superior. La factorización LU se puede utilizar para resolver sistemas de ecuaciones lineales.

Teorema Una matriz \(A\) cuadrada tiene factorización LU única si y solo si todos los menores principales de \(A\) son diferentes de cero.

Supongamos que queremos resolver el sistema de ecuaciones \[Ax = b\] Sea \(A=LU\), así podemos escribir el sistema de ecuaciones como \[LUx = b\] Definamos \(y = Ux\), así podemos escribir el sistema de ecuaciones como \[Ly = b\] Así resolver el sistema original se convierte en solucionar dos sistemas triangulares, lo cual computacionalemente es más eficiente.

Algoritmo de factorización LU

El algoritmo de factorización LU es el siguiente:

Inicializar \(L = I\) y \(U = A\).
Para \(i = 1, 2, \dots, n-1\):
- Para \(j = i+1, i+2, \dots, n\):
  - Calcular \(L_{j,i} = \frac{U_{j,i}}{U_{i,i}}\).
  - Para \(k = i, i+1, \dots, n\):
    - Calcular \(U_{j,k} = U_{j,k} - L_{j,i}U_{i,k}\).

Ejemplo en pyrhon

import numpy as np
## algortimo LU
def LU(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy()
    for i in range(n-1):
        for j in range(i+1, n):
            L[j,i] = U[j,i]/U[i,i]
            for k in range(i, n):
                U[j,k] = U[j,k] - L[j,i]*U[i,k]
    return L, U

Note De forma vectorial podemos escribir el algoritmo de la siguiente forma:


## algortimo LU
import numpy as np

def LU(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy()
    
    for i in range(n - 1):
        L[i+1:, i] = U[i+1:, i] / U[i, i]
        U[i+1:, i:] -= np.outer(L[i+1:, i], U[i, i:])
    
    return L, U

Pivotación parcial

La pivotación parcial es un método para evitar que el algoritmo de factorización LU falle. La pivotación parcial consiste en intercambiar filas de la matriz \(A\) de tal forma que el elemento de la diagonal que se está utilizando para calcular los multiplicadores de la columna sea el mayor posible.

El algoritmo de factorización LU con pivotación parcial es el siguiente:

Inicializar \(L = I\) y \(U = A\).
Para \(i = 1, 2, \dots, n-1\):
- Encontrar \(j\) tal que \(|U_{j,i}| = \max_{k=i, i+1, \dots, n} |U_{k,i}|\).
- Intercambiar las filas \(i\) y \(j\) de \(U\) y \(L\).
- Para \(j = i+1, i+2, \dots, n\):
  - Calcular \(L_{j,i} = \frac{U_{j,i}}{U_{i,i}}\).
  - Para \(k = i, i+1, \dots, n\):
    - Calcular \(U_{j,k} = U_{j,k} - L_{j,i}U_{i,k}\).

algortimo LU con pivotación parcial

import numpy as np

def LU_partial_pivot(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy()
    P = np.eye(n)

    for i in range(n - 1):
        pivot_row = np.argmax(np.abs(U[i:, i])) + i
        
        if pivot_row != i:
            U[[i, pivot_row]] = U[[pivot_row, i]]
            P[[i, pivot_row]] = P[[pivot_row, i]]
            
            if i > 0:
                L[[i, pivot_row], :i] = L[[pivot_row, i], :i]

        L[i+1:, i] = U[i+1:, i] / U[i, i]
        U[i+1:, i:] -= np.outer(L[i+1:, i], U[i, i:])
    
    return P, L, U

Ventajas del pivoteo parcial:

Estabilidad numérica Ayuda a prevenir errores de redondeo y la propagación de errores en el cálculo numérico.
Prevención de divisiones por cero
Mayor robustez
Eficiencia en sistemas grandes Aunque el pivoteo parcial puede aumentar la cantidad de operaciones realizadas en comparación con el pivoteo simple (sin pivoteo), en sistemas grandes, los beneficios en términos de estabilidad numérica superan con creces el costo adicional computacional.

Ejercicios

Encuentre la factorización LU a mano de \[A=\left(\begin{array}{ccc} 1&-3&0\\ -1&0&-9\\ 4&6&0 \end{array}\right), \]
La siguiente matriz no tiene inversa, por tanto no tiene factorizacion LU única, encuentre dos factorizaciones diferentes de la matriz \[A=\left(\begin{array}{ccc} 1&-3&0\\ -2&0&-9\\ 1&3&9 \end{array}\right), \]
Sea el sistema

\[\begin{cases} 10^{-14}x+y&=1\\ x+y&=2 \end{cases} \]

a.  Grafique cada ecuación en un mismo plano cartesiano. ¿Observa algo particular en la gráfica?

b. Resuelva el  sistema de ecuaciones exactamente

c.  Resuelva el sistema de ecuaciones sin usar Pivoteo por medio de LU

d.  Resuelva el  sistema de ecuaciones usando Pivoteo por medio de LU

e.  Observe las soluciones anteriores, realice una comparación entre ellas ¿Que puede observa? De una conclusión del porque hay tienen diferentes respuestas.

Sea la matriz (use format short) \[A=\left(\begin{array}{ccc} 1276929.38 &-930737.61&-617.4\\ −1930737.61 & 3316209.17&1755.3\\ −617.4& 1755.3&2 \end{array}\right), b=\left(\begin{array}{c}2\\1\\3\end{array}\right) \]
1. Resuelva el siguiente sistema de ecuaciones \(Ax=b\) sin usar Pivoteo por medio de LU
2. Resuelva el siguiente sistema de ecuaciones \(Ax=b\) usando Pivoteo por medio de LU
3. Observe las soluciones anteriores, compruebe la solución de la forma \(b-A*x\), realice una comparación entre ellas ¿Que puede observar?
Usando la factorización LU, escriba un pseudocódigo para encontrar la inversa de una matriz. Prográmelo.
Para calcular el esfuerzo de una viga normalmente se tiene que resolver un sistema de ecuaciones de la forma \[\left(\begin{array}{cccccccc} 4&2&1\\ 2&4&2&1\\ 1&2&4&2&1\\ &1&2&4&2&1\\ &&1&2&4&2&1\\ &&&1&2&4&2\\ &&&&1&2&4\\ \end{array}\right)\left(\begin{array}{c} x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7 \end{array}\right)=\left(\begin{array}{c} 4\\1\\1\\1\\1\\1\\4 \end{array}\right) \]

¿Para resolver este sistema de ecuaciones que método usaría?
Usando el programa de LU sin pivoteo encuentre la inversa de la matriz
Usando el programa de LU con pivoteo encuentre la inversa de la matriz
Usando el programa de Cholesky encuentre la inversa de la matriz
¿Cual de los métodos requiere menos esfuerzo computacional? Compare su respuesta con las del punto a.

Sea la siguiente matriz, Note que ella es simétrica y se puede probar que es definida positiva

\[\begin{bmatrix} 1&1/2&1/3&1/4&1/5&1/6&1/7&1/8&1/9&1/10\\ 1/2&1/3&1/4&1/5&1/6&1/7&1/8&1/9&1/10&1/11\\ 1/3&1/4&1/5&1/6&1/7&1/8&1/9&1/10&1/11&1/13\\ 1/4&1/5&1/6&1/7&1/8&1/9&1/10&1/11&1/13&1/14\\ 1/5&1/6&1/7&1/8&1/9&1/10&1/11&1/13&1/14&1/15\\ 1/6&1/7&1/8&1/9&1/10&1/11&1/13&1/14&1/15&1/16\\ 1/7&1/8&1/9&1/10&1/11&1/13&1/14&1/15&1/16&1/17\\ 1/8&1/9&1/10&1/11&1/13&1/14&1/15&1/16&1/17&1/18\\ 1/9&1/10&1/11&1/13&1/14&1/15&1/16&1/17&1/18&1/19\\ 1/10&1/11&1/13&1/14&1/15&1/16&1/17&1/18&1/19&1/20\\ \end{bmatrix}\begin{bmatrix} x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8\\x_9\\x_{10} \end{bmatrix}=\begin{bmatrix} 4\\1\\1\\1\\1\\1\\1\\1\\1\\3 \end{bmatrix}\]

Para obtener esta matriz puede usar desde scipy.linalg la función hilbert

Resuelva el sistema usando Cholescky y compruebe su respuesta
resuelva el sistema usando LU y compruebe su respuesta
compre las dos repuestas. ¿Son diferentes?
Repita esta proceso con hilbert(1000), ¿que pasa? ¿Por que cree que sucede esto?

Número de condición de una matriz

Sea \(A\) una matriz cuadrada de \(n \times n\), notamos como \(\|A\|\) la norma de la matriz \(A\) y \(\|A^{-1}\|\) es la norma de la matriz inversa de \(A\).

Teorema Si \(A\) es una matriz invertible, entonces \(\kappa_2(A) = \frac{\lambda_{\max}(A)}{\lambda_{\min}(A)}\), donde \(\lambda_{\max}(A)\) es el mayor valor propio de \(A\) y \(\lambda_{\min}(A)\) es el menor valor propio de \(A\).

Factor de error de magnificación El factor de error de magnificación es una medida de la magnitud del error relativo en la solución de un sistema de ecuaciones lineales, \(Ax=b\). El factor de error de magnificación se define como:

Sea \(x_a\) la solución aproximada del sistema de ecuaciones lineales \(Ax=b\) y \(x\) la solución exacta del sistema de ecuaciones lineales \(Ax=b\). El residual se define como \(r = b - Ax_a\). El error hacia atrás es la norma del residual \(\|r\|=\|b-Ax\|\). El error hacia adelante es la norma de la diferencia entre la solución aproximada y la solución exacta \(\|x_a - x\|\).

Ejemplo

Encontrar el error hacia adelante y el error hacia atrás para la solución aproximada \(x_a=[1,1]\) del sistema

\[\begin{bmatrix}1&1\\3&-4\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}3\\2\end{bmatrix}\] + Encuentre la solución exacta del sistema de ecuaciones lineales. \([2,1]\) + Encuentre el error hacia atrás. + Encuentre el error hacia adelante. Ejemplo Encuentre los errores hacia adelante y hacia atrás para la solución aproximada \(x_a=[−1, 3.0001]\) de el sistema

\[\begin{bmatrix}1&1\\1.0001&1\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}2\\2.0001\end{bmatrix}\] + Encuentre la solución exacta del sistema de ecuaciones lineales, \([1,1]\) + Encuentre el error hacia atrás. + Encuentre el error hacia adelante. Resolviendo el sistema tenemos \[x =[1,1]\]

Definición Denotamos el error relativo hacia atrás como \(\frac{\|r\|}{\|b\|}\) y el error relativo hacia adelante como \(\frac{\|x_a - x\|}{\|x\|}\).

Usando la anteriro definición podemos introducir el concepto de factor de error de magnificación. El factor de error de magnificación es una medida de la magnitud del error relativo en la solución de un sistema de ecuaciones lineales, \(Ax=b\). El factor de error de magnificación se define como:

\[factor\ magnificación\ error=\frac{error_{realitco\ hacia\ delante}}{error_{realitco\ hacia\ atras}}=\frac{\frac{||x-x_a||_{\infty}}{||x||_\infty}}{\frac{||r||_\infty}{||b||_\infty}}\]

en el sistema anterior tenemos que \[40004.0001\%\]

El número de condición de una matriz es una medida de la sensibilidad de la solución de un sistema de ecuaciones lineales a pequeños cambios en los coeficientes de las ecuaciones. El número de condición de una matriz se define como:

\[\kappa(A) = \|A\| \|A^{-1}\|\]

El numero de condición de una matriz cuadrada se puede interpretar como el máximo factor de ampliación del error posible para resolver \(Ax = b\), sobre todos los lados derechos \(b\).

Ejemplo Calculemos el numero de condición de la matriz \[\begin{bmatrix}1&1\\1.0001&1\end{bmatrix}\]

import numpy as np

# Define the matrix
A = np.array([[1, 1], [1.0001, 1]])

# Calculate the condition number
cond = np.linalg.cond(A)

print("The condition number of the matrix is:", cond)

The condition number of the matrix is: 40002.00007489988

El significado del número de condición es el mismo que en el Capítulo 1. Error de factor de magnificación de una magnitud \(cond( A)\) son posibles. En aritmética de coma flotante, No se puede esperar que el error relativo hacia atrás sea menor que \(\in_{mach}\), ya que almacenar las entradas de \(b\) ya causan errores de ese tamaño. Según el factor de magnificación, la fuerza relativa Son posibles errores de tamaño \(\in_{mach}\cdot cond(A)\) al resolver \(Ax = b\). En otras palabras, si \(cond(A) \aprox 10^k\), deberíamos prepararnos para perder \(k\) dígitos de precisión al calcular \(x\).