Métodos de álgebra linear computacional , Notas de estudo de Matemática Computacional

Baixe Métodos de álgebra linear computacional e outras Notas de estudo em PDF para Matemática Computacional, somente na Docsity!

M´etodos de Algebra Linear Computacional em´

Scilab

Vinicius A. de Souza

Universidade Federal de S˜ao Paulo - UNIFESP

[email protected]

4 de novembro de 2012

Sum´ario

• 1 Introdu¸c˜ao
• 2 Problemas de M´ınimo Quadrado
  • 2.1 Transforma¸c˜oes de Householder
  • 2.2 Rota¸c˜oes de Givens
  • 2.3 Decomposi¸c˜ao QR
• 3 Autovalores e Autovetores
  • 3.1 M´etodo de Leverrier
  • 3.2 M´etodo de Leverrier-Faddeev
  • 3.3 M´etodo das Potˆencias
  • 3.4 M´etodo da Potˆencia Inversa
• 4 Autovalores de Matrizes Sim´etricas
  • 4.1 Decomposi¸c˜ao Sim´etrica de Schur
  • 4.2 M´etodo de Jacobi
  • 4.3 M´etodo de Jacobi C´ıclico
  • 4.4 M´etodo de Rutishauser - LR
  • 4.5 M´etodo de Francis - QR

e1(1,1) = 1

//calcula norma 2 de x norma2 = sqrt(X’*X)

//calcula vetor v v = X + ( (X(1,1)/abs(X(1,1)) * norma2) * e1 )

//calculando a matriz P I = eye(linhas, linhas) P = I - ( (2/(v’v)) * (vv’))

endfunction

Vamos verificar agora um exemplo de utiliza¸c˜ao da rotina. Queremos encon- trar a matriz de Householder que zera os elementos abaixo da primeira posicao do vetor

x =

Podemos fazer isso com o c´odigo

x = [4;4;2] ini = 1 P = Householder(x, ini)

Calculando o produto P x, o resultado ser´a P x =

2.2 Rota¸c˜oes de Givens

A matriz de rota¸c˜ao de Givens para um vetor x ∈ R^2 ´e dada por

R(θ) =

[

cos θ − sin θ sin θ cos θ

]

que rotaciona qualquer vetor x ∈ R^2 no sentido anti-hor´ario por um ˆangulo θ. A matriz de rota¸c˜ao ´e ortogonal, portanto, podemos us´a-la para zerar os elementos abaixo da diagonal.

A fun¸c˜ao abaixo recebe como entrada um vetor coluna x e a posi¸c˜ao ini a partir da qual todos os elementos devem ser zerados. Retorna a matriz de rota¸c˜ao de Givens R.

function R = RotacaoGivens(x,ini)

// extrai as dimensoes [l, c] = size(x)

// verifica qual a posicao para calculo if ini == 1 then X = x end

// se inicio nao for na posicao 1, criamos um novo vetor // comecando da posicao 2 if ini>1 then

j= for i=ini:l

X(j,1) = x(i,1) j=j+

end end

//extrai novas dimensoes volta inicio para 1 [l,c] = size(X) ini = 1

R = eye(l,l)

for j=1:l-

x_i = X(ini,1) y_j = X(j+1,1)

cosseno = x_i / (sqrt( (x_i^2) + (y_j^2) )) seno = -y_j / (sqrt( (x_i^2) + (y_j^2) ))

// monta matriz de givens R(i,j) Rt = eye(l,l)

Rt(ini,ini) = cosseno Rt(ini,j+1) = -seno Rt(j+1,ini) = seno Rt(j+1,j+1) = cosseno

R = Rt*R //realiza produto das matrizes R(i,j) geradas a cada passo

X = Rt*X

end

endfunction

Vamos verificar agora um exemplo de utiliza¸c˜ao da rotina. Queremos encon- trar a matriz de Rota¸c˜ao de Givens que zera os elementos abaixo da primeira

