• Author / Uploaded
• Karla Altamirano Arias

• Categories
• Álgebra lineal
• Análisis matemático
• Conceptos matemáticos
• Objetos matemáticos
• Álgebra

Citation preview

Universidad Politécnica de Madrid–Escuela Técnica Superior de Ingenieros Industriales

Métodos Matemáticos de Especialidad Ingeniería Eléctrica

Sistemas de ecuaciones lineales Métodos directos de solución José Luis de la Fuente O’Connor [email protected] [email protected]

Clase_sisli_12.pdf 1/120

Índice Cuál es el problema; consideraciones teóricas Eliminación de Gauss

 Pivotación  Algoritmo  Número de operaciones  Método de Gauss-Jordan Matlab y los sistemas de ecuaciones lineales Factorización LU

 Métodos explícitos para su obtención ı Método de Crout ı Método de Doolittle 2/120

 Matlab y la factorización LU Solución de sistemas modificados Refinamiento iterativo Sistemas con matrices especiales

 Matrices simétricas ı Factorización LDLT ı Factorización de Cholesky: matrices simétricas definidas positivas ˘ Matlab y la factorización de Cholesky ı Matrices simétricas semidefinidas positivas ı Matrices simétricas indefinidas

3/120

Cuál es el problema; consideraciones teóricas – Se trata de dar solución a sistemas de ecuaciones del tipo a11x1 C a12x2 C


C a1nxn D b1 a21x1 C a22x2 C


C a2nxn D b2 :: :: :: :: am1x1 C am2x2 C


C amnxn D bm; lo que significa determinar los valores de las variables x1; : : : ; xn que hacen que se cumplan todas las igualdades. – A los números aij se les denomina coeficientes del sistema y a los bi términos independientes. 4/120

– Si se introducen las matrices 2 3 a11 a12


a1n 6 a21 a22


a2n 7 AD6 :: :: 7 4 :: 5; am1 am2


amn

2 3 x1 6x2 7 7 xD6 4 :: 5 xn

2

3

b1 6 b2 7 7 y bD6 4 :: 5 ; xm

el sistema se puede representar de forma más compacta por

Ax D b: – En general se supondrá que la matriz de coeficientes A 2 Rmn, x 2 Rn y b 2 Rm.

5/120

– Casos posibles: m=n m=n

·

=

m=n m=n

rango(A) = m = n

·

=

rango(A) = m = n

·

=

m=n

·

=

·

=

·

=

·

=

rango(A) < n < m rango(A) < n < m

·

=

rango(A) < m = n

rango(A) = m = n

rango(A) < m = n

1a

1b

1a

m>n m>n

1b

·

=

·

=

rango(A) = n < m rango(A) = n < m

·

=

m>n

m>n m>n m>n

rango(A) = n < m

rango(A) < n < m

2a

2b

2a

2b

m> b=[2;2;-2;-5]; >> A\b ans = 3.0000 4.0000 -1.0000 -2.0000 62/120

– Utilizando el script Gauss que hemos presentado: >> Gauss(A,b) ans = 3.0000 4.0000 -1.0000 -2.0000 – Utilizando otra utilidad de Matlab muy interesante: >> linsolve(A,b) ans = 3.0000 4.0000 -1.0000 -2.0000 63/120

Índice Cuál es el problema; consideraciones teóricas Eliminación de Gauss Matlab y los sistemas de ecuaciones lineales Factorización LU Solución de sistemas modificados Refinamiento iterativo Sistemas con matrices especiales

64/120

Factorización LU – El cálculo de A D LU se conoce como factorización o descomposición LU . Para resolver un sistema de ecuaciones lineales Ax D b, si A D LU , el problema se convierte en el de resolver LU x D b a través de dos sistemas de ecuaciones triangulares:

Ux D y

y Ly D b:

Esto es muy útil cuando se requiere resolver sistemas de ecuaciones en los que la matriz A es siempre la misma y sólo cambia es el término independiente. 65/120

La factorización LU y la eliminación de Gauss – Una forma indirecta de conseguir esta factorización LU es la propia eliminación de Gauss. En efecto, mediante unas permutaciones y unas transformaciones definidas por matrices elementales triangulares inferiores el método conseguía: Ln 1P n

1


L1 P 1 A

D U:

De este proceso, haciendo P D Pn 1


P1 y L D P.Ln 1P n 1


L2P 2L1P 1/ 1; se puede comprobar que se obtiene la factorización PA D LU : 66/120

