Métodos de Álgebra Linear Computacional em
Scilab
Vinicius A. de Souza
Universidade Federal de São Paulo - UNIFESP
[email protected]
4 de novembro de 2012
Sumário
1 Introdução
2
2 Problemas de Mı́nimo Quadrado
2.1 Transformações de Householder . . . . . . . . . . . . . . . . . . .
2.2 Rotações de Givens . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Decomposição QR . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2
3
5
3 Autovalores e Autovetores
3.1 Método de Leverrier . . . . .
3.2 Método de Leverrier-Faddeev
3.3 Método das Potências . . . .
3.4 Método da Potência Inversa .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
6
6
7
9
11
4 Autovalores de Matrizes Simétricas
4.1 Decomposição Simétrica de Schur .
4.2 Método de Jacobi . . . . . . . . . .
4.3 Método de Jacobi Cı́clico . . . . .
4.4 Método de Rutishauser - LR . . .
4.5 Método de Francis - QR . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
12
15
16
19
.
.
.
.
1
.
.
.
.
1
Introdução
O objetivo deste trabalho é apresentar rotinas desenvolvidas a partir de métodos
de álgebra linear computacional. Em um primeiro momento, serão destacados
apenas os programas computacionais implementados, não dando muita ênfase
na teoria envolvida. A rotinas foram desenvolvidas em Scilab e os fundamentos
teóricos foram extraı́dos das notas de aula do professor Erwin Doescher referentes à disciplina de Álgebra Linear Computacional do curso de Matemática
Computacional da UNIFESP - SJC.
2
2.1
Problemas de Mı́nimo Quadrado
Transformações de Householder
Queremos encontrar matrizes ortogonais Q1 , Q2 , ..., Qn tais que Qn Qn−1 ...Q2 Q1 A
seja uma matriz triangular superior. Assim, a matriz Qk deve zerar os elementos abaixo da diagonal principal na coluna k, sem alterar os zeros das colunas
1, 2, ..., k − 1.
A função abaixo recebe como entrada um vetor coluna x e a posição ini
a partir da qual todos os elementos devem ser zerados. Retorna a matriz de
Householder P .
function P = Householder(x, ini)
// extrai dimensoes
[linhas, colunas] = 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=1
for i=ini:linhas
X(j,1) = x(i,1)
j=j+1
end
end
//extrai novas dimensoes
[linhas, colunas] = size(X)
e1 = zeros(linhas, colunas)
2
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)) * (v*v’))
endfunction
Vamos verificar agora um exemplo de utilização da rotina. Queremos encontrar a matriz de Householder que zera os elementos abaixo da primeira posicao
do vetor


4
x= 4 
2
Podemos fazer isso com o código
x = [4;4;2]
ini = 1
P = Householder(x, ini)


6
Calculando o produto P x, o resultado será P x =  0 .
0
2.2
Rotações de Givens
A matriz de rotação de Givens para um vetor x ∈ R2 é dada por
[
]
cos θ − sin θ
R(θ) =
,
sin θ
cos θ
que rotaciona qualquer vetor x ∈ R2 no sentido anti-horário por um ângulo
θ. A matriz de rotação é ortogonal, portanto, podemos usá-la para zerar os
elementos abaixo da diagonal.
A função abaixo recebe como entrada um vetor coluna x e a posição ini
a partir da qual todos os elementos devem ser zerados. Retorna a matriz de
rotação de Givens R.
function R = RotacaoGivens(x,ini)
// extrai as dimensoes
[l, c] = size(x)
3
// 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=1
for i=ini:l
X(j,1) = x(i,1)
j=j+1
end
end
//extrai novas dimensoes volta inicio para 1
[l,c] = size(X)
ini = 1
R = eye(l,l)
for j=1:l-1
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)
Rt(ini,j+1)
Rt(j+1,ini)
Rt(j+1,j+1)
=
=
=
=
cosseno
-seno
seno
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 utilização da rotina. Queremos encontrar a matriz de Rotação de Givens que zera os elementos abaixo da primeira
4
posicao do vetor


