Capı́tulo 6 Autovalores e Autovetores 6.1 Introdução Neste capı́tulo, apresentaremos alguns dos métodos utilizados para a solução do problema do autovalor, i.e., o sistema de n equações lineares Ax = λx (6.1) para o qual procuramos um vetor solução x tal que xi = 0 para pelo menos algum i, ou seja, uma solução não-trivial. Para que tal seja possı́vel, é necessário que det(A − λI) = 0 (6.2) a qual é uma equação polinomial de grau n na variável λ, chamada de equação caracterı́stica de A; o polinômio det(A − λI) é chamado de polinômio caracterı́stico de A. As n raı́zes de (6.2) são chamadas de autovalores, raı́zes latentes ou valores caracterı́sticos de A. A cada raiz λ corresponde um vetor x ∈ ICn = 0 que satisfaz a equação (6.1), o qual é chamado de autovetor, vetor latente ou vetor caracterı́stico de A. Note que, se x é um autovetor de A, então kx, onde k ∈ IR, também é, pois Akx = kAx = λkx = kλx. Costumeiramente os autovetores são normalizados, i.e. || x || = 1 em alguma norma escolhida (o que pode ser feito pela relação acima). Se todas as raı́zes de (6.2) são distintas entre si, então isso implica em que a matriz A apresenta um conjunto completo de autovetores linearmente independentes (L.I.). No entanto, mesmo para casos em que os autovalores não são todos distintos, podemos encontrar um conjunto completo de autovetores L.I. Podemos também calcular os autovalores da matriz inversa de A, A−1 , a partir dos autovalores de A. Se multiplicarmos a equação (6.1) à esquerda por A−1 , temos x = λA−1 x ou A−1 x = 1 x. λ (6.3) Essa última equação nos diz que λ1 é autovalor de A−1 , onde λ é um autovalor de A, com o autovetor x correspondente. Problemas envolvendo autovalores e autovetores surgem em inúmeras aplicações, como podemos ver nos exemplos que seguem, conforme apresentados em [6]. 110 Introdução ao Cálculo Numérico Autovalores e Autovetores Exemplo 6.1 O estudo das vibrações de sistemas dinâmicos e de estruturas requer a solução de problemas de autovalores e autovetores. Considere, apenas para fins de explanação, o problema de se determinar as vibrações de pequenas partı́culas presas por um fio uniforme, sem peso, ao → − qual é aplicada uma força F nas extremidades (cf. a figura 6.1) e no qual desconsidera-se a ação da gravidade. As partı́culas encontram-se a distâncias iguais entre si e as vibrações das mesmas são consideradas pequenas e perpendiculares à posição de descanso do fio. Escrevendo as equações Figura 6.1: O problema das vibrações. diferenciais para as forças atuantes em cada partı́cula, temos: d2 x1 dt2 d2 x2 m2 2 dt d2 x3 m3 2 dt d2 x4 m4 2 dt m1 = = = = x2 − x1 x1 +F h h x3 − x2 x2 − x1 +F −F h h x3 − x4 x3 − x2 −F −F h h x4 x3 − x4 −F +F h h −F Introduzindo a notação x di = (x1 , x2 , x3 , x4 )T mi h , i = 1, 2, 3, 4 = F podemos escrever o sistema de equações diferenciais acima na forma matricial D onde D é a matriz diagonal D= d2 x = Tx dt2 (6.4) d1 d2 d3 d4 e T é a matriz tridiagonal −2 1 0 0 1 −2 1 0 . T = 0 1 −2 1 0 0 1 −2 A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 111 Introdução ao Cálculo Numérico Autovalores e Autovetores Quando as partı́culas vibram em fase ou em direções opostas, i.e., em modo normal, então a condição d2 x = −w2 x, w ∈ IR (6.5) dt2 é satisfeita. Substituindo a equação (6.5) em (6.4), obtemos o problema de autovalor Dwi2 xi = −T xi , i = 1, 2, 3, 4 (6.6) para as freqüências de vibração w1 , w2 , w3 e w4 e os modos normais correspondentes, i.e., os autovetores x1 , x2 , x3 e x4 . Aparentemente, se isolarmos x no lado direito da equação (6.6), obterı́amos o que se chama de problema generalizado do autovalor, cuja forma geral é (A − λB)x = 0 onde A e B são matrizes de ordem n. Porém, se introduzirmos o vetor y = D1/2 x o que é possı́vel, já que os elementos da diagonal de D são positivos, por definição, então podemos escrever (6.6) como D−1/2 T D−1/2 yi = −wi2 yi o qual recai na forma 6.1. Exemplo 6.2 A teoria de Leontief sobre a compra e a venda de produtos é muito utilizada no campo de estudo da macroeconomia; como exemplo, consideramos as vendas e compras de produtos num setor industrial. Seja bij as vendas da indústria i para a indústria j; bii representa os bens produzidos pela indústria i e retidos por ela própria. As vendas de bens da indústria i para o mercado é denotada por yi e o total de bens produzidos por xi . Então, bij (6.7) xi = yi + j A fim de definirmos bij , assume-se que as vendas da indústria i para a j estão em proporção constante à produção da indústria j, i.e. bij = aij xj onde aij são ditos coeficientes de entrada. Em uma situação estática, podemos escrever, a partir de (6.7), x = y + Ax (6.8) onde x = (x1 , x2 , . . . , xn )T e y = (y1 , y2 , . . . , yn )T e A é matriz de ordem n cujos elementos (i, j) são os coeficientes de entrada aij . Ora, a equação (6.8) pode ser reescrita como (I − A)x = y (6.9) onde I − A é chamada de matriz de Leontief. A equação (6.9) pode ser resolvida calculandose os autovalores e autovetores de A. Sua utilidade reside no fato de que, com ela, é possı́vel determinar-se a quantidade de bens produzidos (x) necessários para satisfazer a uma demanda final (y), pré-estabelecida. Se a produção e a demanda não se encontram em equilı́brio, então devemos considerar um modelo dinâmico, que leve em consideração a taxa de variação da produção. Nesse caso, usualmente considera-se que a produção em cada indústria varia a uma taxa proporcional à diferença entre os nı́veis de venda e de produção. Daı́, dx(t) = D ((A − I)x(t) + y(t)) dt A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha (6.10) 112 Introdução ao Cálculo Numérico Autovalores e Autovetores onde D é uma matriz diagonal de ordem n, cujos elementos dii representam os coeficientes de reação das indústrias. A equação (6.10) substitui nesse caso a equação (6.8) e representa o comportamento dinâmico do sistema econômico em estudo. Uma das questões a serem estudadas, nesse caso, é se o sistema é estável, determinando-se os autovalores e autovetores da matriz D(A − I). Particularmente, procura-se determinar se esses autovalores tem parte real positiva pois, como as soluções do sistema de equações diferenciais (6.10) são da forma eλi t , isso indicaria uma instabilidade, já que a demanda x(t) cresceria exponencialmente com o tempo. A seguir, apresentaremos dois importantes teoremas, os quais nos permitirão desenvolver técnicas de determinação de autovalores e autovetores para um tipo especı́fico de matrizes. 6.2 Teoremas de limites sobre autovalores Teorema 6.2.1 Discos de Gerschgorin: Seja A uma matriz de ordem n, e di , i = 1, 2, . . . , n os discos cujos centros são os elementos aii e cujos raios ri são dados por ri = n |aij |, i = 1, 2, . . . , n. j=1 j=i Seja D a união de todos os discos di . Então, todos os autovalores de A encontram-se contidos em D. Prova: Seja λ um autovalor de A e x um autovetor correspondente, tal que maxi | xi | = 1. Então, λx = Ax de onde (λ − aii )xi = n aij xj , i = 1, 2, . . . n j=1 j=i Supondo que | xk | = 1, então | λ − akk | ≤ n | akj || xj | j=1 j=i ≤ n | akj | = rk j=1 j=i i.e., o autovalor λ está contido no disco dk e, como λ é arbitrário, então todos os autovalores de A devem estar contidos na união de todos os discos, D. ♦ O exemplo a seguir apresenta uma aplicação do teorema 6.2.1. Exemplo 6.3 A matriz −1 0 −1 4 −1 A = −1 −1 −2 10 tem como seus autovalores λ1 = 10, 3863, λ2 = 3, 8037 e λ3 = 0, 8100. Calculando os discos de Gerschgorin, temos: d1 = |z − 1| < |0| + | − 1| = 1 d2 d3 = |z − 4| < | − 1| + | − 1| = 2 = |z − 10| < | − 1| + | − 2| = 3 A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 113 Introdução ao Cálculo Numérico Autovalores e Autovetores Como todos os autovalores de A são reais, e observando (veja a figura 6.2) que em cada disco devemos ter um autovalor, podemos dizer que: • existe um autovalor, λ1 , que está dentro do disco centrado em 10 e raio 3 e, de fato, 7 < 10, 3863 < 13; • existe um autovalor, λ2 , que está dentro do disco centrado em 4 e raio 2 e, realmente, 2 < 3, 8037 < 6; • existe um autovalor, λ3 , que está dentro do disco centrado em 1 e raio 1 e, com efeito, 0 < 0, 81 < 2; A figura 6.2 ilustra esse resultado. Figura 6.2: Discos de Gerschgorin Uma conseqüência do teorema de Gerschgorin é a determinação do maior disco que contém todos os autovalores de A. Podemos obter, a partir dos discos, os extremos ao longo do eixo dos números reais, i.e. o intervalo [α, ω] tal que α = min{aii − ri }, i ω = max{aii + ri }, i i = 1, 2, . . . , n (6.11) e o maior disco é justamente aquele com centro (α + ω)/2 e raio (α + ω)/2. No caso em que todos os autovalores são reais, basta então considerar o intervalo [α, ω]. Teorema 6.2.2 Maior e menor autovalor: Seja A uma matriz real simétrica de ordem n, e x ∈ IC um vetor arbitrário. Então, λ1 = max x=0 xT Ax , xT x λn = min x=0 xT Ax xT x onde os autovalores são ordenados tais que λ1 ≥ λ2 ≥ . . . ≥ λn . A razão xT Ax , x = 0 (6.12) xT x é chamada de quociente de Rayleigh correspondente a x e, juntamente com o teorema 6.2.2, nos permitirá estimar de forma bastante rápida um autovalor de uma matriz simétrica, conforme veremos na seção 6.5. A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 114 Introdução ao Cálculo Numérico 6.3 Autovalores e Autovetores Cálculo de autovalores e autovetores via determinantes Já vimos que, por definição, os autovalores de uma matriz A são as raı́zes do polinômio caracterı́stico de A. Evidentemente, para matrizes de ordem n > 4, não é aconselhável que se utilize a equação (6.2) para se obter o polinômio caracterı́stico, por duas razões: 1. o cálculo de determinantes de ordem superior a 4 envolve considerável custo computacional; 2. o polinômio caracterı́stico de uma matriz grande pode ser instável numericamente. No entanto, algumas aplicações de engenharia, fı́sica e outros campos do conhecimento envolvem a determinação de autovalores de matrizes de ordem n = 2 ou n = 3 e, nesse caso, é possı́vel obterse os autovalores extraindo as raı́zes do polinômio caracterı́stico, conforme mostra o exemplo a seguir. Exemplo 6.4 Seja a matriz 2 5 3 −4 . O seu polinômio caracterı́stico é 2−λ 5 p(λ) = det(A − λI) = 3 −4 − λ = (2 − λ)(−4 − λ) − 15 ou p(λ) = λ2 + 2λ − 23, cujas raı́zes são λ1 = 3, 8990 e λ2 = −5, 8990. Para se determinar os autovetores, utiliza-se a equação (6.1) para cada autovalor λi , na forma (A − λi I)xi = 0, como segue: Exemplo 6.5 Calcule os autovetores do exemplo 6.4. Solução: Para o autovalor λ1 = 3, 8990, escrevemos (A − 3, 8990I)x1 7, 8990 5 (x1 )1 (x1 )2 3 1, 8990 de onde obtemos x1 = k −1, 5798k = = 0 0 0 , k = 0 O autovetor correspondente a λ2 = −5, 8990 é obtido de forma similar: (A + 5, 8990I)x1 = −1, 8990 5 (x2 )1 = (x2 )2 3 −7, 8990 de onde obtemos x2 = k 0, 3798k 0 0 0 , k = 0 Computacionalmente, no entanto, podemos estimar o autovetor correspondente a um autovalor utilizando os métodos da potência com translação da origem (seção 6.5.2) ou da iteração inversa com translação da origem (seção 6.5.3). A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 115 Introdução ao Cálculo Numérico 6.4 Autovalores e Autovetores Autovalores de uma matriz tridiagonal simétrica Em muitas aplicações surgem matrizes tridiagonais simétricas, das quais necessitamos extrair autovalores e/ou autovetores. Por exemplo, ao aproximarmos a equação diferencial parcial ∂2u ∂u = ∂t ∂x2 por diferenças finitas, obtemos uma matriz −2 1 0 0 1 −2 1 0 0 1 −2 1 0 0 1 −2 a qual apresenta aquela caracterı́stica. De forma geral, consideramos uma matriz a1 b 1 b 1 a2 T = 0 ... . . .. . . 0 ... T de ordem n, 0 .. . .. .. . . 0 ... 0 .. .. . . .. . 0 .. . bn−1 bn−1 an e chamamos de Tr a matriz principal de ordem r de T , i.e. a1 a1 b 1 T 1 = a1 , T 2 = , T3 = b1 b 1 a2 0 b1 a2 b2 (6.13) 0 b2 , a3 ... Escrevendo as equações caracterı́sticas p1 (λ), p2 (λ) e p3 (λ) das matrizes T1 , T2 e T3 , obtemos: p1 (λ) p2 (λ) p3 (λ) = det(T1 − λI) = a1 − λ (6.14) b21 b21 = = det(T2 − λI) = (a2 − λ)(a1 − λ) − = (a2 − λ)p1 (λ) − (6.15) 2 2 det(T3 − λI) = (a3 − λ) (a2 − λ)(a1 − λ) − b1 (a3 − λ) − b2 (a1 − λ) = = (a3 − λ)p2 (λ) − b22 p1 (λ) (6.16) de onde podemos escrever, generalizando para r, pr (λ) = (ar − λ)pr−1 (λ) − b2r−1 pr−2 (λ), r = 2, 3, . . . , n, p0 (λ) = 1, (6.17) A equação (6.17) nos permite avaliar o polinômio caracterı́stico da matriz T de forma bastante eficiente; no entanto, estamos preocupados em obter os autovalores de T , ou as raı́zes de pn . O teorema a seguir nos permitirá escrever um algoritmo bastante eficiente para se extrair alguns ou todos os autovalores de T . Teorema 6.4.1 Seqüência de Sturm: Se a matriz tridiagonal (6.13) é não-reduzı́vel 1 , então os r − 1 autovalores µ da matriz Tr−1 separam estritamente os r autovalores λ da matriz Tr : λr < µr−1 < λr−1 < µr−2 < · · · < λ2 < µ1 < λ1 . Mais ainda, se s(λ) representa o número de trocas de sinal na seqüência {p0 (λ), p1 (λ), . . . , pn (λ)} então s(λ) é igual ao número de autovalores de T menores do que λ, onde pr (λ) é dado por (6.17) e assume-se que pr (λ) tem o sinal oposto de pr−1 (λ) se pr (λ) = 0. 1 Uma matriz A é dita não-reduzı́vel se os elementos da diagonal da matriz triangular superior R, resultante de sua fatoração no produto QR, são todos não-nulos, onde Q é uma matriz ortogonal (i.e. QT Q = I). A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 116 Introdução ao Cálculo Numérico Autovalores e Autovetores O teorema 6.4.1 é extremamente importante: ele nos diz que, se tivermos os n autovalores µ de uma matriz triadiagonal Tn (de ordem n), então entre cada par de autovalores consecutivos µ (com exceção do menor e do maior), existe um e apenas um autovalor λ da matriz tridiagonal Tn+1 (de ordem n + 1), obtida acrescentando-se uma linha e uma coluna à matriz Tn . Devido à essa caracterı́stica, podemos utilizar o algoritmo da bissecção (ver algoritmo 2.2.1), juntamente com a equação (6.17), para obtermos rapidamente, e com segurança, um autovalor de Tn+1 , à partir de um intervalo que é um par de autovalores consecutivos de Tn . Para obter-se o menor e o maior autovalores de Tn+1 , utilizamos o teorema de Gerschgorin mais especificamente, calculamos o maior intervalo que engloba todos os autovalores, com a equação (6.11). Assim, o menor autovalor é calculado usando-se como estimativa inicial para o método da bissecção o intervalo [α, µr−1 ]; para o maior autovalor, utiliza-se o intervalo [µ1 , ω]. Os algoritmos 6.4.1, 6.4.2 e 6.4.3 combinam as idéias apresentadas acima. Da maneira como o algoritmo 6.4.3 é apresentado, todos os autovalores são obtidos; no entanto, simples modificações do mesmo nos permitem obter apenas alguns autovalores (por exemplo, o maior e o menor, ou os dois maiores, etc.). O exemplo 6.6 demonstra uma situação tı́pica, resolvido utilizando-se esses algoritmos. Algoritmo 6.4.1 Avalia polinômio caracterı́stico de uma matriz tridiagonal simétrica function pol carac trid(input: x, a, b; output: p) % a e b são os vetores contendo os elementos da % diagonal e subdiagonal, respectivamente, % da matriz tridiagonal p0 ← 1 p 1 ← a1 − x p ← p1 for r ← 2, 3, . . . , n do p ← (ar − x)p1 − b2r−1 p0 p0 ← p1 p1 ← p endfor endfunction A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 117 Introdução ao Cálculo Numérico Autovalores e Autovetores Algoritmo 6.4.2 Método da bissecção com polinômio caracterı́stico proc bissecção trid(input: a, b, α, β, kmax , δ, &; output: χ) % a e b são os vetores contendo os elementos da % diagonal e subdiagonal, respectivamente, % da matriz tridiagonal u ← pol carac trid(α, a, b) v ← pol carac trid(β, a, b) e←β−α if (sign(u) = sign(v)) then “não pode proceder” else k←1 w←1 while ((k ≤ kmax ) AND (| e | ≥ δ) AND (| w | ≥ &)) e ← e/2 χ←α+e w ← pol carac trid(χ, a, b) if (sign(w) = sign(u)) then β←χ v←w else α←χ u←w endif k ←k+1 endwhile endif endproc A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 118 Introdução ao Cálculo Numérico Autovalores e Autovetores Algoritmo 6.4.3 Autovalores simétrica de uma matriz tridiagonal proc autovalores tridiagonal(input: a, b, n; output: λ) % Calcula os raios dos discos de Gerschgorin, cada qual com centro a(i) r1 ← | b1 | for i ← 2, 3, . . . , n − 1 do ri ← | bi−1 | + | bi | endfor rn ← | bn−1 | % Calcula o intervalo [α, ω] na reta dos reais % contendo os autovalores α ← minni=1 (ai − ri ) ω ← maxni=1 (ai + ri ) % Calcula os autovalores, iniciando com o autovalor % de T1 = [a1 ], µ = a1 µ1 = a1 for i ← 2, 3, . . . , n do % Calcula os autovalores de Ti % a. entre α e µ1 call bissecção trid(a, b, α, µ1 , kmax , δ, &, λ1 ) % b. autovalores entre µ1 e µi−1 for j ← 1, 2, . . . , i − 2 do call bissecção trid(a, b, µj , µj+1 , kmax , δ, &, λj+1 ) endfor % c. entre µi−1 e ω call bissecção trid(a, b, µi−1 , ω, kmax , δ, &, λ1 ) µ←λ endfor λ←µ endproc Exemplo 6.6 Seja a matriz tridiagonal 2 −1 −1 2 −1 T = −1 2 −1 −1 2 a qual pode ser representada de forma compacta através dos vetores a = (2, 2, 2, 2)T e b = (−1, −1, −1)T Para se obter os autovalores de T , iniciamos com a matriz T 1 = a1 = 2 a qual tem como seu único autovalor µ1 = a1 = 2. Além disso, calculamos os extremos do intervalo de Gerschgorin, α = 0 e ω = 6, através da equação (6.11). Agora, precisamos calcular os dois autovalores λ da matriz a1 b 1 2 −1 = T2 = b 1 a2 −1 2 os quais, pelos teoremas de Gerschgorin e 6.4.1, satisfazem α < λ2 < µ1 < λ1 < ω A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 119 Introdução ao Cálculo Numérico Autovalores e Autovetores √ Estipulando-se como tolerâncias de convergência para o método da bisseção δ = & = ε e um máximo de 200 iterações, λ2 é obtido em 2 iterações utilizando-se como intervalo de busca [α, µ1 ] = [0, 2], resultando no valor λ2 = 1. O autovalor λ1 também é obtido em 2 iterações, usando-se como intervalo de busca [µ1 , β] = [2, 6], com o qual obtém-se λ1 = 3. Antes de procedermos ao cálculo dos autovalores de T3 , fazemos uma cópia dos λ, armazenandoos em µ; assim, temos µ1 = 3 e µ2 = 1. Procedemos, então, com o cálculo dos autovalores λ de T3 ; para tanto, utilizamos os intervalos de busca • [α, µ2 ] = [0, 1] para calcular o autovalor λ3 ; • [µ2 , µ1 ] = [1, 3] para calcular o autovalor λ2 ; • [µ1 , ω] = [3, 6] para calcular o autovalor λ1 . Os autovalores λ3 , λ2 e λ1 são então obtidos com o método da bissecção, utilizando-se as mesmas tolerâncias especificadas anteriormente, resultando em λ3 = 0, 5858, λ2 = 2 e λ1 = 3, 4142, obtidos em 27, 2 e 27 iterações respectivamente. Finalmente, basta calcularmos os autovalores de T4 . Procedendo de forma similar, fazemos µ3 = 0, 5858, µ2 = 2 e µ1 = 3, 4142 e estipulamos os intervalos de busca • [α, µ3 ] = [0, 0, 5858] para calcular o autovalor λ4 ; • [µ3 , µ2 ] = [0, 5858, 2] para calcular o autovalor λ3 ; • [µ2 , µ1 ] = [2, 3, 4142] para calcular o autovalor λ2 ; • [µ1 , ω] = [3, 4142, 6] para calcular o autovalor λ1 . de onde, após aplicarmos o algoritmo da bissecção a cada um desses intervalos, obtemos os autovalores de T4 ≡ T , λ4 = 0, 3820, λ3 = 1, 3820, λ2 = 2, 6180 e λ1 = 3, 6180, após 27 iterações (para todos os intervalos de busca). Note que não se obtém, com essa técnica, os autovetores correspondentes aos autovalores. Os métodos apresentados na seção a seguir podem ser utilizados para se obter esses autovetores. 6.5 Métodos para aproximação de autovalores e autovetores Em muitas aplicações, não é necessário obter-se todos os autovalores; é comum desejar-se, por exemplo, obter apenas o maior autovalor e seu correspondente autovetor. Os métodos apresentados nessa seção são indicados para o caso em que apenas um dos autovalores (e seu autovetor) necessita ser calculado. Particularmente, tais métodos são iterativos e apresentam boa eficiência quando a matriz em estudo é grande, esparsa e apresenta uma grande separação relativa entre o autovalor desejado e os demais autovalores. 6.5.1 Método da potência Seja A uma matriz de ordem n com autovalores λi tais que |λ1 | = |λ2 | = . . . = |λr | > |λr+1 | ≥ . . . ≥ |λn | (6.18) Nesse caso, diz-se que A apresenta r autovalores dominantes. Por hipótese, assumimos que existem n autovetores linearmente independentes xi , de onde qualquer vetor arbitrário z0 pode ser expresso como combinação linear desses autovetores, i.e. z0 = n αi xi (6.19) i=1 A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 120 Introdução ao Cálculo Numérico Autovalores e Autovetores Considere agora o método de aproximação sucessiva zk = Azk−1 , k = 1, 2, . . . (6.20) onde z0 é um valor inicial, dado. Usando as equações (6.1), (6.19) e escrevendo (6.20) em termos de z0 , temos zk = Azk−1 = A2 zk−2 = . . . = Ak z0 n = αi λki xi (6.21) i=1 Se pelo menos um dos α1 , α2 , . . ., αr não é nulo, então os termos correspondentes a eles, i.e. r k i=1 αi λi xi irão dominar o somatório da equação (6.21). Suponha, por exemplo, que temos um autovalor dominante, λ1 , de A. Considerando que α1 = 0, podemos reescrever (6.21) como k n λi zk = λk1 α1 x1 + αi xi λ 1 i=1 Note agora que, como λ1 > λ2 ≥ . . . ≥ λn , por hipótese, então os termos k λi λ1 tendem a zero à medida que k cresce. Daı́, podemos escrever zk = λk1 (α1 x1 + &k ) (6.22) onde &k é um vetor com elementos próximos a zero. O vetor zk tende, então, a aproximar o autovetor não-normalizado x1 . Essa equação nos permite escrever o assim chamado método da potência. Da equação (6.22), podemos escrever (α1 x1 + &k+1 ) zk+1 = λk+1 1 e, dividindo a i-ésima componente da equação acima pela componente correspondente de (6.22), obtemos α1 x1 + &k+1 (zk+1 )i = λ1 (6.23) → λ1 , quando k → ∞, i = 1, 2, . . . , n (zk )i α1 x1 + &k onde (zk )i indica o elemento i do vetor zk . A equação (6.23) nos diz que a taxa de convergência do método depende não só das constantes αi , mas principalmente das frações λ2 λ3 , , . . . , λn . λ1 λ1 λ1 Quanto menores forem esses frações, mais rápida é a convergência; por isso diz-se que o método da potência é eficiente – converge rapidamente para um autovalor – desde que este autovalor seja dominante, i.e., relativamente distante dos demais. De posse das equações (6.20) e (6.23), podemos escrever um algoritmo para o método da potência. Uma questão que se coloca é: quais valores iniciais, λ0 e z0 , devemos utilizar para o autovalor dominante e seu autovetor? Para z0 , consideraremos um vetor arbitrário, o qual será normalizado antes de se iniciar as iterações. Com essa escolha, valemo-nos da equação (6.1) e escrevemos Az0 z0T Az0 = λ0 z0 ... || z0 || = 1 ... = z0T λ0 z0 = λ0 (z0T z0 ) = λ0 Aplica-se, então, repetidamente a equação (6.20), normalizando o vetor zk a cada iteração, conforme mostrado no algoritmo 6.5.1. O exemplo 6.7 mostra o funcionamento do método. A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 121 Introdução ao Cálculo Numérico Autovalores e Autovetores Algoritmo 6.5.1 Método da potência proc potencia(input: A, z0 , &, kmax ; output: λk , zk ) z0 ← z0 /|| z0 || λ0 ← z0T Az0 for k = 1, 2, . . . , kmax do q ← Azk−1 zk ← q/|| q || λk ← zkT Azk if | λk − λk−1 | < & then break endif endfor endproc Exemplo 6.7 Seja a matriz 8 1 A = −1 5 0 1 2 1 , 90 a qual tem como autovalores e respectivos autovetores, λ1 = 90, 0115, λ2 = 7, 6308, λ3 = 5, 3577 0, 0245 −0, 9353 −0, 0043 x1 = 0, 0115 , x2 = 0, 3539 , x3 = −0, 0111 . 0, 9996 −0, 0043 0, 9996 Utilizando-se o método da potência com um vetor com três elementos escolhidos arbitrariamente e normalizado, z0 = (0, 4394, 0, 6415, 0, 6287)T , obtém-se a seguinte seqüência de valores, com uma tolerância para convergência de 10−5 : k 0 1 2 3 4 5 6 zk (0, 4394, 0, 6415, 0, 6287)T (0, 0940, 0, 0590, 0, 9938)T (0, 0313, 0, 0133, 0, 9994)T (0, 0251, 0, 0115, 0, 9996)T (0, 0246, 0, 0115, 0, 9996)T (0, 0245, 0, 0115, 0, 9996)T (0, 0245, 0, 0115, 0, 9996)T λk 40, 5408 89, 2834 89, 9939 90, 0102 90, 0114 90, 0115 90, 0115 onde pode-se verificar que z6 é uma boa aproximação para x1 , sujeita àquela tolerância. Caso a matriz tenha autovalores dominantes repetidos, i.e. λ1 = λ2 = . . . = λr o método da potência irá obter apenas um autovetor, o qual será combinação linear dos autovetores correspondentes a λ1 . O método da potência diverge se A tiver autovalores diferentes, porém de mesmo valor absoluto como, por exemplo, um par conjugado de autovalores dominantes complexos, λ2 = λ1 ; o exemplo a seguir ilustra tal situação: Exemplo 6.8 Seja a matriz 1 10 2 10 , A = −1 1 10 1 −13 A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 122 Introdução ao Cálculo Numérico Autovalores e Autovetores a qual tem como autovalores e respectivos autovetores, λ1 = −9, 4515 + 3, 8807i, λ2 = −9, 4515 − 3, 8807i, λ3 = 7, 9030 0, 0307 − 0, 4143i 0, 0307 + 0, 4143i 0, 7897 x1 = 0, 2160 + 0, 5518i , x2 = 0, 2160 − 0, 5518i , x3 = 0, 4651 . −0, 4368 − 0, 5343i −0, 4368 + 0, 5343i 0, 4000 Utilizando-se o método da potência com um vetor com três elementos escolhidos arbitrariamente e normalizado, z0 = (0, 4857, 0, 0197, 0, 8739)T , obtém-se a seguinte seqüência de valores, com uma tolerância para convergência de 10−5 : k 0 1 2 3 4 5 6 zk (0, 4857, 0, 0197, 0, 8739)T (0, 2253, 0, 7668, −0, 6010)T (0, 4829, −0, 3946, 0, 7817)T (−0, 2066, 0, 7546, −0, 6229)T (0, 5786, −0, 5001, 0, 6443)T (−0, 4517, 0, 7731, −0, 4454)T (0, 8581, −0, 4338, 0, 2749)T λk −4, 3241 −9, 1975 −8, 1345 −9, 4601 −6, 4876 −6, 2931 −1, 8886 a qual apresenta um comportamento não convergente. 6.5.2 O método da potência com translação da origem Como vimos na seção anterior, a convergência do método da potência depende de | λ2 /λ1 |, para o caso de existir apenas um autovalor dominante. Se essa razão for muito próxima de 1, então a convergência é muito lenta. No entanto, podemos obter a solução de forma mais rápida se procedermos a uma modificação do método da potência. Essa modificação baseia-se no fato de que, se λ é um autovalor de uma matriz A, então λ − σ é o autovalor correspondente da matriz A − σI. Dessa forma, se aplicarmos o método da potência a uma matriz A − σI, tal que λ1 − σ ainda seja dominante, a convergência do método dependerá de λ2 − σ λ1 − σ o que, para um valor adequado de σ, poderá ser menor do que | λ2 /λ1 |. A esse processo, dá-se o nome de translação da origem – é como se os autovalores estivessem distribuı́dos em um novo sistema de referência cuja origem é σ, e não mais zero – e pode ser bastante eficaz, desde que a escolha de σ seja criteriosa. Obviamente, poderı́amos calcular explicitamente a matriz A − σI e utilizar o algoritmo 6.5.1; no entanto, pequenas modificações naquele algoritmo nos permitem utilizar o processo de translação da origem de forma mais eficiente. Novamente, z0 é considerado um vetor unitário, de onde podemos obter a seguinte estimativa para λ0 : (A − σI)z0 z0T Az0 − σz0T z0 λ0 = λ0 z0 ; pré-multiplicando por z0T , = λ0 (z0T z0 ) ... || z0 || = 1 ... = z0T Az0 − σ e, por analogia, escrevemos λk = zkT Azk − σ Além disso, ao invés de calcularmos q = Azk , devemos calcular q = Azk − σzk . Ao final do processo, devemos corrigir λk , adicionando a ele σ. Essas idéias são apresentadas no algoritmo 6.5.2. A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 123 Introdução ao Cálculo Numérico Autovalores e Autovetores Algoritmo 6.5.2 Método da potência (com translação) proc potencia translação(input: A, z0 , σ, &, kmax ; output: λk , zk ) z0 ← z0 /|| z0 || λ0 ← z0T Az0 − σ for k = 1, 2, . . . , kmax do q ← Azk−1 − σzk−1 zk ← q/|| q || λk ← zkT Azk − σ if | λk − λk−1 | < & then break endif endfor λk ← λk + σ endproc O exemplo a seguir ilustra o uso do método da potência com translação de origem. Exemplo 6.9 Seja a matriz 2 −1 2 A = −1 0 −1 0 −1 , 2 a qual tem como autovalores e respectivos autovetores, λ1 = 3, 4142, λ2 = 2, λ3 = 0, 5858 0, 5000 0, 7071 0, 5000 x1 = −0, 7071 , x2 = 0, 0000 , x3 = 0, 7071 . 0, 5000 −0, 7071 0, 5000 Note que | λ2 /λ1 | = 0, 5858. Se utilizarmos o método da potência, a uma tolerância de 10−5 e vetor inicial z0 = (1, 0, 0)T , necessitaremos de 13 iterações para obter a aproximação 3, 4142 para o autovalor λ1 . No entanto, se usarmos o translação da origem, com σ = 1, necessitamos apenas de 9 iterações para obter a mesma aproximação; veja que λ2 − 1 1 = = 0, 4142 < 0, 5858 = λ2 λ1 λ1 − 1 2, 4142 o que sugere o menor número de iterações. 6.5.3 Método da iteração inversa Como vimos, o método da potência nos permite aproximar o autovalor dominante de A; suponha, agora, que desejamos aproximar o menor autovalor (e seu correspondente autovetor) de A. Relembrando que os autovalores de A−1 são o inverso dos autovalores de A (equação (6.3)), então, se utilizarmos o método da potência sobre a matriz A−1 , aproximaremos o menor autovalor de A pois ele é o maior autovalor de A−1 . A essa modificação do método da potência chamamos de método da iteração inversa. O método da iteração inversa procede, basicamente, com o cálculo sucessivo de vetores zk dados por zk = A−1 zk−1 , k = 1, 2, . . . mas já vimos (capı́tulo 4) que, computacionalmente, devemos evitar, se possı́vel, calcular a inversa de uma matriz. Nesse caso, é aconselhado que se resolva o sistema Azk = zk−1 , k = 1, 2, . . . A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 124 Introdução ao Cálculo Numérico Autovalores e Autovetores através da fatoração LU de A (seção 4.3.2), uma vez que várias iterações serão necessárias para se aproximar o menor autovalor e respectivo autovetor. Além disso, o método da iteração inversa é, normalmente, combinado com a translação de origem, o que resulta na fatoração LU da matriz A − σI. O algoritmo 6.5.3 apresenta o método da iteração inversa, incorporando translação de origem. Algoritmo 6.5.3 Método da iteração inversa (com translação) proc iteração inversa translação(input: A, z0 , σ, &, kmax ; output: λk , zk ) Fatore A − σI no produto LU z0 ← z0 /|| z0 || λ0 ← z0T Az0 − σ for k = 0, 1, . . . , kmax do Resolva o sistema Ly = zk Resolva o sistema U q = y zk ← q/|| q || λk ← zkT Azk − σ if | λk − λk−1 | < & then break endif endfor λk ← λk + σ endproc Exemplo 6.10 Seja a matriz do exemplo 6.9, 2 −1 2 A = −1 0 −1 0 −1 , 2 cujo menor autovalor é λ3 = 0, 5858 e o seu correspondente autvetor é x3 = (0, 5000, 0, 7071, 0, 5000)T . Se utilizarmos o algoritmo 6.5.3 com σ = 0, i.e. sem translação da origem, a uma tolerância de 10−5 e vetor inicial z0 = (1, 0, 0)T , necessitaremos de 7 iterações para obter a aproximação 0, 5858 para o autovalor λ3 e (0, 5002, 0, 7071, 0, 4998)T para o correspondente autovetor. No entanto, se usarmos o translação da origem, com σ = 0, 5, necessitamos apenas de 4 iterações para obter a mesma aproximação; veja que −1 −1 λ2 − 0, 5 0 = = 0 < 0, 2929 = λ2 λ−1 λ−1 − 0, 5 0, 0858 3 3 o que sugere o menor número de iterações; na verdade, σ = 0, 5 é a melhor escolha possı́vel, nesse caso. Outro exemplo mostra como usar o método da iteração inversa em conjunto com o teorema de Gerschgorin, a fim de se determinar um autovalor especı́fico. Exemplo 6.11 Seja a matriz 5 2 A= 2 3 1 1 1 1 1 Os discos de Gerschgorin são: d1 : c1 = 5, r1 = 3 d2 d3 : c2 = 3, : c3 = 1, r2 = 3 r3 = 2 A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 125 Introdução ao Cálculo Numérico Autovalores e Autovetores Suponha que desejamos aproximar o menor autovalor; como o disco d3 é aquele que se encontra mais à esquerda em comparação aos demais, podemos utilizar o seu centro como fator de transla√ −1 ção. Então, utilizamos o método da iteração inversa com σ = 1, vetor inicial z0 = 3 (1, 1, 1)T e tolerância 10−5 e obtemos λ = 0, 5764 e z = (−0, 0597, −0, 3380, 0, 9393)T após 10 iterações. Por outro lado, se tivéssemos utilizado σ = c3 + r3 , a convergência para o mesmo autovalor seria obtida em apenas 4 iterações. 6.5.4 O método da iteração inversa e o quociente de Rayleigh Como visto no teorema 6.2.2, o valor do autovalor dominante λ1 de uma matriz real simétrica é o máximo do quociente de Rayleigh, dentre todos os vetores x = 0. Isso nos permite utilizar a expressão (6.12) juntamente com o método da iteração inversa, conforme mostra o algoritmo 6.5.4. Algoritmo 6.5.4 Método da iteração inversa com translação via quociente de Rayleigh) proc iteração inversa translação(input: A, z0 , &, kmax ; output: λk , zk ) z0 ← z0 /|| z0 || for k = 0, 1, . . . , kmax do λk = T zk Azk Tz zk k Resolva o sistema (A − λk I)q = zk zk ← q/|| q || if || zk − zk−1 || < & then break endif endfor endproc Note que, no algoritmo 6.5.4, um sistema de equações diferente é resolvido a cada iteração, já que uma nova estimativa λk é utilizada a cada iteração. É possı́vel, no entanto, que o sistema A − λk I seja singular e, nesse caso, o processo deve ser terminado. 6.6 Exercı́cios Exercı́cio 6.1 Determine os autovalores e autovetores correspondentes da matriz 1 −1 −1 0 2 5 . 0 0 −1 Exercı́cio 6.2 Calcule o autovalor dominante de 10 9 8 3 5 6 . 7 2 −1 Exercı́cio 6.3 Calcule o autovalor dominante 6 2 1 e o autovetor correspondente da matriz 2 1 3 1 , 1 1 usando o método da potência, com z0 = (1, 1, 1)T e tolerância 10−5 . A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 126 Introdução ao Cálculo Numérico Autovalores e Autovetores Exercı́cio 6.4 Calcule o autovalor dominante e o autovetor correspondente da matriz 0 1 1 −1 0 1 , −1 −1 0 usando o método da potência, com z0 = (1, 1, 1)T e tolerância 10−5 . Exercı́cio 6.5 Explique o que acontece com o método da potência para a matriz 2 1 −1 1 2 −1 , 1 −1 2 com z0 = (1, 1, 1)T e tolerância 10−5 , sabendo que os seus autovalores são 1, 2 e 3. Repita para z0 = (1, 0, 10−6)T . Exercı́cio 6.6 Utilize o método da iteração inversa para calcular o menor autovalor da matriz 2 1 −1 1 2 −1 , 1 −1 2 com z0 = (1, 1, 1)T e tolerância 10−5 . Exercı́cio 6.7 Seja a matriz 5 2 A= 2 3 1 1 1 1 1 Explique o que ocorre com o método da iteração inversa utilizado com σ = 3, z0 = (1, 0, 0)T e tolerância 10−5 . Generalize a sua resposta. A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 127