Existencia y unicidad de la factorización LU Teorema Sea A una matriz cuadrada regular de orden n. Existe una matriz de permutación P y dos matrices, una triangular inferior y otra triangular superior, L y U , respectivamente, tales que

PA D LU : La matriz L tiene todos los coeficientes de la diagonal principal igual a 1 (triangular inferior unitaria).

Lema La matriz A admite una factorización LU si y sólo si se cumple que det.A k / ¤ 0; k D 1; : : : ; n:

Teorema Si una matriz regular A de orden n admite una factorización A D LU , donde L es una matriz triangular inferior de coeficientes diagonales 1 y U una triangular superior, esa factorización es única. 67/120

Métodos numéricos directos para la obtención de factorizaciones LU Método de Crout. Versión LU 1 – Supongamos que se desea obtener la factorización en la forma LU 1, donde U 1 designa una matriz triangular superior en la que todos los coeficientes de la diagonal principal son 1.

68/120

– Si la matriz A es de orden 3 y se quiere factorizarla de la forma 2

a11 a12 4a21 a22 a31 a32

3

2 32 3 a13 l11 0 0 1 u12 u13 a23 5 D 4l21 l22 0 5 40 1 u23 5 ; a33 l31 l32 l33 0 0 1

usando las reglas de multiplicación de matrices se obtendrá: 1a col. de L: l11 D a11 l21 D a21 l31 D a31 I 

2a fila de U :

l11 u12 D a12 l11 u13 D a13

2a col. de L:

l21 u12 C l22 D a22 l31 u12 C l32 D a32

2a fila de U :

l21 u13 C l22 u23 D a23

3a col. de L:

! u1j D a1j = l11 ; j D 2; 3I  ! l i 2 D ai 2

li1 u12 ; i D 2; 3I

! u2j D .a2j

l31 u13 C l32 u23 C l33 D a33

l21 u1j /= l22 ; j D 3I

! li 3 D ai 3

i 1 X

lij uj i ; i D 3:

j D1 69/120

– En general, las fórmulas de recurrencia que se pueden deducir de este proceso, denominado factorización LU de Crout, son: li1 D ai1; u1j D a1j = l11; k 1 X lik D aik lip upk ;

i D 1; 2; : : : ; n; j > 1; i  k;

pD1

0 ukj D @akj

k 1 X

1 lkp upj A

 lkk ;

j > k:

pD1

70/120

– Plasmadas en el algoritmo de Crout para factorizar una matriz regular A nn en la forma LU 1 resulta el de la tabla. for k D 1 to n for i D k to n l.i; k/

a.i; k/

k 1 X

l.i; p/u.p; k/

pD1

end for i D k C 1 0 to n u.k; i/

@a.k; i /

k 1 X

1

 l.k; p/u.p; i /A l.k; k/

pD1

end end 71/120

– La versión Matlab de este algoritmo es la que sigue. function [L U]=LUCrout(a) % Factorización LU por Crout n=size(a,1); L=zeros(n); U=eye(n); for k=1:n for i=k:n L(i,k)=a(i,k)-L(i,1:k-1)*U(1:k-1,k); end for i=k+1:n U(k,i)=(a(k,i)-L(k,1:k-1)*U(1:k-1,i))/L(k,k); end end

72/120

– Ahora bien, como apuntábamos en la eliminación de Gauss, se puede aprovechar la estructura de la matriz A para guardar en ella las nuevas matrices L y U . El mismo algoritmo quedaría así.

function [L U]=Crout_1(A) % Factorización LU por Crout n=size(A,1); for k=1:n i=k:n; A(i,k)=A(i,k)-A(i,1:k-1)*A(1:k-1,k); i=k+1:n; A(k,i)=(A(k,i)-A(k,1:k-1)*A(1:k-1,i))/A(k,k); end L=tril(A,0); U=triu(A,1)+eye(n,n);

73/120

– Factorizar

"

10 10 20 20 25 40 30 50 61

#

da como resultado 2 32 3 10 1 1 2 LU D 420 5 5 4 1 05 : 30 20 1 1 >> A=[10 10 20;20 25 40;30 50 61]; >> [L,U]=Crout_1(A) L = 10 0 0 20 5 0 30 20 1 U = 1 1 2 0 1 0 0 0 1 >> 74/120

Ejemplo

