Práctica 7. Sistemas de ecuaciones lineales III (Métodos iterativos)
En esta práctica implementaremos los métodos iterativos de Jacobi y de Gauss-Seidel y veremos un par de criterios de convergencia que se pueden aplicar a estos métodos.
Contents
Matrices D, L y U para una matriz A.
En clase se vió como expresar el método de Jacobi y Gauss-Seidel a partir de unas matrices D, L y U de forma que A=D-L-U donde D es una matriz diagonal cuya diagonal coincide con la de la matriz A, L es una matriz triangular inferior que bajo su diagonal tiene a los opuestos de los coeficientes de A que están bajo la diagonal de A, y U es una matriz triangular superior que sobre su diagonal tiene a los opuestos de los coeficientes de A que están sobre la diagonal de A.
Veamos como construir estas matrices D, L y U usando comandos de Matlab.
A=[
[1 .2 .3 .4]; [.5 6 .7 .8] ; [.9 .10 11 .12]; [.13 .14 .15 16]
];
Si A es una matriz, diag(A) es un vector columna con los elementos de la diagonal de A.
Si v es un vector con n componentes, diag(v) nos da una matriz cuadrada de orden n con las componentes de v en la diagonal. Por tanto diag(diag(A)) es una matriz diagonal cuyos elementos de la diagonal coinciden con los de la matriz A.
D=diag(diag(A))
D =
1 0 0 0
0 6 0 0
0 0 11 0
0 0 0 16
El comando tril nos da la parte inferior de una matriz. Veamos unos ejemplos.
tril(A,0) tril(A,1) tril(A,-1)
ans =
1.0000 0 0 0
0.5000 6.0000 0 0
0.9000 0.1000 11.0000 0
0.1300 0.1400 0.1500 16.0000
ans =
1.0000 0.2000 0 0
0.5000 6.0000 0.7000 0
0.9000 0.1000 11.0000 0.1200
0.1300 0.1400 0.1500 16.0000
ans =
0 0 0 0
0.5000 0 0 0
0.9000 0.1000 0 0
0.1300 0.1400 0.1500 0
La instrucción tril(A,0) se puede escribir como tril(A). Para definir la matriz L usamos la siguiente instrucción:
L=-tril(A,-1)
L =
0 0 0 0
-0.5000 0 0 0
-0.9000 -0.1000 0 0
-0.1300 -0.1400 -0.1500 0
El comando triu nos da la parte superior de una matriz.
triu(A,0) triu(A,1) triu(A,-1)
ans =
1.0000 0.2000 0.3000 0.4000
0 6.0000 0.7000 0.8000
0 0 11.0000 0.1200
0 0 0 16.0000
ans =
0 0.2000 0.3000 0.4000
0 0 0.7000 0.8000
0 0 0 0.1200
0 0 0 0
ans =
1.0000 0.2000 0.3000 0.4000
0.5000 6.0000 0.7000 0.8000
0 0.1000 11.0000 0.1200
0 0 0.1500 16.0000
De igual forma triu(A,0) puede escribirse como triu(A). Para definir la matriz U usamos la siguiente instrucción:
U=-triu(A,1)
U =
0 -0.2000 -0.3000 -0.4000
0 0 -0.7000 -0.8000
0 0 0 -0.1200
0 0 0 0
Con estas definiciones la matriz A se puede expresar como D-L-U. Comprobémoslo con la siguente instrucción, que nos devuelve una matriz llena de ceros indicando que las dos matrices son iguales elemento a elemento.
A-(D-L-U)
ans =
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
Método de Jacobi
Veamos cómo emplear el método de Jacobi para resolver un sistema A*x=b. Ya hemos definido la matriz del sistema A, y definimos ahora el vector término independiente b como un vector columna.
b=[-0.1; -12.6; 33.22; -63.7];
Este método iterativo comienza con una aproximación inicial que podemos llamar x0, y dada una iteración xi calcula la siguiente resolviendo el sistema de ecuaciones lineales con matriz de coeficientes D y vector de términos independientes (L+U)*xi+b. Para resolver dicho sistema con MatLab utilizaremos, por comodidad, la división izquierda.
Vamos a calcular las primeras iteraciones del método de Jacobi para resolver A*x=b tomando como aproximación inicial el vector x0 con cuatro filas, una columna y todas sus componentes iguales a 1. (Se puede comprobar que la solución exacta del sistema es [1; -2; 3; -4] ).
x0=ones(4,1) x1=D\((L+U)*x0+b) x2=D\((L+U)*x1+b) x3=D\((L+U)*x2+b) x4=D\((L+U)*x3+b)
x0 =
1
1
1
1
x1 =
-1.0000
-2.4333
2.9182
-4.0075
x2 =
1.1142
-1.8228
3.1677
-3.9792
x3 =
0.9059
-2.0319
2.9888
-4.0041
x4 =
1.0113
-1.9903
3.0080
-3.9989
Para mostrar los datos en una tabla debemos trasponer los vectores.
disp(' k x(1) x(2) x(3) x(4) ') disp(' -------------------------------------------------') disp([1 transpose(x1)]) disp([2 transpose(x2)]) disp([3 transpose(x3)]) disp([4 transpose(x4)])
k x(1) x(2) x(3) x(4)
-------------------------------------------------
1.0000 -1.0000 -2.4333 2.9182 -4.0075
2.0000 1.1142 -1.8228 3.1677 -3.9792
3.0000 0.9059 -2.0319 2.9888 -4.0041
4.0000 1.0113 -1.9903 3.0080 -3.9989
Ej 1. Halle la quinta iteración en el ejemplo descrito. (Solución: x5=[ 0.9952;-2.0020;2.9990;-4.0003])
Ej 2. Considere el sistema A*x = b donde A=[[3 2 -1 1]; [-1 2 2 0]; [0 2 3 1]; [1 -1 0 2]] y b=[11;2;1;-5] . Halle la vigésima iteración del método de Jacobi comenzando con la aproximación inicial x0=[2; 2; 2; 2]. (Indicación: llame x a todas las iteraciones, y dentro de un bucle for calcule cada iteración x resolviendo, usando la división izquierda, el sistema de ecuaciones lineales con matriz de coeficientes D y vector de términos independientes (L+U)*x+b ). (Solución: x20=[2.0024;2.9997;-0.9982;-2.0022]. Se puede comprobar que la solución exacta del sistema es x=[2;3;-1;-2] )
Ej 3. Considere el sistema A*x = b donde A=[[2 3 4 1]; [-1 2 3 0]; [0 2 3 1]; [1 -1 0 1]] y b=[7;1;1;-3] . Halle la vigésima iteración del método de Jacobi comenzando con la aproximación inicial x0=[2; 2; 2; 2]. (Solución: x20=[20.3593;-35.8838;34.9481;-33.9740]. Se puede comprobar que la solución exacta del sistema es x=[2;3;-1;-2], aunque el método de Jacobi no converge en este caso)
Método de Gauss-Seidel
El método de Gauss-Seidel se puede implementar de forma similar teniendo en cuenta que, dada una iteración xi, dicho método calcula la siguiente resolviendo el sistema de ecuaciones lineales con matriz de coeficientes D-L y vector de términos independientes U*xi+b.
Ej 3b. Modifique el programa realizado en el último ejercicio de forma que se halle la vigésima iteración del método de Gauss-Seidel. (Solución: x20=[88.7700;-686.8669;652.5354;-778.6369]. En este caso el método de Gauss-Seidel tampoco converge)
Convergencia. Radio Espectral.
Dada una aproximación inicial x(0) y un método iterativo de la forma M*x(k)=N*x(k-1)+b con M regular, este método es convergente si y solo si el radio espectral de la matriz B=inv(M)*N es menor que 1.
El radio espectral de una matriz es el máximo de las normas de los valores propios de la matriz.
Dada una matriz, podemos calcular sus valores propios (eigenvalues) con la función de Matlab eig(matriz) y si aplicamos la función abs a este vector obtendremos un vector cuyas componentes son las normas de los valores propios de la matriz; por último, la función max aplicada sobre el vector de las normas proporciona el máximo de las mismas, es decir, el radio espectral de la matriz.
BJ=inv(D)*(L+U) matriz=BJ; r= max(abs(eig(matriz)))
BJ =
0 -0.2000 -0.3000 -0.4000
-0.0833 0 -0.1167 -0.1333
-0.0818 -0.0091 0 -0.0109
-0.0081 -0.0088 -0.0094 0
r =
0.2441
Ej 4. Halle el radio espectral de la matriz B=[[1 0 1]; [0 1 3]; [1 2 3]]. Si M es una matriz regular y B=inv(M)*N ¿es convergente el método iterativo de la forma M*x(k)=N*x(k-1)+b ? ¿Por qué? (Solución: el radio espectral de B es r=4.8284 y, por tanto, el método no es convergente)
Ej 5. Escriba un programa en el que para una matriz A dada se muestren en pantalla los radios espectrales de las matrices B de Jacobi y de Gauss-Seidel así como un mensaje que informe de si el método de Jacobi para resolver un sistema con matriz de coeficientes A converge o no converge, y otro que informe de si el método de Gauss-Seidel para resolver un tal sistema converge o no converge. Pruebe el programa con la matriz A=[[2 1 1]; [2 1 0]; [1 2 3]]. (Sol: el radio espectral de la matriz B del método de Jacobi es rJ=1.2965, el radio espectral de la matriz B del método de Gauss-Seidel es rGS=0.5000, por tanto el método de Jacobi no converge y el método de Gauss-Seidel si converge.)
Convergencia. Matriz estrictamente diagonal dominante.
En algunos casos puede asegurarse la convergencia de los métodos de Jacobi y Gauss-Seidel sin necesidad de hallar el radio espectral de una matriz.
Si en un sistema A*x=b se tiene que para cada fila i de la matriz A el coeficiente A(i,i) en valor absoluto es mayor que la suma de los valores absolutos del resto de los coeficientes de la fila (es decir, A es estrictamente diagonal dominante), entonces los métodos de Jacobi y de Gauss-Seidel convergen.
La función diagdom.m permite saber si una matriz dada es o no estrictamente diagonal dominante. En el código de esta función, la matriz A es el argumento de la función, y el valor que devuelve la función es la variable diagdom.
Ej 6. Estudie la función diagdom.m y explique la condición que aparece en el bucle while.
Ej 7. Construya una matriz de orden 3 por 3 que sea estrictamente diagonal dominante, y otra del mismo orden que no lo sea. Escriba un programa que muestre en pantalla esas dos matrices e informe de si cada una de ellas es estrictamente diagonal dominante o no (dicho programa deberá utilizar la función diagdom.m).
% Matemáticas II. Grado en Ingeniería Química. (A. Palomares)
disp(date)
05-Apr-2017