PDF de programación - Programación: método de Jacobi para resolver sistemas de ecuaciones lineales en forma matricial

Programaci´on: forma matricial

del m´etodo de Jacobi

para resolver sistemas de ecuaciones lineales

Objetivos. Programar una funci´on que resuelva sistemas de ecuaciones lineales con el
m´etodo iterativo de Jacobi 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 Jacobi.

1. Multiplicaci´on y divisi´on de vectores componente a componente. Dados dos
vectores a, b ∈ Rn, denotemos por a (cid:12) b al vector de longitud n cuya j-´esima componente
es

(a (cid:12) b)j = ajbj.

Si adem´as bj (cid:54)= 0 para cada j ∈ {1, . . . , n}, entonces denotemos por a (cid:11) b al vector de
longitud n cuya j-´esima componente es

(a (cid:11) b)j =

aj
bj

.

En el lenguaje de MATLAB (y en GNU Octave, Scilab, FreeMat) la multiplicaci´on y
divisi´on por componentes se denotan por .* y ./:

a = [6; -5; 1]
b = [-2; 2; 3]
a .* b
a ./ b

2. Multiplicaci´on de una matriz diagonal por un vector. Se d = [dj]n
j=1 un vector.
Se denota por diag(d) la matriz diagonal con entradas diagonales d1, . . . , dn. Por ejemplo,
si n = 3, entonces

0
0

 d1
 v1

v2
v3

0
d2
0

0
0
d3

 .
 d1v1

d2v2
d3v3

 =

diag(d) =

 d1

0
d2
0

0
0
d3

 = d (cid:12) v.

Notamos que

diag(d)v =

0
0
En general, si d, v ∈ Rn, entonces

diag(d)v = d (cid:12) v.

Programaci´on: m´etodo de Jacobi en forma matricial, p´agina 1 de 4

3. Multiplicaci´on de una matriz por un vector, prueba num´erica.

d = [5; -1; 2; 3]
v = [-3; 2; 3; 5]
diag(d)
diag(d) * v
d .* v
4. Soluci´on de sistemas de ecuaciones lineales con matrices diagonales. Si d, b ∈
Rn y dj (cid:54)= 0 para cada j ∈ {1, . . . , n}, entonces el sistema de ecuaciones lineales

diag(d)x = b

tiene una ´unica soluci´on la cual se puede calcular como

x = b (cid:11) d.

Prueba num´erica:

d = [-2; 3; 2];
b = [5; -1; 0];
x = b ./ d;
norm(diag(d) * x - b)

5. Sacar el vector de las entradas diagonales de una matriz. Dada una matriz A,
denotemos por d al vector de sus entradas diagonales y consideremos la matriz diagonal
D con entradas diagonales A1,1, . . . , An,n:

(cid:3)n

j=1,

d =(cid:2)Aj,j
 ,

 A1,1

A2,2
A3,3

d =

D = diag(d).

 ,

D =

 A1,1

 .

0
0

0

0 A2,2
0

0 A3,3

Por ejemplo, para n = 3,

 A1,1 A1,2 A1,3

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

A =

En el lenguaje de MATLAB el vector d se puede obtener con el comando diag:

A = rand(3)
d = diag(A)
D = diag(d)

Programaci´on: m´etodo de Jacobi en forma matricial, p´agina 2 de 4

Idea del m´etodo de Jacobi en forma matricial
6. Sea A ∈ Mn(R) tal que Aj,j (cid:54)= 0 para cada j ∈ {1, . . . , n}. Denotemos por d al vector
de las entradas diagonales de A y por D a la matriz diagonal generada por el vector d:

d =(cid:2)Aj,j

(cid:3)n

j=1,

D = diag(d).

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

Ax = b

⇐⇒ b − Ax = 0
⇐⇒ Dx = Dx + (b − Ax)
⇐⇒ x = x + D−1(b − Ax)
⇐⇒ x = x + (b − Ax) (cid:11) d.

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

s=0 mediante

x(s+1) = x(s) + (b − Ax(s)) (cid:11) d.

(1)
7. Residuo y su modificaci´on. Es c´omodo guardar el residuo b − Ax en una variable
r. Notamos que si el vector x se modifica mediante la f´ormula

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

entonces el vector residuo correspondiente se cambia de la siguiente manera:

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

En el algoritmo podemos usar los siguientes comandos:

p = r ./ d;
q = A * p;
x = x + p;
r = r - q;

8. Condici´on de terminaci´on. El ciclo while debe terminarse cuando el n´umero s de
las iteraciones realizadas es mayor o igual al n´umero m´aximo de iteraciones permitidas
smax o cuando la norma del residuo es menor que la “tolerancia” dada. En otras palabras,
la condici´on de terminaci´on es

(s ≥ smax)

∨

((cid:107)r(cid:107) < tol).

Escriba la condici´on de continuaci´on.

Programaci´on: m´etodo de Jacobi en forma matricial, p´agina 3 de 4

9. Esquema del algoritmo.

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

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

p = ???;
q = ???;
x = ???;
r = ???;
s = s + 1;

end

end

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

function smallTestJacobiSolve(),

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

end

11. 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 cputime.

function largeTestJacobiSolve(),

...

end

Programaci´on: m´etodo de Jacobi en forma matricial, p´agina 4 de 4