4
x= 4 
2
Podemos fazer isso com o código
x = [4;4;2]
ini = 1
G = RotacaoGivens(x, ini)


6
Calculando o produto Gx, o resultado será Gx =  0 .
0
2.3
Decomposição QR
O seguinte teorema nos permitirá realizar a decomposição de uma matriz A no
produto A = QR, onde Q e R são matrizes especı́ficas. Vejamos:
Teorema 1: Seja A ∈ Rm,n , m ≥ n, e rank(A) = n. Então, existe uma
matriz ortogonal Q ∈ Rm,n e uma matriz triangular superior R ∈ Rn,n tal que
A = QR.
A ideia é encontrar Q1 , Q2 , ..., Qn matrizes ortogonais, que podemos obter
utilizando as funções Householder ou Givens implementadas anteriormente, que
zeram os elementos abaixo da diagonal principal. Dessa forma, temos R =
Qn Qn−1 ...Q2 Q1 , onde R é uma matriz triangular superior. Assim,
A = Q1 T Q2 T ...Qn T R
Definindo Q = Q1 T Q2 T ...Qn , temos A = Q.R = QR.
A função abaixo recebe uma matriz quadrada A e retorna as matrizes Q e
R da decomposição.
function [Q,R] = DecompQR(A)
[l,c] = size(A)
Q = eye(l,l)
for i=1:c-1
coluna = 0
coluna = A(:,i)
// aqui chamamos a funcao que faz a rotacao de Householder
H = Householder(coluna,i)
// voce tambem poderia usar a rotacao de Givens
5
// 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 utilização da rotina. Queremos encontrar a decomposição QR de uma matriz A dada por


1 5 4
A= 2 2 5 
2 3 6
Podemos fazer isso com o código
A=[1,5,4;2,2,5;2,3,6]
[Q,R] = DecompQR(A)
Calculando o produto QR, o resultado será exatamente QR = A.
3
Autovalores e Autovetores
Vamos primeiramente relembrar algumas definições. Dizemos que λ ∈ R é um
autovalor de A ∈ An,n se existir um vetor v não-nulo tal que Av = λv. O vetor v
é denominado autovetor associado ao autovalor λ. O polinômio caracterı́stico da
matriz A é o polinômio p(λ) = det(A−λI). As raı́zes do polinômio caracterı́stico
são autovalores de A.
3.1
Método de Leverrier
Vamos usar a seguinte notação para o polinômio caracterı́stico:
p(λ) = (−1)n (λn − p1 λn−1 − · · · − pn−1 λ − pn ).
O método de Leverrier determina o polinômio caracterı́stico através das
relações
p1 = S1
6
2p2 = S2 − p1 S1
..
.
npn = Sn − p1 Sn−1 − p2 Sn−2 − ... − pn−1 S1
onde
Sk = λ1 k + ... + λn k = traco(Ak )
k = 1, 2, · · · , n.
Sabemos que a multiplicação de matrizes que teremos que desenvolver no
cálculo de Ak é um procedimento bastante custoso. Por isso, iremos verificar
agora o método de Leverrier-Faddeev.
3.2
Método de Leverrier-Faddeev
Sejam A ∈ Rn,n e tr(A) o traço da matriz A. Considere a sequência de matrizes
e valores
A1 = A
qk = tr(Ak )/k
(k = 1, 2, ..., n)
Bk = Ak − qk I
Ak = ABk−1
(k = 2, 3, ..., n)
Teorema: Os termos qk definidos anteriormente são os coeficientes do polinômio caracterı́stico de A, ou seja, qk = pk , para k = 1, 2, ..., n.
Corolário 1: Bn = 0 (matriz nula).
Corolário 2: Se A é n ao-singular então A−1 = (1/pn )Bn−1 .
Corolário 3: Cada coluna não nula da matriz Qk = λk n−1 I + λk n−2 B1 +
λk Bn−2 + Bn−1 é um autovetor de A correspondente ao autovalor λk .
Apresentamos então a função abaixo. Ela encontra o polinômio caracterı́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=1
7
for i=1:c
A_i = A*B_i
q_i = (trace(A_i) / i)
coef(1,c+1-i) = (-1*q_i) //armazena coeficientes do polinomio
//caracteristico reversamente
B_i = A_i - (q_i*I)
[lB,cB] = size(B_i)
C(1:lB,col:col+cB-1) = B_i //armazena matrizes B_k
col = col+cB
end
//gera polinomio caracteristico
p = poly(coef,"x","coeff")
p = ((-1)^c) * p
//calcula autovalores, que sao as raizes do polinomio caracteristico
autovalores = roots(p)
//daqui para frente, calcula autovetores aasociados
[num_autovalores,lixo] = size(autovalores)
for i=1 : num_autovalores
u0 = eye(l,1)
u_j=0
pulo=1 //selecionar primeira coluna da matriz B correspondente
for j=1:l-1
u_j = (autovalores(i,1)*u0) + C(1:l, pulo:pulo)
u0 = u_j
pulo = pulo+c //pula para primeira coluna da proxima matriz B
end
//armazena na coluna i o autovetor associado ao i-simo autovalor
autovetores(1:l,i:i) = u_j
end
endfunction
Vamos verificar agora um exemplo de utilização da rotina. Queremos encon8
trar os autovalores e autovetores da matriz A dada por