// H = RotacaoGivens(coluna,i)

// completa matriz se necessario [lR, cR] = size(H) I = eye(l,l) I((l-lR+1):c, (l-lR+1):l) = H //completa com identidade

Qp = I // armazena Q da iterao

Q = Q*Qp’ //calculando matriz ortogonal Q

A = Qp*A //calculando matriz triangular R

end

R = A

endfunction

Vamos verificar agora um exemplo de utiliza¸c˜ao da rotina. Queremos encon- trar a decomposi¸c˜ao QR de uma matriz A dada por

A =

Podemos fazer isso com o c´odigo

A=[1,5,4;2,2,5;2,3,6] [Q,R] = DecompQR(A)

Calculando o produto QR, o resultado ser´a exatamente QR = A.

3 Autovalores e Autovetores

Vamos primeiramente relembrar algumas defini¸c˜oes. Dizemos que λ ∈ R ´e um autovalor de A ∈ An,n^ se existir um vetor v n˜ao-nulo tal que Av = λv. O vetor v ´e denominado autovetor associado ao autovalor λ. O polinˆomio caracter´ıstico da matriz A ´e o polinˆomio p(λ) = det(A−λI). As ra´ızes do polinˆomio caracter´ıstico s˜ao autovalores de A.

3.1 M´etodo de Leverrier

Vamos usar a seguinte nota¸c˜ao para o polinˆomio caracter´ıstico:

p(λ) = (−1)n(λn^ − p 1 λn−^1 − · · · − pn− 1 λ − pn).

O m´etodo de Leverrier determina o polinˆomio caracter´ıstico atrav´es das rela¸c˜oes

p 1 = S 1

2 p 2 = S 2 − p 1 S 1 .. . npn = Sn − p 1 Sn− 1 − p 2 Sn− 2 − ... − pn− 1 S 1

onde Sk = λ 1 k^ + ... + λnk^ = traco(Ak) k = 1, 2 , · · · , n. Sabemos que a multiplica¸c˜ao de matrizes que teremos que desenvolver no c´alculo de Ak^ ´e um procedimento bastante custoso. Por isso, iremos verificar agora o m´etodo de Leverrier-Faddeev.

3.2 M´etodo de Leverrier-Faddeev

Sejam A ∈ Rn,n^ e tr(A) o tra¸co da matriz A. Considere a sequˆencia de matrizes e valores

A 1 = A qk = tr(Ak)/k (k = 1, 2 , ..., n) Bk = Ak − qkI Ak = ABk− 1 (k = 2, 3 , ..., n) Teorema: Os termos qk definidos anteriormente s˜ao os coeficientes do po- linˆomio caracter´ıstico de A, ou seja, qk = pk, para k = 1, 2 , ..., n.

Corol´ario 1: Bn = 0 (matriz nula).

Corol´ario 2: Se A ´e n ao-singular ent˜ao A−^1 = (1/pn)Bn− 1.

Corol´ario 3: Cada coluna n˜ao nula da matriz Qk = λkn−^1 I + λkn−^2 B 1 + λkBn− 2 + Bn− 1 ´e um autovetor de A correspondente ao autovalor λk.

Apresentamos ent˜ao a fun¸c˜ao abaixo. Ela encontra o polinˆomio carac- ter´ıstico, os autovalores e os autovetores de uma matriz quadrada A passada como argumento.

function [p,autovetores,autovalores] = LeverrierFaddeev(A)

[l,c] = size(A)

I = eye(l,c) B_i = eye(l,c)

C = 0 //armazena valores das matrizes Bk autovetores = 0 //armazena autovetores associados coef = 0 //armazena valores dos coeficientes p coef(1,c+1) = 1

col=

trar os autovalores e autovetores da matriz A dada por

A =

Podemos fazer isso com o c´odigo

A=[ 1, 1, -1; 0, 0, 1; -1, 1, 0]

[polinomio, autovetores,autovalores] = LeverrierFaddeev(A)

