M` etodes Num` erics Sistemes Lineals.

0.1

Introducci´ on

Consideremos el sistema de ecuaciones lineales  a1,1 x1 + a1,2 x2 +     a2,1 x1 + a2,2 x2 + .. ..  . .    an,1 x1 + an,2 x2 +

··· ··· .. .

+ a1,m xm + a2,m xm .. .

···

+ an,m xm = bn ,

que en lenguaje matricial se escribe como Ax = b, siendo    x1 a1,1 a1,2 · · · a1,m  x2  a2,1 a2,2 · · · a2,m     A= . .. ..  , x =  .. ..   .  .. . . . xm an,1 an,2 · · · an,m



= b1 , = b2 , .. .



     y b=  

(1)

b1 b2 .. .

    

bn

donde A es la matriz de los coeficientes, x es vector de las incognitas y b es el vector independiente. La existencia de soluciones del sistema (1) esta resumida en la siguiente versi´on del Teorema de Rouch´e–Frobenius. • Si Rang(A) = Rang(A|b) (donde por A|b denotamos la matriz ampliada) entonces se dice que (1) es un sistema compatible, equivalentemente tiene soluci´on. – Si Rang(A) = m, entonces la soluci´on es u ´nica y el sistema se denomina sistema determinado. – Si Rang(A) < m, entonces existen infinitas soluciones y el sistema se denomina sistema indeterminado. • Si Rang(A) < Rang(A|b), el sistema (1) no tiene soluci´on y se denomina sistema incompatible. En este cap´ıtulo buscamos algoritmos y estrat´egias para el c´alculo de la soluci´on de un sistema de ecuaciones lineales (1) compatible y determinado que sean f´aciles de implementar en un ordenador y supongan un bajo coste en operaciones. De estos m´etodos los hay de dos tipos. • M´ etodos directos, que permiten obtener la soluci´on exacta (salvo errores de redondeo) en un n´ umero finito de pasos. A esta categor´ıa pertenecen los m´etodos de resoluci´on de matrices especiales (diagonales y triangulares), el m´etodo de Gauss, la factorizaci´on LU, factorizaci´ on QR. • M´ etodos iterativos, en los que se construye una sucesi´on de soluciones que tiende a la soluci´on exacta. A esta categor´ıa pertenecen los m´etodos de Jacobi, de Gauss-Seidel y de sobrerelajaci´ on. Como los sistemas a estudiar son compatibles y determinados, las matrices de coeficientes son cuadradas (n = m) y regulares, esto es, det(A) 6= 0.

1

0.2 0.2.1

M´ etodos directos Sistemas diagonales

Consideremos el caso especial A = (ai,j )1≤i,j≤n con  ai,j = 0 si i 6= j ai,i = ai que se denomina matriz diagonal y se denota por A = diag{a1 , a2 , . . . , an }. Como asumimos que la matriz es regular se sigue que ai 6= 0. • For i = 1, 2, . . . , n xi =

bi ai

Este algoritmo necesita n divisiones para resolver el problema. 0.2.2

Sustituci´ on hacia atras

Consideremos el caso A = (aQ i,j )1≤i,j≤n con ai,j = 0 si i > j que se denomina matriz triangular superior. Como det(A) = ni=1 ai,i , se sigue que ai,i 6= 0 para todo i. • xn =

bn an,n

• For i = n − 1, . . . , 1 P bi − nj=i+1 ai,j xj xi = ai,i Para el calculo de xk se necesitan n − k multiplicaciones, n − k restas y una division, que en total son n(n − 1)/2 multiplicaciones, n(n − 1)/2 restas y n divisiones, equivalentemente n2 operaciones. 0.2.3

Sustituci´ on hacia adelante

Consideremos el caso A = (a )1≤i,j≤n con ai,j = 0 si i < j que se denomina matriz triangular Qi,j n inferior. Como det(A) = i=1 ai,i , se sigue que ai,i 6= 0 para todo i. • x1 =

b1 a1,1

• For i = 2, . . . , n P bi − i−1 j=1 ai,j xj xi = ai,i Para el calculo de xk se necesitan k − 1 multiplicaciones, k − 1 restas y una division, que en total son n(n − 1)/2 multiplicaciones, n(n − 1)/2 restas y n divisiones, equivalentemente n2 operaciones.

2

0.2.4

M´ etodo de Gauss