1 1 −1
1 
A= 0 0
−1 1
0
Podemos fazer isso com o código
A=[ 1, 1, -1;
0, 0, 1;
-1, 1, 0]
[polinomio, autovetores,autovalores] = LeverrierFaddeev(A)
Teremos os seguintes resultados:
autovalores
=
- 1.4142136
1.4142136
1.
autovetores =
1.
- 1.
1.4142136
polinomio =
1.
- 1.
- 1.4142136
1.332D-15
- 1.
- 1.
2
3
- 2 + 2x + x - x
3.3
Método das Potências
Permite encontrar o autovalor de maior módulo, desde que os autovetores sejam
linearmente independentes. Este método é baseado no seguinte teorema:
Teorema: Sejam A ∈ Rn,n , λ1 , λ2 , ..., λn seus autovalores e u1 , u2 , ..., un
seus autovetores. Se os autovetores forem LI e |λ1 | > |λ2 | ≥ |λ3 | ≥ ... ≥ |λn | e
a sequência yk+1 = Ayk , para k = 0, 1, 2, ..., sendo y0 um vetor arbitrário que
pode ser escrito como y0 = c1 u1 + ... + cn un , com ci escalares e c1 ̸= 0, então
lim
x→0
(yk+1 )r
= λ1 ,
(yk )r
onde ()r significa a r-ésima componente do vetor. Além disso,
yk
λ1 k
−→ c1 u1 .
Observação 1 : Caso |λ1 | < 1 ou |λ1 | > 1, podemos ter problemas numéricos.
Assim, o seguinte algoritmo é utilizado:
zk+1 = Ayk
(k = 0, 1, ...)
yk+1 = (1/αk+1 )zk+1
9
(k = 0, 1, ...)
com αk+1 = ||zk+1 ||∞
Observação 2 : Se em algum momento zk ou yk forem iguais a 0, o método
falhou, provavelmente porque os autovetores da matriz não são LI.
A função 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 = y0
z1 = A*y0
z_i = z1
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
10
menorErro = erro(i,1)
indice = i
end
end
//verifica se menor erro encontrado eh, por sua vez,
//menor que precisao
if menorErro < precisao then
flag = 0 //encerra
autovalor = lambda_atual(indice,1)
autovetor = y_i
end
//atualiza erro
lambda_ant = lambda_atual
end
endfunction
Vamos agora utilizar o método das Potências para calcular o maior autovalor
absoluto da matriz A dada por