– Se desea factorizar la matriz 2 3 0,001 2,000 3,000 A D 4-1,000 3,712 4,6235 -2,000 1,072 5,643 en una máquina u ordenador con cuatro dígitos significativos.

75/120

– Las operaciones que se realizan en la máquina son: l11 D 0,001I l21 D -1,000I l31 D -2,000I   2,000 D 2000I u12 D f l 0,001   3,000 u13 D f l D 3000I 0,001 l22 D f l Œ3,712 C .1,000/.2000/ D 2004I l32 D f l Œ1,072 C .2,000/.2000/ D 4001I   4,623 C .1,000/.3000/ D 1,500 y u23 D f l 2004 l33 D f lŒ5,643 C (2,000)(3,000) (4,001)(1,500) D 5,642:

– Obsérvese que el cálculo de l33 conlleva la pérdida de tres dígitos por redondeo: el valor que debería obtenerse es 5,922. 76/120

Pivotación

– El ejemplo pone de manifiesto que, aunque se sepa que una matriz no es singular y que su factorización LU existe teóricamente, los errores de redondeo que se pueden producir al trabajar en una máquina pueden dar al traste con el resultado. – Es aconsejable realizar pivotación. Al final de un proceso con pivotación se obtendría

PA D LU es decir, no la factorización LU de la matriz original sino de PA.

77/120

– El algoritmo de Crout con pivotación parcial es el de la tabla for k D 1 to n for i D k to n l.i; k/

a.i; k/

k X1

l.i; p/u.p; k/

pD1

end Determinar índice p 2 fk; k C 1; : : : ; ng tal que ja.p; i /j D mKaxi j n ja.j; i /j. Intercambiar filas p y k. for i D k C 1 0 to n 1  k X1 @ A u.k; i / a.k; i / l.k; p/u.p; i / l.k; k/ pD1

end end

function [L U p]=CroutP(a) % Factorización LU por Crout con pivotación n=size(a,1); p=1:n; for k=1:n i=k:n; a(i,k)=a(i,k)-a(i,1:k-1)*a(1:k-1,k); [r,m]=max(abs(a(k:n,k))); m=m+k-1; if a(m,k)==0, continue, end if k~=m, a([k m],:)=a([m k],:); p([k m])=p([m k]); end i=k+1:n; a(k,i)=(a(k,i)-a(k,1:k-1)*a(1:k-1,i))/a(k,k); end L=tril(a,0); U=triu(a,1)+eye(n,n); 78/120

– Si se factoriza la matriz 2

3

10 10 20 420 25 405 ; 30 50 61 al final de este proceso, el vector IPVT.
/, que indica las pivotaciones realizadas, es Œ3, 2, 1T . >> A=[10 10 20;20 25 40;30 50 61]; >> [L U p]=CroutP(A) L = 30.0000 0 0 20.0000 -8.3333 0 10.0000 -6.6667 0.2000 U = 1.0000 1.6667 2.0333 0 1.0000 0.0800 p = 3 2 1 >> 79/120

– La matriz PA realmente factorizada es 2 3 2 32 3 30 50 61 30 1 1;6667 2;0333 420 25 405 D 420 8;3333 54 1 0;08005 : 10 10 20 10 6;6667 0; 2 1 – El algoritmo de Crout requiere O.n3=3/ multiplicaciones/divisiones y sumas/restas para la factorización de la matriz.

80/120

Método de Crout. Versión L1U – Si se quiere conseguir la factorización L1U de una matriz 3  3, 2 32 3 2 3 a11 a12 a13 u11 u12 u13 1 0 0 4a21 a22 a23 5 D 4l21 1 05 4 0 u22 u23 5 ; a31 a32 a33 l31 l32 1 0 0 u33

operando: 1a fila de U : u11 D a11 u12 D a12 u13 D a13 I 

1a col. de L:

l21 u11 D a21 l31 u11 D a31

2a fila de U :

l21 u12 C u22 D a22 l21 u13 C u32 D a23

2a col. de L:

l31 u12 C l32 u22 D a32

3a fila de U :

! li1 D ai1 =u11 ; i D 2; 3I  ! u2j D a2j ! li 2 D .ai 2

l31 u13 C l32 u23 C u33 D a33

l21 u1j ; j D 2; 3I li1 u12 /=u22 ; i D 3I

! u3j D a3j

j 1 X

l3i uij ; j D 3:

i D1 81/120