Teremos os seguintes resultados:

autovalores =

autovetores =

1. 1. 1.332D-

• 1.4142136 - 1.4142136 - 1. polinomio =

2 3

• 2 + 2x + x - x

3.3 M´etodo das Potˆencias

Permite encontrar o autovalor de maior m´odulo, desde que os autovetores sejam linearmente independentes. Este m´etodo ´e baseado no seguinte teorema:

Teorema: Sejam A ∈ Rn,n, λ 1 , λ 2 , ..., λn seus autovalores e u 1 , u 2 , ..., un seus autovetores. Se os autovetores forem LI e |λ 1 | > |λ 2 | ≥ |λ 3 | ≥ ... ≥ |λn| e a sequˆencia yk+1 = Ayk, para k = 0, 1 , 2 , ..., sendo y 0 um vetor arbitr´ario que pode ser escrito como y 0 = c 1 u 1 + ... + cnun, com ci escalares e c 1 ̸= 0, ent˜ao

lim x→ 0

(yk+1)r (yk)r

= λ 1 ,

onde ()r significa a r-´esima componente do vetor. Al´em disso, (^) λy 1 kk −→ c 1 u 1.

Observa¸c˜ao 1 : Caso |λ 1 | < 1 ou |λ 1 | > 1, podemos ter problemas num´ericos. Assim, o seguinte algoritmo ´e utilizado:

zk+1 = Ayk (k = 0, 1 , ...) yk+1 = (1/αk+1)zk+1 (k = 0, 1 , ...)

com αk+1 = ||zk+1||∞

Observa¸c˜ao 2 : Se em algum momento zk ou yk forem iguais a 0, o m´etodo falhou, provavelmente porque os autovetores da matriz n˜ao s˜ao LI.

A fun¸c˜ao abaixo retorna o maior autovalor de uma matriz A juntamente com seu autovetor associado.

function [autovalor, autovetor] = Potencia(A)

[linhas,colunas] = size(A)

y0 = zeros(linhas,1) + 1 //aproximacao inicial para autovalor y_i = y

z1 = A*y z_i = z

lambda_ant = [500;500;500] //aproximao inicial para possiveis //autovetores

precisao = 10^-2 //ajusta a preciso requerida

flag = 1 //determina encerramento dos calculos

while flag == 1

alpha_i = max(abs(z_i))

y_i = (1/alpha_i)*z_i

z_i = A*y_i

//aproximacao para lambda for i=1:linhas lambda_atual(i,1) = z_i(i,1)/y_i(i,1) end

//erro relativo for i=1:linhas erro(i,1) = abs(lambda_atual(i,1) - lambda_ant(i,1)) / abs(lambda_atual(i,1)) end

//verifica se algum erro eh menor do que precisao requerida menorErro = erro(1,1) //supoe menor indice = 1

for i=1:linhas //busca menor erro de todos if erro(i,1) < menorErro then

sendo que (^) |λ^1 i| ´e autovalor de A−^1.