3 0 1
A= 2 2 2 
4 2 5
Podemos utilizar a função acima fazendo
A=[3,0,1;
2,2,2;
4,2,5]
[autovalor, autovetor] = Potencia(A)
O resultado será
autovetor
=
0.2771084
0.5060241
1.
autovalor =
7.047619
3.4
Método da Potência Inversa
Observando que se |λ1 | ≥ |λ2 | ≥ · · · ≥ |λn−1 | > |λn | são autovetores de A, então
1
1
1
1
≤
≤ ··· ≤
<
,
|λ1 |
|λ2 |
|λn−1 |
|λn |
11
sendo que
1
|λi |
é autovalor de A−1 .
Assim, basta aplicarmos o método da Potência à matriz inversa de A para
calcularmos o menor autovalor de A.
4
Autovalores de Matrizes Simétricas
Lembrete: Uma matriz A é simétrica se AT = A
4.1
Decomposição Simétrica de Schur
Os métodos que iremos estudar buscam uma aproximação para a decomposição
simétrica de Schur, que se baseia no seguinte teorema:
Teorema: Seja A ∈ R2 uma matriz simétrica e λ1 , λ2 , ..., λn autovalores de
A. Então, existe uma matriz ortogonal Q ∈ Rn,n tal que QT AQ = Λ, sendo
Λ = diag(λ1 , λ2 , · · · , λn ).
4.2
Método de Jacobi
A ideia central do método de Jacobi é aplicar transformações ortogonais na
matriz de modo a torná-la mais próxima de uma matriz diagonal. Daı́, queremos
encontrar uma matriz ortogonal Q de modo que QT AQ diminua o valor dado
por
v
u∑
n
u n ∑
ndiag(A) = t
a2ij
i=1 j=1,i̸=j
A matriz de rotação de Jacobi J(p, q, θ) é igual á matriz de Givens. A função
a seguir recebe uma matriz simétrica A e um valor real tolerancia e retorna os
autovalores e autovetores desta matriz. o Valor dado pela tolerancia indica o
quão próximo queremos estar da decomposição simétrica 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
12
//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)
J(p,q)
J(q,p)
J(q,q)
=
=
=
=
c
s
-s
c
Q = Q*J
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 funções auxiliares:
//funcao para calculo do valor ndiag
function valor = Ndiag(A)
aux = 0
for i=1:linhas
for j=1:linhas
13
if i <> j then
aux = aux + (A(i,j)*A(i,j))
end
end
end
valor = sqrt(aux)
endfunction
//funcao para buscar indices p e q do maior valor absoluto fora
//da diagonal de uma matriz
function [pe,que] = buscaIndicesPeQ(A)
[l,c] = size(A)
maior = abs(A(2,1))
pe = 2
que = 1
for u =1:l
for v =1:l
if u <> v then
if abs(A(u,v)) > maior then
maior = abs(A(u,v))
pe = v
que = u
end
end
end
end
endfunction
Vamos, por exemplo, encontrar os autovalores e autovetores da matriz A
definida abaixo, definindo uma tolerência de 10−1 .


