IT600 – An´ alise matricial de sistemas de energia el´ etrica 1 Semestre 2011 Cap´ıtulo 3 – Elimina¸c˜ ao de Gauss e decomposi¸c˜ ao LDU
3.1
Introdu¸c˜ ao
◮ Os princ´ıpios b´ asicos das chamadas t´ecnicas de esparsidade s˜ ao:
3.2
minimizar a quantidade de dados armazenados
´ ˜ MINIMIZAR O NUMERO DE OPERAC ¸ OES REALIZADAS
PRESERVAR A ESPARSIDADE
Resolu¸c˜ ao de sistemas de equa¸c˜ oes alg´ ebricas lineares
3.2.1
M´ etodo direto – invers˜ ao da matriz A
◮ Considerar o seguinte sistema de equa¸co ˜es alg´ebricas lineares colocado na forma matricial: A·x =b em que:
10 1 2 2 9 2 1 1 8 1 A= 1 9 2 2 2 10 1 1 1 9
–1–
e
b=
1 2 3 2 1 2
Como A tem inversa, a solu¸c˜ ao ´e ´ unica e dada por: x = A−1 · b 1 0,1052 −0,0135 0,0035 −0,0248 0,0064 −0,0011 −0,0250 0,1205 −0,0307 0,0087 −0,0143 0,0050 2 0,0033 −0,0156 0,1308 −0,0015 0,0033 −0,0149 3 = −0,0134 0,0073 −0,0022 0,1196 −0,0250 0,0030 · 2 0,0078 −0,0260 0,0081 −0,0260 0,1091 −0,0130 1 2 −0,0012 0,0046 −0,0154 0,0030 −0,0125 0,1142 0,0433 0,1369 0,3349 = 0,2149 0,0113 0,1838 ◮ Obter a inversa da matriz A e multiplic´ a-la pelo vetor b s˜ ao tarefas simples de serem realizadas, mesmo que a matriz A esteja armazenada de forma compacta.
Problema: a matriz A−1 ´e cheia. Do ponto de vista de espa¸co de mem´ oria, n˜ ao ´e vantajoso armazen´ a-la utilizando os esquemas de armazenamento compacto.
–2–
3.2.2
M´ etodo direto – substitui¸c˜ ao de vari´ aveis
◮ Considerar o seguinte sistema de equa¸co ˜es alg´ebricas lineares:
3x1 − 2x2 − x3 = 0 −2x1 + 2x2 − x3 = −1 −x1 − x2 + 3x3 = 1 ◮ A solu¸c˜ ao pode ser obtida, como mostrado na se¸c˜ ao 3.2.1, atrav´es de: 0 3 −2 −1 x1 −2 2 −1 · x2 = −1 x3 1 −1 −1 3
−1 x1 3 −2 −1 0 x2 = −2 2 −1 · −1 x3 −1 −1 3 1 −5/3 −7/3 −4/3 0 = −7/3 −8/3 −5/3 · −1 −4/3 −5/3 −2/3 1 1 = 1 1
–3–
◮ Outra maneira de resolver o sistema de equa¸co ˜es ´e atrav´es da elimina¸c˜ ao de vari´ aveis nas equa¸co ˜es. Reescrevendo a primeira equa¸c˜ ao: 2 1 x1 = x2 + x3 3 3 e substituindo na segunda e terceira equa¸co ˜es: 3x1 − 2x2 − x3 = 0 5 2 3 x2 − 3 x3 = −1 − 35 x2 + 83 x3 = 1 ◮ As duas ´ ultimas equa¸co ˜es n˜ ao dependem mais de x1 . Reescrevendo a segunda equa¸c˜ ao: 5 3 x2 = x3 − 2 2 e substituindo na terceira equa¸c˜ ao: 3x1 − 2x2 − x3 = 0 5 2 3 x2 − 3 x3 = −1 − 23 x3 = − 23 ◮A´ ultima equa¸c˜ ao agora s´ o depende de x3 . A sequˆencia de c´ alculos envolvendo elimina¸c˜ ao de vari´ aveis e que resulta em uma ´ ultima equa¸c˜ ao contendo apenas uma vari´ avel (cuja resolu¸c˜ ao ´e trivial) ´e chamada ao Forward de Substitui¸c˜
–4–
◮ Colocando o sistema de equa¸co ˜es na forma matricial tem-se:
3 −2 −1 0 x1 0 2/3 −5/3 · x2 = −1 −3/2 x3 0 0 −3/2 angulo inferior da matriz de coeficientes foram tornados ou seja, os elementos do triˆ iguais zero. ◮ Resolvendo a terceira equa¸c˜ ao:
x3 = 1
◮ Substituindo o valor de x3 nas duas primeiras equa¸co ˜es: 3x1 − 2x2 = 1 2 x2 = 23 3 x3 = 1 ◮ Agora a segunda equa¸c˜ ao s´ o depende de x2. Resolvendo a segunda equa¸c˜ ao:
x2 = 1
–5–
◮ Substituindo o valor de x2 na primeira equa¸c˜ ao: 3x1 = 3 x = 1 2 x3 = 1 ◮ Agora a primeira equa¸c˜ ao s´ o depende de x1. Resolvendo a primeira equa¸c˜ ao:
x1 = 1
◮ A sequˆencia de c´ alculos envolvendo a substitui¸c˜ ao das vari´ aveis j´ a calculadas nas demais equa¸co ˜es e que resulta na obten¸c˜ ao da solu¸c˜ ao (uma vari´ avel por vez) ´e ao Backward (ou Back) chamada de Substitui¸c˜ ◮ A solu¸c˜ ao do problema ´e: x1 = 1 x2 = 1 x3 = 1
ou, na forma matricial:
1 0 0 x1 1 0 1 0 · x2 = 1 x3 0 0 1 1 angulo superior da matriz de coeficientes foram tornados ou seja, os elementos do triˆ iguais zero.
–6–
◮ Nota-se que, na solu¸c˜ ao, a matriz de coeficientes A transformou-se na matriz identidade. Al´em disso, o vetor independente transformou-se no vetor solu¸c˜ ao. 3.2.3
M´ etodo direto – elimina¸c˜ ao de Gauss
ao de vari´ avies pode tamb´em ser encarado ◮ O procedimento baseado em substitui¸c˜ sob outro ponto de vista. Repetindo o sistema de equa¸co ˜es original:
3x1 − 2x2 − x3 = 0 −2x1 + 2x2 − x3 = −1 −x1 − x2 + 3x3 = 1
◮ A elimina¸c˜ ao de x1 da segunda e terceira equa¸co ˜es pode ser conseguida atrav´es da defini¸c˜ ao de combina¸co ˜es lineares apropriadas entre as equa¸co ˜es. Neste caso, as combina¸co ˜es s˜ ao as seguintes: 2 Linha 2 ← · Linha 1 + Linha 2 3 1 Linha 3 ← · Linha 1 + Linha 3 3 ◮ As opera¸co ˜es entre as equa¸co ˜es definidas acima resultam em: 3x1 − 2x2 − x3 = 0 2 x2 − 53 x3 = −1 3 − 35 x2 + 83 x3 = 1
que j´ a foi obtido anteriormente.
◮ A elimina¸c˜ ao de x2 da terceira equa¸c˜ ao pode ser conseguida atrav´es da seguinte combina¸c˜ ao linear: 5 · Linha 2 + Linha 3 Linha 3 ← 2 –7–
◮ A opera¸c˜ ao definida acima resulta em: 3x1 − 2x2 − x3 = 0 2 x2 − 53 x3 = −1 3 − 23 x3 = − 23
que tamb´em j´ a foi obtida anteriormente. Novamente, a terceira equa¸c˜ ao s´ o depende de x3 e sua solu¸c˜ ao ´e trivial. At´e aqui foi realizada a substitui¸c˜ ao Forward, e os elementos do triˆ angulo inferior da matriz de coeficientes foram tornados iguais a zero. ◮ Obt´em-se o valor de x3 atrav´es de:
Linha 3 ←
2 − 3
· Linha 3
Esta opera¸c˜ ao corresponde a tornar o coeficiente de x3 igual a 1. Pode-se tamb´em realizar o mesmo tipo de opera¸c˜ ao para a primeira e segunda equa¸co ˜es, tornando iguais a 1 os coeficientes de x1 na primeira equa¸c˜ ao e x2 na segunda equa¸c˜ ao: 1 Linha 1 ← · Linha 1 3 3 Linha 2 ← · Linha 2 2 O sistema de equa¸co ˜es resultante ´e: 1 2 x1 − 3 x2 − 3 x3 = 0 x2 − 52 x3 = − 23 x3 = 1
Em geral, considera-se que as opera¸co ˜es descritas acima (tornar alguns coeficientes iguais a 1) tamb´em fazem parte da substitui¸c˜ ao Forward. –8–
◮ A substitui¸c˜ ao Back ´e realizada fazendo-se as seguintes opera¸co ˜es entre as equa¸co ˜es: 5 Linha 2 ← · Linha 3 + Linha 2 2 1 Linha 1 ← · Linha 3 + Linha 1 3 1 2 x1 − 3 x2 = 3 x2 = 1 x3 = 1 2 Linha 1 ← · Linha 2 + Linha 1 3 x1 = 1 x2 = 1 x3 = 1 ◮ Pode-se notar que, ap´ os a substitui¸c˜ ao Back, os elementos do triˆ angulo superior foram tornados iguais a zero, e a matriz dos coeficientes A transformou-se na matriz identidade. Al´em disso, o vetor independente transformou-se no vetor solu¸c˜ ao.
Conclus˜ ao: os c´ alculos envolvendo a substitui¸c˜ ao de vari´ aveis e elimina¸c˜ ao de Gauss resultaram na solu¸c˜ ao do problema sem que tenha sido necess´ ario inverter a matriz de coeficientes (como feito na resolu¸c˜ ao direta).
–9–
◮ Pode-se resolver um problema do tipo:
A·x =b atrav´es da elimina¸c˜ ao de Gauss partindo-se de uma matriz aumentada contendo A e b e obtendo-se uma nova matriz aumentada que cont´em a solu¸c˜ ao x : [A | b] ⇓ [I | x ] em que I ´e a matriz identidade. ◮ Para a obten¸c˜ ao de I a partir de A deve-se realizar opera¸c˜ oes (combina¸co ˜es lineares) entre as linhas de A, de forma a:
o 1 tornar todos os elementos do triˆ angulo inferior de A iguais a 0; o 2 tornar todos os elementos da diagonal de A iguais a 1; o 3 tornar todos os elementos do triˆ angulo superior de A iguais a 0.
As opera¸co ˜es acima ser˜ ao aplicadas tamb´em sobre o vetor b, que resultar´ a na solu¸c˜ ao x . ◮ Os os 1 e 2 correspondem ` a substitui¸c˜ ao Forward do processo de elimina¸c˜ ao de vari´ aveis descrito anteriormente. O o 3 corresponde ` a substitui¸c˜ ao Back.
– 10 –
3.2.4
Pivoteamento
◮ Considere o seguinte sistema de equa¸co ˜es lineares:
10−10 1 2 1
x1 x2
=
1 5
cuja solu¸c˜ ao, por inspe¸c˜ ao, ´e x ≈ [2 1]. ◮ Resolvendo o problema por elimina¸c˜ ao de Gauss tem-se: ou seja:
1 0 0 1
10−10 0
1 1010 0 1
1 −2 · 1010 + 1
1 −2 · 1010 + 5
10 10 −2 · 1010 + 5 / −2 · 1010 + 1
10 −2 · 1010 + 5 / −2 · 1010 + 1 · −1010 + 10 −2 · 1010 + 5 / −2 · 1010 + 1
x≈
0 1
que n˜ ao corresponde ` a solu¸c˜ ao correta. Problema: o termo 10−10 ´e praticamente igual a zero para a maioria dos computadores.
– 11 –
◮ Considere agora que o sistema de equa¸co ˜es seja rearranjado como:
2 1 −10 10 1
x1 x2
=
5 1
◮ Resolvendo o problema por elimina¸c˜ ao de Gauss tem-se: ou seja:
2 0
1 1/2 0 1
1 −1/2 · 10−10 + 1
5 −5/2 · 10−10 + 1
5/2 −5/2 · 10−10 + 1 / −1/2 · 10−10 + 1
−10 1 0 5/2 − 1/2 · −5/2 · 10−10 + 1 / −1/2 · 10 + 1 0 1 −5/2 · 10−10 + 1 / −1/2 · 10−10 + 1
x≈
2 1
◮ O rearranjo acima ´e chamado pivoteamento, que consiste em manter os maiores elementos da matriz na diagonal durante a substitui¸c˜ ao forward. Note que foi feita uma troca entre as linhas do sistema de equa¸co ˜es. Este tipo de pivoteamento ´e chamado parcial. O pivoteamento completo consiste em realizar trocas de linhas e colunas.
– 12 –
3.2.5
M´ etodos iterativos
◮ Considere um sistema de n equa¸co ˜es alg´ebricas lineares A · x = b Tomando a linha i da equa¸c˜ ao matricial acima: i
xi
Ai i
i
.
n X
Ai i xi +
j=1 n X
=
x
A
bi
b
A i j xj = bi
i = 1, . . . , n
A i j xj = bi
i = 1, . . . , n
j=1,j6=i
Resolvendo para xi tem-se: 1 xi = · Ai i
bi −
n X
Ai j xj
j=1,j6=i
– 13 –
!
i = 1, . . . , n
◮ Nota-se que:
Cada inc´ ognita xi pode ser escrita em fun¸c˜ ao das demais.
Isto pode sugerir um processo iterativo do tipo: • Arbitre valores iniciais para xi0, i = 1, . . . , n. ao e calcule os • Substitua os valores de xi0 no termo do lado direito da equa¸c˜ 1 ao). novos valores de xi (lado esquerdo da equa¸c˜ ao e calcule os • Substitua os valores de xi1 no termo do lado direito da equa¸c˜ 2 ao). novos valores de xi (lado esquerdo da equa¸c˜ • ··· que p´ ara quando os valores de xi convergirem para a solu¸c˜ ao
◮ O processo iterativo implica na existˆencia de uma seq¨ uˆencia de valores para as inc´ ognitas, come¸cando por valores iniciais arbitrados: xi0 → xi1 → xi2 · · · enquanto que no m´etodo de Elimina¸c˜ ao de Gauss mostrado anteriormente o vetor solu¸c˜ ao ´e obtido em um s´ o o.
– 14 –
◮ Para uma itera¸c˜ ao (m + 1), o processo iterativo pode ser definido como:
(m+1)
xi
=
1 · Ai i
bi −
n X
(m)
Ai j xj
j=1,j6=i
!
i = 1, . . . , n
M´etodo de Gauss-Jacobi
Para a obten¸c˜ ao de xi(m+1) s˜ ao utilizados os valores de xi(m) (todos os valores da itera¸c˜ ao anterior) Uma forma alternativa para o processo iterativo ´e:
xi(m+1) =
1 · Ai i
bi −
i −1 X
Ai j xj(m+1) −
j=1
n X
Ai j xj(m)
j=i +1
!
i = 1, . . . , n
M´etodo de Gauss-Seidel
(m+1)
Para a obten¸c˜ ao de xi elementos do vetor x
s˜ ao utilizados os valores mais recentes dispon´ıveis dos
– 15 –
Exemplo Processo iterativo utilizando os m´etodos de Gauss-Jacobi e Gauss-Seidel para n = 3: Gauss-Jacobi (m+1) x1 (m+1)
x2
x3(m+1)
1 = A11 1 = A22 1 = A33
h
(m) A12x2
(m) A13x3
i
+ · b1 − i h (m) (m) · b2 − A21x1 + A23x3 i h (m) (m) · b3 − A31x1 + A32x2
Gauss-Seidel x1(m+1) (m+1)
x2
(m+1)
x3
1 = A11 1 = A22 1 = A33
h
A12x2(m)
A13x3(m)
i
+ · b1 − i h (m) (m+1) + A23x3 · b2 − A21x1 i h (m+1) (m+1) + A32x2 · b3 − A31x1
– 16 –
◮ A ideia geral desses m´etodos iterativos ´e converter o sistema de equa¸co ˜es:
A·x =b em:
x = C · x + g = ϕ (x ) em que C ´e uma matriz (n × n), g ´e um vetor (n × 1) e ϕ (x) ´e uma fun¸c˜ ao de itera¸c˜ ao matricial. ◮ A formula¸c˜ ao acima sugere um processo iterativo. No caso do m´etodo de Gauss-Jacobi, tem-se a seguinte express˜ ao geral para o processo iterativo:
x
(k+1)
=C·x
(k)
(k) +g =ϕ x
◮ A matriz A pode ser colocada na forma:
A=L+D+U em que L ´e uma matriz que cont´em os elementos abaixo da diagonal (Lower triangle) de A, D ´e uma matriz que cont´em os elementos da diagonal de A, e U ´e uma matriz que cont´em os elementos acima da diagonal (Upper triangle) de A. Ent˜ ao, a equa¸c˜ ao de Gauss-Jacobi pode ser escrita na forma matricial como: x (k+1) = D−1 b − D−1 (L + U) x (k)
– 17 –
◮ A convergˆencia pode ser acelerada atrav´es da utiliza¸c˜ ao de parˆ ametros de acelera¸c˜ ao. O m´etodo mais popular ´e o chamado m´etodo SOR (successive overrelaxation). Para o m´etodo Gauss-Seidel:
(m+1)
zi
(m+1) xi
=
1 · Ai i
=
(m+1) ωzi
bi −
i −1 X
(m+1)
Ai j xj
−
j=1
+ (1 −
n X
(m)
Ai j xj
j=i +1 (m) ω)xi
!
i = 1, . . . , n
◮ Para ω = 1, tem-se o m´etodo Gauss-Seidel sem acelera¸c˜ ao. ◮ A acelera¸c˜ ao corresponde a uma extrapola¸c˜ ao: x
(m+1)
=x
(m)
(m+1) (m) +ω· z −x
= x (m) + ω · ∆x
x (m+1) z (m+1) x (m) ∆x ω · ∆x
– 18 –
i = 1, . . . , n
◮ Para cada problema existe um ω ´ otimo, mas para a maioria dos problemas pr´ aticos, valores aproximados (emp´ıricos) s˜ ao escolhidos.
◮ Outro m´etodo iterativo comumente usado ´e o do gradiente conjugado. Este m´etodo consiste na minimiza¸c˜ ao da fun¸c˜ ao: E (x) = kA x − b k2
◮ Uma caracter´ıstica interessante do m´etodo ´e que garante-se que a convergˆencia ocorrer´ a em no m´ aximo n itera¸co ˜es (n ´e a dimens˜ ao de A), caso A seja definida positiva. ◮ Algoritmo:
Inicializa¸c˜ ao: k =1 r0 = A x0 − b ρ0 = −AT r 0
Enquanto kr k k ≥ ε: αk+1 = kAT r k k2/kAρk k2 x k+1 = x k + αk+1 ρk r k = A x k+1 − b Bk+1 = kAT r k+1 k2 /kAT r k k2 ρk+1 = −AT r k+1 + Bk+1 ρk k =k +1
– 19 –
3.3
Elimina¸c˜ ao de Gauss e decomposi¸c˜ ao LDU
◮ Considerar novamente a matriz A do exemplo anterior:
10 1 2 2 9 2 1 1 8 1 A= 1 9 2 2 2 10 1 1 1 9
A matriz aumentada ´e: 10 1 2 2 9 2 1 1 8 1 [A | b] = 1 9 2 2 2 10 1 1 1 9
1 2 3 2 1 2
◮ O processo de elimina¸c˜ ao de Gauss ´e iniciado com o o 1 (tornar todos os elementos do triˆ angulo inferior de A iguais a 0)
– 20 –
A primeira posi¸c˜ ao a ser zerada ´e a posi¸c˜ ao (2,1) (A2,1 = 2)
10 1 2 2 9 2 1 1 8 1 1 9 2 2 2 10 1 1 1 9
10,0000 0 0 1,0000 0 0
1 2 3 2 1 2
linha 2 ← linha 1 × (−2/10) + linha 2
1,0000 0 2,0000 0 0 8,8000 2,0000 -0,4000 1,0000 0 1,0000 8,0000 0 0 1,0000 0 0 9,0000 2,0000 0 2,0000 0 2,0000 10,0000 1,0000 0 1,0000 0 1,0000 9,0000
1,0000 1,8000 3,0000 2,0000 1,0000 2,0000
◮ Observa¸co ˜es:
A2,1 = 0
A2,4 6= 0 (essa posi¸c˜ ao era originalmente igual a zero) Estes novos elementos n˜ ao nulos que aparecem em posi¸co ˜es originalmente nulas devido ` as opera¸co ˜es realizadas s˜ ao chamados de fill-ins
todos os esquemas de armazenamento compacto de matrizes devem prever um espa¸co de mem´ oria adicional para levar em conta o poss´ıvel aparecimento de fill-ins: MAX = NB + 2 · NR + folga
– 21 –
a opera¸c˜ ao realizada equivale a premultiplicar a matriz aumentada por uma matriz L1 tal que a matriz resultante tenha A2,1 = 0:
1
α2,1 1 1 L1 = 1 1
1
α2,1 = −A2,1 /A1,1 = −2/10
[A1 | b1 ] = L1 · [A | b]
α2,1 ´e chamado de fator triangular
Observar que a coluna da matriz aumentada correspondente ao vetor b foi alterada na posi¸c˜ ao 2, resultando no vetor b1 :
b2 ← b1 · α2,1 + b2
– 22 –
◮ O pr´ oximo o ´e zerar a posi¸c˜ ao (4,1) (A4,1 = 1). Para isso, a seguinte opera¸c˜ ao ´e realizada:
linha 4 ← linha 1 × (−1/10) + linha 4 e o resultado ´e:
10,0000 1,0000 0 2,0000 0 0 0 8,8000 2,0000 −0,4000 1,0000 0 0 1,0000 8,0000 0 0 1,0000 0 -0,1000 0 8,8000 2,0000 0 0 2,0000 0 2,0000 10,0000 1,0000 0 0 1,0000 0 1,0000 9,0000
1,0000 1,8000 3,0000 1,9000 1,0000 2,0000
Nota-se que um novo fill-in apareceu na posi¸c˜ ao correspondente a A4,2 e que a posi¸c˜ ao correspondente a b4 foi alterada ◮ O processo continua at´e que todo o triˆ angulo inferior de A tenha sido zerado:
10,0000 1,0000 0 2,0000 0 0 1,0000 0 8,8000 2,0000 -0,4000 1,0000 0 1,8000 1,0000 2,7955 0 0 7,7727 0,0455 -0,1136 0 0 0 8,7953 2,0117 -0,0029 1,9123 0 0 0 0 9,2872 1,0592 0,2992 0 0 0 0 0 8,7555 1,6089
Os elementos marcados correspondem a posi¸co ˜es que eram originalmente iguais a zero.
– 23 –
◮ O o 2 corresponde a tornar todos os elementos da diagonal de A iguais a 1. Para tornar o elemento A1,1 igual a 1, deve-se realizar a seguinte opera¸c˜ ao:
linha 1 ← linha 1 × (1/10) ou, da mesma forma, premultiplicar a matriz aumentada por uma matriz D1 tal que a matriz resultante tenha A1,1 = 1:
D1 =
α1,1 1 1 1 1 1
α1,1 = 1/A1,1 = 1/10
◮ Realizando essas opera¸co ˜es para todas as linhas, obt´em-se finalmente:
1,0000 0,1000 0 0,2000 0 0 0,1000 0 1,0000 0,2273 −0,0455 0,1136 0 0,2045 0 0 1,0000 0,0058 −0,0146 0,1287 0,3596 0 0 0 1,0000 0,2287 −0,0003 0,2174 0 0 0 0 1,0000 0,1140 0,0322 0 0 0 0 0 1,0000 0,1838
– 24 –
◮ O o 3 consiste em zerar todos os elementos do triˆ angulo superior de A. Come¸cando pelo elemento A5,6 = 0,1140, deve-se realizar a seguinte opera¸c˜ ao:
linha 5 ← linha 6 × (−0,1140) + linha 5 ou, da mesma forma, premultiplicar a matriz aumentada por uma matriz U1 tal que a matriz resultante tenha A5,6 = 0:
U1 =
1 1 1 1 1 α5,6 1
α5,6 = −A5,6 = −0,1140
O resultado da opera¸c˜ ao ´e:
1,0000 0,1000 0 0,2000 0 0 0,1000 0 1,0000 0,2273 −0,0455 0,1136 0 0,2045 0 0 1,0000 0,0058 −0,0146 0,1287 0,3596 0 0 0 1,0000 0,2287 −0,0003 0,2174 0 0 0 0 1,0000 0 0,0113 0 0 0 0 0 1,0000 0,1838
◮ Repetindo o procedimento para todos os elementos n˜ ao nulos do triˆ angulo superior de A chega-se a:
1,0000 0 0 0 0 0 0 1,0000 0 0 0 0 0 0 1,0000 0 0 0 0 0 0 1,0000 0 0 0 0 0 0 1,0000 0 0 0 0 0 0 1,0000
– 25 –
0,0433 0,1369 0,3349 0,2149 0,0113 0,1838 | {z x
}
◮ Todos os valores de α ao longo do processo de elimina¸c˜ ao de Gauss s˜ ao chamados de fatores triangulares. O processo de obten¸c˜ ao dos valores de α tamb´em ´e chamado de fatora¸c˜ ao triangular ao resultou em 8 fill-ins (oito ◮ No caso da matriz A do exemplo, o processo de fatora¸c˜ elementos adicionais que devem ser armazenados) ◮ Id´eia: armazenar os fatores triangulares (α) nas posi¸co ˜es correspondentes aos elementos que se desejou zerar ou tornar iguais a 1:
A=
0,1000 −0,1000 0 −0,2000 0 0 −0,2000 0,1136 −0,2273 0,0455 −0,1136 0 0 −0,1136 0,1287 −0,0058 0,0146 −0,1287 −0,1000 0,0114 −0,0029 0,1137 −0,2287 0,0003 0 −0,2273 0,0585 −0,2380 0,1077 −0,1140 0 0 −0,1287 0,0007 −0,1094 0,1142
Conhecendo-se os fatores triangulares e a ordem em que eles devem ser utilizados, pode-se repetir as opera¸co ˜es para v´ arios vetores b
– 26 –
◮ Existem dois esquemas b´ asicos de fatora¸c˜ ao triangular para se zerar os elementos dos triˆ angulos inferior e superior da matriz:
Fatora¸c˜ ao por colunas:
Zerar triˆ angulo inferior superior
Colunas
Linhas
esquerda → direita cima → baixo direita → esquerda baixo → cima
– 27 –
Fatora¸c˜ ao por linhas:
Zerar triˆ angulo inferior superior
Linhas
Colunas
cima → baixo esquerda → direita baixo → cima direita → esquerda
A sequˆencia das opera¸co ˜es ´e importante pois uma opera¸c˜ ao n˜ ao deve criar elementos n˜ ao nulos em posi¸co ˜es que foram zeradas em opera¸co ˜es anteriores. Tipicamente a substitui¸c˜ ao forward (zerar triˆ angulo inferior) ´e feita por colunas e a substitui¸c˜ ao back (zerar triˆ angulo superior) ´e feita por linhas.
– 28 –
◮ Resumo das opera¸co ˜es (combina¸co ˜es lineares) realizadas para transformar A em I na matriz aumentada: (U11 · U10 · . . . · U2 · U1) · (D6 · . . . · D1 ) · (L11 · L10 · . . . · L2 · L1 ) · A = I As matrizes L1 → U11 s˜ ao chamadas de matrizes elementares. ◮ As mesmas opera¸co ˜es s˜ ao realizadas para transformar b em x: (U11 · U10 · . . . · U2 · U1 ) · (D6 · . . . · D1 ) · (L11 · L10 · . . . · L2 · L1 ) · b = x Logo: A = [(U11 · U10 · . . . · U2 · U1) · (D6 · . . . · D1 ) · (L11 · L10 · . . . · L2 · L1 )]−1 ou: A = L1−1 · L2−1 · . . . · L10−1 · L11−1 · D1 −1 · . . . · D6 −1 · U1−1 · U2−1 · . . . · U10−1 · U11−1 =L·D·U Este processo tamb´em ´e chamado de decomposi¸c˜ ao (ou fatora¸c˜ ao) LDU da matriz A:
L ´e uma matriz triangular inferior (Lower triangular) – o triˆ angulo inferior tem elementos n˜ ao nulos e os elementos da diagonal s˜ ao iguais a 1
D ´e uma matriz diagonal – somente a diagonal tem elementos n˜ ao nulos
U ´e uma matriz triangular superior (Upper triangular) – o triˆ angulo superior tem elementos n˜ ao nulos e os elementos da diagonal s˜ ao iguais a 1
– 29 –
◮ A matriz L ´e igual a: L = L1 −1 · L2 −1 · L3 −1 · L4−1 · L5−1 · L6−1 · L7−1 · L8 −1 · L9 −1 · L10−1 · L11−1
1,0000 0 0 0 0 0 0,2000 1,0000 0 0 0 0 0 0,1136 1,0000 0 0 0 = 0,1000 −0,0114 0,0029 1,0000 0 0 0 0,2273 −0,0585 0,2380 1,0000 0 0 0 0,1287 −0,0007 0,1094 1,0000
Os elementos do triˆ angulo inferior correspondem ao negativo dos fatores triangulares α obtidos anteriormente. ◮ Notar que:
Lk
−1
◮ A matriz D ´e igual a:
−1
1
1 0 = ... α
1
1
1 0 = ... −α
1
D = D1 −1 · D2 −1 · D3 −1 · D4 −1 · D5−1 · D6 −1 10,0000 0 0 0 0 0 0 8,8000 0 0 0 0 0 0 7,7727 0 0 0 = 0 0 0 8,7953 0 0 0 0 0 0 9,2872 0 0 0 0 0 0 8,7555
Os elementos da diagonal correspondem ao inverso dos fatores triangulares α obtidos anteriormente.
– 30 –
◮ Notar que:
Dk
−1
◮ A matriz U ´e igual a:
−1
α
=
1 0 0 ... 1
1/α
=
1 0 0 ... 1
U = U1 −1 · U2−1 · U3−1 · U4−1 · U5−1 · U6−1 · U7−1 · U8−1 · U9 −1 · U10−1 · U11−1
1,0000 0,1000 0 0,2000 0 0 0 1,0000 0,2273 −0,0455 0,1136 0 0 0 1,0000 0,0058 −0,0146 0,1287 = 0 0 0 1,0000 0,2287 −0,0003 0 0 0 0 1,0000 0,1140 0 0 0 0 0 1,0000
Os elementos do triˆ angulo superior correspondem ao negativo dos fatores triangulares α obtidos anteriormente. ◮ Notar que:
Uk
−1
=
1
−1
1 α . . . 0 1
– 31 –
=
1
1 −α . . . 0 1
◮ Normalmente a matriz D ´e incorporada a L (L ← L · D): A=L·U No caso da matriz A do exemplo: L=L·D 10,0000 0 0 0 0 0 0,2000 8,8000 0 0 0 0 0 0,1136 7,7727 0 0 0 = 0,1000 −0,0114 0,0029 8,7953 0 0 0 0,2273 −0,0585 0,2380 9,2872 0 0 0 0,1287 −0,0007 0,1094 8,7555
◮ A solu¸c˜ ao ´e dada por: x = A−1 · b = (L · U)−1 · b = U−1 · L−1 · b ou w = L−1 · b x = U−1 · w
Substitui¸c˜ ao Forward Substitui¸c˜ ao Back
◮ A fatora¸c˜ ao LDU n˜ ao preserva a esparsidade original da matriz, j´ a que fill-ins s˜ ao criados ao longo do processo. Existem v´ arias t´ecnicas propostas para minimizar o n´ umero de fill-ins gerados (algumas ser˜ ao vistas adiante). Elas se baseiam na reordena¸c˜ ao das linhas e colunas da matriz.
– 32 –
◮ Procedimento geral:
armazenar a matriz A utilizando algum esquema de armazenamento compacto
reordenar linhas e colunas de A de forma a minimizar o n´ umero de fill-ins gerados
obter os fatores triangulares
realizar substitui¸c˜ ao Forward
realizar substitui¸c˜ ao Back
→
– 33 –
˜ SOLUC ¸ AO
3.4
Bifatora¸c˜ ao
◮ Procedimento desenvolvido por Zollenkopf ´ particularmente interessante no caso de matrizes estruturalmente sim´etricas1 ◮E ◮ Considerando um sistema de n equa¸co ˜es alg´ebricas lineares dado por A · x = b, a bifatora¸c˜ ao pode ser representada por: L(n) L(n−1) · · · L(2) L(1) AR(1)R(2) · · · R(n−1)R(n) = I em que: L s˜ ao as matrizes de fatores ` a esquerda, R s˜ ao as matrizes de fatores ` a direita e I ´e a matriz identidade ◮ Logo: A−1 = R(1)R(2) · · · R(n−1)R(n) L(n) L(n−1) · · · L(2)L(1)
◮ O processo de obten¸c˜ ao das matrizes de fatores ` a esquerda e ` a direita ´e dado por: A(0) = A A(1) = L(1)A(0)R(1) A(2) = L(2)A(1)R(2) ··· ··· A(n) = L(n) A(n−1)R(n) = I
1 Denominamos
aqui uma matriz A de estruturalmente sim´ etrica caso Ai j e Aj i sejam simultaneamente iguais a zero ou diferentes de zero.
– 34 –
Exemplo Considere a seguinte matriz:
10 2 3 A= 3 8 1 1 3 15 Inicialmente, ser˜ ao obtidos os fatores triangulares atrav´es da decomposi¸c˜ ao LDU:
1 0 0 L1 = −3/10 1 0 −1/10 0 1 10 2 3 1/10 A = L1 · A = 0 37/5 0 14/5 147/10 1 0 0 1 0 L2 = 0 0 −14/37 1 10 2 3 1/10 A = L2 · A = 0 37/5 0 0 1085/74 1/10 0 0 0 5/37 0 D= 0 0 74/1085 1 1/5 3/10 1 1/74 A=D·A= 0 0 0 1
– 35 –
1 0 −3/10 U3 = 0 1 −1/74 0 0 1 1 1/5 0 1 0 A = U3 · A = 0 0 0 1 1 −1/5 0 1 0 U2 = 0 0 0 1 1 0 0 A = U3 · A = 0 1 0 0 0 1 Agrupando os fatores triangulares: 1/10 −1/5 −3/10 5/37 −1/74 FT = −3/10 −1/10 −14/37 74/1085
O processo de bifatora¸c˜ ao ´e: 1/10 0 0 L1 = −3/10 1 0 −1/10 0 1
1 −1/5 −3/10 1 0 R1 = 0 0 0 1 1 0 0 1/10 A = L1 · A · R1 = 0 37/5 0 14/5 147/10
– 36 –
1 0 0 5/37 0 L2 = 0 0 −14/37 1
1 0 0 R2 = 0 1 −1/74 0 0 1 1 0 0 0 A = L2 · A · R2 = 0 1 0 0 1085/74 1 0 0 1 0 0 0 L3 = 0 1 R3 = 0 1 0 0 0 74/1085 0 0 1 1 0 0 A = L3 · A · R3 = 0 1 0 0 0 1
Agrupando os fatores triangulares: 1/10 −1/5 −3/10 5/37 −1/74 FT = −3/10 −1/10 −14/37 74/1085
Portanto, os fatores triangulares s˜ ao os mesmos, por´em, obtidos em uma sequˆencia diferente.
– 37 –
◮ Na bifatora¸c˜ ao, a obten¸c˜ ao dos fatores triangulares ´e feita por colunas para o triˆ angulo inferior e por linhas para o triˆ angulo superior.
◮ Toda a metodologia proposta por Zollenkopf tem por base a bifatora¸c˜ ao, portanto, a estrutura de armazenamento compacto leva em conta esse processo. Al´em disso, a substitui¸c˜ ao forward ´e realizada por colunas, enquanto que a substitui¸c˜ ao back ´e realizada por linhas.
– 38 –
3.5
Elimina¸c˜ ao de Gauss e redu¸c˜ ao de circuitos
◮ Considerar o seguinte circuito de 6 barras e 7 ramos: 3 2
1
6
4 5
◮ Considerando que todos os ramos do circuito tenham itˆ ancias s´erie iguais a 1 pu, a matriz itˆ ancia nodal fica:
Y=
2 0 −1 −1 0 0 0 2 −1 0 0 −1 −1 −1 3 0 −1 0 −1 0 0 2 −1 0 0 0 −1 −1 3 −1 0 −1 0 0 −1 ∞
A matriz Y tem dimens˜ ao [6 × 6] (rede tem 6 n´ os). Como o n´ o 6 foi definido como sendo o n´ o de referˆencia, atribuiu-se ao elemento da diagonal Y66 um valor grande.
– 39 –
◮ Considerar o seguinte vetor das inje¸co ˜es de corrente nodais:
I=
1 1 1 1 1 X
em que a posi¸c˜ ao 6 se refere ` a inje¸c˜ ao de corrente no n´ o de referˆencia. ◮ O processo de elimina¸c˜ ao de Gauss come¸ca tornando nulos os elementos do triˆ angulo inferior de Y. Para zerar os elementos n˜ ao nulos da primeira coluna, deve-se premultiplicar a matriz Y por:
Logo:
L1 =
1 0 1/2 1/2 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
Y 1 = L1 · Y 2 0 −1 −1 0 0 0 2 −1 0 0 −1 0 −1 5/2 - 1/2 −1 0 = 0 0 - 1/2 3/2 −1 0 0 0 −1 −1 3 −1 0 −1 0 0 −1 ∞
Nota-se o aparecimento de dois fill-ins nas posi¸co ˜es Y34 e Y43.
– 40 –
◮ Como a primeira coluna (correspondente ao n´ o 1) foi zerada, considerar agora somente as colunas 2 → 6 da matriz Y1:
Y1′ =
2 0 −1 −1 0 0 0 2 −1 0 0 −1 0 −1 5/2 −1/2 −1 0 0 0 −1/2 3/2 −1 0 0 0 −1 −1 3 −1 0 −1 0 0 −1 ∞
◮ A constru¸c˜ ao da rede a partir de Y1′ resulta em: n´ o eliminado 3 2
1 ramo fict´ıcio criado 4
6 5 ◮ Considera¸co ˜es:
Zerar a primeira coluna de Y corresponde a eliminar o n´ o 1 do circuito
A elimina¸c˜ ao da coluna 1 resulta no aparecimento de 2 fill-ins
Os fill-ins correspondem ao aparecimento de um ramo fict´ıcio. Como o n´ o 1 era conectado aos n´ os 3 e 4, o ramo fict´ıcio conecta os n´ os 3 e 4 → fill-ins nas posi¸co ˜es (3,4) e (4,3) de Y
A itˆ ancia do ramo fict´ıcio do exemplo ´e igual a 1/2 pu, que corresponde ` a associa¸c˜ ao em s´erie dos ramos 1-3 e 1-4: y34 =
1 y13 · y14 = pu y13 + y14 2
– 41 –
◮ Aplicando a matriz L1 ao vetor de correntes tem-se: I1 = L1 · I 1 1 3/2 = 3/2 1 X
ou seja, a inje¸c˜ ao do n´ o 1 ´e distribu´ıda entre os n´ os 3 e 4, que s˜ ao diretamente conectados a ele. A forma como se d´ a essa distribui¸c˜ ao ´e fun¸c˜ ao dos valores das itˆ ancias conectadas ao n´ o 1. Pode-se tamb´em eliminar a linha de I1 correspondente ao n´ o 1, resultando em:
I1′ =
1 1 3/2 3/2 1 X
◮ Para a elimina¸c˜ ao do n´ o 2, a seguinte matriz elementar ´e definida:
L2 =
1 0 0 0 0 0 1 0 0 0 0 1/2 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1/2 0 0 0
– 42 –
0 0 0 0 0 1
A nova matriz itˆ ancia nodal ´e: Y2 = L2 · Y1′ 2 0 −1 −1 0 0 0 2 −1 0 0 −1 0 0 2 −1/2 −1 -1/2 = 0 0 0 −1/2 3/2 −1 0 0 −1 −1 3 −1 0 0 -1/2 0 −1 ∞ Nota-se o aparecimento de dois fill-ins nas posi¸co ˜es Y26 e Y62. ◮ Eliminando a linha e a coluna 2 de Y2 (que corresponde ao n´ o 2):
′ Y2 =
2 0 0 0 0 0
0 −1 −1 0 0 2 −1 0 0 −1 0 2 −1/2 −1 −1/2 0 −1/2 3/2 −1 0 0 −1 −1 3 −1 0 −1/2 0 −1 ∞
◮ Reconstruindo o circuito a partir de Y2′ : n´ o eliminado
n´ o eliminado 3
2
1
ramo fict´ıcio criado
ramo fict´ıcio criado 4
6 5
– 43 –
◮ Considera¸co ˜es:
Zerar a segunda coluna de Y corresponde a eliminar o n´ o 2 do circuito
A elimina¸c˜ ao da coluna 2 resulta no aparecimento de 2 fill-ins
Novamente aparece um ramo fict´ıcio. Como o n´ o 2 era conectado aos n´ os 3 e 6, o ramo fict´ıcio conecta os n´ os 3 e 6
A itˆ ancia do ramo fict´ıcio do exemplo ´e igual a 1/2 pu, que corresponde ` a associa¸c˜ ao em s´erie dos ramos 2-3 e 2-6: 1 y23 · y26 = pu y23 + y26 2
y36 =
Se a linha e a coluna correspondentes ao n´ o de referˆencia n˜ ao estivessem representadas na matriz Y, o aparecimento de certos fill-ins s´ o ser´ a notado atrav´es da mudan¸ca dos valores de elementos da diagonal.
◮ Aplicando a matriz L2 ao vetor de correntes tem-se: I2 = L2 · I1′
1 1 2 = 3/2 1 (X + 1/2)
ou seja, a inje¸c˜ ao do n´ o 2 ´e distribu´ıda entre os n´ os 3 e 6, que s˜ ao diretamente conectados a ele. Pode-se tamb´em eliminar a linha de I2 correspondente ao n´ o 2, resultando em:
I2′ =
1 1 2 3/2 1 (X + 1/2)
– 44 –
◮ Ap´ os a elimina¸c˜ ao do n´ o 3:
Y3′ =
2 0 0 0 0 0
0 −1 −1 0 0 2 −1 0 0 −1 0 2 −1/2 −1 −1/2 0 0 11/8 −5/4 −1/8 0 0 −5/4 5/2 −5/4 0 0 −1/8 −5/4 ∞
I3′ =
1 1 2 2 2 (X + 1)
n´ o eliminado n´ o eliminado
n´ o eliminado 3
2
1
6
4 5
ramo fict´ıcio criado O ramo fict´ıcio tem uma itˆ ancia igual a 1/8 pu. Os elementos da matriz que correspondem aos ramos 4-5 e 5-6 (que j´ a existiam) tiveram seus valores alterados. Ap´ os a elimina¸c˜ ao do n´ o 3, todos os seus ramos vizinhos (4, 5 e 6) aram a estar conectados entre si.
– 45 –
◮ Ap´ os a elimina¸c˜ ao do n´ o 4:
Y4′ =
2 0 0 0 0 0
0 −1 −1 0 0 2 −1 0 0 −1 0 2 −1/2 −1 −1/2 0 0 11/8 −5/4 −1/8 0 0 0 15/11 −15/11 0 0 0 −15/11 ∞
I4′ =
1 1 2 2 42/11 (X + 13/11)
n´ o eliminado n´ o eliminado
n´ o eliminado 3
2
1
6
4 5 n´ o eliminado
Os elementos da matriz que correspondem ao ramo 5-6 (que j´ a existia) tiveram seus valores alterados.
– 46 –
◮ Ap´ os a elimina¸c˜ ao do n´ o 5:
Y5′ =
2 0 0 0 0 0
0 −1 −1 0 0 2 −1 0 0 −1 0 2 −1/2 −1 −1/2 0 0 11/8 −5/4 −1/8 0 0 0 15/11 −15/11 0 0 0 0 ∞
I5′ =
1 1 2 2 42/11 (X + 5)
n´ o eliminado n´ o eliminado
n´ o eliminado 3
2
1
6
4 5 n´ o eliminado
n´ o eliminado
◮ O processo de elimina¸c˜ ao de n´ os at´e que reste somente um n´ o na rede corresponde ` a substitui¸c˜ ao Forward
– 47 –
Exerc´ıcio Obtenha a matriz dos fatores triangulares para a rede exemplo de 6 barras e 7 ramos, conforme mostrado na p´ agina 26.
◮ Pode-se tamb´em restaurar a rede a partir do conhecimento das matrizes elementares L. Por exemplo, pode-se restaurar o n´ o 2 atrav´es de: Y1′ = L−1 · Y2 2 2 0 −1 −1 0 0 0 2 −1 0 0 −1 0 −1 5/2 −1/2 −1 0 = 0 0 −1/2 3/2 −1 0 0 0 −1 −1 3 −1 0 −1 0 0 −1 ∞
Importante: a restaura¸c˜ ao de n´ os em princ´ıpio deve seguir exatamente a ordem inversa de elimina¸c˜ ao.
– 48 –
◮ Em cada o do processo de elimina¸c˜ ao de Gauss os seguintes os s˜ ao executados (do ponto de vista da redu¸c˜ ao de circuitos):
eliminar o n´ o
eliminar os ramos diretamente conectados ao n´ o eliminado
interligar todos os n´ os vizinhos do n´ o eliminado. Neste ponto, duas situa¸co ˜es podem ocorrer: • dois n´ os vizinhos j´ a conectados continuar˜ ao assim e os elementos da matriz correspondentes a esse ramo ter˜ ao seus valores alterados • dois n´ os vizinhos ainda n˜ ao conectados o ser˜ ao atrav´es de um ramo fict´ıcio. Consequentemente, aparecer˜ ao elementos correspondentes a esse ramo fict´ıcio na matriz, em posi¸co ˜es que anteriormente eram nulas (fill-ins)
Exerc´ıcio Realize a elimina¸c˜ ao de Gauss para a rede a seguir, mostrando as redes reduzidas em cada o. N˜ ao utilize matrizes na resolu¸c˜ ao do problema. 5 1 2
4 7
6
10
3 8
9
– 49 –
Exemplo A elimina¸c˜ ao de Gauss sob o ponto de vista de elimina¸c˜ ao de n´ os ´e utilizada na obten¸c˜ ao de equivalentes externos. Fronteira (F)
Rede Interna (I)
Rede Externa (E)
.. .
´ Area de Interesse Equivalente externo Fronteira (F)
Rede Interna (I)
shunt equivalente
liga¸c˜ ao equivalente .. .
inje¸c˜ ao equivalente
A id´eia ´e substituir a rede externa (E) por uma rede equivalente que reaja de maneira semelhante a dist´ urbios na rede interna (I). Os n´ os da rede externa (as inje¸co ˜es de potˆencia) s˜ ao eliminados (distribu´ıdas) atrav´es da elimina¸c˜ ao de Gauss.
– 50 –
Considerando que os n´ os s˜ ao numerados de maneira que aqueles a serem eliminados apare¸cam primeiro, tem-se:
1 .. . k k +1 0
.. .
Yr
n
Exemplo A transforma¸c˜ ao Y-∆ ´e um caso particular da elimina¸c˜ ao de Gauss sob o ponto de vista de elimina¸c˜ ao de n´ os.
4
4 eliminar n´ o1
1 2
3
2
– 51 –
3
Considerando o caso particular em que todas as itˆ ancias dos ramos do circuito Y s˜ ao iguais a 1 pu, tem-se:
3 −1 −1 −1 −1 1 0 0 Y= −1 0 1 0 −1 0 0 1
Eliminando o n´ o 1:
Y′ = L · Y 1000 3 −1 −1 −1 1/3 1 0 0 −1 1 0 0 = 1/3 0 1 0 · −1 0 1 0 1/3 0 0 1 −1 0 0 1 3 −1 −1 −1 0 2/3 −1/3 −1/3 = 0 −1/3 2/3 −1/3 0 −1/3 −1/3 2/3 Eliminando a linha e a coluna correspondente ao n´ o 1:
Yred
3 −1 −1 −1 0 2/3 −1/3 −1/3 = 0 −1/3 2/3 −1/3 0 −1/3 −1/3 2/3
que corresponde ` a matriz itˆ ancia do circuito ∆ reduzido. Atrav´es de Yred nota-se que as itˆ ancias dos ramos do circuito reduzido valem 1/3 pu, ou seja:
Y∆ =
YY 3
– 52 –
3.6
Algoritmo para elimina¸c˜ ao de Gauss
◮ Algoritmo para matrizes sim´etricas e assim´etricas: do k = 1 , (n-1)
alfa = a(k,k) a(k,k) = - 1.0 do i = k , n
(D e U)
a(k,i) = -a(k,i)/alfa end do do i = (k+1) , n
beta = a(i,k) a(i,k) = -beta/alfa
(L)
do j = (k+1) , n a(i,j) = a(i,j) + beta * a(k,j) end do end do
ao das atualiza¸c˜ demais colunas da linha i
end do → invers˜ ao do ´ ultimo elemento
a(n,n) = 1.0/a(n,n)
– 53 –
Exemplo Considere a matriz (3 × 3) a seguir.
a11
a12
a13
a21
a22
a23
a31
a32
a33
Executando o algoritmo para k = 1:
(U)
1/a11
−a12/a11 −a13/a11
−a21/a11
′ a22
′ a23
−a31/a11
′ a32
′ a33
(L)
atualiza¸c˜ ao
– 54 –
◮ Uma vez conhecida a matriz dos fatores triangulares e dado um vetor independente b, pode-se obter a solu¸c˜ ao do problema A · x = b por:
Substitui¸c˜ ao forward (por coluna) do k = 1 , n alfa = b(k) b(k) = b(k) * a(k,k) do i = (k+1) , n b(i) = b(i) + a(i,k) * alfa end do end do Substitui¸c˜ ao back (por coluna) do k = n , 2 , -1 do i = (k-1) , 1 , -1 b(i) = b(i) + a(i,k) * b(k) end do end do
– 55 –
Exerc´ıcio Considere a rede de 4 n´ os e 4 ramos a seguir. As itˆ ancias dos ramos em pu s˜ ao fornecidas na pr´ opria figura. O n´ o 4 ´e tomado como o n´ o de referˆencia.
1
2
1
4
1
1/2
1
3
(a) Obtenha a matriz itˆ ancia da rede. (b) Implemente a rotina de obten¸c˜ ao dos fatores triangulares apresentada anteriormente, e determine a matriz resultante de sua aplica¸c˜ ao ` a matriz do item (a). (c) Considere o seguinte vetor de correntes nodais: I=
1 1 1 X
T
Implemente as rotinas das substitui¸co ˜es Forward e Back apresentadas anteriormente e obtenha a solu¸c˜ ao do problema I = Y · V , em que V ´e o vetor das tens˜ oes nodais. Mostre o resultado intermedi´ ario ao final da substitui¸c˜ ao Forward.
– 56 –
◮ Referˆ encias
A.J. Monticelli, “Fluxo de Carga em Redes de Energia El´etrica”, Edgard Bl¨ ucher, 1983.
K. Zollenkopf, “Bi-factorization – basic computational algorithms and programming techniques”, In: J.K. Reid, Large sparse sets of linear equations, Academic Press, 1971.
M. Crow, “Computational Methods for Electric Power Systems”, CRC Press, 2003.
W.F. Tinney, J.W. Walker, “Direct solutions of sparse network equations by optimally ordered triangular factorization”, Proceedings of the IEEE, vol.55, n.11, 1967.
◮ Materiais interessantes
www.labspot.ufsc.br/˜campagno/Pos.html – Cap´ıtulo II – Partes 1 e 2
www.dm.ufscar.br/profs/silvia/AulaGradientesConjugados.pdf
www.eee.ufg.br/˜colemar/Espars01.pdf
www.ufmt.br/icet/matematica/geraldo/sela2.pdf
– 57 –