– Las fórmulas de recurrencia que se pueden deducir de este proceso son: u1j D a1j ; li1 D ai1=u11; k 1 X ukj D akj lkp upj ;

j D 1; 2; : : : ; n; j > 1; j  k;

pD1

0 lik D @aik

k 1 X

1 lip upk A

 ukk ;

i > k:

pD1

82/120

– El algoritmo de Crout para factorizar una matriz regular A nn en la forma L1U es el que sigue. for k D 1 to n for j D k to n u.k; j /

a.k; j /

k 1 X

l.k; p/u.p; j /

pD1

end for i D k C 10to n l.i; k/

@a.i; k/

k 1 X

1

 l.i; p/u.p; k/A u.k; k/

pD1

end end 83/120

– Su implementación en Matlab:

function [L,U]=Croutl1u(a) % Factorización L1U por Crout n=size(a,1); for k=1:n-1 i=k+1:n; a(i,k)=a(i,k)/a(k,k); a(i,i)=a(i,i)-a(i,k)*a(k,i); end L=tril(a,-1)+eye(n,n); U=triu(a);

El resultado con la matriz precedente es: >> [L U]=Croutl1u(A) L = 1 0 0 2 1 0 3 4 1 U = 10 10 20 0 5 0 0 0 1 >> L*U ans = 10 10 20 20 25 40 30 50 61 84/120

– La versión del algoritmo con pivotación en Matlab es esta. function [L U p]=CroutP1(a) % Factorización L1U por Crout con pivotación n=size(a,1); p=1:n; for k=1:n-1 [r,m]=max(abs(a(k:n,k))); m=m+k-1; if a(m,k)==0, continue, end if k~=m, a([k m],:)=a([m k],:); p([k m])=p([m k]); end i=k+1:n; a(i,k)=a(i,k)/a(k,k); j=k+1:n; a(i,j)=a(i,j)-a(i,k)*a(k,j); end L=tril(a,-1)+eye(n,n); U=triu(a);

85/120

– El resultado con este script para el último ejemplo: >> [L U p]=CroutP1(A) L = 1.0000 0 0.6667 1.0000 0.3333 0.8000 U = 30.0000 50.0000 0 -8.3333 0 0 p = 3 2 1 >> L(p,:)*U ans = 10 10 20 20 25 40 30 50 61

0 0 1.0000 61.0000 -0.6667 0.2000

86/120

– Con los recursos de Matlab: >> [L U P]=lu(A) L = 1.0000 0 0.6667 1.0000 0.3333 0.8000 U = 30.0000 50.0000 0 -8.3333 0 0 P = 0 0 1 0 1 0 1 0 0 >> P’*L*U ans = 10 10 20 20 25 40 30 50 61

0 0 1.0000 61.0000 -0.6667 0.2000

87/120

Obtención de la matriz inversa a partir de la factorización LU – Si se designa por X la matriz inversa de A 2 Rnn, los n vectores columna de X son los vectores solución de los sistemas Ax i D e i , i D 1; : : : ; n. – Si suponemos que tenemos la factorización PA D LU , donde P es una matriz de permutación, para obtener la inversa de A hay que resolver los 2n sistemas siguientes: Ly i D Pe i ;

U xi D y i ;

i D 1; : : : ; n:

Es decir 2n sistemas de ecuaciones lineales con matrices triangulares en los que sólo cambian los términos independientes. 88/120

Matlab y la factorización LU – Para resolver Ax D b con Matlab, usando la factorización LU de A, se utiliza [L U P]=lu(A) y luego se obtiene la solución del sistema haciendo

x=U\(L\(P*b))

89/120

– Apliquemos esta idea a uno de los ejemplos que manejamos: >> A=[2 1 0 4;0 -3 -12 -1;0 -1 -2 0;0 0 3 1]; >> b=[2;2;-2;-5]; >> [L U P]=lu(A) L = 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0.3333 0.6667 1.0000 U = 2.0000 1.0000 0 4.0000 0 -3.0000 -12.0000 -1.0000 0 0 3.0000 1.0000 0 0 0 -0.3333 P = 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 >> x=U\(L\(P*b)) x = 3.0000 4.0000 -1.0000 -2.0000 90/120

– Otra cuestión a tener muy en cuenta: % Ensayo tiempos LU: Tiemp_LU.m A=rand(200,200); tic for i=1:1000 b=rand(200,1); x=A\b; end toc tic [L U P] = lu(A); for i=1:1000 b=rand(200,1); x=U\(L\(P*b)); end toc 91/120