El m´etodo de Gauss es el primero de los m´etodos que estudiaremos para la resoluci´on de un sistema de ecuaciones lineales con matriz de coeficientes general y no singular. El m´etodo de Gauss se basa en transformar el sistema de ecuaciones conservando invariante el conjunto de las soluciones, hasta obtener un sistema triangular superior. Las trasnformaciones que permiten hacer esto de una manera sencilla son: • Permutar la posici´ on de dos ecuaciones. • Multimplicar por un n´ umero real no nulo una de las ecuaciones. • Restar a una ecuaci´ on una combinaci´ on lineal del resto de las ecuaciones. Otra forma de enterder el m´etodo de Gauss es pensar que en cada paso despejamos una variable en una ecuaci´ on y la sustituimos en el resto de las ecuaciones. Partiendo de la matriz ampliada (A1 |b1 ) = (A|b) pasamos a (A2 |b2 ) seg´ un:  1   1  a1,1 a11,2 · · · a11,n b11 , a1,1 a11,2 · · · a11,n b11 , F11 → −−1−−→  0 a2 1 1 1 1  −−−−−−−−−− a22,n b22 ,  F21 →  = a1 2,2 · · · i,1 /a 1, 1  a2,1 a2,2 · · · a2,n b2 ,  mi,1   2 1 1 F = F − m F    i,1 1 .. .. .. .. .. .. .. .. ..  . i i .. .. i = 2, 3, . . . , n    . . . . . . . . . . .  Fn1 →

a1n,1 a1n,2 · · ·

a1n,n b1n ,

0

a2n,2 · · ·

a2n,n b2n ,

En este nuevo sistema tenemos que la variable x1 u ´nicamente interviene en la primera ecuaci´on. An´alogos razonamientos nos llevan del sistema (A2 |b2 ) al sistema (A3 |b3 ) que para el que la segunda variable x2 intervine u ´nicamente en las dos primeras ecuaciones. El procedimiento concluye al obtener en n − 1 pasos una matriz An triangular superior, que corresponde a un sistema de ecuaciones lineales con el mismo conjunto de soluciones que el original, pero ahora facil de resolver por el m´etodo de sustituci´on hacia atras. • For k = 1, 2, . . . , n − 1 For i = k + 1, k + 2, . . . , n ai,k mi,k = ak,k ai,k = 0 For j = k + 1, k + 2, . . . , n ai,j = ai,j − mi,k ak,j bi = bi − mi,k bk En el algoritmo anterior aparece una operaci´on que no siempre est´a definida, y a´ un cuando lo est´a, puede no ser num´ericamente aceptable, no referimos a la divisi´on. Recordamos que la divisi´on por n´ umeros peque¨ı¿ 21 os amplifica el error absoluto. Para evitar este problema existen dos variantes del m´etodo de Gauss conocidas como m´etodo de Gauss con pivote m´aximo por columnas y m´etodo de Gauss con pivote maximal. Pivotage m´ aximo por columnas. Para este m´etodo, en el paso k antes de calcular mi,k buscamos en la columna k la fila m tal que |akm,k | ≥ |aki,k | para k ≤ i ≤ n e intercambiamos la fila k por la fila m. Supongamos que akm,k = 0, en ese caso aki,k = 0 para k ≤ i ≤ n y det(A) = 0, en contradicci´ on con las hip´ otesis. Por tanto akm,k 6= 0. Como hemos visto, si det(A) = 0, entonces para alg´ un k el pivote akm,k podr´ıa ser cero. En tal caso, podemos dejar esta fila como est´a y continuar el proceso con la fila siguiente, no siendo necesaria la detenci´ on del algoritmo. Llegamos al siguiente resultado. 3

Teorema 1. El m´etodo de Gauss con pivotage m´ aximo por columnas se puede aplicar a toda matriz cuadrada. En el caso de matrices regulares, esta alternativa disminuye la cota del error en la soluci´ on en comparaci´ on con el m´etodo de Gauss general. • For k = 1, 2, . . . , n − 1 For i = k + 1, k + 2, . . . , n m=k For r = k + 1, k + 2, . . . , n If |ar,k | > |am,k | then m = r If am,k = 0 then Alerta! → det(A) = 0 mi,k = 1 else intercabiar filas k y m ai,k mi,k = ak,k ai,k = 0 For j = k + 1, k + 2, . . . , n ai,j = ai,j − mi,k ak,j bi = bi − mi,k bk Pivotage maximal. En esta variante del m´etodo de Gauss el pivote no se busca u ´nicamente entre los elementos de la columna, si no que se busca de entre todos los elementos de la submatriz Ak , formada por las filas k, . . . , n y las columnas k, . . . , n. Ejemplo: Resolver el siguiente sistema de ecuaciones usando aritmetica flotante con mantisa de 4 d´ıgitos y los algoritmos de Gauss sin y con pivote.  x1 +x3 = 1  x1 +0.0001x2 +2x3 = 2  x1 +x2 +x3 = 0 En ambos casos:    .1 × 101 0 .1 × 101 .1 × 101 .1 × 101 0 .1 × 101 .1 × 101  .1 × 101 .1 × 10−3 .2 × 101 .2 × 101  →  0 .1 × 10−3 .1 × 101 .1 × 101  .1 × 101 .1 × 101 .1 × 101 0 0 .1 × 101 0 −.1 × 101 

