Métodos iterativos para sistemas lineales provenientes de problemas de punto de silla
Abstract
Grado en Matemáticas
Full text
FacultaddeCiencias Trabajo Fin de Grado Grado en Matemáticas Curso 2018-2019 Métodositerativosparasistemaslinealesprovenientesdeproblemas depuntodesilla Autor:AlbaCrespoMartínez Tutor:LuisMªAbiaLlera Dpto.MatemáticaAplicada
ii
Prefacio El objetivo de este trabajo Fin de Grado es generalizar algunos resultados sobre el m´etodo iterativo de sobrerrelajaci´on (SOR) de par´ametro ω, para sistemas lineales aumentados. El inter´es de estos sistemas estar´ıa ya justificado porque surgen de forma natural en la discretizaci´on del problema de Stokes en mec´anica de fluidos, de gran inter´es para las aplicaciones. Este tipo de sistemas aparecen tambi´en en problemas punto de silla en optimizaci´on. El punto de partida son los resultados b´asicos de convergencia de los m´etodos iterativos basados en escisiones de la matriz Aestudiados en el Grado de Matem´aticas. En el cap´ıtulo primero ampliamos la teor´ıa de convergencia del m´etodo SOR para sistemas lineales generales con vistas a proporcionar resultados relativos a la elecci´on del par´ametro ´optimo. Esto nos lleva a introducir clases de matrices para las que es factible esta optimizaci´on del par´ametro. El cap´ıtulo 2 aporta extensiones del m´etodo SOR para tratar los sistemas aumentados. Constatamos en los art´ıculos manejados una diversidad de extensiones del m´etodo SOR para sistemas aumentados, cada una con su propia teor´ıa de convergencia y de optimizaci´on del par´ametro. Hemos seleccionado dos an´alisis: uno se refiere al m´etodo SOR con un par´ametro (con y sin preacondicionamiento), y el segundo al que se denomina GSOR (SOR generalizado), que depende de dos par´ametros. El cap´ıtulo 3 considera m´etodos iterativos en los que la escisi´on de la matriz del sistema es diferente a la de los m´etodos SOR. Se basan en escisiones de la matriz del sistema en una parte Herm´ıtica y otra antiherm´ıtica, con desplazamiento de la diagonal. Desarrollamos una teor´ıa de convergencia para sistemas generales, se˜nalando las modificaciones que permiten aplicar los resultados a sistemas aumentados. El ´enfasis de la memoria es en el an´alisis antes que en un estudio comparativo de la eficiencia de las distintas propuestas. En Valladolid, a 15 de julio de 2019 iii
iv
´ Indice general 1. M´etodos iterativos para sistemas lineales 1 1.1. Introducci´on............................ 1 1.2. Teor´ıa general de convergencia . . . . . . . . . . . . . . . . . 1 1.3. El m´etodo SOR de par´ametro ω................ 4 1.4. Convergencia del m´etodo SOR . . . . . . . . . . . . . . . . . 6 1.5. Optimizaci´on del par´ametro para el m´etodo SOR . . . . . . . 9 2. M´etodos SOR para sistemas aumentados 19 2.1. Los sistemas aumentados y su inter´es . . . . . . . . . . . . . . 19 2.2. El m´etodo SOR para sistemas aumentados . . . . . . . . . . . 20 2.3. Convergencia del m´etodo SOR para sistemas aumentados . . 21 2.4. Optimizaci´on del par´ametro del m´etodo SOR para sistemas aumentados............................ 24 2.5. El m´etodo SOR generalizado para sistemas aumentados . . . 27 2.6. Teoremas de convergencia para el m´etodo SOR generalizado . 29 2.7. Optimizaci´on de los par´ametros en el m´etodo SOR generalizado 31 3. M´etodos HSS para sistemas generales y sistemas aumentados 41 3.1. Introducci´on............................ 41 3.2. El m´etodo HSS para sistemas aumentados . . . . . . . . . . . 42 3.3. Convergencia del m´etodo HSS . . . . . . . . . . . . . . . . . . 43 3.4. Optimizaci´on del par´ametro para el m´etodo HSS . . . . . . . 48 v
vi
Cap´ıtulo 1 M´etodos iterativos para sistemas lineales 1.1. Introducci´on Sean A∈Cn×nuna matriz cuadrada compleja y b∈Cnun vector complejo dado. Consideramos el sistema lineal Ax =b. (1.1) La soluci´on x∈Cnexiste y es ´unica si y solo si Aes no singular. Entonces, x=A−1b. Supodremos a partir de aqu´ı que Aes no singular y todas las entradas de la diagonal son n´umeros complejos distintos del cero. Estudiaremos ahora distintos m´etodos iterativos. Estos m´etodos iterativos no calculan, en general, la soluci´on exacta del sistema, sino que a partir de una aproximaci´on inicial x(0) van obteniendo sucesivas aproximaciones a la soluci´on hasta que se alcanza la precisi´on deseada. La convergencia de las aproximaciones as´ı obtenidas a la soluci´on exacta puede garantizarse bajo diferentes condiciones que estudiaremos. El n´umero de iteraciones necesarias para alcanzar una precisi´on prefijada depender´a, en general, del m´etodo utilizado. Aunque el punto de partida de este trabajo es el an´alisis del m´etodo SOR y la optimizaci´on de su par´ametro, en la siguiente secci´on revisamos sucintamente la teor´ıa general de convergencia de m´etodos iterativos para (1.1) ligados a escisiones A=M−Nde la matriz, para lo que nos hemos basado principalmente en el libro de Varga [6] y el de Young [7]. 1.2. Teor´ıa general de convergencia Consideramos la siguiente escisi´on para la matriz A A=M−N 1
y podemos reescribir (1.1) como Mx =Nx +b(1.2) Entonces, para un iterante inicial x(0), el m´etodo iterativo genera una sucesi´on definida recurrentemente por Mx(m+1) =Nx(m)+b, m = 0,1,2, ... (1.3) Para que el m´etodo sea ´util en la pr´actica necesitaremos que resolver el sistema lineal (1.3) con matriz Msea mucho m´as barato computacionalmente que resolver el sistema inicial (1.1) de matriz A. Adem´as para que (1.3) defina el iterante x(m+1) de forma ´unica es necesario que Msea regular. En ese caso, podemos escribir (1.3) como x(m+1) =Gx(m)+d, m = 0,1,2, ... (1.4) donde G=M−1Nse llama la matriz de iteraci´on del m´etodo iterativo y d=M−1b. En particular, si x∗es la soluci´on de (1.1), restando de (1.4) la identidad x∗=Gx∗+d, se obtiene la siguiente recurrencia para el error e(m)=x(m)−x∗, e(m+1) =Ge(m), m = 0,1,2, . . . Obviamente, e(m)=Gme(0), y se tiene que la sucesi´on {e(m)}∞ m=0 converge a cero, cualquiera que sea e(0), si y solo si la matriz Gmconverge a cero. 1.1 Definici´on. Sea G∈Cn×n. Se dice que Ges convergente (a cero) si la secuencia de matrices G, G2, G3, . . .converge a la matriz nula O, y se dice que Ges divergente en caso contrario. 1.1 Teorema. Sea G∈Cn×n. Entonces, G es convergente si y solo si el radio espectral (recordemos que el radio espectral es el mayor de entre los valores absolutos de los autovalores de G,ρ(G) := m´ax i(|λi|)) satisface ρ(G)<1. Demostraci´on. Recordemos que para una matriz Gdada, existe una matriz P∈Cn×ncon la cual podemos reducir Ga su forma can´onica de Jordan, es decir, PGP−1=J= J1 J2 ... Jr ,(1.5) donde cada matriz Ji∈Cni×nies de la forma 2
Ji= λi1 λi1 λi1 ...... 1 λi . Como cada matriz Jies triangular superior, tambi´en lo es J. Luego, {λi}r i=1 incluye todos los autovalores de GyJ(que tienen los mismos autovalores por ser matrices semejantes). De (1.5) se tiene Jm= Jm 1 Jm 2... Jm n , m ≥1 (1.6) y por la forma que tienen las Ji, J2 i= λ2 i2λi1 λ2 i2λi1 ......... λ2 i2λi1 λ2 i2λi λ2 i . En general, Jm i=d(m) k,l (i)para 1 ≤k, l ≤ni, donde d(m) k,l (i) = 0 si l < k m l−kλm−l+k isi k≤l≤m´ın(ni, m +k) 0 si m+k < l ≤ni (1.7) es decir, Jm i= λm im 1λm−1 i. . . m ni−1λm−(ni−1) i 0λm i. . . m ni−2λm−(ni−2) i . . ..... . . 0 0 . . . λm i . Supongamos entonces que Ges convergente, entonces Gm→Ocuando m→ ∞. Como Jm=PGmP−1, entonces Jm→Ocuando m→ ∞. Entonces, cada Jm i→Ocuando m→ ∞, luego las entradas diagonales λide Ji 3
donde los bloques diagonales Ai,i para 1 ≤i≤Nson matrices cuadradas y no vac´ıas. Cada matriz Ai,i es de dimensi´on ni×nicon ni≥1. Asumimos que cada submatriz diagonal es no singular, de forma que D= A1,1 A2,2 ... AN,N (1.21) es no singular. La matriz B∈Cn×n B:= −D−1A+I(1.22) es la matriz de Jacobi por bloques asociada a la partici´on (1.20) de A. 1.2 Definici´on. Sea A∈Cn×n(no necesariamente no negativa o irreducible). Se dice que Aes d´ebilmente c´ıclica de ´ındice k > 1 si existe una matriz de permutaci´on Ptal que PAPTes de la forma PAPT= 0 0 . . . 0A1,k A2,10. . . 0 0 0A3,2. . . 0 0 . . ..... . . 0 0 . . . Ak,k−10 ,(1.23) donde las submatrices diagonales nulas son cuadradas. La teor´ıa de matrices no negativas de Perron-Frobenius establece la invariancia frente a rotaciones de ´angulos 2π/k de matrices no negatives irreducibles que tienen la forma (1.23). Esta invariancia se puede establecer directamente a partir de la forma (1.23) para matrices generales, un resultado que se debe a Romanovsky [6] 1.6 Teorema. (1936, [6]) Sea A= (aij)una matriz n×nd´ebilmente k-c´ıclica. Entonces φ(t) = det (tI −A) = tm r Y i=1 (tk−σk i), donde m+rk =n. 1.3 Definici´on. Si la matriz de Jacobi por bloques, B, de (1.22) para la matriz Ade (1.20) es d´ebilmente c´ıclica de ´ındice p≥2, entonces Aes p-c´ıclica relativa a la partici´on de (1.20). 10
Las siguientes matrices tienen un especial inter´es: A1= A1,10 0 . . . 0A1,p A2,1A2,20. . . 0 0 . . ........ . . . . ........ . . . . ........ . . 0 0 0 . . . Ap,p−1Ap,p , p ≥2 (1.24) A2= A1,1A1,2 A2,1A2,2A2,3 ......... ......AN−1,N AN,N−1AN,N .(1.25) A las matrices de la forma de A2se les llama matrices tridiagonales por bloques. De acuerdo con (1.22), las matrices A1yA2dan lugar, respectivamente, a las siguientes matrices de Jacobi: B1= 0 0 0 . . . 0B1,p B2,10 0 . . . 0 0 0B3,20. . . 0 0 . . ........ . . . . ........ . . 0 0 0 . . . Bp,p−10 , p ≥2 (1.26) B2= 0B1,2 B2,10B2,3 ......... ......BN−1,N BN,N−10 .(1.27) Observemos que la matriz B1es una matriz d´ebilmente c´ıclica de ´ındice p. Por tanto, por definici´on, se tiene que la matriz A1es p-c´ıclica. Por otro lado, los bloques de la matriz B2se pueden permutar para ver que B2es una matriz d´ebilmente c´ıclica de ´ındice 2. Por tanto, la matriz A2es una matriz 2-c´ıclica, es decir, las matrices tridiagonales por bloques son un caso importante de matrices 2-c´ıclicas. 11
1.4 Definici´on. Si la matriz Ade (1.20) es p-c´ıclica, entonces se dice que A es consistentemente ordenada si todos los autovalores de la matriz B(α) := αL +α−(p−1)U, que deriva de la matriz de Jacobi B=L+U, son independientes de α, para α6= 0. En ese caso tambi´en diremos que Bes consistentemente ordenada. En caso contrario diremos que AyBson inconsistentemente ordenadas. Veamos que las matrices A1de (1.24) y A2de (1.25) son matrices consistentemente ordenadas y, por tanto, lo son tambi´en B1de (1.26) y B2de (1.27). Para B1consideremos B1(α) := 0 0 0 . . . 0α−(p−1)B1,p αB2,10 0 . . . 0 0 0αB3,20. . . 0 0 . . ........ . . . . ........ . . 0 0 0 . . . αBp,p−10 y es f´acil comprobar que Bp 1(α) = Bp 1, para cualquier α6= 0. Por tanto, los autovalores de B1(α) son independientes de α. As´ı pues, las matrices A1y B1son consistentemente ordenadas. Centr´emonos ahora en la matriz B2y consideremos los autovalores λde la matriz B2(α), es decir, B2(α)x=λx, donde x6= 0. De acuerdo con la partici´on (1.27), podemos tomar una partici´on tambi´en para x, de forma que se tiene αBj,j−1Xj−1+1 αBj,j+1Xj+1 =λXj,1≤j≤N, donde tanto B1,0como BN,N+1 se definen como matrices nulas. Definamos entonces Zj:= (1/αj−1)Xj, para 1 ≤j≤N, de forma que la ecuaci´on anterior es Bj,j−1Zj−1+Bj,j+1Zj+1 =λZj,1≤j≤N. Luego todo autovalor de B2(α), para α6= 0 es tambi´en un autovalor de B2. Luego, B2yA2son ambas matrices consistentemente ordenadas. Esto es, cualquier matriz tridiagonal por bloques es una matriz consistentemente ordenada. 12
1.5 Definici´on. Se dice que una matriz Ade orden ntiene la propiedad A de Young si existen dos subconjuntos disjuntos S1yS2que particionan W={1,2, . . . , n}y tales que si i6=jy si o bien aij 6= 0 o aji 6= 0 entonces i∈S1yj∈S2o bien i∈S2yj∈S1. Cuando para un par (i, j) se tiene que o bien aij 6= 0 o bien aji 6= 0 se dice que los ´ındices iyjest´an asociados. En esta definici´on puede suceder que S1oS2sea vac´ıo, en cuyo caso la matriz Aes diagonal. La definici´on anterior es equivalente a la siguiente: 1.1 Proposici´on. Una matriz Atiene la propiedad A si y s´olo si Aes una matriz diagonal o existe una matriz de permutaci´on Ptal que PAPTtiene la forma PTAP=D1B C D2(1.28) donde D1yD2son matrices diagonales cuadradas y no necesariamente del mismo orden. Demostraci´on. Si Atiene la propiedad A, sean S1yS2los conjuntos que especifican la definici´on. Si S1oS2son vac´ıos entonces Aes diagonal. En otro caso, denotemos con s1 y s2 el n´umero de elementos de S1yS2, respectivamente, y denotemos los ´ındices de Skpor ik 1< ik 2< . . . < ik sk,k= 1,2. Construyamos la permutaci´on σdefinida mediante σ(i1 1) = 1, σ(i1 2) = 2, . . . , σ(i1 s1) = s1, σ(i2 1) = s1+ 1, σ(i2 2) = s1+ 2, . . . , σ(i2 s2) = s1+s2, Demostraremos que si Pes la matriz de permutaci´on asociada a σ, entonces A0=PTAP tiene la forma (1.28). Sea T1={1,2, . . . , s1}yT2={s1+ 1, . . . , s1+s2}. Basta demostrar que si a0 ij 6= 0 e i6=jentonces j∈T2si i∈T1, y j∈T1si i∈T2. Si a0 ij 6= 0 e i6=j, entonces aσ−1(i),σ−1(j)6= 0 y por tanto, σ−1(i) y σ−1(j) est´an asociados. Puesto que Atiene la propiedad A y puesto que σ−1(i)6=σ−1(j) se tiene que o bien σ−1(i)∈S1yσ−1(j)∈S2 oσ−1(i)∈S2yσ−1(j)∈S1. Por la construcci´on de σentonces o bien i=σ(σ−1(i)) ∈T1yj=σ(σ−1(j)) ∈T2o bien i∈T2yj∈T1. Por tanto A0tiene la forma que postula el teorema si T1yT2son ambos no vac´ıos o A0es diagonal. Rec´ıprocamente, si para alguna permutaci´on P, la matriz A0=PTAP tiene la forma (1.28), entonces A0tiene la propiedad A puesto que A0es una matriz tridiagonal por bloques y, por tanto, est´a consistentemente ordenada. 1.7 Teorema. Sea Ade la forma (1.20) una matriz p-c´ıclica consistentemente ordenada cuyas submatrices diagonales Ai,i para 1≤i≤Nson regulares. 13
Si ω6= 0, si λes un autovalor no nulo de la matriz Lω= (I−ωE)−1{ωF + (1 −ω)I}y si µsatisface (λ+ω−1)p=λp−1ωpµp,(1.29) entonces µes un autovalor de la matriz por bloques de Jacobi B=L+U. Rec´ıprocamente, si µes un autovalor de Byλsatisface (1.29), entonces λ es un autovalor de Lω. Demostraci´on. Los autovalores de Lωson las ra´ıces del polinomio caracter´ıstico det(λI −Lω)=0.(1.30) Como I−ωL es no singular y det(I−ωL) = 1, det(λI −Lω) = det(I−ωL) det(λI −Lω) = det{(λ+ω−1)I−ωλL −ωU}. Entonces, llamando φ(λ) = det{(λ+ω−1)I−ωλL −ωU}, (1.30) es equivalente a φ(λ) = det{(λ+ω−1)I−ωλL −ωU}= 0. Se puede probar que φ(λ) = det{(λ+ω−1)I−λ(p−1)/pωB}.(1.31) Como Aes p-c´ıclica, entonces B es d´ebilmente c´ıclica de ´ındice p y, por tanto, λ(p−1)/pωB es d´ebilmente c´ıclica de ´ındice p. Aplicando el teorema de Romanovsky 1.6, φ(λ)=(λ+ω−1)m r Y i=1{(λ+ω−1)p−λp−1ωpµp i},(1.32) donde los µison no nulos si r≥1. Veamos en primer lugar la segunda implicaci´on del teorema. Supongamos que µes un autovalor de By que λsatisface (1.29). Entonces, uno de los factores de (1.32) desaparece y tendr´ıamos que φ(λ) = 0 lo que implica directamente que λes autovalor de Lω. Veamos ahora la primera parte del teorema. Sea ω6= 0 y sea λun autovalor no nulo de Lω. Entonces, φ(λ) = 0, luego al menos uno de los factores de (1.32) ha de ser cero. Si µ6= 0 y µsatisface (1.29), entonces λ+ω−16= 0. Luego debe ser (λ+ω−1)p=λp−1ωpµp i, para alg´un i, 1≤i≤r, donde µies no nulo. Combinando esto con (1.29) tenemos que λp−1ωp(µp−µp i) = 0 y como λyµson no nulos, debe ser µp=µp i. Tomando ra´ıces p-´esimas, µ=µie2πir/p con rentero que satisface 0 ≤r < p. Como B es d´ebilmente c´ıclica, entonces µ6= 0 es autovalor de B. Si µ= 0, ω6= 0 y λes un autovalor no nulo de Lω, como µsatisface (1.29) entonces de (1.31) se deduce que φ(λ) = det B= 0 , luego µ= 0 es autovalor de B. 14
1.2 Corolario. Sea la matriz Acon la partici´on (1.20) una matriz p-c´ıclica y consistentemente ordenada cuyas submatrices diagonales son no singulares. Si µes un autovalor de la matriz de Jacobi por bloques B=L+U, entonces µpes un autovalor de L1. Rec´ıprocamente, si λes un autovalor no nulo de L1yµp=λ, entonces µes un autovalor de B. Por tanto, el m´etodo iterativo de Jacobi por bloques converge si y solo si el m´etodo iterativo de Gauss-Seidel por bloques converge, y si ambos convergen, entonces ρ(L1)=(ρ(B))p<1. Supongamos que la matriz Aes una matriz p-c´ıclica y consistentemente ordenada de la forma (1.20) y cuyas submatrices diagonales son no singulares. Supongamos tambi´en que la matriz de Jacobi B=L+Ues convergente. Por el corolario anterior, tambi´en la matriz L1es convergente y, por continuidad, tambi´en es convergente en un intervalo de ωque contenga al 1. Buscamos ωbtal que ρ(Lωb) = m´ınω∈Rρ(Lω). En concreto, si los autovalores de Bpson reales y no negativos, el valor de ωbque minimiza ρ(Lω) de forma ´unica es la ´unica raiz real y positiva de la ecuaci´on (ρ(B)ωb)p= [pp(p−1)1−p](ωb−1) donde ρ(B) denota el radio espectral de la matriz de Jacobi por bloques. 1.8 Teorema. Sea la matriz Ade la forma (1.20) una matriz p-c´ıclica consistentemente ordenada cuyas submatrices diagonales, Ai,i para 1≤i≤N son regulares. Si todos los autovalores de la p-´esima potencia de la matriz Bde Jacobi por bloques son reales y no negativos y 0≤ρ(B)<1entonces para ωb, la ´unica raiz real , positiva y menor que p/(p−1) de (ρ(B)ωb)p= [pp(p−1)1−p](ωb−1) se tiene que 1. ρ(Lωb) = (ωb−1)(p−1); 2. ρ(Lω)> ρ(Lωb), ω 6=ωb. Adem´as, la matriz Lωdel m´etodo SOR es convergente para todo ωcon 0< ω < p/(p−1). Demostraci´on. Para p= 2 la ra´ız de la ecuaci´on anterior que nos interesa se puede expresar como ωb=2 1 + p1−ρ2(B)= 1 + ρ(B) 1 + p1−ρ2(B)!2 .(1.33) Veamos que, efectivamente, la expresi´on (1.33) minimiza ρ(Lω) para p= 2 (la misma prueba, con las correspondientes modificaciones, ser´ıa v´alida para el caso general). Si los autovalores de B2son n´umeros reales no 15
negativos, entonces como Bes d´ebilmente c´ıclica de ´ındice 2, los autovalores no nulos de B,µivienen dados en pares de autovalores con signos opuestos. Por tanto, −ρ(B)≤µi≤ρ(B)<1. Por el teorema anterior sabemos que los autovalores de la matriz para el m´etodo SOR, λ, y los autovalores de la matriz de Jacobi, µ, satisfacen (λ+ω−1)2=λω2µ2. Operando tenemos λ+ω−1 ω=±λ1/2µ. Definamos entonces gω(λ) := λ+ω−1 ω, ω 6= 0 y m(λ) := λ1/2µ, 0≤µ≤ρ(B)<1. Figura 1.1: Optimizaci´on del par´ametro para el m´etodo SOR Como se puede observar en la Figura 1.1 , gω(λ) es una recta que pasa por el punto (1,1) y cuya pendiente decrece al aumentar el valor de ω. Como ten´ıamos gω(λ) = m(λ), se puede interpretar como la intersecci´on de estas dos funciones. Adem´as, el valor ´optimo para ωse obtendr´a cuando gω(λ) sea tangente a m(λ), ya que el mayor de los dos puntos de intersecci´on alcanzar´a 16
su valor m´ınimo cuando esto ocurra. En dicho caso ser´a ωb=2 1 + p1−µ2. Entonces, el valor de la abscisa en este punto es ˜ λ=ωb−1.Por tanto, m´ın ω∈Rρ(Lω) = ρ(Lωb) = ωb−1. Figura 1.2: Optimizaci´on del par´ametro para el m´etodo SOR Como podemos apreciar en la Figura 1.2, donde el pico corresponde al valor ´optimo ωb, es mejor sobreestimar el valor ´optimo en una peque˜na cantidad que subestimar en esta misma cantidad. 17
18
Cap´ıtulo 2 M´etodos SOR para sistemas aumentados 2.1. Los sistemas aumentados y su inter´es Sean A∈Rm×muna matriz sim´etrica y definida positiva y B∈Rm×n. Consideramos el siguiente sistema lineal aumentado A B BT0x y=b q(2.1) donde b∈Rmyq∈Rnson dos vectores dados. En ocasiones, denotamos a la matriz del sistema (2.1) con la letra A. Cuando en (2.1) las matrices AyBson grandes y dispersas, es cuando los m´etodos iterativos que vamos a tratar toman mayor relevancia, puesto que resolver´an el sistema con mayor eficiencia que los m´etodos directos. El sistema aumentado (2.1) aparece en diversos problemas como pueden ser problemas de optimizaci´on sujetos a restricciones, el m´etodo de los elementos finitos para resolver las ecuaciones de Navier-Stokes, problemas de elasticidad, problemas de m´ınimos cuadrados o problemas de puntos de silla. Muchas veces, por simplicidad, reescribiremos el sistema (2.1) como A B −BT0x y=b −q(2.2) Para sistemas aumentados, la matriz del segundo bloque diagonal es nula, por lo que el m´etodo SOR est´andar no es aplicable. Desarrollaremos a continuaci´on el m´etodo SOR para sistemas aumentados como (2.1), bas´andonos en el art´ıculo [5]. Adem´as, estudiaremos en este cap´ıtulo la teor´ıa del m´etodo SOR generalizado para sistemas aumentados utilizando [2] como referencia. 19
y fω(λ) = 1 −4λµ. Observemos que gω(1) = 1 y fω(0) = 1, esto es,gω(λ) pasa por el punto (1,1) yfω(λ) pasa por (0,1). Entonces, la recta fωcruza la curva parab´olica gω. El valor ´optimo para el par´ametro ωse obtendr´a cuando fωsea tangente a gω. Esto sucede cuando ωb=2√ρ−1 ρ ya que g0 ω(λ) = 8(λ−1 + ω/2)/ω2yf0 ω(λ) = −4µ. Adem´as, ρ(Mωb) = √1−ωb=|√ρ−1| √ρ. Figura 2.1: Curva del radio espectral del m´etodo SOR aumentado para ρ= 2 26
2.3 Teorema. Supongamos que ρ0≤1 4. Entonces, ρ(Mω) = h|2−ω−ω2µ0|+ωp(ωµ0+ 1)2−4µ0i 2,0< ω ≤ωb h|2−ω−ω2ρ|+ωp(ωρ + 1)2−4ρi 2, ωb< ω < 4 √4ρ+1+1 (2.12) donde ωb<2es la ra´ız positiva de la ecuaci´on |2−ω−ω2µ0|+ωp(ωµ0+ 1)2−4µ0=|2−ω−ω2ρ|+ωp(ωρ + 1)2−4ρ. Demostraci´on. De (2.11) deducimos que para todo µ0< µ ≤1 4, ocurre |λ(µ)| ≤ |λ(µ0)|. Tenemos |λ(µ0)|= 0,5h|2−ω−ω2µ0|+ωp(ωµ0+ 1)2−4µ0i. Adem´as, |λ(µ0)| ≥ √1−ωpara µ > 1 4. Para ω > (2√ρ−1) ρ,se tiene |λ(µ0)|≤|λ(ρ)|, donde |λ(ρ)|= 0,5h|2−ω−ω2ρ|+ωp(ωρ + 1)2−4ρi. Combinando las dos igualdades que hemos obtenido, llegamos al resultado que quer´ıamos probar. 2.5. El m´etodo SOR generalizado para sistemas aumentados Con el objetivo de mejorar la velocidad de convergencia del m´etodo SOR para sistemas aumentados, se desarrolla, mediante la introducci´on de un nuevo par´ametro, el m´etodo generalizado SOR (GSOR). Consideramos, para el sistema aumentado (2.2), la misma escisi´on que en el m´etodo SOR, es decir, A B −BT0=D−L−U donde D=A0 0QL=0 0 BT0U=0−B 0Q yQ∈Rn×nes una matriz sim´etrica y no singular. 27
Sean ahora ωyτdos reales no nulos. Llamemos Ω = ωIm0 0τIn donde Im∈Rm×meIn∈Rn×nson, respectivamente, la matriz identidad m×myn×n. Generalizamos, para el sistema aumentado (2.2), el m´etodo SOR como sigue x(k+1) y(k+1) = (D−ΩL)−1[(I−Ω)D+ΩU]x(k) y(k)+(D−ΩL)−1Ωb −q. Equivalentemente, si definimos las matrices H(ω, τ)=(D−ΩL)−1[(I−Ω)D+ ΩU] =(1 −ω)I−ωA−1B (1 −ω)τQ−1BTI−ωτQ−1BTA−1B y M(ω, τ)=Ω−1(D−ΩL) = 1 ωA0 −BT1 τQ podemos escribir la versi´on matricial del m´etodo GSOR para sistemas ampliados x(k+1) y(k+1) =H(ω, τ)x(k) y(k)+M(ω, τ)−1b −q. De esta forma, el m´etodo generalizado SOR en su versi´on componente a componente se escribe como x(k+1) = (1 −ω)x(k)+ωA−1(b−By(k)) y(k+1) =y(k)+τQ−1(BTx(k+1) −q)(2.13) donde Qes una matriz que aproxima al complemento de Schur, BTA−1B, y sirve como preacondicionamiento. Observemos que de forma trivial cuando τ=ω, el m´etodo generalizado SOR se reduce al m´etodo SOR expuesto en el apartado anterior. Denotemos ahora N(ω, τ) = M(ω, τ)−A =(1 ω−1)A−B 01 τQ Entonces, podemos considerar la siguiente escisi´on para la matriz A, matriz de coeficientes del sistema aumentado (2.2), A=M(ω, τ)− N(ω, τ). Entonces, es f´acil ver que H(ω, τ) = M(ω, τ)−1N(ω, τ) es la matriz de iteraci´on del m´etodo SOR generalizado. 28
2.6. Teoremas de convergencia para el m´etodo SOR generalizado Hemos visto que la matriz H(ω, τ) es la matriz de iteraci´on del m´etodo GSOR. Por tanto, el m´etodo ser´a convergente si y solo si el radio espectral de la matriz H(ω, τ) es estrictamente menor que 1, esto es, ρ(H(ω, τ)) <1. 2.4 Teorema. Sean A∈Rm×myQ∈Rn×nmatrices sim´etricas y definidas positivas. Sea tambi´en B∈Rm×nde rango m´aximo. Denotemos J=Q−1BTA−1Byµm´ın,µm´ax al menor y mayor autovalor de J, respectivamente. Entonces, el m´etodo GSOR converge si ωverifica que 0< ω < 2 yτsatisfacela condici´on: 0< τ < 2(2 −ω) ωµm´ax . Demostraci´on. Por las condiciones establecidas en el enunciado del teorema, se tiene que todos los autovalores, µ, de la matriz J=Q−1BTA−1Bson reales y no nulos. Sea λun autovalor no nulo de la matriz de iteraci´on H(ω, τ), con autovector asociado (xT, yT)T∈Rm+n. Recordemos que H(ω, τ) = M(ω, τ)−1N(ω, τ) = 1 ωA0 −BT1 τQ−1(1 ω−1)A−B 01 τQ Entonces, la relaci´on H(ω, τ)x y=λx yse puede escribir como: (1 ω−1)A−B 01 τQx y=λ1 ωA0 −BT1 τQx y Equivalentemente, (1 −ω)A−ωB 0Qx y=λA0 −τBTQx y Entonces tenemos (1 −ω)Ax −ωBy =λAx Qy =−λτBTx+λQy luego (1 −ω−λ)Ax =ωBy λτBTx= (λ−1)Qy (2.14) De la primera ecuaci´on de (2.14) obtenemos, (1 −ω−λ)x=ωA−1By. Luego cuando λ6= 1 −ω, se tiene que λτ(1 −ω−λ)BTx=λτωBTA−1By. 29
De aqu´ı y de la segunda ecuaci´on de (2.14), se tiene que (λ−1)(1 −ω−λ)Qy =λτωBTA−1By =⇒ (λ−1)(1 −ω−λ)y=λτωQ−1BTA−1By =⇒ (λ−1)(1 −ω−λ)y=λτωJy Si λ= 1 −ω6= 0, entonces de la primera ecuaci´on de (2.14) tenemos 0 = By y de la segunda ecuaci´on de (2.14), −ωQy =λτBTx. Por tanto, y= 0 y x∈ker(BT), donde ker(BT) es el n´ucleo de BT. Luego, λ= 1−ωes un autovalor de H(ω, τ) cuyo autovector correspondiente es (xT,0)T, donde x∈ker(BT). Entonces, los autovalores λ(excepto λ= 1 −ω)de la matriz H(ω, τ) y los autovalores µde la matriz Jsatisfacen la relaci´on siguiente (1 −ω−λ)(λ−1) = λτωµ es decir, λsatisface la ecuaci´on cuadr´atica λ2+ (τωµ +ω−2)λ+ 1 −ω= 0 (2.15) Entonces, de acuerdo con el lema de Young, sabemos que tanto λ= 1 −ω como las dos ra´ıces de la ecuaci´on cuadr´atica anterior, satisfacen |λ|<1 si y solo si |1−ω|<1 y |τωµ +ω−2|<2−ω. De la primera ecuaci´on se obtiene directamente que 0 < ω < 2 y de la segunda, −2< τωµ +ω−2<2−ω=⇒0< τωµ < 4−2ω=⇒ 0< τ < 2(2 −ω) ωµ =⇒0< τ < 2(2 −ω) ωµm´ax . De forma an´aloga se puede probar que si Q∈Cn×nes sim´etrica y definida negativa, entonces el m´etodo GSOR converge si ωsatisface 0 < ω < 2 y τ satisface 2(2 −ω) ωµm´ın < τ < 0. 2.1 Corolario. Sean A∈Rm×muna matriz sim´etrica y positiva definida,B∈ Rm×nde rango m´aximo y Q∈Rn×nno singular y sim´etrica. Sea tambi´en J=Q−1BTA−1B. Si µes un autovalor de la matriz J, entonces el λ determinado por la ecuaci´on cuadr´atica (2.15) es un autovalor de la matriz H(ω, τ). Rec´ıprocamente, si λes un autovalor de la matriz H(ω, τ), entonces el µdeterminado por (2.15) es un autovalor de la matriz J. Por tanto, los autovalores no nulos de la matriz H(ω, τ)vienen dados por λ= 1 −ωo λ=1 22−ω−τωµ ±p(2 −ω−τωµ)2−4(1 −ω). 30
2.2 Corolario. Sean A∈Rm×muna matriz sim´etrica y positiva definida,B∈ Rm×nde rango m´aximo y Q∈Rn×nno singular y sim´etrica. Sea tambi´en J=Q−1BTA−1B. Si µes un autovalor de la matriz J, entonces el λ determinado por la ecuaci´on cuadr´atica (2.15) con τ=ωes un autovalor de la matriz L(ω).Rec´ıprocamente, si λes un autovalor de la matriz L(ω), entonces el µdeterminado por (2.15) con τ=ωes un autovalor de la matriz J. Por tanto, los autovalores no nulos de la matriz L(ω)vienen dados por λ= 1 −ωo λ=1 22−ω−ω2µ±p(2 −ω−ω2µ)2−4(1 −ω). 2.7. Optimizaci´on de los par´ametros en el m´etodo SOR generalizado Determinamos en esta secci´on los valores ´optimos de los par´ametros del m´etodo GSOR y el correspondiente factor ´optimo de convergencia. Para aligerar la notaci´on llamemos: µk(k= 1,2, . . . , n) son los autovalores de la matriz J=Q−1BTA−1B µm´ın = m´ın 1≤k≤nµk µm´ax = m´ax 1≤k≤nµk Asumimos, sin p´erdida de generalidad, que estos autovalores est´an ordenados, es decir 0< µm´ın =µ1≤µ2≤. . . ≤µn−1≤µn=µm´ax Hemos probado en la secci´on anterior que los par´ametros del m´etodo deben satisfacer 0< ω < 2 y 0 < τ < 2(2 −ω) ωµm´ax Denotemos tambi´en f1(ω, τ, µ) = 1 22−ω−τωµ +p(2 −ω−τωµ)2−4(1 −ω) para 4τµ (1+τµ)2≤ω < 2 1+τµ yτµ < 1 f2(ω, τ, µ) = 1 2τωµ +ω−2 + p(τωµ +ω−2)2−4(1 −ω) para ω≥4τµ (1+τµ)2yτµ > 1 oω > 2 1+τµ yτµ < 1 f3(ω, τ, µ) = g(ω) = √1−ω para ω≤4τµ (1+τµ)2 31
Estas funciones se obtienen de calcular el valor absoluto de los autovalores no nulos de la matriz de iteraci´on H(ω, τ) del m´etodo GSOR, donde hemos distinguidos tres casos seg´un sea positivo o negativo el discriminante de la ecuaci´on cuadr´atica y seg´un el signo del t´ermino 2 −ω−τωµ. Se puede ver que fj(ω, τ, µ)≥√1−ω≥1−ω, j = 1,2. Estudiemos ahora la monoton´ıa de estas funciones. ∂f1(ω,τ,µ) ∂µ =−τω 21 + 2−ω−τωµ √(2−ω−τωµ)2−4(1−ω) ∂f2(ω,τ,µ) ∂µ =τω 21 + τωµ+ω−2 √(τωµ+ω−2)2−4(1−ω) ∂f3(ω,τ,µ) ∂µ =dg(ω) dµ = 0 Luego, ∂f1(ω,τ,µ) ∂µ <0 para 4τµ (1+τµ)2< ω < 2 1+τµ yτµ < 1 ∂f2(ω,τ,µ) ∂µ >0 para ω > 4τµ (1+τµ)2yτµ > 1 oω > 2 1+τµ yτµ < 1 Es decir, la funci´on f1(ω, τ, µ) decrece con respecto a µpara 4τµ (1+τµ)2< ω < 2 1+τµ yτµ < 1. La funci´on f2(ω, τ, µ) crece con respecto a µpara ω > 4τµ (1+τµ)2yτµ > 1 oω > 2 1+τµ yτµ < 1. Por otro lado, ∂f1(ω,τ,µ) ∂ω =−1 21 + τµ +(2−ω−τωµ)(1+τµ)−2 √(2−ω−τωµ)2−4(1−ω) ∂f2(ω,τ,µ) ∂ω =1 21 + τµ +(τωµ+ω−2)(1+τµ)+2 √(τωµ+ω−2)2−4(1−ω) ∂f3(ω,τ,µ) ∂ω =dg(ω) dµ =−1 2√1−ω Luego, ∂f1(ω,τ,µ) ∂ω >0 para 4τµ (1+τµ)2< ω < 2 1+τµ yτµ < 1 ∂f2(ω,τ,µ) ∂ω >0 para ω > 4τµ (1+τµ)2yτµ > 1 oω > 2 1+τµ yτµ < 1 ∂f3(ω,τ,µ) ∂ω <0 para ω≤4τµ (1+τµ)2 Es decir, la funci´on f1(ω, τ, µ) crece con respecto a ωpara 4τµ (1+τµ)2< ω < 2 1+τµ yτµ < 1. La funci´on f2(ω, τ, µ) crece con respecto a ωpara ω > 4τµ (1+τµ)2yτµ > 1 oω > 2 1+τµ yτµ < 1. La funci´on f3(ω, τ, µ) decrece con respecto a ωpara ω≤4τµ (1+τµ)2. 32
Adem´as, para cualesquiera ˜µy ˜τreales positivos para los que las funciones est´en bien definidas, f1(ω, ˜τ, ˜µ) = f3(ω, ˜τ, ˜µ) para ω=4˜τ˜µ (1+˜τ˜µ)2 f2(ω, ˜τ, ˜µ) = f3(ω, ˜τ, ˜µ) para ω=4˜τ˜µ (1+˜τ˜µ)2 f1(ω, ˜τ, ˜µ) = f2(ω, ˜τ, ˜µ) para ω=2 1+˜τ˜µ y para dos autovalores de J,µα, µβ∈ { µk|1≤k≤n}para los que f1 yf2est´en bien definidas, tenemos f1(ω, ˜τ, µα) = f2(ω, ˜τ, µβ) si ω=4 ˜τ(µα+µβ)+2 Por otra parte, definimos las siguientes funciones ω−(τ) = 4τµm´ın (1+τµm´ın )2 ω+(τ) = 4τµm´ax (1+τµm´ax)2 ω0(τ) = 4 τ(µm´ın+µm´ax)+2 Usando el primer corolario de la secci´on anterior, podemos expresar el valor absoluto de los autovalores λde la matriz H(ω, τ) como: Cuando µτ < 1 o bien |λ|=|1−ω|o bien |λ|= f1(ω, τ, µ) para 4τµ (1+τµ)2< ω < 2 1+τµ f2(ω, τ, µ) para ω≥2 1+τµ g(ω) para ω≤4τµ (1+τµ)2 y cuando µτ ≥1 entonces o bien |λ|=|1−ω|o bien |λ|=(f2(ω, τ, µ) para ω > 4τµ (1+τµ)2 g(ω) para ω≤4τµ (1+τµ)2 2.5 Teorema. Consideremos el m´etodo GSOR. Sea A∈Rm×myQ∈Rn×n sim´etrica y definida positiva, y sea B∈Rm×nde rango m´aximo. Denotemos µm´ın,µm´ax al menor y mayor autovalor de J=Q−1BTA−1B, respectivamente. Entonces, (i) cuando τ≤1 √µm´ınµm´ax , se tiene que ρ(H(ω, τ)) = √1−ω para 0< ω < ω−(τ). 1 22−ω−τωµm´ın +p(2 −ω−τωµm´ın)2−4(1 −ω) para ω−(τ)≤ω≤ω0(τ). 1 2τωµm´ax +ω−2 + p(τωµm´ax +ω−2)2−4(1 −ω) para ω0(τ)< ω < 2. (2.16) 33
(ii) cuando 1 √µm´ınµm´ax < τ < 2(2−ω) ωµm´ax ρ(H(ω, τ)) = √1−ω para 0< ω < ω+(τ). 1 2τωµm´ax +ω−2 + p(τωµm´ax +ω−2)2−4(1 −ω) para ω+(τ)≤ω < 2. (2.17) Adem´as, los par´ametros ´optimos ωopt yτopt vienen dados por ωopt =4√µm´ınµm´ax (√µm´ın+√µm´ax)2 τopt =1 √µm´ınµm´ax (2.18) y por tanto,el factor de convergencia ´optimo para el m´etodo GSOR es ρ(H(ωopt, τopt)) = √µm´ax −√µm´ın √µm´ax +√µm´ın . Demostraci´on. Para simplificar la prueba vamos a distinguir tres casos en funci´on del par´ametro τ. Caso (a): τ≤1 µm´ax . Caso (b): τ≥1 µm´ın . Caso (c): 1 µm´ax < τ < 1 µm´ın Caso (a): τ≤1 µm´ax . Como τµm´ax ≤1, entonces para todo µautovalor de Jse tiene τµ < 1. Adem´as, ω−(τ)≤ω+(τ)<1< ω0(τ)<2. Recordemos que para λ, autovalor de H(ω, τ) o bien |λ|=|1−ω|o bien |λ|= f1(ω, τ, µ) para 4τµ (1+τµ)2< ω < 2 1+τµ f2(ω, τ, µ) para ω≥2 1+τµ g(ω) para ω≤4τµ (1+τµ)2 Consideremos fijos ω, τ > 0, entonces por la monoton´ıa estudiada para las funciones fi(ω, τ, µ), i = 1,2,3 con respecto a µse tiene que m´axµnf1(ω, τ, µ)|4τµ (1+τµ)2< ω < 2 1+τµ o=f1(ω, τ, µm´ın) m´axµnf2(ω, τ, µ)|ω≥2 1+τµ o=f2(ω, τ, µm´ax) m´axµnf3(ω, τ, µ)|ω≤4τµ (1+τµ)2o=g(ω) Adem´as, sabemos que el punto de intersecci´on de f1(ω, τ, µm´ın) con f2(ω, τ, µm´ax) es ω0(τ). El punto de intersecci´on de f1(ω, τ, µm´ın) y de g(ω) es ω−(τ). 34
Entonces, por la monoton´ıa estudiada de las funciones fi(ω, τ, µ), i = 1,2,3 con respecto a ωse tiene que ρ(H(ω, τ)) = g(ω) para 0 < ω < ω−(τ) f1(ω, τ, µm´ın) para ω−(τ)≤ω≤ω0(τ) f2(ω, τ, µm´ax) para ω0(τ)< ω < 2. Por tanto, para cualquier τfijo argminωρ(H(ω, τ)) = ω−(τ) Como ρ(H(ω−(τ), τ)) = p1−ω−(τ), el valor de τpara el cual se minimiza ρ(H(ω−(τ), τ)) es aquel que maximiza ω−(τ). En conclusi´on, para este caso, los par´ametros ´optimos son τ(1) opt =1 µm´ax yω(1) opt =4µm´ınµm´ax (µm´ın +µm´ax)2 y el factor de convergencia ´optimo asociado es ρ(H(ω(1) opt, τ(1) opt)) = µm´ax −µm´ın µm´ax +µm´ın . Caso (b): τ≥1 µm´ın . Como τµm´ın ≥1, entonces para todo µautovalor de Jse tiene τµ > 1. Recordemos que para λ, autovalor de H(ω, τ) o bien |λ|=|1−ω|o bien |λ|=(f2(ω, τ, µ) para ω > 4τµ (1+τµ)2 g(ω) para ω≤4τµ (1+τµ)2 Consideremos fijos ω, τ > 0, entonces por la monoton´ıa estudiada para las funciones fi(ω, τ, µ), i = 2,3 con respecto a µse tiene que m´axµnf2(ω, τ, µ)|ω > 4τµ (1+τµ)2o=f2(ω, τ, µm´ax) m´axµnf3(ω, τ, µ)|ω≤4τµ (1+τµ)2o=g(ω) Adem´as, sabemos que el punto de intersecci´on de f2(ω, τ, µm´ax) y de g(ω) es ω+(τ). Entonces, por la monoton´ıa estudiada de las funciones fi(ω, τ, µ), i = 2,3 con respecto a ωse tiene que ρ(H(ω, τ)) = g(ω) para 0 < ω < ω+(τ) f2(ω, τ, µm´ax) para ω+(τ)< ω < 2. Por tanto, para cualquier τfijo argminωρ(H(ω, τ)) = ω+(τ) 35
3.2. El m´etodo HSS para sistemas aumentados Consideremos el sistema Ax=bque proviene de un problema de punto de silla, donde A=A BH −B0, x =u p, b =f −g. donde A∈Cn×nes una matriz herm´ıtica y semidefinida positiva, B∈ Cm×n tiene rango m´aximo (m≤n), f∈ Cnyg∈ Cm. Denotemos ahora, H=1 2(A+AH), a la parte herm´ıtica de A. Entonces, H=A0 0 0 . Denotemos tambi´en, S=1 2(A − AH), a la parte anti-herm´ıtica de A. Entonces, S=0BH −B0. Consideremos las dos siguientes escisiones de la matriz A, donde α > 0 es un par´ametro e Idenota la matriz indentidad (n+m)×(n+m): A= (H+αI)−(αI −S) A= (S+αI)−(αI −H) Observemos que tanto H+αI=A+αIn0 0αIm como S+αI=αInBH −B αIm son matrices regulares. La iteraci´on herm´ıtica/anti-herm´ıtica (HSS) se obtiene alternando las dos escisiones anteriores en la definici´on de los iterantes. Para un vector inicial x(0) = (u(0), p(0)), la secuencia de vectores {x(m)}para el m´etodo HSS viene dada por: (H+αI)x(m+1/2) = (αI −S)x(m)+b, (S+αI)x(m+1) = (αI −H)x(m+1/2) +b. (3.2) La primera ecuaci´on se traduce en A+αIn0 0αImu(m+1/2) p(m+1/2) =αIn−BH B αImu(m) p(m)+f −g 42
Entonces, (A+αIn)u(m+1/2) =αu(m)−BHp(m)+f p(m+1/2) =p(m)+1 α(Bu(m)−g)(3.3) Como la matriz (A+αIn) es no singular, el sistema tiene soluci´on ´unica que podemos determinar. Adem´as, la matriz es herm´ıtica y definida positiva, por lo que para determinar la soluci´on se puede emplear cualquiera de los m´etodos para resolver sistemas con matriz herm´ıtica y definida positiva, como pueden ser la factorizaci´on de Cholesky o el m´etodo del gradiente conjugado. Para resolver la segunda parte de (3.2), nos encontramos con el siguiente sistema αInBH −B αImu(m+1) p(m+1) =αIn−A0 0αImu(m+1/2) p(m+1/2) +f −g es decir, αu(m+1) +BHp(m+1) = (αIn−A)u(m+1/2+f −Bu(m+1) +αp(m+1) =αp(m+1/2) −g(3.4) Para resolver este sistema, denotemos f(m):= (αIn−A)u(m+1/2) +f y g(m):= αp(m+1/2) −g. Entonces podemos reducir el sistema anterior eliminando u(m+1) obteniendo (BBH+α2Im)p(m+1) =Bf(m)+αg(m) y despu´es, despejando de la primera ecuaci´on tenemos que u(m+1) =1 α(f(m)−BHp(m+1)). Sin embargo, observemos que el m´etodo HSS se puede usar como preacondicionamiento, sin necesidad de resolver los sistemas (3.3) y (3.4) de forma exacta. Se pueden resolver de forma inexacta (IHSS) para hacer este m´etodo m´as competitivo. 3.3. Convergencia del m´etodo HSS Para analizar la convergencia del m´etodo HSS (3.2) podemos eliminar el vector intermedio x(m+1/2) de (3.2), para ello utlizamos el siguiente lema que describe un criterio de convergencia para m´etodos con iteraci´on en dos pasos. 43
3.1 Lema. Sea A ∈ Cn×n. Consideremos dos escisiones para A,A=Mi−Ni para i= 1,2, donde las matrices M1yM2son no singulares. Sea entonces x(0) ∈Cnun vector inicial dado. Si la secuencia de iterantes {x(m)}viene descrita por: M1x(m+1/2) =N1x(m)+b, M2x(m+1) =N2x(m+1/2) +b, para m= 0,1,2, . . .. Entonces, se puede escribir x(m+1) =M−1 2N2M−1 1N1x(m)+M−1 2(I+N2M−1 1)b, m = 0,1,2, . . . Adem´as, si el radio espectral de la matriz de iteraci´on M−1 2N2M−1 1N1es ρ(M−1 2N2M−1 1N1)<1, entonces la secuencia de iterantes {x(m)}converge a la soluci´on ´unica x∗∈Cndel sistema de ecuaciones lineales anterior para todo vector inicial x(0) ∈Cn. Tomemos entonces en el lema anterior, M1=αI+H,N1=αI − S, M2=αI+SyN2=αI −H. As´ı, podemos escribir el m´etodo como x(m+1) =Tαx(m)+c, (3.5) donde la matriz de iteraci´on del m´etodo es Tα:= (S+αI)−1(αI −H)(H+αI)−1(αI −S) (3.6) y c:= (S+αI)−1[I+ (αI −H)(H+αI)−1]b. Por otro lado, este m´etodo HSS tambi´en puede escribirse en ”forma corregida”: x(m+1) =x(m)+M−1 αr(m), donde la matriz de iteraci´on es M−1 αpara Mα:= 1 2α(H+αI)(S+αI), y r(m)=b−Ax(m). Observemos que Tα=I −M−1 αA. 3.1 Teorema. Sea A ∈ Cn×nuna matriz definida positiva y αuna constante positiva. Consideremos el m´etodo HSS para esta matriz, donde la matriz de iteraci´on del m´etodo viene dada por Tα= (S+αI)−1(αI −H)(H+αI)−1(αI −S), 44
para S=1 2(A−AH)yH=1 2(A+AH). Entonces, para el radio espectral de esta matriz se tiene que ρ(Tα)≤σ(α)<1,∀α > 0, es decir, el m´etodo HSS converge a la soluci´on ´unica x∗∈Cndel sistema de ecuaciones lineales Ax=b, para σ(α) = m´ax λi∈λ(H) α−λi α+λi donde λ(H)es el espectro de la matriz H. Demostraci´on. Observemos que las matrices (αI+S)−1(αI − H)(αI+ H)−1(αI − S) y (αI − H)(αI+H)−1(αI − S)(αI+S)−1son matrices semejantes. Por tanto, como el espectro de una matriz es invariante por semejanza, se tiene que ρ(Tα) = ρ((αI+S)−1(αI −H)(αI+H)−1(αI −S)) =ρ((αI −H)(αI+H)−1(αI −S)(αI+S)−1) ≤ k(αI −H)(αI+H)−1(αI −S)(αI+S)−1k2 ≤ k(αI −H)(αI+H)−1k2k(αI −S)(αI+S)−1k2. Veamos que k(αI −S)(αI+S)−1k2= 1. Denotemos para ello Q(α) = (αI−S)(αI+S)−1y como sabemos que Ses anti-herm´ıtica, es decir SH= −S, tenemos que Q(α)HQ(α) = {(αI −S)(αI+S)−1}H(αI −S)(αI+S)−1 ={(αI+S)−1}H(αI −S)H(αI −S)(αI+S)−1 ={(αI+S)H}−1{(αI)H−(S)H}(αI −S)(αI+S)−1 = (αI −S)−1(αI+S)−1(αI −S)(αI+S)−1 = (αI −S)−1(αI −S)−1(αI+S)(αI+S)−1 =I Por tanto, kQ(α)k2= 1. De aqu´ı se deduce que ρ(Tα)≤ k(αI −H)(αI+H)−1k2= m´ax λi∈λH α−λi α+λi. Como Aes definida positiva, tambi´en lo es Hy, por tanto λi>0 para todo i= 1,2, . . . , n. Adem´as, α > 0, luego se deduce que ρ(Tα)≤σ(α)<1. Como acabamos de ver, para el caso en que la matriz Aes definida positiva, el m´etodo HSS es convergente para cualquier valor de α. Sin embargo, para el caso de problemas generales de punto de silla, la matriz Hes solo semidefinida positiva y es, en general, singular. Por tanto, no es aplicable a este caso el an´alisis de convergencia anterior. 45
3.2 Teorema. Para el problema (3.1) supongamos que la matriz Aes una matriz real y positiva, Ces una matriz sim´etrica y semidefinida positiva yBtiene rango m´aximo. Entonces, para la iteraci´on (3.5) se tiene que ρ(Tα)<1para todo α > 0. Demostraci´on. La matriz de iteraci´on del m´etodo, Tα:= (S+αI)−1(αI − H)(H+αI)−1(αI − S) es semejante a ˆ Tα:= (αI − H)(H+αI)−1(αI − S)(S+αI)−1=RU, para R= (αI −H)(H+αI)−1 y U= (αI −S)(S+αI)−1. Observemos que la matriz Res sim´etrica por serlo Hy la matriz Ues ortogonal, es decir, su inversa coincide con su traspuesta. La matriz Res ortogonalmete semejante a la matriz diagonal (n+m)×(n+m), D, es decir, existe una matriz ortogonal Ptal que R=PTDP, donde la matriz Dviene dada por D= α−µ1 α+µ1α−µ2 α+µ2 ... α−µn α+µnα−ν1 α+ν1α−ν2 α+ν2... α−νm α+νm =D10 0D2 para µi,i= 1,...n autovalores de Hyνi,i= 1,...n son los autovalores no negativos de C.D1es la matriz diagonal n×ncuyas entradas son los autovalores de HyD2es la matriz diagonal m×mcuyas entradas son los autovalores no negativos de C. Observemos que α−µi α+µi<1 para i= 1, . . . , n yα−νi α+νi<1 para i= 1, . . . , m. Como para Rtenemos PTRP=D, entonces RU es ortogonalmente semejante a PTRUP= (PTRP)(PTUP) = DQ para Q=PTUP. Adem´as, por ser producto de matrices ortogonales es una matriz ortogonal. Entonces, la matriz Tαes semejante a DQ, luego ρ(Tα) = ρ(DQ) = ρ(QD). 46
Veamos ahora que ρ(QD)<1 para todo α > 0. Consideremos la partici´on de Q, Q=Q1,1Q1,2 Q2,1Q2,2. Entonces, QD =Q1,1D1Q1,2D2 Q2,1D1Q2,2D2. Sea λ∈Cun autovalor de QD y sea x∈Cn+mun autovector asociado. Podemos suponer, sin p´erdida de generalidad que kxk2= 1. Si λ= 0, entonces es obvio que |λ|<1, luego no hay nada que probar. Supongamos entonces que λ6= 0. Como λes autovalor, tenemos QDx= λx, luego por ser Qortogonal se tiene que Dx=λQtx. Luego kDxk2=|λ|kQTxk2=|λ|. Luego, |λ|2=kDxk2 2= n X i=1 α−µi α+µi2 xi˜xi+ n+m X i=n+1 α−νi α+νi2 xi˜xi≤ kxk2 2= 1. (3.7) Por tanto, |λ| ≤ 1. Veamos que, de hecho, esta desigualdad es estricta. Para ello vamos a probar que existe, al menos, un i(1 ≤i≤n) tal que xi6= 0. Si fuera xi= 0 para todo i= 1, . . . , n, es decir, si fuera x=0 ˆx, entonces QDx=λx ser´ıa QDx=Q1,1D1Q1,2D2 Q2,1D1Q2,2D20 ˆx=Q1,2D2ˆx Q2,2D2ˆx=0 λˆx. De forma que ser´ıa Q1,2D2ˆx= 0. Veamos ahora que Q1,2tiene rango m´aximo. Recordemos que Q= PTUP, con P=P1,10 0P2,2. Donde P1,1∈Rnes la matriz ortogonal que diagonaliza (αIn−H)(αIn+H)−1yP2,2∈Rmes la matriz ortogonal que diagonaliza (αIm−C)(αIm+C)−1. Recordemos adem´as que la matriz Uviene dada por U= (αI −S)(αI+S)−1=αIn−S−BT B αImαIn+S BT −B αIm−1 =U1,1U1,2 U2,1U2,2. Operando, U1,2=−(αIn−S)(αIn+S)−1+InBTαIm+B(αIn+S)−1BT−1. 47
Observemos que (αIn−S)(αIn+S)−1no puede tener a −1 como autovalor, por tanto, (αIn−S)(αIn+S)−1+Ines no singular. Adem´as, como B(αIn+ S)−1BTes real y positiva, la matriz αIn+B(αIn+S)−1BTes no singular. Tambi´en, Q=PTUP=PT 1,1U1,1P1,1PT 1,1U1,2P2,2 PT 2,2U2,1P1,1PT 2,2U2,2P2,2. Luego, Q1,2=PT 1,1U1,2P2,2 =−PT 1,1(αIn−S)(αIn+S)−1+InBTαIm+B(αIn+S)−1BT−1P2,2. Esto implica que Q1,2tiene rango m´aximo puesto que tanto PT 1,1como P2,2son ortogonales y BTtiene rango m´aximo. Recordemos que ten´ıamos que Q1,2D2ˆx= 0 y como acabamos de probar que Q1,2tiene rango m´aximo, entonces ha de ser D2ˆx= 0. Pero por otro lado ten´ıamos que Q2,2D2ˆx=λˆx, luego tendr´ıamos que λˆx= 0 y como hab´ıamos asumido que λ6= 0, entonces deber´ıa ser ˆx= 0, y por tanto, x= 0. Sin embargo, esto contradice el hecho de que kxk2= 1, por tanto, debe existir al menos un i(1 ≤i≤n) tal que xi6= 0. Entonces de (3.7) y como α−µi α+µi<1 para 1 ≤i≤n, se deduce que |λ|<1. Por tanto, el m´etodo es convergente. 3.4. Optimizaci´on del par´ametro para el m´etodo HSS Usando el primer teorema de la secci´on anterior y conociendo el m´aximo autovalor y el m´ınimo autovalor de la matriz Hpodemos obtener el valor de αque optimiza σ(α). 3.1 Proposici´on. Sea A∈Cn×nuna matriz definida positiva y sean γm´ın yγm´ax los autovalores m´ınimos y m´aximos, respectivamente, de la matriz H=1 2(A+AH)y sea αuna constante positiva. Entonces, α∗= arg m´ın αm´ax γm´ın≤λ≤γm´ax α−λ α+λ=√γm´ınγm´ax y σ(α∗) = √γm´ax −√γm´ın √γm´ax +√γm´ın =pκ(H)−1 pκ(H)+1 donde κ(H)es el n´umero de condici´on espectral de H. Demostraci´on. En este caso, σ(α) = m´ax γm´ın≤λ≤γm´ax α−λ α+λ= m´ax α−γm´ın α+γm´ın , α−γm´ax α+γm´ax . 48
Para obtener una aproximaci´on al factor αque optimice ρ(Tα), podemos minimizar la cota superior que tenemos para ρ(Tα), es decir, minimizamos σ(α). Sea α∗ese valor que minimiza σ(α). Entonces, α∗debe satisfacer, α∗−γm´ın >0, α∗−γm´ax <0 y α∗−γm´ın α∗+γm´ın =γm´ax −α∗ γm´ax +α∗. Operando, (α∗−γm´ın)(γm´ax +α∗)=(γm´ax −α∗)(α∗+γm´ın) α∗γm´ax + (α∗)2−γm´ınγm´ax −α∗γm´ın =α∗γm´ax +γm´axγm´ın −(α∗)2−α∗γm´ın 2(α∗)2= 2γm´axγm´ın α∗=√γm´axγm´ın Entonces, σ(α) = √γm´axγm´ın −γm´ın √γm´axγm´ın+γm´ın =√γm´axγm´ın−√γ2 m´ın √γm´axγm´ın+√γ2 m´ın =√γm´ın √γm´ın (√γm´ax−√γm´ın) (√γm´ax+√γm´ın) =√γm´ax−√γm´ın √γm´ax+√γm´ın . N´otese que este corolario nos da el valor ´optimo α∗que minimiza la cota superior del radio espectral de la matriz de iteraci´on, pero no minimiza el radio espectral en s´ı. 49
50
Bibliograf´ıa [1] Z. Z. Bai, G. H. Golub, M. K. Ng: Hermitian and Skew-Hermitian splitting methods for non-Hermitian positive definite linear systems. Siam J. Matrix Appl. Vol 24, No. 3, pp. 603-626 (2003) (Cited on p. 41) [2] Z.Z. Bai, B. N. Parlett, Z.Q. Wang: On generalized successive overrelaxation methods for augmented linear systems . Nummer. Math. 102; 1-38 (2005) (Cited on p. 19) [3] M. Benzi, M. J. Gander, G. H. Golub: Optimization of the Hermitian and Skew-Hermitian splitting iteration for saddle-point problems. BIT Numerical Mathematics 43: 881-900 (2003) (Not cited) [4] M. Benzi, G. H. Golub: A preconditioner for generalized saddle point problems. Siam J. Matrixx Anal. Appl. Vol. 26, No. 1, pp. 20-41 (2004) (Cited on p. 41) [5] G. H. Golub, X. Wu, J.Y. Yuan:SOR-like methods for augmented systems. BIT 41, 71-85 (2001) (Cited on p. 19) [6] R. S. Varga: Matrix Iterative Analysis. Prentice-Hall, Englewood Cliffs, N.J.(1962) (Cited on pp. 1, 6, 8, 9, 10) [7] D. M. Young, Iterative Solution of Large Linear Systems, Academic Press, New York, USA, 1971. (Cited on p. 1) 51