Índice Cuál es el problema; consideraciones teóricas Eliminación de Gauss Matlab y los sistemas de ecuaciones lineales Factorización LU Solución de sistemas modificados Refinamiento iterativo Sistemas con matrices especiales

92/120

Solución de sistemas modificados – Si en Ax D b se modifica el vector b pero no A, no es necesario rehacer la factorización LU para resolver el nuevo sistema. – Si se modifica ligeramente la matriz A, por ejemplo el coeficiente (j; k), con lo que A D A ˛ej e Tk , puede que no sea necesario tampoco recalcular la factorización en su totalidad. La fórmula de Sherman-Morrison-Woodbury proporciona la inversa de una matriz en términos de los vectores de una modificación a la misma de rango uno del tipo uvT :



A

uv

 1 T

D 

A 1 C A 1u 1

 1 vT A 1u vT A 1: 93/120

– Para resolver el nuevo sistema .A fórmula, se obtendría



x D A

 1 T

b  D A 1b C A 1u 1 uv

uvT /x D b, usando la

 1 vT A 1u vT A 1b;

operación que podría hacerse por partes: 1. Resolviendo Az D u, obteniendo z. 2. Resolviendo Ay D b, obteniendo y. 3. Calculando x D y C ..vT y/=.1

vT z//z.

– Como A ya está factorizada, este procedimiento requiere solo sustituciones inversas y productos interiores; es decir O.n2/ operaciones frente a las O.n3=3/ de la factorización.

94/120

Ejemplo – Consideremos la matriz 2 2 4 AD4 4 9 2 3

2 3 2 3 1 0 0 2 4 2 35 D 4 2 1 05 40 1 1 1 1 0 0 7

3 2 15 4

›› L

U

a la que se le efectúa una modificación consistente en cambiar el coeficiente (3,2) de -3 a -1. – En este caso

2

3 0 u D 4 05 2

2 3 0 y v D 415 ; 0

con lo que la matriz resultante es A

uvT . 95/120

– Con la factorización LU de A, se resuelve Az D u y Ay D b, dando 2 3 2 3 3=2 1 z D 4 1=25 y y D 4 25 : 1=2 2 – Por último, 2

3 2 3 2 3 1 3=2 7 vT y 2 4 1=25 D 4 45 : xDyC z D 4 25 C T 1 v z 1 1=2 2 1=2 0

96/120

Índice Cuál es el problema; consideraciones teóricas Eliminación de Gauss Matlab y los sistemas de ecuaciones lineales Factorización LU Solución de sistemas modificados Refinamiento iterativo Sistemas con matrices especiales

97/120

Refinamiento iterativo – Supongamos que se tiene una solución aproximada, x 0, del sistema de ecuaciones lineales Ax D b, y sea y una corrección o mejora de la misma tal que la solución exacta, x, cumple que x D x 0 C y: Sustituyendo esta expresión en Ax D b se tiene que Ay D b

Ax 0 D r 0;

donde r 0 es el vector de residuos.

98/120

– Si este vector no cumple unos requisitos de precisión que nos interesen, se pude resolver el sistema Ay D r 0 y hacer x 1 D x 0 C y; lo que hará que la solución se aproxime más a x que x 0. – Si es necesario, se calcula un nuevo vector de residuos, r 1 D b Ax 1 y se continua el proceso hasta que la solución se aproxime tanto como se quiera a la esperada. – El vector de residuos debe calcularse con más precisión que la usada para calcular la solución inicial.

99/120

– Este script de Matlab lleva a cabo el proceso a mano. % Script_Ref.m - Script de Refinamiento Iterativo n=6; format short A=hilb(n); b=A*ones(n,1); pause x=A\b

% Matriz de Hilbert (muy mal condicionada) % Elegimos término independiente para sol. x=1. % Solución, evidentemente, =1

B=A; % En B está A perturbada un poquito B(6,1)=B(6,1)+1.e-06; pause x1=B\b % Veamos la nueva solución; difiere bastante pause xex=ones(n,1); norm(xex-x1,2) norm(xex-x,2) pause

% Calculemos cuánto % Como magnitud calculemos la norma 2 de la desviaci.

res=b-A*x1; x1=x1+B\res norm(xex-x1,2) pause

% Hagamos una iteración del Refinamiento iterativo

res=b-A*x1; format long x1=x1+B\res norm(xex-x1,2) pause