Sin pivotaje se sigue:   .1 × 101 0 .1 × 101 .1 × 101 x ¯1 = 0 −3 1 1   0 .1 × 10 .1 × 10 .1 × 10 ¯2 = 0 → x 0 0 −.1 × 105 −.1 × 105 x ¯3 = 1 Con pivotaje:  .1 × 101 0 .1 × 101 .1 × 101 x ˜1 = 0  0 .1 × 101 0 −.1 × 101  → x ˜2 = −1 0 0 .1 × 101 .1 × 101 x ˜3 = 1 

¯, x ˜ cometemos menos error. Si x es la soluci´on correcta Veamos en cual de las dos soluciones x ¯ entonces el error es la soluci´ de un sistema Ax = b y x ¯ es la aproximada, esto es, A¯ x = b, on 4



del sistema Aex = b − ¯b. En nuestro caso ex¯ = −10−4 , −1, 10−4 y ex˜ = −10−4 , 0, 10−4 . De donde se sigue que la soluci´ on x ˜ es m´ as aproximada. Contemos el n´ umero de operaciones que necesitamos para resolver un sistema lineal n×n por el m´etodo de eliminaci´ on Gaussiana sin pivotaje. Para para k = 1, . . . , n − 1 y i = k + 1, . . . , n hacemos 1 divisi´ on, n − k productos, n − k restas, 1 producto y 1 resta. Para k = 1, . . . , n − 1 tenemos n − k divisiones, (n − k)2 + (n − k) productos y (n − k)2 + (n − k) restas. En total: 3

n−1 X

k+2

k=1

n−1 X k=1

k2 = 3

n(n − 1)(2n − 1) 2n3 n(n − 1) +2 = + O(n2 ). 2 6 3

Supongamos que queremos utilizar el m´etodo de Gauss para el c´alculo de la inversa de A. Habriamos de resolver los n sistemas Ax = ek para k = 1, 2, . . . , n. Como cada uno necesita del orden de n3 operaciones, en total necesitamos del orden de n4 operaciones. Si miramos el algoritmo de Gauss vemos que cada nuevo sistema hemos de volver a calcular los multiplicadores y todos los pasos de la triangulaci´ on cuando en realidad la matriz de coeficientes no cambia. Parece pues necesario almacenar los multiplicadores y la matriz triangularizada una vez obtenida y no volver a calcularlos. El m´etodo de factorizaci´on LU trata de obtener beneficios de esta idea. 0.2.5

M´ etodo LU

En este apartado presentamos otra m´etodo para resolver el sistema de ecuaciones lineales Ax = b, donde la matriz A es regular. En este caso se expresa la matriz de los coeficientes A como producto de dos matrizes A = LU tal que L = (li,j )1≤i,j≤n es triangular inferior y U = (ui,j )1≤i,j≤n triangular superior. Veamos que siempre que podemos aplicar el m´etodo de Gauss, existe una posible factorizaci´ on para la matriz A. Llamemos U a la matriz triangular superior obtenida al aplicar el m´etodo de Gauss a A y sean mi,k los factores que construimos en el m´etodo de Gauss. Entonces:  1      F1 F11 F11  F1   F 2 = F 1 − m2,1 F 1    F22 2 1   2   2    F1   F 2 = F 1 − m3,1 F 1   F 3 = F 2 − m3,2 F 2  A = A1 =  3  → A2 =  3 3 1  → A3 =  3 3 2  → . . . U.  ..      .. ..  .      . . Fn1

Fn2 = Fn1 − mn,1 F11    1 0 0 ... 0 1 0 0 ...  m2,1 1 0 . . . 0   0 1 0 ...       . ..  ←  0 m ← 3,2 1  m3,1 0 1    ..  .. .. .. . . ..  .. .. . .  .  . . .  . . . . . mn,1 0 0 . . . 1 0 mn,2 0 . . . De donde se sigue que si li,i = 1 para 1 ≤ i ≤ n y li,j = mi,j A = LU. Sabiendo que la diagonal de L est´a formada por unos, almacenar en la misma matriz A. • For k = 1, 2, . . . , n − 1 For i = k + 1, k + 2, . . . , n ai,k ai,k = ak,k For j = k + 1, k + 2, . . . , n ai,j = ai,j − ai,k ak,j 5

