Métodos iterativos para resolver sistemas de ecuaciones lineales¶
Para poner en práctica la teoría que hemos visto sobre sistemas de ecuaciones, usaremos dos funciones que crean una matriz $A$ de dimensiones $N\times N$ y un vector $\mathbf{b}$ de longitud $N$, para estudiar el sistema $A\mathbf{x}=\mathbf{b}$.
$$ A = \left(\begin{array}{rrrrr} 2 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \end{array}\right)\quad \mathbf{b}= \left(\begin{array}{r} 1 \\0 \\0 \\0 \\0 \end{array}\right) $$function A=An(N)
A = zeros(N);
for i=1:N
A(i,i)=2;
end
for i=1:N-1
A(i,i+1)=-1;
A(i+1,i)=-1;
end
end
function b=bn(N)
b = zeros(N,1);
b(1) = 1;
end
An(5)
An(5)\bn(5)
Normas vectoriales¶
Hemos estudiado dos normas vectoriales muy importantes: $$ \|\mathbf{x}\|_2=\left\{ \sum_{i=1}^{n}x_i^2 \right\}^{1/2} \qquad \|\mathbf{x}\|_{\infty}=\max_{1\leq i\leq n}|x_i| $$
N = 5;
A = An(5);
b = bn(5);
function n = norma2(x)
%Recuerda que x.^2 es el array que tiene cada elemento al cuadrado,
% es igual a x.*x, pero es distinto de x*x', que es el producto escalar
n = sqrt(sum(x.^2))
end
matlab ya tiene una función para calcular la norma 2.
x = randn(1,4)
norma2(x)
norm(x)
sqrt(x*x')
Ejercicio¶
- Define una función que calcule la norma infinito
- Compárala con la función "norma infinito" de
matlab:norm(b, inf)
x = randn(1,4)
norm(x, inf)
Normas matriciales¶
Hemos estudiado dos normas matriciales:
- La norma $\|A\|_2$ inducida por la norma vectorial $\|\mathbf{x}\|_2$.
- La norma $\|A\|_\infty$ inducida por la norma vectorial $\|\mathbf{x}\|_\infty$. La definición de norma inducida es: $$ \|A\|_2=\max_{\mathbf{z}\neq\mathbf{0}}\frac{\|A\mathbf{z}\|_2}{\|\mathbf{z}\|_2} $$ $$ \|A\|_\infty=\max_{\mathbf{z}\neq\mathbf{0}}\frac{\|A\mathbf{z}\|_\infty}{\|\mathbf{z}\|_\infty} $$ pero la forma de calcularlas es muy diferente, porque estas definiciones no son muy constructivas.
La norma infinito es sencilla de calcular $$ \|A\|_{\infty}=\max_{1\leq i\leq n}\sum_{j=1}^{n}|a_{ij}|. $$ La norma $2$ es más complicada:
- El radio espectral de una matriz $A$ es el máximo de los valores absolutos de los autovalores de $A$. $$ \rho(A)=\max|\lambda| $$ Si $\lambda$ es complejo, entonces $|\lambda|=\sqrt{\alpha^2+\beta^2}$.
- $\|A\|_2=\left[\rho(A^tA)\right]^{1/2}$ (si $A$ es simétrica: $\|A\|_2=\rho(A)$)
function r=radio_espectral(A)
r = max(abs(eigs(A)));
end
function n = matrix_norm_inf(A)
[N, M] = size(A);
n = -inf;
for j=[1:N] %máximo sobre todas las filas
sumafila = sum(abs(A(j,:))); %suma de los elementos de la fila j
n = max(n, sumafila);
end
end
De nuevo, matlab ya tiene implementadas estas normas: comparamos el resultado con nuestras funciones.
% Una matriz con entradas aleatorias
M = randn(5,5);
matrix_norm_inf(M)
norm(M, inf)
sqrt(radio_espectral(M*M'))
norm(M, 2)
Método de Jacobi¶
Para escribir el método de Jacobi de la manera anterior es necesario partir la matriz $A$ en tres.
% el primer diag(A) extrae los elementos de la diagonal
% el segundo diag(diag(A)) crea una matriz cuadrada
% que tiene ceros fuera de la diagonal, y en su diagonal
% tiene diag(A)
D = diag(diag(A))
%tril(A) extra la parte de A que está debajo de la diagonal
% incluyendo la diagonal
% np.tril(A,k=-1) extra la parte de A que está debajo de la diagonal
% sin incluir la diagonal
L = -tril(A,-1)
U = -triu(A,1)
D - L - U
Cada iteración del método de Jacobi consiste en resolver un sistema de ecuaciones donde la matriz de coeficientes es diagonal: $$ D\mathbf{x}^{(k)}=(L+U)\mathbf{x}^{(k-1)}+\mathbf{b} $$ Multiplicando por la inversa de $D$: $$ \mathbf{x}^{(k)}=D^{-1}(L+U)\mathbf{x}^{(k-1)}+D^{-1}\mathbf{b}\qquad k=1,2,\ldots. $$ Con la notación $$ T_J=D^{-1}(L+U)\quad \mathbf{c}_J=D^{-1}\mathbf{b} $$ se puede escribir de forma equivalente: $$ \mathbf{x}^{(k)}=T_J\mathbf{x}^{(k-1)}+\mathbf{c}_J. $$
%También vale Dinv = inv(D), pero es menos eficiente
Dinv = diag(1./diag(A))
T = Dinv*(L+U)
c = Dinv*b
Criterio de parada¶
Finalmente, hay que iterar $$ \mathbf{x}^{(k)}=T_J\mathbf{x}^{(k-1)}+\mathbf{c}_J. $$ hasta que se satisfaga el criterio de parada. Dos criterios de parada habituales:
- [tolerancia relativa en la x] Detener cuando $$ \frac{\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\|}{\|\mathbf{x}^{(k)}\|} < \text{rtol} $$
- [tolerancia absoluta en la y] Detener cuando $$ \|A\cdot\mathbf{x}^{(k)}-\mathbf{b}\| < \text{ytol} $$
rtol = 1e-6;
%nuestro código del método de Jacobi
% con criterio de parada de tolerancia relativa en la x
xk = ones(N,1);
xk1 = T*xk + c;
k = 0;
while norm(xk1-xk, inf )>rtol*norm(xk, inf)
xk = xk1;
xk1 = T*xk + c;
k = k + 1;
end
xk1
k
% Comprobamos el error usando el otro criterio: Ax-b debe ser pequeño
norm(A*xk1 - b, inf )
Ejercicio¶
"Empaqueta" la función que resuelve un sistema de ecuaciones lineales de forma aproximada por el método de Jacobi como una función:
- inputs: matriz de coeficientes, vector de términos independientes, tolerancia
- outputs: solución aproximada, número de iteraciones
Ejercicio¶
Modifica el código anterior para cambiar el criterio de parada: [tolerancia absoluta en la y]: detener cuando: $$ \|A\cdot\mathbf{x}^{(k)}-\mathbf{b}\| < \text{ytol} $$
Ejercicio¶
"Empaqueta" la función que resuelve un sistema de ecuaciones lineales de forma aproximada por el método de Gauss-Seidel como una función:
- inputs: matriz de coeficientes, vector de términos independientes, tolerancia
- outputs: solución aproximada Usa el criterio de parada que prefieras.
Cada iteración del método de Gauss-Seidel consiste en resolver un sistema de ecuaciones donde la matriz de coeficientes es triangular inferior: $$ (D-L)\mathbf{x}^{(k)}=U\mathbf{x}^{(k-1)}+\mathbf{b} $$
Ejercicio¶
Fija una matriz, por ejemplo la matriz An(5) que hemos usado en los ejemplos.
- Genera un vector de tolerancias, y otro de número de iteraciones que necesita el método de Jacobi para terminar su ejecución con esa tolerancia.
- Repite el ejercicio anterior para el método de Gauss-Seidel.
- Haz una gráfica con tolerancias en las "x" (en escala logarítmica) y número de iteraciones en las "y" (en escala normal). Dibuja el número de iteraciones para alcanzar una tolerancia dada para cada uno de los métodos.
help semilogx