% Hagamos otra iteración del Refinamiento iterativo

res=b-A*x1; x1=x1+B\res norm(xex-x1,2)

% Hagamos otra iteración del Refinamiento iterativo

100/120

Índice Cuál es el problema; consideraciones teóricas Eliminación de Gauss Matlab y los sistemas de ecuaciones lineales Factorización LU Solución de sistemas modificados Refinamiento iterativo Sistemas con matrices especiales

101/120

Sistemas con matrices especiales Matrices simétricas Factorización LDLT Lema Si todas las submatrices principales de una matriz A 2 Rnn son regulares, existen dos matrices triangulares inferiores unitarias únicas, L y M , y otra diagonal también única, D D diag.d1; : : : ; dn/, tales que A D LDM T .

Teorema Si A admite una factorización LDM T y es simétrica, L D M.

102/120

– Para derivar unas fórmulas simbólico de orden 3, 2 3 2 a11 a12 a13 1 4a21 a22 a235 D 4l21 l31 a31 a32 a33

de recurrencia, a partir de un ejemplo 32

32

3

0 0 d11 1 l21 l31 5 40 1 l325 ; 1 05 4 d22 l32 1 d33 0 0 1

operando de acuerdo con las reglas de multiplicación matricial se obtiene que a11 D d11 a21 D l21d11 a31 D l31d11 2 a22 D l21 d11 C d22 a32 D l31l21d11 C l32d22 2 2 a33 D l31 d11 C l32 d22 C d33: 103/120

– Generalizando se obtiene el algoritmo de la tabla. for k D 1 to n

d.k/

a.k; k/

k 1 X

a2 .k; p/d.p/

pD1

if d.k/ D 0 then stop for i D k C 1 0 to n

a.i; k/

@a.i; k/

k 1 X

1 a.i; p/a.k; p/d.p/A

 d.k/

pD1

end end

– Requiere O.n3=6/ multiplicaciones y divisiones y sumas y restas. – Si no se efectúan pivotaciones, el procedimiento numérico puede fallar por la posible presencia de coeficientes pivote muy pequeños, o por la acumulación de errores de redondeo o de cancelación importantes. 104/120

Factorización de Cholesky – Una matriz es definida positiva si para todo x ¤ 0 se cumple que x T Ax > 0: Todos sus valores propios son positivos. – Las matrices simétricas definidas positivas admiten una descomposición de la forma A D GT G; donde G es una matriz triangular superior. – Esta descomposición fue formulada por André Louis Cholesky (1875-1918), comandante del ejército francés de la época, durante la ocupación internacional de Creta en 1906–09. 105/120

– Las matrices simétricas definidas positivas se presentan habitualmente en: Problemas relacionados con el análisis de sistemas eléctricos de generación y transporte de energía. Ajuste de funciones por mínimos cuadrados. Análisis de estructuras mecánicas. En muchos procedimientos de optimización lineal y no lineal. – En general, en todas aquellas aplicaciones donde al modelizar un sistema, la expresión x T Ax mide la energía presente, o disponible, o cualquier otra magnitud física que sólo admite cantidades positivas en un entorno determinado. 106/120

Lema Las submatrices principales de una matriz definida positiva son definidas positivas.

Teorema Si A es una matriz definida positiva de orden n, tiene una descomposición de la forma LDM T , siendo todos los coeficientes de la matriz diagonal D positivos.

Teorema Si A es una matriz simétrica definida positiva de orden n, existe una única matriz triangular superior, G , con todos sus coeficientes diagonales positivos, tal que A D G T G .

107/120

– Procedamos a simular el algoritmo con la descomposición simbólica de una matriz 3  3. – Si se desea obtener la factorización 2 32 3 2 3 a11 a12 a13 g11 0 0 g11 g12 g13 4a12 a22 a235 D 4g12 g22 0 5 4 0 g22 g235 ; g13 g23 g33 a13 a23 a33 0 0 g33 operando de acuerdo con las obtiene que: a11 D a12 D a13 D a22 D a23 D a33 D

reglas de multiplicación matricial se 2 g11 g11g12 g11g13 2 2 g12 C g22 g12g13 C g22g23 2 2 2 g13 C g23 C g33 : 108/120

– Generalizando este proceso se obtiene el algoritmo que describe la tabla. for i D 1 to nv g.i; i /

u u u ta.i; i /

i 1 X

g 2.k; i /

kD1

for j D i C 1 0 to n g.i; j /