Fn3 = Fn2 − mn,2 F22  0 0   ..  .   ← ... ..  .  1 con 1 ≤ j < i ≤ n, entonces las matrices L y U se pueden

En este algoritmo guardamos los multiplicadores en debajo de la diagonal y la matriz U en la diagonal y arriba de ella. Teorema 2.

(a) La factorizaci´ on de A existe si se le puede aplicar el m´etodo de Gauss.

(b) Si existe, la factorizaci´ on LU de una matriz A es u ´nica. Demostraci´ on: (a) Es consecuencia del algoritmo. (b) Consideremos las siguientes afirmaciones sobre matrices triangulares. • El producto de dos matrices triangulares inferiores (superiores) es una matriz triangular inferior (superior). • La inversa de una matriz triangular inferior (superior) es una matriz triangular inferior (superior). Supongamos que A = L1 U1 = L2 U2 . Entonces U1 U2−1 = L−1 1 L2 . Aplicando las afirmaciones anteriores se sigue que la matriz de la izquierda de la igualdad es triangular superior y la de la derecha es triangular inferior, por tanto ambas son diagonales. Como la de la derecha tienes unos en la diagonal ambas son la identidad, esto es U1 U2−1 = L−1 1 L2 = I, de donde se sigue que U1 = U2 y L1 = L2 .  Sabemos que la descomposici´ on LU esta basada en el m´etodo de Gauss. Sin embargo este no puede aplicarse directamente a todas las matrices, ni tan siquiera a las regulares. Para ello introduciamos una modificaci´ on de m´etodo de Gauss que consiste en permutar convenientemente las filas. Como la filas elegidas como pivote u ´nicamente se mueven una vez, si sabemos de antemano las transformaciones a realizar las podemos hacer antes y aplicar el m´etodo de Gauss sin pivotar. De esta forma tenemos el siguiente resultado. Teorema 3. Sea A una matriz regular. Existen P una matriz de permutaci´ on, L una matriz triangular inferior con unos en la diagonal y U una matriz triangular superior, de forma que P A = LU. Ejemplo: 

  1 2 3 −−−−→ 7 A =  4 5 6  (3, 2, 1)  4 7 8 1 1  7 8 −−−−→ (3, 1, 2)  1/7 6/7 4/7 3/7

Y por tanto, las matrices L, U  1 0  1/7 1 L= 4/7 1/2 0.2.6

  8 1 −−−−−−−→ = −4/7  5 6  ll2,1 3,2 = −1/7 2 3   1 −−−−−−−→ 20/7  l2,1 = −1/2  38/7

 7 8 1 4/7 3/7 38/7  1/7 6/7 20/7  7 8 1 1/7 6/7 20/7  4/7 1/2 4

y P son las siguientes.      0 7 8 1 0 0 1 0  , U =  0 6/7 20/7  y P =  1 0 0  1 0 0 4 0 1 0

M´ etodo QR

Antes de describir el m´etodo, introducimos algunas definiciones. Dada una matriz A = (ai,j )1≤i,j≤n definimos matriz transpuesta y la denotamos por AT a la matriz (˜ ai,j )1≤i,j≤n donde ai,j = a ˜j,i . T Diremos que A es una matriz sim´ etrica si A = A . En el caso de matrices reales diremos que

6

A es una matriz ortogonal si A−1 = AT . De la definici´on se sigue que las filas (columnas) de una matriz ortogonal son vectores ortonormales. El objeto del m´etodo QR es factorizar la matriz A en dos matrices una triangular superior que denotamos por R y otra ortogonal que denotamos por Q. En consecuencia, para resolver el sistema Ax = b es suficiente resolver el sistema triangular Rx = QT b. Las matrices ortogonales conservan la norma eucl´ıdea (kAxk2 = xT AT Ax = kxk2 ). El determinante de una matriz ortogonal es ±1 (ejercicio). El producto de matrices ortogonales es una matriz ortogonal (ejercicio). Por tanto, multiplicaremos la matriz A por matrices ortogonales convenientemente elegidas hasta conseguir la mariz R. Pasemos a introducir las matrices ortogonales de Householder. Matrices de Householder. Dado un vector u de Rn denotamos por uT = (u1 , . . . , un ) al vector transpuesto, por uT u al producto escalar, por uuT la matriz sim´etrica (ui uj ))1≤i,j≤n y por P (u) la matriz de Householder 2 P (u) = I − T uuT u u que tiene las siguientes propiedades: • Para todo c ∈ R se sigue que P (u)cu = cu −

