Métodos Numéricos. 2009.
Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
Sistemas de ecuaciones.
Primera práctica Computacional.
Método de Descomposición LU.
Se resolvió un sistema Ax=b mediante el método de descomposición LU, que consiste en descomponer a la matriz A del
sistema en un producto de una triangular inferior (L) y una triangular superior (U).
L se forma eliminando los elementos debajo de la diagonal, siendo L, una matriz tal que si se la multiplica su inversa por
izquierda a A, elimina los elementos de una columna dada. Luego se obtiene U por medio de invL*A, operación que se
A es una matriz tridiagonal, y b es un vector segundo miembro.
realiza al obtener L.
2 1 0 0 0 0.027778
Para n=5 la matriz A y el vector b tienen la siguiente pinta:
1 2 1 0 0
0.027778
0 1 2 1 0 0.027778
0 0 1 2 1 0.027778
0 0 0 1 2 0.027778
La descomposición LU resulto:
1 0 0 0 0 2 1 0 0 0
0,5 1 0 0 0
0 1,5 1 0 0
0 0,66667 1 0 0 0 0 1,33333 1 0
0 0 0,75 1 0 0 0 0 1,25 1
0 0 0 0,80 1 0 0 0 0 1,2
El código LU usado es
#Resolucion de un sistema de ecuacion nxn
n=5
#Creacion de la Matriz A con formato lleno
S1 = sparse ([1:n],[1:n],2*ones(1,n),n,n,0);
S2 = sparse ([1:n-1],[2:n],-1*ones(1,n-1),n,n,0);
S3 = sparse ([2:n],[1:n-1],-1*ones(1,n-1),n,n,0);
A=full(S1+S2+S3);#Si se quiere la matriz en formato lleno
b=ones(1,n)/(n+1)**2;#Para generar el vector segundo miembro
l=eye(n,n);#Inicio a l como una matriz identidad nxn
u=A;#u=U comienza como A
t1=time();#empiezo “contar” el tiempo que tardo
for h=1:n-1
z=eye(n,n);#z se iniciara como identidad en cada iteracion
#Estoy haciendo a l
for i=h+1:n
z(i,h)=-u(i,h)/(u(h,h));
l(i,h)=-z(i,h);#escribo a l columna por columna con los z, debajo de la diagonal
#z tiene el signo contrario a l este me sirve para hacer z*aux(invl*A)
endfor
u=z*u;
endfor
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
Pivoteo Parcial.
A continuación se agrega al la descomposición LU, el pivoteo parcial, consiste en aplicar permutaciones a la matriz L*A
(L inicial= identidad) de modo que cuando haga la eliminación gaussiana en una determinada columna no cometa errores
de truncamiento o redondeo cuando el divisor es muy cercano a cero, para esto se intercambia las filas de tal manera que
en el vector por debajo de la diagonal el primer elemento sea el mayor en valor absoluto.
#Resolucion de un sistema de ecuacion nxn
n=4
#Creacion de la Matriz A con formato lleno o ralo
S1 = sparse ([1:n],[1:n],2*ones(1,n),n,n,0);
S2 = sparse ([1:n-1],[2:n],-1*ones(1,n-1),n,n,0);
S3 = sparse ([2:n],[1:n-1],-1*ones(1,n-1),n,n,0);
# A=full(S1+S2+S3);#Si se quiere la matriz en formato lleno
b=ones(1,n)/(n+1)**2;#Para generar el vector segundo miembro
# A=[2,1,1,0;4,3,3,1;8,7,9,5;6,7,9,8];#se puede probar con esta matriz nxn tambien
#Cargando
Perm=eye(n,n);#matriz de permutacion
l=eye(n,n);
u=A;
t1=time();#comienzo a “contar” el tiempo
for h=1:n-1
z=eye(n,n);
#Busco el maximo
[M,o]=max(abs(u(h:n,h)));
# p=z;
o=o+h-1;#corrijo la ubicación de la fila con el mayor elemento
#Genero un matriz que permutara, las permutaciones las puedo efectuar con matrices(#)
usando p mediante las operaciones sin #.
# if(o>h)
# p(o,h)=p(h,h);p(h,h)=0;p(h,o)=p(o,o);p(o,o)=0;
# endif
# p;
#Aplicando permutaciones
# u=p*u;
u2=u(h,1:n);u(h,1:n)=u(o,1:n);u(o,1:n)=u2;#u2 es un vector auxiliar que me ayudara en
distintas permutaciones
# l=p*l*p;
u2=l(h,1:n);l(h,1:n)=l(o,1:n);l(o,1:n)=u2;u2=l(1:n,h);l(1:n,h)=l(1:n,o);l(1:n,o)=u2;
# Perm=p*Perm;
u2=Perm(h,1:n);Perm(h,1:n)=Perm(o,1:n);Perm(o,1:n)=u2;
#Estoy haciendo a l
for i=h+1:n
z(i,h)=-u(i,h)/(u(h,h));
l(i,h)=-z(i,h);
endfor
u=z*u;
endfor
La descomposición LU quedo igual que sin la permutaciones debido al tipo de matriz A con la que se trabajo, tridiagonal.
El método LU es útil para resolver sistema de ecuaciones del tipo Ax=b, ya que te permite resolver el sistema mediante la
siguiente descomposición.
L*y=b
U*x=y
Donde las matrices de los dos sitemas son matrices triangulares fáciles de resolver.
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
Resolución del sistema.
La resolución del sistema Ly=b, se resuelve mediante el método conocido como Backward Sustitution, y se debe tener
presente que si consiguió la matriz L aplicando el pivoteo parcial, se le deben aplicar las mismas permutaciones al vector
b.
La resolución del sistema Ux=y, que involucra una matriz triangular inferior se puede resolver fácilmente con el método
de Forward Sustitution.
0,069444
El resultado para n=5 con la matriz A dada al inicio fue:
0,111111
0,125
0,111111
0,069444
Se modifico el tamaño de la matriz como lo muestra el siguiente cuadro, y se comparo los tiempos que tarda el programa
en resolver el sistema:
Metodo LU Metodo LU
tamaño de la matriz A Tiempo(seg) tamaño de la matriz A Tiempo(seg)
10 0,04497 160 2,72490
20 0,05228 170 3,05270
30 0,11282 180 3,42080
40 0,18348 190 4,80220
50 0,34082 200 4,86526
60 0,42939 210 4,83560
70 0,65442 220 5,22540
80 0,88536 230 7,42480
90 1,5217 240 12,313
100 1,39860 250 13,84000
110 1,32340 260 15,54100
120 1,52620 270 16,48900
130 1,82640 280 18,31900
140 2,08630 290 18,97200
150 2,39750 300 20,84500
Se obtuvo el siguiente grafico:
21,00000
t 16,00000
11,00000
6,00000
1,00000
0 50 100 150 200 250 300 350 n
-4,00000
Figura1. Tiempos en contra del tamaño de la matriz del sistema, método LU.
El comportamiento aproxima a una exponencial.
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
El código usado para la resolución es
#Resolucion de un sistema de ecuacion nxn
n=5
#Creacion de la Matriz A con formato lleno o ralo
S1 = sparse ([1:n],[1:n],2*ones(1,n),n,n,0);
S2 = sparse ([1:n-1],[2:n],-1*ones(1,n-1),n,n,0);
S3 = sparse ([2:n],[1:n-1],-1*ones(1,n-1),n,n,0);
A=full(S1+S2+S3);#Si se quiere la matriz en formato lleno
b=ones(1,n)/(n+1)**2;#Para generar el vector segundo miembro
#Cargando
Perm=eye(n,n);l=eye(n,n);u=A;
t1=time();
for h=1:n-1
z=eye(n,n);
#Busco el maximo
[M,o]=max(abs(u(h:n,h)));
# p=z;
o=o+h-1;
#Genero un matriz que permutara
# if(o>h)
# p(o,h)=p(h,h);p(h,h)=0;p(h,o)=p(o,o);p(o,o)=0;
# endif
#Aplicando permutaciones
# u=p*u;
u2=u(h,1:n);u(h,1:n)=u(o,1:n);u(o,1:n)=u2;
# l=p*l*p;
u2=l(h,1:n);l(h,1:n)=l(o,1:n);l(o,1:n)=u2;u2=l(1:n,h);l(1:n,h)=l(1:n,o);l(1:n,o)=u2;
# Perm=p*Perm;
u2=Perm(h,1:n);Perm(h,1:n)=Perm(o,1:n);Perm(o,1:n)=u2;
#Estoy haciendo a l
for i=h+1:n
z(i,h)=-u(i,h)/(u(h,h));l(i,h)=-z(i,h);
endfor
u=z*u;
endfor
b=Perm'*b';
#Forward sustitution
for k=1:n
y(1)=b(1);
q=0;
for i=1:k-1
q=(l(k,i)*y(i))+q;
endfor
y(k)=b(k)-q;
endfor
#Backward sustitution
x(n)=y(n)/u(n,n);
for k=1:n-1
d=0;
for i=1:k
d=d+u(n-k,n-k+i)*x(n-k+i);
endfor
x(n-k)=(y(n-k)-d)/u(n-k,n-k);
endfor
t2=time()-t1
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
Método Gradiente Conjugado
Se planteo resolver el sistema ya planteado Ax=b por medio del método Gradiente Conjugado. Este método pide que A
sea una matriz simétrica.
Se plantea el siguiente funcional " # , el gradiente de este funcional es $ y la solución x que hace
!
cero al gradiente es la solución exacta de sistema de ecuaciones. Se busca la solución empezando por una solución x
cualquiera, luego de esta nos moveremos en dirección del valor del gradiente evaluado en x, que forma un vector residuo.
El siguiente paso es resolver hasta donde moverse, la Figura2 muestra una idea de esto, ej. en dos dimensiones llegara a
la solución en dos pasos o dos iteraciones. Nos movemos sobre el semiespacio restringido por el gradiente de manera de
llegar a la solución en el menor número de pasos. Este método es conocido como Steepest Descent
Figura 2. Método del Gradiente.
Los residuos se calculan con $ % &% y los x se moverán de la siguiente manera %'( % ) *&% , r es la
dirección para donde se moverá x, y alfa indica hasta donde se moverá, y se la obtiene con la siguiente condición
, - % ) *&% . &%# &%
+ / 0 1 *
,* 0 &%# &%
El método converge en n pasos teóricamente, pero por redondeos del computador no logra llegar a la solución exacta, lo
que hacemos es poner un límite de tolerancia, de modo que converja en algún numero de iteraciones.
Un método de la elección de la dirección de descenso que hace más eficiente al método es hacer cada dirección calculada
ortogonal a las demás direcciones, por lo que indirectamente se calcula una base en la que puede ser escrita la solución,
siendo la dimensión de esta igual al tamaño de la matriz.
x se mueve en una dirección p, % %2( ) *3% y * p se mueve 3% &%2( ) 93%2( donde 9 4567 4567 , cada
4567 4567 4 4
85 !85 56: 56:
,
residuo irá cambiando de la siguiente forma & &%2( *3% .
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
El código usado en Gradientes conjugados es
#Resolucion de un sistema de ecuacion nxn
n=5
#Creacion de la Matriz A con formato lleno o ralo
S1 = sparse ([1:n],[1:n],2*ones(1,n),n,n,0);
S2 = sparse ([1:n-1],[2:n],-1*ones(1,n-1),n,n,0);
S3 = sparse ([2:n],[1:n-1],-1*ones(1,n-1),n,n,0);
# A=full(S1+S2+S3);#Si se quiere la matriz en formato lleno
A=S1+S2+S3#Si se quiere la matriz en formato RALO
b=ones(1,n)/(n+1)**2;#Para generar el vector segundo miembro
# Metodo del gradiente conjugado
x=rand(n,1);#Semilla
r0=b'-A*x;
tol=0.000001;#tolerancia
k=0;#iteracion
r=r0;
t3=time();#comienzo a “contar” el tiempo
while(((norm(r,2))>tol))#impongo mi condición, si no está al cuadrado converge mas rápido
k=k+1;
if(k==1)3condicion trivial
p=r0;
t=r0;
else
beta=(r'*r)/(t'*t);
p=r+(beta*p);
endif
alfa=(r'*r)/(p'*A*p);
x=x+alfa*p;
t=r;
r=r-alfa*A*p;
endwhile
x=x;
condicion=cond(A)
normr=norm(r,2)
t4=time()-t3
k
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
Se comparo para diferentes tamaños los tiempos de resolución del problema para matrices A con el formato lleno y ralas,
las del formato ralo tienen menores valores de tiempo porque evitan el realizar cálculos triviales.
25
t
20
15
Formato Lleno
10 Formato Ralo
0
-50 150 350 550 750 950 1150 n
Figura3. Comparación de tiempos de ejecución para matrices con formatos llenos y ralos.
Precondicionamiento con Escaleo Diagonal.
Cuando un sistema de ecuaciones Ax=b tiene una matriz con un numero de condición muy alto se puede proponer
precondicionar a la matriz A del sistema, la teoría pide que la matriz del sistema sea simétrica y definida positiva.
Lo que se hace es multiplica a A por una matriz C-1 de tal forma que se pase a resolver un sistema con una matriz mejor
condicionada (C-1A C-T) (Cx)= (C-1b).
Se puede demostrar que en realidad se estaría resolviendo el siguiente problema C2x=b, donde se define
convenientemente a M=C2, ya que luego ayudara en los cálculos. M es la matriz de precondicionamiento.
Usaremos el precondicionador de Jacobi o escaleo diagonal que consiste en M=diagonal(A).
La matriz A del sistema a resolver es una matriz mal condicionada, que se obtiene con las siguientes líneas.
rand("seed",1.0);
v=5/(n-1) * (0:n-1);
D=(10.**v);
A=eye(n)+0.01*rand(n,n);
[Q R]=qr(A);
A=Q'*diag(D)*Q;
Para n=5 la matriz A tiene la pinta
6.4102; ) 000 1.3118; ) 000 7.4171; 001 1.5980; ) 001 7.3425; ) 002
1.3118; ) 000 1.8436; ) 001 6.7005; 001 5.0926; ) 001 1.4622; ) 002
7.4171; 001 6.7005; 001 3.1631; ) 002 1.6258; ) 001 5.7095; ) 001
1.5980; ) 001 5.0926; ) 001 1.6258; ) 001 5.6308; ) 003 8.6142; ) 002
7.3425; ) 002 1.4622; ) 002 5.7095; ) 001 8.6142; ) 002 9.9986; ) 004
1 0.1206681 0.0164720 0.0841102 0.9171398
La matriz Ã= C-1A C-T
0.1206681 1 0.0087745 0.1580612 0.1076980
à 0.0164720 0.0087745 1 0.0121820 0.0101525
0.0841102 0.1580612 0.0121820 1 0.0363045
0.9171398 0.1076980 0.0101525 0.0363045 1
Sus números de condición son <-Ã. 24.005 > <-. 1.0000; ) 005
Donde evidentemente se mejora la condición de la matriz del sistema de ecuaciones.
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
El código del gradiente conjugado con precondicionamiento es
#Resolucion de un sistema de ecuacion nxn
n=5
#Creacion de la Matriz A con formato lleno o ralo
S1 = sparse ([1:n],[1:n],2*ones(1,n),n,n,0);
S2 = sparse ([1:n-1],[2:n],-1*ones(1,n-1),n,n,0);
S3 = sparse ([2:n],[1:n-1],-1*ones(1,n-1),n,n,0);
# A=full(S1+S2+S3);#Si se quiere la matriz en formato lleno
A=S1+S2+S3;#Si se quiere la matriz en formato RALO
b=ones(1,n)/(n+1)**2;#Para generar el vector segundo miembro
#Genera una matriz A simetrica, definida positiva y con autovalores entre 1 y 10^5 equiespaciados
rand("seed",1.0);
v=5/(n-1) * (0:n-1);
D=(10.**v);
A=eye(n)+0.01*rand(n,n);
[Q R]=qr(A);
A=Q'*diag(D)*Q;
#Precondicionador de Jacobi
M=diag(diag(A));C=sqrt(M);Cin=inv(C);
#Gradiente precondicionado
t5=time();
x=rand(n,1);
r=b'-A*x;
tol=0.000001;
e=0;
w=r'*(inv(M)*r);
while((norm(r,2))>tol)#si no está al cuadrado me converge mas rápido
z=inv(M)*r;
wold=w;w=r'*z;
e=e+1;
if(e==1)
d=z;
z=d;
else
beta=w/wold;
d=z+(beta*d);
endif
alfa=(r'*z)/(d'*A*d);
x=x+alfa*d;
r=r-alfa*A*d;
endwhile
t6=time()-t5
e
Ac=Cin*A*Cin';#nueva matriz del sistema
condicion=cond(Ac);
normr=norm(r,2);
Universidad de Cuyo. Instituto Balseiro
Métodos Numéricos. 2009. Burgos Fernando Ezequiel
Practica Nº4. Ingeniería Mecánica
En la siguiente grafica se muestra para diferentes tamaños de matrices los números de condición para un matriz bien
condicionada y otra mal condicionada.
1,20E+05
κ
1,00E+05
8,00E+04
6,00E+04 Matriz mal condicionada
Matriz precondicionada
4,00E+04
2,00E+04
0,00E+00 n
0 100 200 300 400 500 600
Figura 4. Números de condiciones contra tamaños de matrices del sistema, método gradiente conjugado.
Vemos que a medida que se agranda la matriz, el precondicionamiento, se acerca al número de condición de la matriz mal
condicionada.
Universidad de Cuyo. Instituto Balseiro