Assim, basta aplicarmos o m´etodo da Potˆencia `a matriz inversa de A para calcularmos o menor autovalor de A.

4 Autovalores de Matrizes Sim´etricas

Lembrete: Uma matriz A ´e sim´etrica se AT^ = A

4.1 Decomposi¸c˜ao Sim´etrica de Schur

Os m´etodos que iremos estudar buscam uma aproxima¸c˜ao para a decomposi¸c˜ao sim´etrica de Schur, que se baseia no seguinte teorema:

Teorema: Seja A ∈ R^2 uma matriz sim´etrica e λ 1 , λ 2 , ..., λn autovalores de A. Ent˜ao, existe uma matriz ortogonal Q ∈ Rn,n^ tal que QT^ AQ = Λ, sendo Λ = diag(λ 1 , λ 2 , · · · , λn).

4.2 M´etodo de Jacobi

A ideia central do m´etodo de Jacobi ´e aplicar transforma¸c˜oes ortogonais na matriz de modo a torn´a-la mais pr´oxima de uma matriz diagonal. Da´ı, queremos encontrar uma matriz ortogonal Q de modo que QT^ AQ diminua o valor dado por

ndiag(A) =

v u u t

∑^ n

i=

∑^ n

j=1,i̸ =j

a^2 ij

A matriz de rota¸c˜ao de Jacobi J(p, q, θ) ´e igual ´a matriz de Givens. A fun¸c˜ao a seguir recebe uma matriz sim´etrica A e um valor real tolerancia e retorna os autovalores e autovetores desta matriz. o Valor dado pela tolerancia indica o qu˜ao pr´oximo queremos estar da decomposi¸c˜ao sim´etrica de Schur.

function [autovalores, autovetores] = JacobiClassico(A, tolerancia)

[linhas, colunas] = size(A)

Q = eye(linhas,linhas) J = eye(linhas,linhas)

//norma F da matriz A NF = norm(A,’fro’)

eps = tolerancia*NF

ndiag = Ndiag(A)

while ndiag>eps

//escolhe p,q [p,q] = buscaIndicesPeQ(A)

//calcula ’c’ e ’s’

if A(p,q) == 0 then mi = 1 else mi = (A(q,q)-A(p,p)) / (2*A(p,q)) end

t1 = -mi + sqrt(1+(mi^2)) t2 = -mi - sqrt(1+(mi^2))

//recebe menor raiz m modulo t = min( list(abs(t1), abs(t2)) )

c = 1/sqrt(1+t^2) s = t*c

J(p,p) = c J(p,q) = s J(q,p) = -s J(q,q) = c

Q = QJ A = J’A*J

//verifica valor novamente ndiag = Ndiag(A)

J = eye(linhas,linhas)

end //while

autovalores = diag(A) autovetores = Q

endfunction

Ainda temos duas fun¸c˜oes auxiliares:

//funcao para calculo do valor ndiag function valor = Ndiag(A)

aux = 0

for i=1:linhas for j=1:linhas

[autovalores, autovetores] = JacobiClassico(A, tolerancia)

Obtemos como resposta

autovetores =

• 0.7116322 0.4081016 - 0. 0.7025375 0.4080968 - 0.
• 0.0045490 - 0.8166456 - 0. autovalores =

Os resultados nos indicam que temos os autovalores λ 1 = 6, λ 2 = − 11. 99 e λ 3 = 17.99, com os respectivos autovalores v 1 = (− 0. 7 0. 7 0)T^ , v 2 = (0. 4 0. 4 − 0 .8)T^ e v 3 = (− 0. 5 − 0. 5 − 0 .5)T^.

4.3 M´etodo de Jacobi C´ıclico

O m´etodo de Jacobi cl´assico possui o problema de que a busca pelo elemento de maior valor absoluto fora da diagonal (func˜ao auxiliar buscaIndicesP eQ(A)) requer O(n^2 ) opera¸c˜oes. Por isso, o algoritmo C´ıclico de Jacobi elimina esta busca, diminuindo assim o total de opera¸c˜oes. Assim, temos a fun¸c˜ao abaixo:

function [autovalores, autovetores] = JacobiCiclico(A, tolerancia)

[linhas, colunas] = size(A)

Q = eye(linhas,linhas) J = eye(linhas,linhas)

//norma F da matriz A NF = norm(A,’fro’)

eps = tolerancia*NF ndiag = Ndiag(A)

while ndiag>eps

for p=1:(linhas-1)

for q = (p+1):linhas

//calcula ’c’ e ’s’ if A(p,q) == 0 then mi = 1 else mi = (A(q,q)-A(p,p)) / (2*A(p,q)) end

t1 = -mi + sqrt(1+(mi^2)) t2 = -mi - sqrt(1+(mi^2))

//recebe menor raiz em modulo t = min( list(abs(t1), abs(t2)) )

c = 1/sqrt(1+t^2) s = t*c

J(p,p) = c J(p,q) = s J(q,p) = -s J(q,q) = c

Q = QJ A = J’A*J

//verifica ndiag novamente ndiag = Ndiag(A) J = eye(linhas,linhas)

end end

end //while

autovalores = diag(A) autovetores = Q

endfunction

Veja que esta fun¸c˜ao agora somente utiliza a fun¸c˜ao auxiliar N diag(A). N˜ao se esque¸ca de inclu´ı-la em seu c´odigo.

4.4 M´etodo de Rutishauser - LR

O m´etodo de Rutishauser ´e baseado no seguinte teorema:

Teorema: Seja uma matriz A ∈ Rn,n^ sim´etrica com autovalores λ 1 , · · · , λn tais que |λ 1 | > |λ 2 | > · · · > |λn| e a sequˆencia Ak definida por A 1 = A e Ak+1 = RkLk, para k = 1, 2 , ..., sendo que Ak = LkRk ´e a decomposi¸c˜ao de Ak no produto de matrizes triangular inferior Lk e triangular superior Rk. Ent˜ao,

lim k→∞

Ak = S,

sendo S uma matriz triangular superior cuja diagonal ´e formada por auto- valores de A.

Note que o teorema acima faz uso da decomposi¸c˜ao LU de uma matriz. A fun¸c˜ao abaixo recebe uma matriz A sim´etrica e uma precisao real e retorna os

end

cont = cont+

end

U = A(:,1:colunas)

endfunction

//funcao para verificar se elementos abaixo da diagonal //principal de uma matriz sao menores que // uma dada precisao ’p’ function indicador = verifica(A,p)

[l,c] = size(A)

indicador = 0 //supoe verdadeira

for i=1:l for j=1:c if i>j then //para elementos abaixo da diag

if A(i,j) > p then

indicador = 1 //nao atingiu precisao

end end end end

endfunction

Como exemplo, vamos usar o m´etodo de Rutishauser para encontrar os au- tovalores da matriz A indicada abaixo, com precis˜ao de 10−^2.

A =

Fazemos isto com o c´odigo

A = [2, 0, 1; 0, 1, 0; 1, 0, 1]

precisao = 0.

autovalores = Rutishauser(A, precisao)

Obtemos o seguinte resultado

autovalores =

O resultado nos mostra que os autovalores da matriz A s˜ao λ 1 = 2.6, λ 2 = 1 e λ 3 = 0.3.

4.5 M´etodo de Francis - QR

O m´etodo de Francis trabalha de forma parecida ao m´etodo de Rutishauser, com a diferen¸ca de que ´e utilizada a decomposic¸c˜ao QR, ao inv´es da decomposi¸c˜ao LU.

Teorema: Seja uma matriz A ∈ Rn,n^ sim´etrica com autovalores λ 1 , · · · , λn tais que |λ 1 | > |λ 2 | > · · · > |λn| e a sequˆencia Ak definida por A 1 = A e Ak+1 = RkQk, para k = 1, 2 , ..., sendo que Ak = QkRk ´e a decomposi¸c˜ao de Ak no produto de matrizes ortogonal Qk e triangular superior Rk. Ent˜ao,

lim k→∞ Ak = S,

sendo S uma matriz triangular superior cuja diagonal ´e formada por auto- valores de A.

A fun¸c˜ao abaixo calcula os autovalores de uma matriz A pelo m´etodo de Francis.

function autovalores = Francis(A, precisao)

flag = 1 A1 = A

while flag == 1

[Q,R] = qr(A) //decomposicao QR da matriz A

A = R*Q

//verificamos se elementos abaixo da diagonal //sao menores quea precisao flag = verifica(A,precisao)

end

autovalores = diag(A)

endfunction

Esta fun¸c˜ao utiliza a fun¸c˜ao auxiliar verif ica() descrita anteriormente e a fun¸c˜ao qr() do pr´oprio Scilab.