B @a.i; j /

i 1 i

1  C g.k; i /g.k; j /A

g.i; i /

kD1

end end

– El algoritmo requiere O.n3=6/ operaciones de multiplicación+división y de suma+resta. 109/120

– Este algoritmo en Matlab sería como sigue. function G=Chols_1(A) % Factorización de Cholesky n=size(A,1); for i=1:n, j=i+1:n; A(i,i)=sqrt(A(i,i)); A(i,j)=A(i,j)/A(i,i); A(j,j)=A(j,j)-A(i,j)’*A(i,j); end G=triu(A);

110/120

– La factorización de

2

5 6 1 6 4 2 0

1 2 0 0

2 0 4 1

3

0 07 7W 15 3

>> A=[5 1 -2 0;1 2 0 0;-2 0 4 1;0 >> G=Chols_1(A) G = 2.2361 0.4472 -0.8944 0 1.3416 0.2981 0 0 1.7638 0 0 0 >> G=chol(A) G = 2.2361 0.4472 -0.8944 0 1.3416 0.2981 0 0 1.7638 0 0 0 >>

0 1 3];

0 0 0.5669 1.6366

0 0 0.5669 1.6366 111/120

Matlab y la factorización de Cholesky – Para resolver un sistema lineal de ecuaciones Ax D b con Matlab utilizando la factorización de Cholesky hay que utilizar la función G=chol(A). – La solución del sistema correspondiente se puede obtener, teniendo en cuenta que se realiza A D G T G , haciendo

x=G\(G’\b)

112/120

– Para utilizar esta operación con un ejemplo de los que estamos manejando, habría que hacer algo parecido a lo que sigue. >> A=[5 1 -2 0;1 2 0 0;-2 0 4 1;0 0 1 3]; >> b=[1;5;14;15]; >> G=chol(A) G = 2.2361 0.4472 -0.8944 0 0 1.3416 0.2981 0 0 0 1.7638 0.5669 0 0 0 1.6366 >> x=G\(G’\b) x = 1.0000 2.0000 3.0000 4.0000 113/120

Matrices simétricas semidefinidas positivas – Una matriz A se dice semidefinida positiva, si para todo x ¤ 0, x T Ax  0. Teorema Si A 2 Rnn es simétrica semidefinida positiva: jaij j  .ai i C ajj /=2 p ai i ajj .i ¤ j / jaij j  mKax jaij j D mKax ai i i;j

i

ai i D 0 ) aij D aj i D 0; j D 1; : : : ; n:

– Si el algoritmo de Cholesky se aplica a una matriz semidefinida positiva, y en un paso akk es cero, entonces aj k D 0; j D k; : : : n, por lo que no habría que hacer nada más en la columna k. En la práctica, los errores de redondeo internos impiden los ceros exactos por lo que se recurre a la pivotación. 114/120

Pivotación

– Para mantener la simetría, las pivotaciones han de ser simétricas: si se intercambian dos filas, también hay que intercambiar las columnas simétricas: A PAP T . – La pivotación en Cholesky se lleva adelante así: En cada etapa k del proceso se determina el elemento de mayor valor de la diagonal principal, mKaxkj n ajj : Si no es cero, se intercambian las filas/columnas p y k, siempre y cuando k ¤ p; si es cero, el resto de la matriz a factorizar sería nula y no se haría nada más. 115/120

– Este es el algoritmo de Cholesky con pivotación para matrices semidefinidas positivas. for i D 1 to n Determinar índice p 2 fi; i C 1; ng tal que ja.p; p/j D maK xi j n fja.j; j /jg if a.p; p/ > 0 Intercambiar sfilas/columnas p y i .

g.i; i /

a.i; i /

i 1 X

g 2 .k; i /

kD1

for j D i C 1 0 to n

g.i; j /

B Ba.i; j / @

i 1 i

1 C g.k; i /g.k; j /C A

 g.i; i /

kD1

end end end

116/120

Matrices simétricas indefinidas – Una matriz A se dice indefinida si para algún vector x ¤ 0 la forma cuadrática x T Ax es positiva y para otros negativa. – Para factorizar este tipo de matrices se recurre a descomposiciones de pivotación diagonal en bloques de la forma PAP T D LBLT donde la matriz L es triangular inferior unitaria y la matriz B es tridiagonal, o diagonal en bloques, con bloques de dimensión 1  1 ó 2  2, bidiagonal en este caso. – Casi todos los códigos modernos utilizan alguna variedad de matriz bidiagonal B en bloques, aunque todavía se usan mucho rutinas que implementan algún método en el que esa matriz es tridiagonal (T ). 117/120

