Story Transcript
Lecci´ on E
Sistemas de ecuaciones lineales. M´ etodos directos En esta lecci´on estudiamos la soluci´on de un sistema de Cramer Ax = B, lo que significa que A es regular o invertible o que verifica det(A) 6= 0, mediante un m´etodo directo. Siendo ´este cualesquiera que permita, en ausencia de errores, mediante un n´ umero finito de pasos obtener la soluci´on exacta. En propiedad, esto no ocurre en general debido a los inevitables errores de redondeo.
E.1.
Sistema triangular inferior. Sustituci´ on progresiva
Consideremos el sistema Lx = B con L una matriz triangular inferior, es decir que verifica lij = 0 para j > i lo que significa que nuestro sistema toma la forma: l11 x1
= b1
l21 x1 + l22 x2 .. .. . .
.
= b2 .. .
ln1 x1 + ln2 x2 + . . .
+ lnn xn = bn
..
La regla de sustituci´on progresiva consiste en despejar x1 de la primera ecuaci´on (es inmediato poniendo x1 = b1 /l11 ) y continuar el proceso mediante la sustituci´on de x1 en la segunda ecuaci´on para obtener el valor de x2 y as´ı sucesivamente. En definitiva, la u ´nica soluci´on (x1 , x2 , . . . , xn ) est´a caracterizada por ser: xk =
1 lkk
bk −
k−1 X j=1
xj
lkj lkk
Pr´ actica a Con la orden type e_sp obtenemos el listado de la funci´on Mt\e_sp.m definida en el libro [6], p´ag 163, y que permite resolver un sistema cuadrado, regular y triangular inferior en el que se dan como datos la matriz regular L y el vector B. Como ejemplo de aplicaci´on de dicha funci´on resolvemos el sistema de ecuaciones x1
= 1
2x1 + 3x2
= 2
x1 + 2x2 + 3x3 = 4 31
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
32
para lo que ejecutamos el gui´on flops(0); %comenzamos a contar el n´ umero de operaciones A=[1 0 0;2 3 0;1 2 3]; b=[1 2 4].’; x=e_sp(A,b).’, flops con ´el obtenemos la u ´nica soluci´on del sistema, x = (1, 0, 1), y el n´ umero de operaciones “flops” necesarias para obtenerlo, 12. En general, para un sistema regular y triangular inferior de n ecuaciones reales el n´ umero de operaciones necesarias para resolverlo es del orden de 2 n(n + 1) = O(n ). Resolvemos tambi´en los sistemas triangulares inferiores siguientes: 7x1
= 7
x1 − x2
= 4.1
;
0.5x1 + 7.82x2 − 24.97x3 = −151.114
E.2.
28x1
= 7
34x1 − 26.2x2
= 4.1
3x1 + 7.82x2 − 2.497x3 = −15.1
Sistema triangular superior. Sustituci´ on regresiva
Consideramos ahora el sistema U x = B con U una matriz triangular superior, es decir que verifica sij = 0 si i > j. La regla de sustituci´on regresiva consiste en despejar de abajo hacia arriba de forma an´aloga a la anterior con lo que obtenemos que la u ´nica soluci´on x = (x1 , . . . , xn ) del sistema U x = B est´a dado por la regla: n X ukj 1 xj xk = bk − u kk ukk j=k+1
Pr´ actica b Con la orden type e_sr obtener y estudiar el listado de la funci´on Mt\e_sr.m definida en el libro [6], p´ag 165, que resuelve mediante la regla de sustituci´on regresiva un sistema cuadrado, regular y triangular superior, a partir de los datos: U la matriz del sistema y B el vector de los t´erminos independientes. Probar la eficacia del programa resolviendo los sistemas: x1 + 4x2 + 5x3 = 1 3x2 + 2x3 = 2 ; 3x3 = 6
E.3.
3x1 − x2 + 2x3 = 12 7x2 + 7x3 = 21 −21x3 = −42
M´ etodos exactos para sistemas generales
Antes de comenzar a describir los m´etodos de resoluci´on de sistemas generales vamos a recordar ´ algunos elementos de Algebra Lineal Elemental Sea E : Mm×n (K) → Mm×n (K) una transformaci´on en el conjunto de matrices de orden m × n. Diremos que E es una transformaci´ on elemental si E(A) se obtiene a partir de A intercambiando la columna (resp. fila) i-´esima por la j-´esima, i 6= j.
´ E.3. METODOS EXACTOS PARA SISTEMAS GENERALES
33
E(A) se obtiene a partir de A multiplicando la columna (resp. fila) i-´esima por λ 6= 0. E(A) se obtiene a partir de A sumando a la i-´esima columna (resp. fila) la j-´esima columna (resp. fila) multiplicada por λ = 6 0, i 6= j. Es importante observar: 1. Si denotamos por C(A) el resultado de una transformaci´on elemental en columnas entonces C(A) = A · C(Idn ) y que si F (A) denota el resultado de una transformaci´on en filas entonces F (A) = F (Idm ) · A Adem´as si C = Cn ◦· · ·◦C1 y F = Fr ◦· · ·◦F1 entonces tambi´en C(A) = M ·Cn (. . . (C1 (Idn )) . . . ) y F (A) = Fr (. . . (F1 (Idm )) . . . ) · M . Las dos igualdades encuadradas anteriores justifican que las matrices C(Idn ) y F (Idm ) se denominen matrices de las transformaciones elementales C y F respectivamente. 2. Cada transformaci´on elemental posee determinante distinto de cero y su inversa es otra transformaci´on elemental. En consecuencia, si E es una transformaci´on cualesquiera en filas o columnas y F es una transformaci´on en filas entonces se verifica: a) Rango E(A) = Rango A. b) Si A∗ = [A, B] es la matriz ampliada del sistema Ax = B entonces F (A∗ ) es la matriz ampliada de un sistema equivalente, en el sentido de que ambos poseen el mismo conjunto de soluciones. Como consecuencia del apartado (b) anterior se deduce: 1. Todo sistema Ax = B regular es equivalente a un sistema U x = B 0 donde U es una matriz triangular superior. 2. Existe una permutaci´on en columnas P de la matriz regular A de modo que A · P factoriza en el producto de una matriz triangular inferior L y otra triangular superior U . Es decir, LU = AP . El anterior enunciado se puede tambi´en poner: Existe una permutaci´ on en filas P de la matriz regular A de modo que P · A factoriza en el producto de una matriz triangular inferior L y otra triangular superior U . Es decir, LU = P A. En efecto, la mec´anica que convierte A∗ puede detallar con un ejemplo: Consideramos la matriz a 11 a21 a31
en una matriz triangular superior equivalente U ∗ se a12 a13 b1
a22 a23 b2 a32 a33 b1
34
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
Supuesto a11 6= 0 podemos considerar la transformaci´on elemental F1 que en filas suma a la segunda el resultado de multiplicar por −a21 /a11 la primera y que a la tercera fila suma el resultado de multiplicar por −a31 /a11 la primera fila. Obtenemos:
a a a b 11 12 13 1 F1 (A ) = 0 a022 a023 b02 ; 0 a032 a033 b01
1 0 0 F1 = −a21 /a11 1 0 −a31 /a11 0 1
∗
El proceso realizado corresponde a lo que se denomina utilizar a11 como pivote. En el siguiente paso, supuesto a22 6= 0, utilizando a22 como pivote conseguimos:
a a a b 11 12 13 1 0 0 F2 (A ) = 0 a22 a23 b02 ; 0 0 a0033 b001 ∗
1 0 0 F2 = a 1 0 b c 1
Con lo que habr´ıamos terminado ya que F2 · A∗ = U ∗ y tambi´en F2 · A = U de donde Ax = B es equivalente a un sistema triangular superior y A = F2−1 U = LU . Lo descrito anteriormente puede fallar si alguno de los pivotes aii que hemos ido eligiendo es 0. En dicho caso existe un elemento en su misma fila distinto de cero, aik 6= 0 para alg´ un k > i, ya que en caso contrario la matriz A ser´ıa singular, con lo cual si hacemos sobre A la transformaci´on elemental Tik que corresponde a intercambiar las columnas i, k tendr´ıamos que a la matriz A · Tij (Id) se le puede aplicar un paso m´as en el proceso de factorizar como producto de dos matrices triangulares y en definitiva eso nos prueba que existe una permutaci´on P de modo que A · P es factorizable de la forma A · P = LU con L una matriz triangular inferior y U una matriz triangular superior. Si aplicamos esto a la matriz At concluimos que existe P = P t de modo que L0 U 0 = At P t pero entonces U 0t L0t = P A luego LU = P A. Describimos ahora los algoritmos m´as importantes de resoluci´on directa de sistemas regulares generales. Eliminaci´ on gaussiana La eliminaci´on gaussiana es la que empleamos para reducir un sistema regular Ax = B a un sistema triangular superior U x = B 0 . El coste operativo de resolver un sistema utilizando la eliminaci´on gaussiana es O( 32 n3 ). Factorizaci´ on LU Sea P una permutaci´on en filas de la matriz A para la que existe la factorizaci´on P · A = LU , claramente para P se verifica P t = P −1 . El sistema Ax = B se puede expresar LU x = P · B de modo que se puede reducir a la resoluci´on de los sistemas: Lz = P · B
Triangular inferior
Ux = z
Triangular superior
La factorizaci´on LU est´a implementada en Matlab mediante la orden [L,U,P]=lu(A)
´ E.3. METODOS EXACTOS PARA SISTEMAS GENERALES
35
y donde P · A = LU . El coste operativo de resolver un sistema utilizando la factorizaci´on LU es el mismo que en la eliminaci´on gaussiana pero tiene la ventaja de que una vez conseguida la factorizaci´on P · A = LU la podemos emplear en todos aquellos sistemas Ax = B que tengan la misma matriz A. Un ejemplo de esta situaci´on es el c´alculo de la inversa. Factorizaci´ on de Cholesky ´ Tal como se puede ver en el u ´ltimo problema de Algebra Lineal de primero toda matriz sim´etrica t A definida positiva posee una u ´nica descomposici´on A = L L con L una matriz triangular inferior. En consecuencia, el sistema Ax = B es equivalente a resolver: Lt z = B
Triangular inferior
Lx = z
Triangular superior
La factorizaci´on de Cholesky est´a implementada en Matlab mediante la orden L=chol(A) y donde A = Lt L. El coste operativo de la resoluci´on por el m´etodo de Cholesky es = O( 13 n3 ). Factorizaci´ on QR Se factoriza A en la forma A = QR, con Q una matriz ortogonal, Qt Q = Id, y R triangular superior. De este modo el sistema Ax = B queda equivalente a la resoluci´on de los sistemas: Qz = B
De soluci´on z = Qt · B
Rx = z
Triangular superior
La factorizaci´on QR est´a implementada en Matlab mediante la orden [Q,R]=qr(A) y donde A = QR. El coste operativo de la factorizaci´on QR es O(2n3 ). M´ etodo empleado por Matlab El m´etodo que emplea Matlab para resolver el sistema regular Ax = B cuando ejecutamos la orden A\B es el ´optimo de entre los anteriores. Con la precauci´on de que no siempre Matlab resuelve el sistema AX = B. Esto que decimos lo observamos en las siguientes pr´acticas: Pr´ actica c Consideramos el sistema x+2y = 3 del cual observamos que posee infinitas soluciones. Si ejecutamos en Matlab el listado A=[1 2], b=[3] x=A\b, x=x.’
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
36
obtenemos como respuesta x = (0, 1.5) de lo que podemos observar que Matlab s´olo nos da una de las soluciones. Pr´ actica d
Consideramos el sistema r1 : x − y
= 0
r2 : x + y
= 0
r3
: 2x + y = 1
(E.1)
Con el listado x=-0.5:0.05:1.5; y1=x; y2=-x; y3=1-2*x; a1= 0.2857; a2= 0.1429; plot(x,y1,’r’,x,y2,’b’,x,y3,’g’,a1,a2,’+’) obtenemos la representaci´on gr´afica de las rectas r1 , r2 y r3 que definen el sistema (E.1). De esta representaci´on observamos que el sistema es incompatible ya que no hay intersecci´on com´ un a las tres rectas. Sin embargo, si resolvemos el sistema por Matlab, es decir, si ejecutamos A=[1 -1; 1 1; 2 1], b=[0 0 1].’ a=A\b; obtenemos como u ´nica soluci´on a = (a1,a2) = (0.2857; 0.1429). ¿Qu´e significa esto?: Del sistema AX = b dado en (E.1) obtenemos At Ax = At b multiplicando ambos miembros por la traspuesta de A. Para este nuevo sistema se verifica que su matriz asociada, At A, es regular. En efecto, At A es claramente cuadrada y verifica que Ker At A = Ker A = 0. La igualdad Ker At A = Ker A es general y f´acil de probar: Claramente Ker At A ⊃ t t t Ker √ A y ahora si x ∈ Ker A A entonces x (A Ax) = 0 lo que implica que ||Ax||2 = xt At · Ax = 0 y por tanto que Ax = 0 y x ∈ Ker A. En consecuencia, At Ax = At b posee una u ´nica soluci´on de la que podemos comprobar con el listado c=(A.’*A)\(A.’*b) que coincide con a. Adem´as, para todo vector z se tiene z t (At Ac−At b) = 0 y por tanto que A(z) es perpendicular a Ac − b para todo z de donde ||Ay − b||2 = ||Ac − b + A(y − c)||2 = ||Ac − b||2 + ||A(y − c)||2
(E.2)
lo que implica que la expresi´on ||Ay − b||2 es m´ınima precisamente cuando y = c. Podemos resumir todo lo anterior diciendo: Para A una matriz inyectiva, es decir, tal que su rango sea m´aximo, se verifica que es equivalente: Resolver el sistema regular At Ax = At b. Resolver en el sentido de los m´ınimos cuadrados el sistema Ax = b. Por ello entendemos obtener el vector c con la propiedad de ser ||Ac − b||2 m´ınimo.
´ E.4. METODO DIRECTO PARA SISTEMAS TRIDIAGONALES
37
Calcular c = (At A)−1 At b. La aplicaci´on cuya matriz es (At A)−1 At recibe el nombre de pseudo-inversa de A y est´a implementada en Matlab con el comando pinv(A). Ejecutar A\b ´o pinv(A)*b en Matlab. Como observaci´on final a esta pr´actica y a esta exposici´on decir que si la transformaci´on TA : E → F no es inyectiva entonces la transformaci´on que TA induce en el cociente de forma natural, T A : E/ Ker TA → F , es inyectiva con lo que todo lo anterior es directamente aplicable a T A y por tanto podemos afirmar: El sistema At Ax = At b es siempre compatible y su conjunto de soluciones es una variedad af´ın de direcci´on Ker A. El sistema Ax = b siempre es compatible en el sentido de los m´ınimos cuadrados y su conjunto de soluciones es una variedad af´ın de direcci´on Ker A. Pr´ actica dd Matlab tambi´en cuenta con la posibilidad de resolver de forma ‘exacta’ o simb´ olica un sistema de ecuaciones. Como ejemplo resolvemos simb´ olicamente y num´ericamente los siguientes sistemas lineales
4 −3 3 x 0 1 6 1 −9 x2 = 9 2 −5 −6 x3 5
;
2 1 3 x 0 1 4 2 −1 x2 = 0 6 3 2 x3 0
con el listado
A1=[4 -3 3;6 1 -9;2 -5 -6]; As1=sym(A1) b1=[0 9 5].’ x1=A1\b1 %soluci´ on num´ erica xs1=As1\b1 %soluci´ on simb´ olica A2=[2 1 3;4 2 -1;6 3 2]; As2=sym(A2) b2=[0 0 0].’ x2=A2\b2 %soluci´ on num´ erica xs2=As2\b2 %soluci´ on simb´ olica
E.4.
M´ etodo directo para sistemas tridiagonales
Los sistemas tridiagonales aparecen en la resoluci´on de numerosos problemas num´ericos, de aqu´ı su importancia. Un sistema tridiagonal tiene la forma:
38
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
B1 C1
A2 B2 C2 A3 B3 C3 .. .. . . Ai
..
.
Bi .. .
Ci .. . An
.. . Bn
x1
x2 x3 .. . = xi .. . xn
d1
d2 d3 .. . di .. . dn
Es f´acil observar lo siguiente: 1. Que utilizando como pivotes los sucesivos elementos de la diagonal desde i = 1 hasta i = n la matriz ampliada del sistema tridiagonal anterior es equivalente a la matriz:
B10 C1
d01
B20
C2
d02
B30 C3 .. .
d03 .. .
..
.
Bi0
Ci .. .
..
.
d0i .. .
Bn0 d0n donde los elementos Bi0 , d0i y, por tanto, la soluci´on es f´acilmente calculable utilizando el siguiente algoritmo: a) Definimos las variables Bi0 , d0i , i = 1, . . . , n por igualdades B10 = B1 ;
d01 = d1 ;
Bi0 = Bi − rCi−1 ;
d0i = di − rd0i−1 ;
i = 2, . . . , n
0 . donde r es la variable auxiliar Ai /Bi−1
´ltima componente de la soluci´on con b) Calculamos la u xn =
d0n Bn0
c) Calculamos las siguientes componentes en orden decreciente por la igualdad: xi =
d0i − Ci xi+1 Bi0
para i = n − 1, . . . , 2, 1. 2. Que lo anterior no es directamente aplicable si uno de los elementos diagonales, pongamos Bi0 es cero, en cuyo caso la ecuaci´on i-´esima tomar´ıa la forma Ci xi+1 = d0i . Pero entonces, la
´ E.4. METODO DIRECTO PARA SISTEMAS TRIDIAGONALES
39
inc´ognita xi+1 se obtendr´ıa de forma autom´atica por la igualdad xi+1 = d0i /Ci con lo que el sistema quedar´ıa reducido a un sistema de orden n − 1 en las variables {x1 , . . . , xd i+1 , . . . , xn } que es tridiagonal y por tanto resoluble mediante el anterior algoritmo o reducible a un sistema tridiagonal de orden estrictamente inferior y por tanto tambi´en resoluble mediante una ligera modificaci´on del algoritmo descrito en el punto anterior (ver la pr´actica eg). Pr´ actica e Ejecutar el comando type e_trid con el fin de estudiar el listado de dicho archivo y comprobar que corresponde a la implementaci´ on del algoritmo descrito en el punto 1 anterior. Vamos ahora a utilizar la funci´on citada para resolver el sistema tridiagonal 2.1x1 − x1
x2
=
− 2.2x2 + x2
x3
= 1.6
− 3.1x3 + x3
1
x4
= 2.1
− 2x4 = 1.2
para lo que ejecutamos el siguiente gui´on: b=[2.1 -2.2 -3.1 -2], a=[1 1 1], c=[-1 1 1], d=[1 1.6 2.1 1.2], e_trid(a,b,c,d) obteniendo como respuesta la soluci´on (−0.2927; −1.6146; −1.6595; −1.4297). Pr´ actica f Las matrices tridiagonales son un tipo particular de lo que se denominan matrices huecas o dispersas. Matlab cuenta con un comando especial para tratar dichas matrices, se trata del comando sparse. En esta pr´actica vemos su utilidad: Consideramos una matriz tridiagonal con la diagonal principal id´enticamente igual a 4, y las diagonales superior e inferior igual a 2, y estudiamos la variaci´ on en el n´ umero de operaciones a realizar al tratar la matriz como llena o como hueca: D=sparse(4*diag(ones(1,32))); DInf=sparse(2:32,1:31,2*ones(1,31),32,32); A=D+DInf+DInf.’; Allena=full(A); b=[1:32]’; flops(0), y=A\b; oper_hueca=flops, flops(0), oper_llena=flops, norm(y-x)
%Se define la diagonal, %los elementos infradiagonales, %la matriz del sistema, %la matriz llena A %los t´ erminos independientes x=Allena\b;
Obtenemos, oper_hueca=645, oper_llena=15400 y ky − xk = 1.7381e − 14 y observamos que para un sistema de orden 32 el n´ umero de operaciones se reduce en un 95 % respecto al sistema pleno. Podemos tambi´en comparar el comando sparse con la funci´on e_trid con el listado: D=sparse(4*diag(ones(1,32))); DInf=sparse(2:32,1:31,2*ones(1,31),32,32); A=D+DInf+DInf.’; b=[1:32]’; flops(0), y=A\b; oper_hueca=flops, a=2*ones(1,31); diago=4*ones(1,32); flops(0), x=e_trid(a,diago,a,b); oper_trid=flops, norm(y-x)
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
40
y comprobar que la funci´on e_trid necesita aproximadamente la mitad de operaciones. Pr´ actica puntuable g (0.4 puntos) Modificar el fichero e_trid de la forma necesaria para que pueda resolver sistemas tridiagonales generales y no s´olo aquellos en los que todos los elementos de la diagonal Bi0 son distintos de cero. La modificaci´on que hay que introducir es la especificada en el punto 2 de la enumeraci´ on del comienzo de esta secci´on.
E.5.
Estabilidad
Pr´ actica h Consideramos la matriz y los vectores siguientes (el ejemplo presente de mal condicionamiento o de inestabilidad es debido a R. S. Wilson): 10 7
8
5
6
7 A= 8 7
7
5 ; 6 10 9 5 9 10
32
23 ; b= 33 31
32.1
22.9 c= 33.1 30.9
Con el listado A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; det(A) b=[32 23 33 31].’; c=[32.1 22.9 33.1 30.9].’; x=A\b, y=A\c, obtenemos que A es una matriz sim´etrica de determinante 1 para la que las soluciones a los sistema Ax = b, Ay = c son respectivamente x = (1, 1, 1, 1) y y = (9.2; −12.6; 4.5; −1.1). Con ello observamos que a pesar de estar b y c separados u ´nicamente por una d´ecima en cada una de sus componentes se verifica que x e y est´ an separados en 13.6 unidades en una de sus componentes lo que significa que los errores relativos se han multiplicado, con la norma infinito, por un factor del orden de ¡4300!. Con la norma 2 los errores relativos se han multiplicado por un factor 1352. El condicionamiento de sistemas de ecuaciones lineales trata de explicar por qu´e se produce este fen´omeno. Las relaciones kbk = kAxk ≤ kAkkxk kx − yk ≤ kxk
kA−1 b
A−1 ck
− kbk/kAk
≤ kA−1 kkAk
(E.3) kb − ck kbk
(E.4)
nos dicen que el par´ametro que mide el factor de propagaci´on de los errores relativos de la soluci´on viene dado por el n´ umero cond(A) = kA−1 kkAk,
donde k · k es una norma inducida
a dicho n´ umero se le denomina n´ umero de condici´ on de la matriz invertible A. Matlab permite el c´alculo directo del n´ umero de condici´on, con la norma 2, de una matriz invertible A con el comando cond(A). En nuestro ejemplo resulta cond(A)=2984. Una estimaci´on del n´ umero 1/cond(A) se consigue con el comando rcond(A), este comando tiene la ventaja de que se computa con muchas menos operaciones.
E.5. ESTABILIDAD
41
La desigualdad (E.3) es una igualdad para los vectores x y b − c que verifican kA(x)k = kAkkxk,
kA−1 (b − c)k = kA−1 kkb − ck
(E.5)
De ello se deduce que si k · k es una norma inducida entonces el n´ umero de condici´on es una cota para el refuerzo de la propagaci´on de los errores relativos que siempre se alcanza. Es decir, siempre existen sistemas Ax = b, Ax = c de modo que kx − yk kb − ck = cond(A) · kxk kbk
(E.6)
´ltima afirmaci´on vamos a encontrar vectores x y b − Pr´ actica i Como ejemplo de nuestra u c que verifican la igualdad (E.5) anterior. Ello equivale, claramente, a encontrar vectores kA−1 (x)k v1 = x y v2 = b − c que sean m´aximos para las funciones f (x) = kAxk kxk y g(x) = kxk respectivamente. Denotemos entonces M (A) al conjunto de m´aximos de la funci´on f y m(A) a su conjunto de m´ınimos. Por ser A una matriz sim´etrica es ortogonalmente diagonalizable, de hecho ejecutando [V,D]=eig(A) comprobamos que A = V DV −1 para D = diag(0.0102, 0.8431, 3.8581, 30.2887) y V la matriz ortogonal dada por la igualdad
0.5016 −0.3017 −0.6149 0.5286 −0.8304 0.0933 −0.3963 0.3803 V = 0.2086 0.7603 0.2716 0.5520 −0.1237 −0.5676 0.6254 0.5209 La igualdad A = V DV −1 con V ortogonal implica: 1. kAk = kDk = m´axi |dii |, kA−1 k = kD−1 k =
1 , para D = (dij ). m´ıni |dii |
2. M (A) = V (M (D)). 3. M (A−1 ) = V (m(D)). m´axi |dii | 4. cond(A) = m´ıni |dii | Ahora, es claro que (0, 0, 0, 1)t ∈ M (D) y (1, 0, 0, 0)t ∈ m(D) de modo que v1 = V (:, 4) = (0.5286, 0.3803, 0.5520, 0.5209)t y v2 = 10−1 V (:, 1) = 10−1 (0.5016, −0.8304, 0.2086, −0.1237)t son elementos de M (A) y M (A−1 ) respectivamente. Todo lo anterior lo podemos llevar al ordenador con el listado: A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10] [V,D]=eig(A); x=V(:,4); v2=0.1*V(:,1); b=A*x; c=b-v2; rb=norm(b-c)/norm(b); y=A\c; rx=norm(x-y)/norm(x); factor=rx/rb; resultados=[x,y,b,c] fprintf(’....x....\t ....y....\t ....b....\t ....c....\t\n’) fprintf(’%f\t %f\t %f\t %f\t\n’,resultados.’), disp(’factor de propagaci´ on y n´ umero de condici´ on’) factor, condic=cond(A)
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
42
Los vectores b = (16.0096, 11.5176, 16.7180, 15.7781)t y c = (15.9595, 11.6007, 16.6971, 15.7905)t tienen la propiedad de que las soluciones x e y de los sistemas Ax = b y Ay = c verifican la igualdad (E.6). Pr´ actica j Veamos ahora c´omo se procede al c´alculo de vectores v1 ∈ M (A) y v2 ∈ M (A−1 ) cuando A es una matriz regular general que para fijar ideas supondremos que es la matriz de Vandermonde del vector (1, 2, 3.5, 5). Ahora A no es sim´etrica pero en cualquier caso s´ı lo es At A; en consecuencia, At A es ortogonalmente diagonalizable y de esto se deduce que existe una descomposici´on A = U DV t con D una matriz diagonal de valores no negativos y U, V ortogonales (ver [3], p´ag. 264) Se llama descomposici´ on en valores singulares de la matriz A a cualesquiera de las descomt posiciones A = U DV en que D es una matriz diagonal y U, V ortogonales. Es importante observar: La matriz D es u ´nica, salvo permutaci´ on en los elementos de la diagonal, ya que los 2 elementos diagonales de D coinciden con los valores propios de At A. Las afirmaciones 1–2) de la pr´actica anterior siguen siendo ciertas ahora para A = U DV t una descomposici´on en valores singulares. Adem´as del hecho de ser A−1 = V D−1 U t una descomposici´on de A−1 en valores singulares se deduce que la afirmaci´on 3) de la pr´actica anterior es cierta a condici´on de poner M (A−1 ) = U (m(D)). Matlab cuenta con un comando, svd, que permite el c´alculo de una descomposici´on en valores singulares. Lo aplicamos ahora al c´alculo de los vectores v1 , v2 con el listado: A=vander([1 2 3.5 5]); %A es la matriz de vandermonde de (1 2 3.5 5) [U,D,V]=svd(A), x=V(:,1); v2=0.1*U(:,4); b=A*x; c=b-v2; rb=norm(b-c)/norm(b); y=A\c; rx=norm(x-y)/norm(x); factor=rx/rb; resultados=[x,y,b,c] fprintf(’....x....\t ....y....\t ....b....\t ....c....\t\n’) fprintf(’%f\t %f\t %f\t %f\t\n’,resultados.’), display(’factor de propagaci´ on y n´ umero de condici´ on’) factor, condic=cond(A) Los vectores b = (1.2375, 8.7417, 44.5976, 127.5692)t y c = (1.1859, 8.8174, 44.5586, 127.5781)t tienen la propiedad de que las soluciones x e y de los sistemas Ax = b y Ay = verifican la igualdad (E.6).
E.6.
Ejercicios
Pr´ actica p Sean los n´ umeros complejos p = 4 + 1.5i, q = −2 − 0.5i, r = 6 + 2i y s = 8 + 2.5i. Se pide resolver el sistema
E.6. EJERCICIOS
43
px1 + qx2
= 5
qx1 + rx2 + qx3 = 0 qx2 + rx3 + qx4 = 0 qx3 + sx4 Pr´ actica q
= 0
Resolver los sistemas de ecuaciones siguientes: (−1 + i)z1 + (2 − 2i)z2 + 4z3 = 1
1)
(2 − 2i)z1 + 3z2 + iz3 3z1 + (1 − 5i)z2 + (1 + 2i)z3 x + 3x2 + 2x3 = 0 1 3)
Pr´ actica r
= 0 ;
2x1 − x3
3x1 + 3x2 + x3 = 0
= 4
x + 2x2 + 4x3 1 ;
2)
= −1
4x1 − x2 + x3 = 8 2x1 + 5x2 + 2x3 = 3
x1 + x2 + x4 2x + x − x + x 1 2 3 4 4) 4x1 − x2 − 3x3 + 2x4 3x1 + x2 − x3 + 7x4
;
= 2 = 1 = 0
;
= −3
Resolver el sistema tridiagonal
x1
+
x2
−x1 +
x2
= −
x3 x3
2
= −3
2x2 − 4x3 + 3x4
Pr´ actica s
= 11
−
=
2
x4
+ x5 =
0
x4
− x5 =
1
Considerar el siguiente sistema x2 + x3 x2 + 3x3 − x4 x −x 3 5 x4 + x5 −5x1 + x6 x1 + 2x2
= 6 = 1 = 3 = −1 = 1 = 0
Se pide resolverlo como sistema hueco y pleno comparando en cada caso el n´ umero de operaciones Pr´ actica puntuable t
(0.3 puntos) Consideramos el sistema de orden n definido por la igualdad
´ E. SISTEMAS DE ECUACIONES LINEALES. METODOS ´ LECCION DIRECTOS
44
2
−1
−1 2 −1 −1 2 −1 .. .. .. . . . −1 2
x1
x2 x3 = .. . xn
0
0 0 .. . 1
Se pide: 1. Comprobar que el vector x =
1 n+1 (1, 2, . . . , n)
es la u ´nica soluci´on del sistema anterior.
2. Resolver el sistema anterior para los casos n = 10, n = 20 y n = 30 utilizando el comando ‘\’, el comando sparse y la funci´on e_trid. Adem´as, se pide comparar la eficacia entre ambos m´etodos obteniendo para ello una tabla donde se refleje la distancia de la soluci´on exacta y el n´ umero de operaciones necesarias para obtener la soluci´on en cada caso y para cada uno de los m´etodos. 3. Resolver el sistema anterior para los casos n = 400, n = 600 y n = 800 utilizando la funci´on e_trid. Adem´as, se pide calcular la distancia con la soluci´on exacta y el n´ umero de operaciones necesarias para obtener la soluci´on. Pr´ actica u Sea H10 la matriz de Hilber de orden 10; se pide encontrar valores para b y c de modo que los sistemas H10 x = b y H10 y = c verifiquen: 1. kb − ck < 0.05. 2.
kx−yk kxk
= cond(H10 ) ·
kb−ck kbk
Pr´ actica v Para n = 3 hasta n = 20 se pide calcular el n´ umero de condici´on de la matriz de Hilbert Hn de orden n y representar en una gr´afica el resultado obtenido. Pr´ actica w Para n = 3, 7, 11, 15, 19, 23 resolver el sistema de ecuaciones Hn x = b para b = Hn (1, 1, . . . , 1) y Hn la matriz de Hilbert. La soluci´on debe ser el vector cuyas componentes son todas iguales a 1; sin embargo, comprobar c´omo para valores peque˜ nos de n la soluci´on en exacta y como se va perdiendo precisi´on para dimensiones m´as grandes. (Apartado puntuable con 0.2 puntos) Obtener una gr´afica que represente para cada n = 1, . . . , 20 la distancia con la soluci´on exacta. Pr´ actica puntuable x (0.4 puntos) Demostrar que una matriz tridiagonal con diagonal dominante (ver la pr´actica fu) puede escribirse como el producto de dos matrices del tipo (factorizaci´on de Crout):
L=
l11
0
...
...
l21
l22
0
...
0
l32
l33
...
...
...
...
...
0
0
0
lnn−1
0 0 0 ; ... lnn
U=
u11 u12
0
...
0
u22 u23 . . .
0
u33 . . .
...
...
...
...
0
0
0
0
0 0 0 ... unn