7 1 10
A =  1 7 10 
10 10 −2
Utilizamos o código
A=[7, 1, 10;
1, 7, 10;
10, 10, -2]
tolerancia = 0.02
14
[autovalores, autovetores] = JacobiClassico(A, tolerancia)
Obtemos como resposta
autovetores
=
- 0.7116322
0.4081016
0.7025375
0.4080968
- 0.0045490 - 0.8166456
autovalores =
- 0.5718677
- 0.5830078
- 0.5771215
6.0007446
- 11.999998
17.999253
Os resultados nos indicam que temos os autovalores λ1 = 6, λ2 = −11.99
e λ3 = 17.99, com os respectivos autovalores v1 = (−0.7 0.7 0)T , v2 =
(0.4 0.4 − 0.8)T e v3 = (−0.5 − 0.5 − 0.5)T .
4.3
Método de Jacobi Cı́clico
O método de Jacobi clássico possui o problema de que a busca pelo elemento de
maior valor absoluto fora da diagonal (funcão auxiliar buscaIndicesP eQ(A))
requer O(n2 ) operações. Por isso, o algoritmo Cı́clico de Jacobi elimina esta
busca, diminuindo assim o total de operações. Assim, temos a função 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
15
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)
J(p,q)
J(q,p)
J(q,q)
=
=
=
=
c
s
-s
c
Q = Q*J
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 função agora somente utiliza a função auxiliar N diag(A). Não
se esqueça de incluı́-la em seu código.
4.4
Método de Rutishauser - LR
O método de Rutishauser é baseado no seguinte teorema:
Teorema: Seja uma matriz A ∈ Rn,n simétrica com autovalores λ1 , · · · , λn
tais que |λ1 | > |λ2 | > · · · > |λn | e a sequência Ak definida por A1 = A e
Ak+1 = Rk Lk , para k = 1, 2, ..., sendo que Ak = Lk Rk é a decomposição de Ak
no produto de matrizes triangular inferior Lk e triangular superior Rk . Então,
lim Ak = S,
k→∞
sendo S uma matriz triangular superior cuja diagonal é formada por autovalores de A.
Note que o teorema acima faz uso da decomposição LU de uma matriz. A
função abaixo recebe uma matriz A simétrica e uma precisao real e retorna os
16
seus autovalores. Aqui, o valor dado pela precisao nos indica quão próximo de
uma matriz triangular superior queremos estar.
function autovalores = Rutishauser(A, precisao)
flag = 1
A1 = A
while flag == 1
[L,U] = decomposicaoLU(A)
A = U*L
//verificamos se elementos abaixo da diagonal sao menores
//que a precisao
flag = verifica(A,precisao)
end
autovalores = diag(A)
endfunction
Esta função faz uso de duas outras funções auxiliares:
//funcao para decomposicao LU
function [L,U] = decomposicaoLU(A)
[linhas, colunas] = size(A)
cont = 1
L = eye(linhas,colunas)
U = zeros(linhas,colunas)
for i=2:linhas
cont2 = i;
while cont2<=linhas
aux = A(cont2,cont)/A(i-1,i-1)
L(cont2,cont) = aux
for j=1:colunas
A(cont2,j) = A(cont2,j)-(aux*A(i-1,j))
end
cont2 = cont2+1
17
end
cont = cont+1
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 método de Rutishauser para encontrar os autovalores da matriz A indicada abaixo, com precisão de 10−2 .


2 0 1
A= 0 1 0 
1 0 1
Fazemos isto com o código
A = [2, 0, 1;
0, 1, 0;
1, 0, 1]
precisao = 0.02
autovalores = Rutishauser(A, precisao)
18
Obtemos o seguinte resultado
autovalores =
2.6153846
1.
0.3846154
O resultado nos mostra que os autovalores da matriz A são λ1 = 2.6, λ2 = 1
e λ3 = 0.3.
4.5
Método de Francis - QR
O método de Francis trabalha de forma parecida ao método de Rutishauser, com
a diferença de que é utilizada a decomposicção QR, ao invés da decomposição
LU.
Teorema: Seja uma matriz A ∈ Rn,n simétrica com autovalores λ1 , · · · , λn
tais que |λ1 | > |λ2 | > · · · > |λn | e a sequência Ak definida por A1 = A e
Ak+1 = Rk Qk , para k = 1, 2, ..., sendo que Ak = Qk Rk é a decomposição de Ak
no produto de matrizes ortogonal Qk e triangular superior Rk . Então,
lim Ak = S,
k→∞
sendo S uma matriz triangular superior cuja diagonal é formada por autovalores de A.
A função abaixo calcula os autovalores de uma matriz A pelo método 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 função utiliza a função auxiliar verif ica() descrita anteriormente e a
função qr() do próprio Scilab.
19