41.2: Resolver sistemas lineales

import numpy as np
import sympy as sym

Cuando tenemos el mismo número de ecuaciones que incógnitas, tenemos el siguiente sistema\(Ax=b\) con una matriz cuadrada\(A\). En esta sección, asumimos que la matriz\(A\) tiene rango completo. Es decir, el sistema tiene una solución única. Hablamos de muchas formas de resolver este sistema de ecuaciones. Algunos ejemplos son:

• Iteración jacobia/Iteración de Gauss-Seidel
• \(x=A^{−1}b\)
• Eliminación gaussiana
• Descomposición LU

En esta tarea, mostraremos que podemos resolverlo por QR decomposion.

Considere el siguiente sistema de ecuaciones:

\[\begin{split}\begin{bmatrix}5&-2&2 \\ 4 & -3 &4 \\ 4& -6 &7 \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}1\\2\\3\end{bmatrix}\end{split} \nonumber \]

A = np.matrix([[5, -2, 2], [4, -3, 4], [4,-6,7]])
b = np.matrix([[1],[2],[3]])
display(sym.Matrix(A))
display(sym.Matrix(b))

Sustitución de espalda

Primero implementemos la sustitución de espalda en Python. La función de sustitución inversa back_subst resuelve el sistema\(Rx=b\) donde\(R\) se encuentra una matriz triangular superior.

Cuando resolvemos para\(x\), comenzamos con\(x_n\):\(x_n = {b_n\over R_{n,n}}\). Entonces resolvemos para\(x_{n-1}\) como\(x_{n-1} = {b_{n-1}-R_{n-1,n}x_n\over R_{n-1,n-1}}\). Entonces podemos encontrar\(x_{n-2},\cdots,x_1\):\(x_{n-2} = {b_{n-2}-R_{n-2,n-1}x_{n-1}-R_{n-2,n}x_n\over R_{n-2,n-2}}\). Podemos resolver para todos los componentes de\(x\) en el orden inverso. Entonces a esto se le llama sustitución posterior.

Hacer esto

Complete el siguiente código para la sustitución inversa.

def back_subst(R,b):
    n = R.shape[0]; x = np.zeros(n);
    for i in reversed(range(n)):
        x[i] = b[i]
        for j in range(i+1,n):            
## Your code starts ##            
            x[i] =   # Complet this line of code.
## Your code ends ##            
        x[i] = x[i]/R[i,i]
    return np.matrix(x).T

Hacer esto

Cuando resolvemos para\(Ax=b\) con descomposición QR. Tenemos los siguientes pasos:

• Encuentra la descomposición QR de\(A\)
• De\(QRx=b\), obtenemos\(Rx=Q^{\top}b\)
• Resuelve por\(x\) usar la sustitución de espalda.

Utilice estos pasos para resolver\(Ax=b\) con lo dado\(A\) y\(b\). Compara el resultado con np.linalg.solve.

## Your code starts ##
x =  
## Your code ends ##
print(type(x))   # x is a column vector in the np.matrix type
np.allclose(x, np.linalg.solve(A,b))