2 2 uuT cu = cu − T c(uT u)u = −cu uT u u u

• Si u⊥ = {v ∈ Rn : uT v = 0} es el hiperplano ortogonal al vector u y v ∈ u⊥ , entonces P (u)v = v −

2 uT u

u(uT v) = v

• Si w ∈ Rn , entonces existe c ∈ R tal que w = cu + v. Pot tanto P (u)w = −cu + v. Por lo que P (u) es una simetr´ıa respecto del hiperplano u⊥ . Adem´as, uT w = cuT u, de donde c=

uT w . uT u

De esta forma podemos calcular P (u)w sin construir la matriz de Householder y con un m´ınimo coste de operaciones P (u)w = w −

2 uT u

(uT w)u.

• Como uuT es un matriz sim´etrica, entonces P (u) es una matriz sim´etrica. • Adem´as P (u)P (u) = P (u)(I −

2 2 2 uuT ) = P (u) − T P (u)uuT = P (u) + T uuT = I. uT u u u u u

Por tanto P (u) = P (u)−1 . Como adem´as P (u) es sim´etrica, se sigue que P (u) es ortogonal. • Para todo α 6= 0 se sigue que P (αu) = I −

2 α 2 uT u

7

α2 uuT = P (u).

• Si a, b ∈ Rn con a 6= b y kak2 = kbk2 , entonces P (a − b)a = b. Tomando u = a − b, como la norma eucl´ıdea de ambos vectores es igual si definimos v = 21 (a + b), entonces v ∈ u⊥ y a = 12 u+v. De donde P (a − b)a = − 21 u + v = b. Como las matrices ortogonales conservan la norma eucl´ıdea, no tiene sentido preguntarse que pasa si kak2 6= kbk2 . • Finalmente, si aT = (a1 , . . . , an ) ∈ Rn con ai 6= 0 para i = 2, . . . , n, entonces P (a + kak2 e1 )a = −kak2 e1 y P (a − kak2 e1 )a = kak2 e1 . Veamos c´ omo utilizar las matrices de Householder para triangularizar una matriz A = (ai,j )1≤i,j≤n . i Denotemos A0 = A y por aj = (a1,j , . . . , an,j )T el vector de la j–esima columna de A0 . ii Si a22,1 + . . . + a2n,1 = 0, entonces pasamos al punto . iii Si a22,1 + . . . + a2n,1 6= 0, entonces definimos s = ka1 k2 si a1,1 < 0 o s = −ka1 k2 si a1,1 > 0. iv Calculamos u = a1 − se1 . Por tanto P (u)a1 = se1 . v Para calcular P (u)A0 = (se1 , P (u)a2 , . . . , P (u)an ) hacemos α = P (u)aj 0.2.7

2 uT u

=

ka1 k

−1 −1 , = 1 su1 2 (ka k2 + |a1,1 |)

= aj − α(uT aj )u

M´ etodos iterativos

Los m´etodos iterativos para resolver un sistema de ecuaciones lineales Ax = b consisten en encontrar una matriz B, y dos vectores c, x0 ∈ Rn de forma que la sucesi´on {xk }∞ k=1 construida seg´ un xk+1 = Bxk + c sea convergente a la soluci´on del sistema original. Para ello necesitamos de los siguientes conceptos previos. Preliminares: • Dado un K–espacio vectorial V, donde K = R, C, llamamos norma vectorial a una aplicaci´ on kk : V → R tal que: kxk = 0 si y solo si x = 0, kλxk = |λ|kxk, kx + yk ≤ kxk + kyk. P ´ • Para todo p ∈ N la aplicaci´ on kxkp = ( ni=1 |xi |p )1/p es una norma vectorial. Unicamente es necesario probar la desigualdad triangular, el resto son inmediatas. Desigualdad de Minkowski: n X k=1

!1

p

|xk + yk |p

≤

n X k=1

!1

p

|xk |p

+

n X

!1

p

|yk |p

.

k=1

Elevamos a p ambos lados de la desigualdad. Desarrollando el lado izquierdo obten  Pn P P P p p−1 n r p−r . Desarrollando el lado dereemos k=1 |xk |p + nk=1 |yk |p + r=1 k=1 |xk | |yk | r   p−r r P Pn Pn Pp−1 p P p p cho obtenemos k=1 |xk | + k=1 |yk | + r=1 ( nk=1 |xk |p ) p ( nk=1 |yk |p ) p . De r 8

