Métodos de factorización#
Factorización LU#
Si se necesita resolver el sistema de la forma \(\bf AX = B\), se puede utilizar el método de descomposición triangular \(\bf A=LU\) y reescribir el sistema como
(22)#\[\begin{equation}
\bf LUX = B.
\end{equation}\]
donde \(\bf U\) es una matriz triangular superior y \(\bf L\) es una matriz triangular inferior.
Considere un sistema de ecuaciones de dimensión 3, tenemos
(23)#\[\begin{equation}
\begin{bmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{bmatrix}
\begin{bmatrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33} \\
\end{bmatrix}=
\begin{bmatrix}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33} \\
\end{bmatrix}
\end{equation}\]
Introduciendo una variable auxiliar \(\bf Y\), tenemos
(24)#\[\begin{equation}
\bf UX = Y, \quad LY = B.
\end{equation}\]
Dado un vector \(\bf B\), se despeja \(\bf Y\) del sistema \(\bf LY = B\) y luego se despeja \(\bf X\) del sistema \(\bf UX = Y\).
Los valores diagonales de \(\bf L\) son iguales a 1, por lo tanto, el sistema de ecuaciones resultante es
(25)#\[\begin{equation}
l_{21}u_{11} = a_{21} \quad \rightarrow \quad l_{21} = \frac{a_{21}}{u_{11}},
\end{equation}\]
(26)#\[\begin{equation}
l_{31}u_{11} = a_{31} \quad \rightarrow \quad l_{31} = \frac{a_{31}}{u_{11}},
\end{equation}\]
(27)#\[\begin{equation}
l_{31}u_{12} + l_{32}u_{22} = a_{32} \quad \rightarrow \quad l_{32} = \frac{\left(a_{32} - l_{31}u_{12}\right)}{u_{22}}.
\end{equation}\]
Truco
Este algoritmo realiza una cantidad de operaciones igual que el método de eliminación gaussiana. Por tal razón, se recomienda para uso general donde sea necesario algún método directo.
import numpy as np
import matplotlib.pyplot as plt
# Suprimir advertencias menores de integración durante los bucles de optimización
import warnings
warnings.filterwarnings("ignore")
# Configurar el estilo de graficación para una mejor visualización
plt.style.use('ggplot')
# Matriz de coeficientes del sistema de ecuaciones lineales.
A = np.array([
[9, 3, 3, 0, 1],
[4, 4, 5, 5, 8],
[3, 5, 8, 3, 6],
[8, 4, 4, 6, 1],
[5, 6, 8, 2, 1]
], dtype=float)
# Vector de términos independientes (entradas) del sistema.
B = np.array([4, 0, 6, 8, 5], dtype=float)
# Número de incógnitas del sistema de ecuaciones (calculado como el rango de A).
N = np.linalg.matrix_rank(A)
# Proceso de eliminación gaussiana para obtener U
U = A.copy() # Inicialmente U es una copia de la matriz A.
for k in range(N - 1):
for m in range(k + 1, N):
MT = -U[m, k] / U[k, k] # Multiplicadores
U[m, :] = U[m, :] + MT * U[k, :] # Modificación de la matriz U
# Proceso para obtener L (Matriz Triangular Inferior)
L = np.eye(N) # Se inicializa con una matriz identidad.
# Se resuelve la primera columna de L (índice 0 en Python).
L[1:N, 0] = A[1:N, 0] / U[0, 0]
for k in range(1, N - 1):
for m in range(k + 1, N):
MT = 0
for l in range(k):
# Sumatoria de elementos conocidos (adaptada a indexación de 0)
MT = MT + L[m, l] * U[l, k]
# Se calcula el elemento (m, k) de la matriz L.
L[m, k] = (A[m, k] - MT) / U[k, k]
# Proceso de sustitución progresiva (LY = B)
Y = np.zeros(N)
Y[0] = B[0] # El primer elemento del vector auxiliar Y.
for k in range(1, N):
suma_y = 0
for m in range(k):
suma_y = suma_y - L[k, m] * Y[m]
Y[k] = B[k] + suma_y
# Proceso de sustitución regresiva (Ux = y)
x = np.zeros(N)
x[N - 1] = Y[N - 1] / U[N - 1, N - 1] # Se calcula la última variable.
for k in range(N - 2, -1, -1):
suma_x = 0
ind = N - 1 - k
for m in range(ind):
suma_x = suma_x - U[k, k + 1 + m] * x[k + 1 + m]
x[k] = (Y[k] + suma_x) / U[k, k]
# Mostrar resultados
print("Matriz L:")
print(L)
print("\nMatriz U:")
print(U)
print("\nVector solución x:")
print(x)
Matriz L:
[[ 1. 0. 0. 0. 0. ]
[ 0.44444444 1. 0. 0. 0. ]
[ 0.33333333 1.5 1. 0. 0. ]
[ 0.88888889 0.5 -0.33333333 1. 0. ]
[ 0.55555556 1.625 0.25 -2.5 1. ]]
Matriz U:
[[ 9. 3. 3. 0. 1. ]
[ 0. 2.66666667 3.66666667 5. 7.55555556]
[ 0. 0. 1.5 -4.5 -5.66666667]
[ 0. 0. 0. 2. -5.55555556]
[ 0. 0. 0. 0. -24.30555556]]
Vector solución x:
[ 1.256 -7.072 4.95657143 1.22857143 -0.95771429]
Factorización de Doolittle-Crout#
La matriz \(\bf A\) se factoriza en una triangular inferior \(\bf D\) y en una triangular superior \(\bf C\), con unos en la diagonal que se obtiene haciendo una eliminación Gauss-Jordan, quedando como \(\bf DC = A\).
Considere la siguiente matriz
(28)#\[\begin{equation}
A = \begin{bmatrix}
4 & 1 & 1 & 1 \\
1 & 3 & -1 & 1 \\
1 & -1 & 2 & 0 \\
1 & 1 & 0 & 2 \\
\end{bmatrix}.
\end{equation}\]
Si aplicamos Gauss-Jordan hasta poner unos en la diagonal, tenemos
(29)#\[\begin{equation}
C = \begin{bmatrix}
1 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} \\ \\
0 & 1 & -\frac{5}{11} & \frac{3}{11} \\ \\
0 & 0 & 1 & \frac{1}{13} \\ \\
0 & 0 & 0 & 1 \\
\end{bmatrix}.
\end{equation}\]
Dado que \(\bf DC = A\), tenemos
(30)#\[\begin{equation}
\begin{bmatrix}
D_{11} & 0 & 0 & 0 \\ \\
D_{21} & D_{22} & 0 & 0 \\ \\
D_{31} & D_{32} & D_{33} & 0 \\ \\
D_{41} & D_{42} & D_{43} & D_{44} \\
\end{bmatrix}
\begin{bmatrix}
1 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} \\ \\
0 & 1 & -\frac{5}{11} & \frac{3}{11} \\ \\
0 & 0 & 1 & \frac{1}{13} \\ \\
0 & 0 & 0 & 1 \\
\end{bmatrix}=
\begin{bmatrix}
4 & 1 & 1 & 1 \\ \\
1 & 3 & -1 & 1 \\ \\
1 & -1 & 2 & 0 \\ \\
1 & 1 & 0 & 2 \\
\end{bmatrix}.
\end{equation}\]
Multiplicando los cuatro renglones por la primera columna se obtienen las siguientes ecuaciones
(31)#\[\begin{equation}
{\bf R_{1}}{\bf C_{1}} = A_{11}, \quad D_{11}(1) = 4, \quad D_{11} = 4,
\end{equation}\]
(32)#\[\begin{equation}
{\bf R_{2}}{\bf C_{1}} = A_{21}, \quad D_{21}(1) = 1, \quad D_{21} = 1,
\end{equation}\]
(33)#\[\begin{equation}
{\bf R_{3}}{\bf C_{1}} = A_{31}, \quad D_{13}(1) = 1, \quad D_{31} = 1,
\end{equation}\]
(34)#\[\begin{equation}
{\bf R_{4}}{\bf C_{1}} = A_{41}, \quad D_{41}(1) = 1, \quad D_{41} = 1.
\end{equation}\]
Para la columna dos sólo se necesitan calcular tres valores
(35)#\[\begin{equation}
{\bf R_{2}}{\bf C_{2}} = A_{22}, \quad D_{21}\left(\frac{1}{4} \right) + D_{22}\left(1 \right) = 3, \quad D_{22} = \frac{11}{4},
\end{equation}\]
(36)#\[\begin{equation}
{\bf R_{3}}{\bf C_{2}} = A_{32}, \quad D_{31}\left(\frac{1}{4} \right) + D_{32}\left(1 \right) = -1, \quad D_{32} = -\frac{5}{4},
\end{equation}\]
(37)#\[\begin{equation}
{\bf R_{4}}{\bf C_{2}} = A_{42}, \quad D_{41}\left(\frac{1}{4} \right) + D_{42}\left(1 \right) = 1, \quad D_{42} = \frac{3}{4}.
\end{equation}\]
En la columna tres sólo hay dos incógnitas
(38)#\[\begin{equation}
{\bf R_{3}}{\bf C_{3}} = A_{33}, \quad D_{31}\left(\frac{1}{4} \right) + D_{32}\left( -\frac{5}{11} \right) + D_{33}(1) = 2, \quad D_{33} = \frac{13}{11},
\end{equation}\]
(39)#\[\begin{equation}
{\bf R_{4}}{\bf C_{3}} = A_{43}, \quad D_{41}\left(\frac{1}{4} \right) + D_{42}\left(-\frac{5}{11} \right) + D_{43}(1) = 0, \quad D_{43} = \frac{1}{11}.
\end{equation}\]
Finalmente, para la columna cuatro, se tiene una incógnita en función de los valores de las columnas anteriores
(40)#\[\begin{equation}
{\bf R_{4}}{\bf C_{4}} = A_{44}, \quad D_{41}\left(\frac{1}{4} \right) + D_{42}\left(\frac{3}{11} \right) + D_{43}\left(\frac{1}{13}\right) + D_{44}(1) = 2, \quad D_{44} = \frac{20}{13}.
\end{equation}\]
Así, se tiene la factorización final \(\bf DC = A\)
(41)#\[\begin{equation}
\begin{bmatrix}
4 & 0 & 0 & 0 \\ \\
1 & \frac{11}{4} & 0 & 0 \\ \\
1 & -\frac{5}{4} & \frac{13}{11} & 0 \\ \\
1 & \frac{3}{4} & \frac{1}{11} & \frac{20}{13} \\
\end{bmatrix}
\begin{bmatrix}
1 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} \\ \\
0 & 1 & -\frac{5}{11} & \frac{3}{11} \\ \\
0 & 0 & 1 & \frac{1}{13} \\ \\
0 & 0 & 0 & 1 \\
\end{bmatrix}=
\begin{bmatrix}
4 & 1 & 1 & 1 \\ \\
1 & 3 & -1 & 1 \\ \\
1 & -1 & 2 & 0 \\ \\
1 & 1 & 0 & 2 \\
\end{bmatrix}.
\end{equation}\]
De modo que el sistema a resolver está dado como sigue
(42)#\[\begin{equation}
\begin{bmatrix}
4 & 0 & 0 & 0 \\ \\
1 & \frac{11}{4} & 0 & 0 \\ \\
1 & -\frac{5}{4} & \frac{13}{11} & 0 \\ \\
1 & \frac{3}{4} & \frac{1}{11} & \frac{20}{13} \\
\end{bmatrix}
\begin{bmatrix}
1 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} \\ \\
0 & 1 & -\frac{5}{11} & \frac{3}{11} \\ \\
0 & 0 & 1 & \frac{1}{13} \\ \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
x_{1} \\ \\ x_{2} \\ \\ x_{3} \\ \\ x_{4}
\end{bmatrix} =
\begin{bmatrix}
2 \\ \\
4 \\ \\
6 \\ \\
8
\end{bmatrix}.
\end{equation}\]
Utilizando un paso intermedio de la forma
(43)#\[\begin{equation}
\begin{bmatrix}
4 & 0 & 0 & 0 \\ \\
1 & \frac{11}{4} & 0 & 0 \\ \\
1 & -\frac{5}{4} & \frac{13}{11} & 0 \\ \\
1 & \frac{3}{4} & \frac{1}{11} & \frac{20}{13} \\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\ \\ y_{2} \\ \\ y_{3} \\ \\ y_{4}
\end{bmatrix} =
\begin{bmatrix}
2 \\ \\
4 \\ \\
6 \\ \\
8
\end{bmatrix}.
\end{equation}\]
Por sustitución hacia adelante se obtiene la solución
(44)#\[\begin{equation}
y_{1} = \frac{2}{4},
\end{equation}\]
(45)#\[\begin{equation}
y_{2} = \frac{14}{11},
\end{equation}\]
(46)#\[\begin{equation}
y_{3} = \frac{78}{13},
\end{equation}\]
(47)#\[\begin{equation}
y_{4} = \frac{78}{20}.
\end{equation}\]
De esta forma se llega a
(48)#\[\begin{equation}
\begin{bmatrix}
1 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} \\ \\
0 & 1 & -\frac{5}{11} & \frac{3}{11} \\ \\
0 & 0 & 1 & \frac{1}{13} \\ \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
x_{1} \\ \\ x_{2} \\ \\ x_{3} \\ \\ x_{4}
\end{bmatrix} =
\begin{bmatrix}
\frac{2}{4} \\ \\
\frac{14}{11} \\ \\
\frac{78}{13} \\ \\
\frac{39}{10}
\end{bmatrix}.
\end{equation}\]
Finalmente, por sustitución regresiva obtenemos la solución del sistema
(49)#\[\begin{equation}
x_{4} = \frac{39}{10},
\end{equation}\]
(50)#\[\begin{equation}
x_{3} = \frac{57}{10},
\end{equation}\]
(51)#\[\begin{equation}
x_{2} = \frac{14}{5},
\end{equation}\]
(52)#\[\begin{equation}
x_{1} = -\frac{13}{5}.
\end{equation}\]
# Matriz de coeficientes del sistema de ecuaciones lineales.
A = np.array([
[1, 1, 4, 2, 2],
[2, 5, 3, 6, 6],
[4, 3, 9, 1, 7],
[6, 5, 2, 7, 1],
[6, 4, 9, 0, 2]
], dtype=float)
# Vector de términos independientes (entradas) del sistema.
B = np.array([5, 2, 7, 6, 6], dtype=float)
# Número de incógnitas del sistema de ecuaciones (calculado como el rango de A).
N = np.linalg.matrix_rank(A)
# Proceso para obtener la matriz triangular superior (C)
C = A.copy() # Inicialmente C es una copia de la matriz A.
for k in range(N):
# Normalización del pivote en la fila k.
C[k, :] = C[k, :] / C[k, k]
for m in range(k + 1, N):
MT = -C[m, k] / C[k, k] # Multiplicadores
C[m, :] = C[m, :] + MT * C[k, :] # Modificación de la matriz C
# Proceso para obtener la matriz triangular inferior (D)
D = np.zeros((N, N)) # Se inicializa con una matriz de ceros.
# Se resuelve la primera columna de D (índice 0 en Python).
D[0:N, 0] = A[0:N, 0] / C[0, 0]
for k in range(1, N):
for m in range(k, N):
MT = 0
for l in range(k):
# Sumatoria de elementos conocidos (adaptada a indexación de 0)
MT = MT + D[m, l] * C[l, k]
# Se calcula el elemento (m, k) de la matriz D.
D[m, k] = (A[m, k] - MT) / C[k, k]
# Proceso de sustitución progresiva (Dy = B)
y = np.zeros(N)
y[0] = B[0] / D[0, 0] # El primer elemento del vector auxiliar y.
for k in range(1, N):
suma_y = 0
for m in range(k):
suma_y = suma_y - D[k, m] * y[m]
y[k] = (B[k] + suma_y) / D[k, k]
# Proceso de sustitución regresiva (Cx = y)
x = np.zeros(N)
x[N - 1] = y[N - 1] / C[N - 1, N - 1] # Se calcula la última variable.
for k in range(N - 2, -1, -1):
suma_x = 0
ind = N - 1 - k
for m in range(ind):
suma_x = suma_x - C[k, k + 1 + m] * x[k + 1 + m]
x[k] = y[k] + suma_x
# Mostrar resultados
print("Matriz D (Inferior):")
print(D)
print("\nMatriz C (Superior):")
print(C)
print("\nVector solución x:")
print(x)
Matriz D (Inferior):
[[ 1. 0. 0. 0. 0. ]
[ 2. 3. 0. 0. 0. ]
[ 4. -1. -8.66666667 0. 0. ]
[ 6. -1. -23.66666667 12.96153846 0. ]
[ 6. -2. -18.33333333 2.73076923 -5.97626113]]
Matriz C (Superior):
[[ 1. 1. 4. 2. 2. ]
[ 0. 1. -1.66666667 0.66666667 0.66666667]
[-0. -0. 1. 0.73076923 0.03846154]
[ 0. 0. 0. 1. -0.72700297]
[-0. -0. -0. -0. 1. ]]
Vector solución x:
[ 0.85004965 -1.85203575 0.93843098 1.19314796 -0.06901688]
Método de Cholesky#
La factorización de Choleski calcula y almacena los elementos de \(\bf L\) sólo si los elementos de \(\bf L\) y \(\bf U\) son iguales, entonces \(\bf U = L^{T}\). Entonces, para el primer pivote se tiene
(53)#\[\begin{equation}
C_{h}(1,1) = \sqrt{A(1,1)},
\end{equation}\]
y para la primera columna
(54)#\[\begin{equation}
C_{h}(2:n,1) = \frac{A(2:n,1)}{A(1,1)}.
\end{equation}\]
Los siguientes pivotes son
(55)#\[\begin{equation}
C_{h}(i,i) = \left[ A(i,i) - \sum_{k=1}^{i-1} \left\{ C_{h}(i,k) \right\}^{2} \right], \quad \text{para} \quad i=2,3,4,\dots,n.
\end{equation}\]
Las siguientes columnas se calculan con la siguiente expresión
(56)#\[\begin{equation}
C_{h}(i,j) = \frac{1}{C_{h}(i,i)} \left[ A(i,j) - \sum_{k=1}^{j-1} C_{h}(i,k)\, C_{h}(j,k) \right].
\end{equation}\]
Nota
El método de Cholesky calcula los factores de la forma \(\bf A = CC^{T}\) por lo que sólo se calcula la matriz triangular inferior.
Definición 8
Una matriz \(\bf A\) se dice matriz definida positiva si para cada vector no nulo \(\bf X\) se tiene que \(\bf X^{T}AX>0\).
# Matriz de coeficientes del sistema de ecuaciones lineales.
# (Nota: Para Cholesky, la matriz debe ser simétrica y definida positiva).
A = np.array([
[9, 1, 3, 2, 1],
[1, 7, 1, 3, 4],
[3, 1, 9, 1, 2],
[2, 3, 1, 7, 1],
[1, 4, 2, 1, 8]
], dtype=float)
# Vector de términos independientes (entradas) del sistema.
B = np.array([2, 1, 4, 1, 3], dtype=float)
# Número de incógnitas del sistema de ecuaciones (calculado como el rango de A).
N = np.linalg.matrix_rank(A)
# Proceso para obtener la matriz triangular inferior (CInf)
CInf = np.zeros((N, N))
# Se calcula el primer pivote (índice 0 en Python).
CInf[0, 0] = np.sqrt(A[0, 0])
# Se calcula la primera columna de CInf.
CInf[1:N, 0] = A[1:N, 0] / CInf[0, 0]
# Se calculan los siguientes pivotes y columnas.
for k in range(1, N):
# Cálculo del pivote en la posición diagonal (k, k).
CInf[k, k] = np.sqrt(A[k, k] - np.sum(CInf[k, 0:k] ** 2))
for m in range(k + 1, N):
# Sumatoria de los elementos conocidos.
suma = np.sum(CInf[m, 0:k] * CInf[k, 0:k])
# Cálculo del elemento (m, k).
CInf[m, k] = (1 / CInf[k, k]) * (A[m, k] - suma)
# Matriz triangular superior (CSup)
# En Cholesky, la matriz superior es la transpuesta de la inferior.
CSup = CInf.T
# Proceso de sustitución progresiva (CInf * y = B)
y = np.zeros(N)
y[0] = B[0] / CInf[0, 0] # Primer elemento del vector auxiliar y.
for k in range(1, N):
suma_y = 0
for m in range(k):
suma_y = suma_y - CInf[k, k - 1 - m] * y[k - 1 - m]
y[k] = (B[k] + suma_y) / CInf[k, k]
# Proceso de sustitución regresiva (CSup * x = y)
x = np.zeros(N)
x[N - 1] = y[N - 1] / CSup[N - 1, N - 1] # Última variable.
for k in range(N - 2, -1, -1):
suma_x = 0
ind = N - 1 - k
for m in range(ind):
suma_x = suma_x - CSup[k, k + 1 + m] * x[k + 1 + m]
x[k] = (y[k] + suma_x) / CSup[k, k]
# Mostrar resultados
print("Matriz Triangular Inferior (CInf):")
print(CInf)
print("\nMatriz Triangular Superior (CSup):")
print(CSup)
print("\nVector solución x:")
print(x)
Matriz Triangular Inferior (CInf):
[[ 3. 0. 0. 0. 0. ]
[ 0.33333333 2.62466929 0. 0. 0. ]
[ 1. 0.25400025 2.81699909 0. 0. ]
[ 0.66666667 1.05833439 0.02290243 2.33129993 0. ]
[ 0.33333333 1.48166815 0.45804863 -0.34350493 2.31640761]]
Matriz Triangular Superior (CSup):
[[ 3. 0.33333333 1. 0.66666667 0.33333333]
[ 0. 2.62466929 0.25400025 1.05833439 1.48166815]
[ 0. 0. 2.81699909 0.02290243 0.45804863]
[ 0. 0. 0. 2.33129993 -0.34350493]
[ 0. 0. 0. 0. 2.31640761]]
Vector solución x:
[ 0.0628659 -0.15137998 0.35405631 0.09018679 0.34304433]
Factorización QR#
Este método consiste en realizar transformaciones ortogonales dentro del proceso para proporcionar condiciones de estabilidad. Es decir, las matrices se descomponen en un producto \(\bf Q_{r}R_{r}\), donde \(\bf R_{r}\) es triangular superior
(57)#\[\begin{equation}
\bf A = A_{1} = Q_{1}R_{1},
\end{equation}\]
(58)#\[\begin{equation}
\bf A_{r} = R_{r-1}Q_{r-1} = Q_{r}R_{r}, \quad r = 2,3,4,\dots .
\end{equation}\]
Se puede observar que las matrices son similares debido a que
(59)#\[\begin{equation}
\bf A_{r} = R_{r-1}Q_{r-1} = Q_{r}R_{r} = Q_{r-1}^{-1}A_{r-1}Q_{r-1}.
\end{equation}\]
La descomposición básica \(\bf A = Q_{r}R_{r}\) se puede realizar para cualquier matriz, lo que no pasa en el caso de la descomposición \(\bf LU\).
# Matriz de coeficientes del sistema de ecuaciones lineales.
Aqr = np.array([
[9, 5, 8, 9, 1],
[3, 6, 8, 2, 4],
[5, 7, 3, 6, 2],
[3, 1, 4, 5, 2],
[1, 5, 6, 3, 4]
], dtype=float)
# Vector de términos independientes (entradas) del sistema.
B = np.array([2, 3, 1, 3, 1], dtype=float)
# Número de incógnitas del sistema de ecuaciones (calculado como el rango de Aqr).
N = np.linalg.matrix_rank(Aqr)
A1 = Aqr.copy() # Asignación de la matriz original a A1.
Qi = np.eye(N) # Se inicializa Qi como matriz unitaria.
Qd = np.eye(N) # Se inicializa Qd como matriz unitaria.
# Ciclo iterativo para calcular Qd y Qi para la triangulación (Algoritmo QR)
# En Python, 'range(1, 31, 2)' emula el incremento de 1 a 30 con paso de 2 (1, 3, 5, ..., 29).
for kr in range(1, 31, 2):
Q1, R1 = np.linalg.qr(A1)
A2 = R1 @ Q1
Q2, R2 = np.linalg.qr(A2)
A1 = R2 @ Q2
# Actualización de Qd y Qi.
Qd = Qd @ Q1 @ Q2
Qi = np.linalg.inv(Q2) @ np.linalg.inv(Q1) @ Qi
# Triangular superior a la que se llega utilizando el método QR.
As = Qi @ Aqr @ Qd
# Solución del sistema de ecuaciones.
x = Qd @ np.linalg.inv(As) @ Qi @ B
# Mostrar resultados
print("Matriz Triangular Superior (As):")
print(As)
print("\nVector solución x:")
print(x)
Matriz Triangular Superior (As):
[[ 2.25104624e+01 7.55156334e-01 2.37447855e+00 3.59185254e+00
-6.60024070e+00]
[-3.55357804e-15 6.45865053e+00 1.94422204e+00 -1.91070019e+00
9.68106235e-01]
[-2.11665787e-15 4.22387329e-06 -4.08553885e+00 -4.52429044e-01
-1.31832197e+00]
[ 1.12047556e-15 -3.87898657e-12 1.35775224e-05 2.73655361e+00
-1.27841962e-01]
[ 1.01581394e-15 -4.78723422e-16 -3.16477837e-16 4.77579321e-16
-6.20127696e-01]]
Vector solución x:
[ 1.70634921 -0.88690476 -0.61706349 -0.70833333 2.38888889]