PDF de programación - Programación: Gauss-Seidel en forma matricial

Programaci´on: el m´etodo de Gauss–Seidel

en forma matricial

Objetivos. Programar una funci´on que resuelva sistemas de ecuaciones lineales con el
m´etodo iterativo de Gauss–Seidel en la forma matricial.

Requisitos. Operaciones con vectores (operaciones lineales y el producto punto), matri-
ces triangulares, soluci´on de sistemas con matrices triangulares, ciclo while, la idea del
m´etodo de Gauss–Seidel.

1. Soluci´on de sistemas triangulares inferiores. Escribir una funci´on que resuelva
sistemas de ecuaciones lineales de la forma T x = b, donde T es una matriz cuadrada
triangular inferior con entradas diagonales no nulas.

Entrada: una matriz T , un vector b.

Propiedades de la entrada: suponer (y no verificar) que T es una matriz
cuadrada triangular inferior con entradas diagonales no nulas y que la longitud de
b coincide con el orden de T .

Salida: un vector x tal que T x = b.

Idea: aplicar la sustituci´on hacia adelante. Se recomienda usar operaciones de nivel 1 (el
producto punto de algunos fragmentos de columnas o las operaciones lineales con algunos
fragmentos de columnas).

function x = solvelt(T, b),

n = length(b); x = zeros(n, 1);
...

end

2. Descomposici´on de una matriz en la suma de una matriz triangular infe-
rior y una matriz estrictamente triangular superior. Cada matriz cuadrada A se
puede escribir de manera ´unica como A = T + U , donde T es triangular inferior y U es
estrictamente triangular superior. Por ejemplo, para n = 3,


(cid:125)

 0 A1,2 A1,3
(cid:123)(cid:122)
(cid:124)

0 A2,3
0

0
0

0

U


(cid:125)

.

 A1,1 A1,2 A1,3
(cid:124)

A2,1 A2,2 A2,3
A3,1 A3,2 A3,3

(cid:123)(cid:122)

A


(cid:125)

 A1,1
(cid:124)

=

A2,1 A2,2
A3,1 A3,2 A3,3

+

0
0

0

(cid:123)(cid:122)

T

En el lenguaje de MATLAB las matrices T y U se pueden obtener con los siguientes
comandos:

T = tril(A);
U = triu(A, 1);

Programaci´on: Gauss–Seidel en forma matricial, p´agina 1 de 4

M´etodo de Gauss–Seidel en forma matricial 1,
con dos matrices triangulares

3. Usando la notaci´on anterior T y U , escribimos el sistema original Ax = b en la forma

y luego en la forma

T x = b − U x

x = T −1(b − U x).

Empezando con alg´un vector x(0), construimos la sucesi´on de vectores (x(s))∞
la f´ormula recursiva

s=0 mediante

x(s+1) = T −1(b − U x(s)).

(1)

En el programa se recomienda trabajar solamente con el vector actual y el vector previo,
esto es, en cada paso del m´etodo de Gauss–Seidel guardar una copia del vector actual x
en la variable xprev, luego calcular el vector actual por la f´ormula (1):

xprev = x;
x = solvelt(T, b - U * xprev);

4. Esquema del algoritmo.

function [x, s] = GaussSeidelMatrixForm1(A, b, tol, smax),

T = ???; U = ???;
n = length(b); x = zeros(n, 1);
er = tol + 1; s = ???;
while (er ??? tol) && (s ??? smax),

xprev = x;
...
d = norm(x - xprev); s = s + 1;

end

end

Programaci´on: Gauss–Seidel en forma matricial, p´agina 2 de 4

M´etodo de Gauss–Seidel en forma matricial 2,
con una matriz triangular y con el vector residuo

5. M´etodo de Gauss escrito en t´erminos del vector residuo. Como antes, denota-
mos por T a la parte triangular inferior de la matriz A, incluso la diagonal. En cada paso
s denotamos por r(s) al vector residuo:

r(s) := b − Ax(s).

Notamos que el sistema original de ecuaciones lineales Ax = b se puede escribir en las
siguientes formas equivalentes:

0n = b − Ax

⇐⇒

0n = T −1(b − Ax)

⇐⇒

x = x + T −1(b − Ax).

La f´ormula recursiva ser´ıa

x(s+1) = x(s) + T −1r(s).

Por supuesto, en vez de calcular la matriz inversa T −1, hay que calcular T −1r(s) como el
vector soluci´on de un sistema de ecuaciones lineales con la matriz triangular inferior T .

6. Cambio del residuo al cambiar la aproximaci´on de la soluci´on. Si

x(s+1) = x(s) + p(s),

entonces

r(s+1) = b − Ax(s+1) = b − A(x(s) + p(s)) = b − Ax(s) − Ap(s) = r(s) − Ap(s).

Denotaremos el vector Ap(s) por q(s).

7. Esquema del algoritmo. Ahora es c´omodo usar el vector r(s) en la condici´on del
ciclo while.

function [x, s] = GaussSeidelMatrixForm2(A, b, tol, smax),

T = ???;
n = length(b); x = zeros(n, 1); r = b;
s = ???;
while (s ??? smax) && (norm(r) ??? tol),

p = solvelt(???, r);
q = A * p;
x = ???;
r = r - q;
s = ???;

end

end

Programaci´on: Gauss–Seidel en forma matricial, p´agina 3 de 4

Pruebas

8. Prueba con matrices peque˜nas. Para garantizar la convergencia generamos una
matriz estrictamente diagonal dominante. En vez de generar una matriz pseudoaleatoria
puede construir el ejemplo a mano.

function [] = testGaussSeidel1(),

n = 4;
A = 2 * rand(n) - ones(n) + n * eye(n);
b = 2 * rand(n, 1) - ones(n, 1);
[x, s] = GaussSeidelMatrixForm1(A, b, 1.0E-5, 100);
display(s); display(norm(A * x - b));

end

Tambi´en hacer pruebas del otro algoritmo, GaussSeidelMatrixForm2.

9. Pruebas con matrices grandes. Elegir n de tal manera que el tiempo de ejecuci´on
sea de 1 a 100 segundos, medir el tiempo de ejecuci´on con el comando comando cputime.

function [] = testGaussSeidel2(),

...

end

Programaci´on: Gauss–Seidel en forma matricial, p´agina 4 de 4