– Los métodos más conocidos se citan a continuación. Método

Estrategia

Operaciones

Parlett y Reid

PAP T D LT LT

Aasen

PAP T D LT LT

O.n3=3/

Bunch y Parlett

PAP T D LBLT O.n3 =6/ C O .n3 =6/ compara.

O.n3=6/

Bunch y Kaufman PAP T D LBLT O.n3=6/ C .n2

1/ compara.

– El del Bunch y Kaufman (1977), en alguna de sus variantes, es el más utilizado y el que emplean los códigos profesionales para factorizar matrices simétricas.

118/120

function [L D P rho] = diagpiv(A) %DIAGPIV Diagonal pivoting factorization with pivoting of a symetric A. % P*A*P’=L*D*L’; L is triangular and D a block diagonal D 1x1 or 2x2. % Rho is the growth factor. This routine does not exploit symmetry. % Bunch and Kaufman (1977), Some stable methods for calculating inertia % and solving symmetric linear systems, Math. Comp. 31(137):163-179. if norm(triu(A,1)’-tril(A,-1),1), error(’Matrix must be symmetric.’), end n = max(size(A)); k = 1; D = eye(n); L = eye(n); pp = 1:n; normA = norm(A(:),inf); rho = normA; alpha = (1 + sqrt(17))/8; while k < n [lambda r] = max(abs(A(k+1:n,k))); r = r(1) + k; if lambda > 0 swap = 0; if abs(A(k,k)) >= alpha*lambda s = 1; else temp = A(k:n,r); temp(r-k+1) = 0; sigma = norm(temp, inf); if alpha*lambda^2 =alpha*sigma swap = 1; m1 = k; m2 = r; s = 1; else swap = 1; m1 = k+1; m2 = r; s = 2; end end if swap A([m1 m2],:) = A([m2 m1],:); L([m1 m2],:) = L([m2 m1],:); A(:,[m1 m2]) = A(:,[m2 m1]); L(:,[m1 m2]) = L(:,[m2 m1]); pp([m1 m2]) = pp([m2 m1]); end if s == 1 % s = 1 D(k,k) = A(k,k); A(k+1:n,k) = A(k+1:n,k)/A(k,k); L(k+1:n,k) = A(k+1:n,k); i = k+1:n; A(i,i) = A(i,i) - A(i,k)*A(k,i); else % s = 2 E = A(k:k+1,k:k+1); D(k:k+1,k:k+1) = E; C = A(k+2:n,k:k+1); temp = C/E; L(k+2:n,k:k+1) = temp; A(k+2:n,k+2:n) = A(k+2:n,k+2:n) - temp*C’; end

if k+s = 3, P = eye(n); P = P(pp,:); end rho = rho/normA;

119/120

– Con una matriz de Hankel, el programa funciona así: >> A = gallery(’ris’,6) A = 0.0909 0.1111 0.1429 0.1111 0.1429 0.2000 0.1429 0.2000 0.3333 0.2000 0.3333 1.0000 0.3333 1.0000 -1.0000 1.0000 -1.0000 -0.3333 >> cond(A) ans = 2.2185 >> [L D P rho]=diagpiv(A) L = 1.0000 0 0 0 1.0000 0 -0.1760 0.2160 1.0000 -0.3143 0.1714 -1.1905 -0.1048 0.3429 0.2646 -0.9778 0.2000 -0.6173 D = 0.0909 1.0000 0 1.0000 -0.1111 0 0 0 -0.9216 0 0 0 0 0 0 0 0 0 P = 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 rho = 1.9264

0.2000 0.3333 1.0000 -1.0000 -0.3333 -0.2000

0.3333 1.0000 -1.0000 -0.3333 -0.2000 -0.1429

1.0000 -1.0000 -0.3333 -0.2000 -0.1429 -0.1111

0 0 0 1.0000 -0.6667 0.6222

0 0 0 0 1.0000 0

0 0 0 0 0 1.0000

0 0 0 1.7415 0 0

0 0 0 0 -0.8256 1.9264

0 0 0 0 1.9264 0.1284

0 1 0 0 0 0

>> eig(A) ans = -1.5708 -1.5705 -1.4438 0.7080 1.5622 1.5708 >>

120/120