Pn r p−r ≤ donde resta por probar que para todo 1 ≤ r ≤ p tenemos que k=1 |xk | |yk | p−r r P P ( nk=1 |xk |p ) p ( nk=1 |yk |p ) p . Para ello probamos la desigualdad de H¨older !1/p n !1/q n n X X X |xk yk | ≤ |xk |p |yk |q k=1

k=1

k=1

para cualesquiera p > 1 y q > 1 tales que 1/p + 1/q = 1. Notemos que si la desigualdad de H¨older es cierta para vectores lo es para vectores λx Pn x e y,pentonces Pn tambien Pynµy. Por tanto p podemos considerar que k=1 |xk | = k=1 |yk | = 1. Veamos pues que k=1 |xk yk | ≤ 1. Ra Consideremos la funci´ on y = xp−1 y su inversa x = y q−1 . Si S1 = 0 xp−1 dx = ap /p y Rb S2 = 0 y q−1 dy = bq /q, mirando la representaci´ on gr´afica inmediato que P de S1 y S2 esP P ab ≤ S1 + S2 = ap /p + bq /q. De donde nk=1 |xk yk | ≤ 1/p nk=1 |xk |p + 1/q nk=1 |yk |q = 1. A las normas p se las denomina normas H¨ older. El caso particular p = 2 se denomina norma euclideana. • La aplicaci´ on kxk∞ = maxi=1,...,n {|xi |} es una norma vectorial y se denomina norma del m´ aximo. • Diremos que una aplicaci´ on kk : Mn (K) → R es una norma matricial, si es una norma vectorial y adem´ as satisface que kABk ≤ kAkkBk. • Una norma matricial es una norma consistente con la norma vectorial del espacio si kAxk ≤ kAkkxk. Aunque la notaci´on puede inducir a confusi´on las normas de la desigualdad anterior significan cosas diferentes, y es bueno tenerlo presente. La aplicaci´ on kAk = max{kAxk/kxk} es una norma matricial consistente con la norma vectorial que se denomina norma subordinada. • Si A ∈ Mn (Rn ) diremos que un vector v 6= 0 es un vector propio si existe un valor λ ∈ R tal que Av = λv. Al valor λ se le llama valor propio asociado al vector propio v. Denominamos polinomio caracter´ıstico de A al polinomio det(A − λI). Los valores propios son las raices del polinomio caracter´ıstico. • La matriz A y su transpuesta AT tienen los mismos valores propios. • Una matriz es inversible si y solo si tiene todos sus valores propios diferentes de cero. Como det(A−1 − λI) = (−λ)n det(A−1 ) det(A − λ1 I), los valores propios de la inversa son los inversos de los valores propios de A. Teorema 4 (Teorema de Gerschgorin). Los valores propios de una matriz A = (ai,j )1≤i,j≤n est¨ı¿ 12 n contenidos dentro del plano complejo en la uni´ on F de los discos       n   X |ai,j | para todo i = 1, 2, . . . , n, Fi = z : |z − ai,i | ≤     j =1   j 6= i

y tambien en la uni´ on C de los discos       n   X Cj = z : |z − aj,j | ≤ |ai,j | para todo j = 1, 2, . . . , n.     i=1   i 6= j

Adem´ as, si una componente conexa de F o de C esta formada por k discos, entonces esta componente contiene ex´ actamente k valores propios teniendo en cuenta la multiplicidad. 9

Demostraci´ on: Sea v un vector propio de A y λ el valor propio asociado, esto es (A − |vi | . + (ak,k − λI)v = 0. Sea |vk | = kvk∞ 6= 0, entonces |v | ≤ 1 para todo i. De ak,1 v1 + . .P k P vi λ)vk +. . .+ak,n vn = 0 se sigue que ak,k −λ = i6=k ai,k vk , de donde |ak,k −λ| ≤ i6=k |ai,k |. En consecuencia λ ∈ Fk ⊂ F. Para ver que λ ∈ C tomamos la matriz transpuesta. m n Supongamos que ∪m i=1 Fi es una componente conexa de F, esto es ∪i=1 Fi ∩ ∪i=m+1 Fi = ∅ y veamos que A tiene m valores propios en la componente conexa. Descomponemos A como suma de la matriz diagonal y del resto, A = D + R. Definimos la familia de matrices At = D + tR en el intervalo [0, 1] de forma que A0 = D y A1 = A. Los valores propios de A0 son los centros de los discos Fi por lo que A0 tiene exactamente m valores propios en la componente conexa. Como los valores propios de At varian continuamente con t, existen exactamente m en la componente conexa. Lo que prueba el teorema. 

• Si λi para i = 1, 2, . . . , m son los valores propios de la matriz A llamamos radio espectral de A y lo denotamos por ρ(A) al maxi {|λi |}. • Las normas matriciales subordinadas a las normas vectoriales kk1 , kk2 y kk∞ son: P kAk1 = maxj=1,...,n { ni=1 |ai,j |} , P kAk1 = maxkxk1 =1 {kAxk1 } ≥ kAek k1 = ni=1 |ai,k | parta todo k por tanto, es mayor que el m¨ı¿ 21 ximo de la suma por columnas, SupongamosPque el m´ aximo de la sumaPpor columnas se alcanza en la columna n n |ai,j |} . Para j0 , esto es i=1 |ai,j0 | = maxj=1,...,n Pn vector x de Pncualquier Pn { Pi=1 n ( a x | ≤ | norma 1 se tiene que kAxk = i,j j 1 i=1 |ai,j |) |xj | ≤ j=1 j=1 i=1 Pn Pn Pn |a |) . |a |) |x | ≤ ( ( j i=1 i,j0 i=1 i,j0 j=1 p kAk2 = ρ(AT A). ∗ En primer lugar definimos el cociente de Rayleigh (1842-1919) de una matriz M como la aplicaci¨ı¿ 21 n RM : V \ {0} → C definida seg´ un RM (v) =

vT M v . vT v

∗ Si v es una valor propio de M de valor propio λ, entonces RM (v) = λ. ∗ Si M es sim´etrica, tiene todos los valores propios reales (λ1 ≤ . . . ≤ λn ) y existe una base de vectores propios ortogonal {p1 , . . . pn }(ver ”Introduction `a l’analyse num´erique matricielle et `a l’optimisation” [Philippe G. Ciarlet] Dunod). Si U es la matriz con columnas los vectores propios, entonces U T M U es una matriz diagonal D con los valores propios en la diagonal. Por tanto si v = U w, entonces RM (v) = Pk 2 Pk i=1 λi |αi | RD (w). Tomando v = i=1 αi pi , se sigue que RM (v) = Pk . De donde 2 i=1

|αi |

es sencillo ver que λk = max{RM (v) : v ∈ Vk }, siendo Vk =< p1 , . . . pk > . ∗ De lo anterior se deduce que cuando M es sim´etrica, ρ(M ) = max{RM (v) : v ∈ T T V }. Por otra parte, kAk22 = sup v vAT vAv = max{RAT A (v) : v ∈ V } = ρ(AT A). nP o n kAk∞ = maxi=1,...,n |a | . j=1 i,j Pn Supongamos que el m´ a ximo por filas se da en la fila i , esto es 0 j=1 |ai0 ,j | = P maxi { nj=1 |ai,j |}. Sea v tal que P vj = 1 si ai0 ,j >P 0 y vj = −1 en otro caso. Por tanto, kAk∞ ≥ kAvk∞ = maxi {| nj=1 ai,j vj |} ≥ nj=1 |ai0 ,j |. Pn Para cualquier vector x de norma 1 se sigue que kAxk = max {| ∞ i j=1 ai,j xj |} ≤ Pn Pn maxi { j=1 |ai,j ||xj |} ≤ maxi { j=1 |ai,j |}.

10

• El radio espectral es menor que cualquier norma matricial de la matriz: ρ(A) ≤ kAk. En efecto, si λ es el valor propio tal que |λ| = ρ(A) y v es el vector propio asociado, entonces para cualquier vector w tal que la matriz vwT no tenga norma nula tenemos que ρ(A)kvwT k = kλvwT k = kAvwT k ≤ kAkkvwT k. De donde se sigue el resultado• Dada una matriz A y un valor ε > 0 existe una norma matricial subordinada a una norma vectorial tal que kAk ≤ ρ(A) + ε. Una demostraci´on se puede encontrar en ”Introduction `a l’analyse ....” [Philippe G. Ciarlet] Dunod. Criterios de convergencia de los m´ etodos iterativos. Para ver que la sucesi´ on {xk } generada por un m´etodo iterativo es convergente veamos que es de Cauchy (Rn es un espacio de Banach) kxk+1 − xk k ≤ kB(xk − xk−1 )k ≤ kBkkxk − xk−1 k ≤ kBkk kx1 − x0 k. Luego una condici´ on suficiente para la convergencia es que kBk < 1 para alguna norma matricial consistente con la norma vectorial. En particular, si ρ(B) < 1 entonces el m´etodo es convergente para alguna norma; mientras que si ρ(A) ≥ 1 el m´etodo es divergente en cualquier norma. Finalmente, si kBk < 1 y la sucesi´on de iterados converge a z, entonces kz − xk k ≤ kz − xk+n k +

n−1 X

kxk+1+i − xk+i k ≤ . . . ≤

i=0

kBkk kx1 − x0 k 1 − kBk

M´ etodo de Jacobi. Si los elementos de la diagonal de la matriz A son diferentes de 0 podemos descomponer A como A = D(L + I + U ) donde D es la matriz de la diagonal de A, L es una matriz triangular inferior con ceros en la diagonal y U es una matriz triangular superior con ceros en la diagonal. El sistema Ax = b es por tanto equivalente al sistema x = −(L + U )x + D−1 b. El m´etodo de Jacobi consiste en tomar BJ = −(L + U ) y cJ = D−1 b, que en lenguaje algor´ıtmico se puede escribir como:   1 − kBJ k ln ε 1 kx − x0 k • Sea n ≥ ln kBJ k • For k = 1, 2, . . . , n  P 1  xk+1 bi − j6=i ai,j xkj = i ai,i M´ etodo de Gauss–Seidel. Al igual que el m´etodo de Jacobi, se aplica a sistemas lineales donde la matriz A no tiene ceros en la diagonal. Notar que con transformaciones elementales podemos eliminar los ceros de la diagonal sin alterar el conjunto de soluciones del sistema. Partiendo de A = D(L+I +U ) y de que D y L + I son matrices invertibles, el sistema es equivalente al sistema x = −(L + I)−1 U x + (L+I)−1 D−1 b. El m´etodo de Gauss–Seidel consiste pues en tomar BGS = −(L+I)−1 U y cGS = (L+I)−1 D−1 b. Sin embargo consideraremos la siguiente versi´on: xk+1 = −Lxk+1 −U xk +D−1 b, donde se aprecia que el m´etodo de Gauss–Seidel aprovecha la informaci¨ı¿ 21 n que va opteniendo en el paso k + 1. Este m´etodo en leguaje algoritmico se puede escribir seg´ un:   1 − kBGS k ln ε 1 kx − x0 k • Sea n ≥ ln kBGS k • For k = 1, 2, . . . , n  Pi−1 Pn 1  k+1 k xk+1 = b − a x − a x i i j=1 i,j j j=i+1 i,j j ai,i 11

M´ etodo de sobrerelajaci´ on. Se trata de una familia de m´etodos que generaliza el m´etodo de Gauss–Seidel. De nuevo exige que A no tenga ceros en la diagonal, y partimos del sistema equivalente x = x−ω(Lx+(I + U )x − D−1 x), que podemos escribir como x = Bω x + cω donde Bω = (I + ωL)−1 [(1 − ω)I − ωU ] y cω = ω(I + ωL)−1 D−1 b. A partir de aqui podemos construir el siguiente algoritmo:   1 − kBω k ln ε 1 kx − x0 k • Sea n ≥ ln kBω k • For k = 1, 2, . . . , n xk+1 = xk+1 + i i

 Pi−1 Pn ω  k bi − j=1 ai,j xk+1 − a x j j=i i,j j ai,i

Converg´ encia de Jacobi, Gauss–Seidel y sobrerelajaci´ on. Consideraciones sobre la convergencia. • El m´etodo de Jacobi, Gauss–Seidel o sobrerelajaci´on converge si el r´adio espectral de BJ , BGS o Bω es menor que 1. • Si A es estr´ıctamente diagonal, los m´etodos de Jacobi y Gauss–Seidel convergen. • Si ω 6∈ (0, 2) el r´ adio espectral de la matriz de sobrerelajaci´on es mayor que 1. Q ) Como det(Bω ) = det((1−ω)I−ωU = (1 − ω)n , entonces ni=1 λi = (1 − ω)n . Por otra det(I+ωL) parte ρ(Bω ) ≥ |λi | para todo i. En consecuencia, ρ(Bω ) ≥ |1 − ω|. • Si A sim´etrica, con valores positivos en la diagonal y definida positiva, el m´etodo de sobrerelajaci´ on converge siempre que ω ∈ (0, 2). Consideraciones sobre la velocidad de convergencia. • Para matrices A definidas positivas, tridiagonal por bloques y tales que los bloques de la diagonal sean matrices diagonales tenemos que ρ(BGS ) = ρ(BJ )2 . • Si adem´ as los valores propios de BJ est¨ı¿ 12 n en el intervalo (−1, 1), entonces ρ(Bω ) alcanza su m´ınimo en 2 p 1