Wilson Lopes Di Mambro
DOS MÉTODOS DE REFATORIZAÇÃO DE MATRIZES NA
RESOLUÇÃO DE SISTEMAS MATRICIAIS
Pontifícia Universidade Católica de Minas Gerais
Belo Horizonte
2004
Wilson Lopes Di Mambro
DOS MÉTODOS DE REFATORIZAÇÃO DE MATRIZES NA
RESOLUÇÃO DE SISTEMAS MATRICIAIS
Dissertação apresentada ao Programa de Pós-Graduação em
Engenharia Elétrica da Pontifícia Universidade Católica de
Minas Gerais, Campus Coração Eucarístico, na Área de
Concentração em Sistemas de Potências, como requisito
parcial para obtenção do grau de mestre em Engenharia
Elétrica sob a orientação do Prof. Dr. Roberto de Maria Nunes
Mendes e a co-orientação do Prof. Dr. Luiz Danilo Barbosa
Terra.
Pontifícia Universidade Católica de Minas Gerais
Belo Horizonte
2004
Agradecimentos
À Deus por sempre me levar no colo quando o caminho é áspero.
Aos orientadores Prof. Dr. Roberto de Maria Nunes e Prof. Dr. Luiz Danilo
Barbosa Terra que tiveram paciência e dedicação em orientar-me nesta tarefa.
Ao Geraldo Magela dos Anjos Silva que com sua ajuda consegui digitar este
trabalho, tarefa que da qual tive muito que aprender.
Aos professores do mestrado que com dedicação, no transcorrer do curso, tiveram
papel importante para a realização deste trabalho.
3
Para meus filhos
Ana, Moisés e Mateus
4
RESUMO
Este trabalho apresenta métodos diretos para solução da equação matricial A x = b com A
esparsa e quadrada. O procedimento predominante para a solução será a decomposição de A
em um produto L D U evitando, sempre que possível, alterar a esparsidade da matriz A,
ordenando-a. Inicia-se o estudo pelo método de Gauss até, finalmente, obter-se os métodos de
refatorização parcial, os quais atualizam os fatores L D U da matriz A, para refletirem as
alterações em alguns de seus elementos. Discutem-se também, com mais detalhes, os métodos
de compensação e o método de Bennett, antes de se chegar aos métodos de refatorização
parcial I e II para matrizes esparsas. O resultado dos testes de refatorização parcial I e II
indica que eles são significativamente mais eficientes que as técnicas inicialmente
apresentadas e, que têm o potencial de acelerar a solução de diversos problemas.
ABSTRACT
This work presents direct methods to solve the matricial equation A x = b with A square
and sparse. The main way to solve it is supported by the decomposition of A in a L D U
product, avoiding to change the sparsity of the matrix, through the ordenation of A. At the
beginning of this work was used the Gauss method in order to get partial refatorization, that
updates the factories L D U of matrix A, showing the transformations of some elements of it.
Also, with details, was discussed both methods of compensation and the method of Bennett,
before arriving to the parcials I and II refatorization methods, for sparse matrices. The
conclusions for the parcial refactorization tests I and II show that they are more eficients than
the methods presented before and they can also accelerate the solutions of the problems.
LISTA DE FIGURAS
FIGURA 1 – Grafo de uma rede de energia elétrica com cinco nós .................................. 28
FIGURA 2 – Matriz associada ao grafo da figura 1............................................................ 28
FIGURA 3 – Árvore associada à matriz da Figura 2 partindo dos nós 1 e 2.................... 28
FIGURA 4 – Rede de energia elétrica com modificação de um ramo (1 - 4) ................... 31
FIGURA 5 – Rede de energia elétrica com modificações de dois ramos (1-4) e (3-5) ...... 32
FIGURA 6 – Estrutura da matriz admitância nodal: rede com cinco nós........................ 37
FIGURA 7 – Estrutura da matriz admitância nodal: rede de cinco nós, com a conexão
dos nós 4 e 5............................................................................................................................. 37
FIGURA 8 – Rede referente à matriz A da rede nodal modificada................................... 43
FIGURA 9 – Esquema de fatorização linha superior /coluna inferior..............................64
FIGURA 10 – Rede nodal com 12 nós .................................................................................. 65
FIGURA 11 – Estrutura esparsa da matriz L ..................................................................... 66
FIGURA 12 – Tabela de caminhos referentes à rede nodal da FIG. 10............................66
FIGURA 13 – Gráfico do caminho com 12 nós (árvore).....................................................66
FIGURA 14 – Rede de energia elétrica com 11 nós onde não foi feita uma ordenação
prévia dos mesmos .................................................................................................................. 68
FIGURA 15 – Estrutura da matriz L conseguida ao se fatorar a matriz A da rede nodal
acima, apresentando 3 fill-ins ................................................................................................ 69
FIGURA 16 – Grafo do caminho sem ordenação (árvore de 11 nós) ................................ 69
FIGURA 17 – Rede de energia elétrica com 11 nós onde foi feita a ordenação de grau
mínimo ..................................................................................................................................... 70
FIGURA 18 – Estrutura da matriz L ao se fatorar a matriz A, após a ordenação
observando o grau mínimo. Neste caso aparece apenas um fill-in..................................... 70
FIGURA 19 – Grafo após se fazer uma ordenação de grau mínimo (árvore c/ 11 nós) .. 71
SUMÁRIO
INTRODUÇÃO ........................................................................................................................ 9
CAPÍTULO 1 - FATORAÇÃO COMPLETA – MÉTODO DE GAUSS.......................... 13
1.1
Fatoração completa – método de Gauss...................................................................... 13
1.2
Algoritmo – Método de Gauss...................................................................................... 16
1.3
Exemplo.......................................................................................................................... 17
1.4
Decomposição L U e L D U ....................................................................................... 18
1.5
Exemplo.......................................................................................................................... 22
1.6
Variante do método de Gauss apresentada por W. F. Tinney e seus
∗
colaboradores ................................................................................................................ 23
1.7
Exemplo.......................................................................................................................... 25
1.8
Grafos............................................................................................................................. 26
CAPÍTULO 2 - MÉTODO DE COMPENSAÇÃO ............................................................. 30
2.1
Método de Compensação.............................................................................................. 30
2.2
Modificação por ramo orientado ................................................................................. 31
2.3
Modificação por nó orientado...................................................................................... 33
2.4
Exemplo
...................................................................................................................... 38
CAPÍTULO 3 - MÉTODOS DE REFATORIZAÇÃO PARCIAL ....................................43
3.1
Método de refatorização parcial simples .................................................................... 43
3.2
Decomposição L U em blocos ....................................................................................... 44
3.3
Exemplo.......................................................................................................................... 46
3.4
Método de refatorização parcial simples com arranjo especial................................ 48
3.5
Observações do que ocorre quando se decompõe A em L D U .............................. 48
∗
CAPÍTULO 4 - MÉTODO DE BENNETT.......................................................................... 50
4.1
Método de Bennett ....................................................................................................... 50
4.2
Prova do processo.......................................................................................................... 55
4.3
Exemplo.......................................................................................................................... 56
4.4
Algoritmo ....................................................................................................................... 58
4.5
Exemplo.......................................................................................................................... 59
CAPÍTULO 5 - MÉTODOS DE REFATORIZAÇÃO PARCIAL I E II.......................... 63
5.1
Métodos de refatorização parcial I e II ....................................................................... 63
5.2
Ordenação dos nós por grau mínimo .......................................................................... 64
5.3
Método de refatorização parcial I ............................................................................... 67
5.4
Método de refatorização parcial II.............................................................................. 68
5.5
Exemplo: (método de refatorização parcial I)............................................................ 68
CONCLUSÃO......................................................................................................................... 76
REFERÊNCIAS BIBLIOGRÁFICAS ................................................................................77
APÊNDICE I........................................................................................................................... 80
APÊNDICE II ........................................................................................................................ 78
9
INTRODUÇÃO
As áreas de estudo cujos problemas envolvem a solução de sistemas lineares são
várias, tais como programação linear, análise de estruturas, teoria dos grafos, solução
numérica de equações diferenciais, sistemas de transmissão e distribuição de energia elétrica
entre outras.
Este trabalho trata dos sistemas lineares da forma A x = b, onde são conhecidas as
matrizes A e b de tipos n x n e n x m, respectivamente. A matriz x, do tipo n x 1, é a solução
procurada.
Volta-se o interesse para sistemas de grande porte, onde em geral, a matriz A é
esparsa.
Concentra-se a atenção na solução de problemas de rede de energia elétrica onde A é a
matriz das admitâncias, x é o vetor das tensões e b o vetor das correntes. A pode ser também a
∆V
∆P
matriz jacobiana, J, do modelo linear de rede de energia elétrica = J , derivado da
∆θ
∆Q
formulação não linear de potências versus tensão. Sendo f(P,Q,V, θ ) = 0, onde f é uma função
vetorial.
Quando são projetadas futuras ampliações na rede, ou mesmo com a adição ou
remoção de linhas ou transformadores, haverá alterações das admitâncias; logo, a matriz A
será modificada. É importante observar que estas alterações levam a novos cálculos, que
acarretam mudanças no tempo de computação e no espaço de memória exigidos, obrigando a
procura de métodos eficientes para a solução de problemas com matrizes modificadas.
10
Como este assunto não se encontra disseminado nos livros, procura-se descrever
didaticamente os vários métodos e compará-los entre si. Com o desenvolvimento científico,
alguns métodos ficaram obsoletos, porém, têm caráter seminal, donde o interesse no seu
estudo.
Neste trabalho são tratados apenas os métodos diretos.
Foram feitos estudos comparando o desempenho do métodos diretos e iterativo na
solução de equação Ax = b, empregando-se as mesmas condições de teste geral. A conclusão
final é que os métodos diretos ainda são mais rápidos, embora o método iterativo seja mais
fácil de se implementar [Monticelli,99].
O objetivo geral deste trabalho é organizar didaticamente num único texto os métodos
diretos de refatorização de matrizes na solução de sistemas matriciais.
Os objetivos específicos são:
• Dar condições para usar o método adequado verificando as vantagens que
existem em cada um.
• Verificar as vantagens de se usar L D U em vez de usar a inversa A −1 .
• Usar o método de Bennett, [Bennett,65] empregando o programa
implementado computacionalmente.
A metodologia adotada consiste na pesquisa bibliográfica, no seu estudo e seleção. A
pesquisa bibliográfica concentra-se em publicações do IEEE, como também em livros de
Álgebra Linear.
A contribuição deste trabalho inclui a interpretação e registro sistemático de
bibliografia pertinente e aplicações dos métodos de refatorização de matrizes para
planejamento de sistemas de energia elétrica. Além do mais, o trabalho pode contribuir como
incentivo ao desenvolvimento de novas pesquisas.
11
Organização do texto
A presente Dissertação está organizada em cinco capítulos:
Capítulo 1:
O capítulo apresenta o método de Gauss para solução de sistema A x = b como,
também, a decomposição de A em L D U e uma variante do método de Gauss feita por W.F.
Tinney. Ele introduz o conceito de grafos associados à matriz A.
Capítulo 2:
0 capítulo considera as modificações da matriz A, modificações estas feitas por nó ou
ramos orientados. A partir da equação matricial (A+ ∆A ) x = b, chega-se à formula
x = [ A −1 − A −1 ( I + ∆A. A −1 ) −1 ∆A. A −1 ] b através de deduções matemáticas. A idéia principal é
tornar o cálculo de x rápido e eficiente. Mostra-se que basta armazenar A −1 uma única vez e
usar este fato a cada modificação efetuada.Verifica-se que pode-se armazenar também as
inversas de L e U para calcular o valor de x. Finaliza-se com exemplos.
Capítulo 3.
O capítulo mostra que o método de refatorização simples é mais conveniente, pois, o
método de refatorização com arranjo especial degrada a esparsidade da matriz A. Mostra-se
que para a solução do sistema deve-se decompor a matriz A em L D U por blocos. É feita a
decomposição da matriz em blocos, finalizando-se com um exemplo, aplicando o método.
12
Capítulo 4.
O capítulo mostra que John M. Bennett substituindo ∆A por X C Y t
na matriz
( A + ∆A) −1 chegou à identidade ( A + X C Y t ) = A −1 X [C −1 + Y t A −1 X ] −1 Y t A −1 .
Verifica-se que não é necessário armazenar a inversa da matriz a cada modificação. É
suficiente armazenar a inversa de A uma única vez, o que é fundamental para a computação.
É desenvolvido um algoritmo para chegar à decomposição da matriz A + X C Y t em
L• D • U • usando o L D U da matriz A, e é mostrado que a inversa completa da matriz
modificada pode ser evitada a cada alteração. Isto permite armazenar apenas os fatores L D U
da matriz A.
Capítulo 5.
O capítulo mostra que ao se alterar uma ou mais linhas/colunas de uma matriz é
afetado somente um subgrupo de linhas/colunas subseqüentes. Mostra-se a vantagem de
ordenar uma matriz antes de fatorá-la em LU, pois proporciona menos preenchimentos
nas posições, onde a matriz A possía inicialmente elementos nulos, ao se fazer a
decomposição. A associação que se faz de uma matriz com um grafo mostra as
vantagens que existem ordenando a matriz antes de fatorá-la, pois o número de* “fillins“ é menor com a ordenação prévia.
Finalmente, apresenta-se a conclusão, as referências bibliográficas e o apêndice.
* “fill-ins” – elementos não nulos.
13
CAPÍTULO 1
FATORAÇÃO COMPLETA – MÉTODO DE GAUSS
Neste capítulo serão discutidos os métodos diretos para solução de um sistema Ax = b,
entre eles o método de Gauss como também uma variante deste método apresentada por W.F.
Tinney [Tinney,67], o qual decompõe a matriz A em L U. Introduz-se também a associação
de uma figura, isto é, um grafo linear (GA) com uma matriz A. Tal associação se faz
necessária para os estudos dos capítulos posteriores.
1.1 Fatoração completa – método de Gauss
Para resolver o sistema linear A x = b, onde A é n x n, x é n x 1 e b é n x 1, usa-se o
método de eliminação de Gauss [Leon, 98]. Este método consiste em aplicar à matriz
completa C = [A / b] as chamadas transformações elementares sobre as linhas, a saber:
a) Tij , trocar de posição as linhas i e j;
b) Tij (λ ) , somar à linha i a linha j multiplicada por
λ ∈ℜ.
14
Supõe-se que a matriz A seja tal que apenas a transformação (b) seja necessária. Uma
condição suficiente para isso, é que A seja simétrica definida positiva ou que A = ( aij ) seja
diagonalmente dominante ( a jj >
n
∑a
ij
, j = 1,..., n ) .
i =1, i ≠ j
Neste método, resolve-se a equação linear A x = b através de operações elementares
efetuadas sobre as linhas da matriz A e b, até que seja obtido um sistema equivalente ao dado,
de solução imediata, ou pelo menos mais simples.
A decomposição triangular da matriz por Gauss [Leon,98]é feita pela eliminação dos
elementos abaixo da diagonal principal em todas as colunas. O método é direto, porém, é
lento. O tempo de computação é grande e ocupa também grande espaço de armazenamento.
Como as matrizes de admitância são esparsas e simétricas, pode-se obter algumas
vantagens na solução dos problemas. A primeira, e mais evidente, é a grande economia de
memória possível de ser obtida quando armazenam-se apenas os elementos não nulos. A
segunda vantagem é a redução do tempo de computação. Uma terceira é a redução do erro de
arredondamento.
O método de Gauss tende a destruir a esparsidade da matriz a ele submetida.
Entretanto, um controle pode ser exercido sobre o processo de triangularização, o que não
pode ser feito no caso de matrizes cheias.
Antes que se dê um exemplo de aplicação do método de Gauss, deve-se comentar três
esquemas que levam à solução da equação A x = b, para justificar a opção por este método.
a) Inverter explicitamente a matriz dos coeficientes, no caso a matriz A.
A x = b ⇒ x = A −1 . b
Este esquema é inconveniente porque:
15
a.1) a inversão da matriz A exige um número de operações bastante grande em relação a
outros esquemas;
a.2) a inversão destrói a esparsidade das matrizes;
a.3) o acúmulo de erros de arredondamentos, durante o processo de inversão, pode
inutilizar os resultados obtidos, sobretudo quando se trata de sistemas de grande
porte.
b) Admitir uma solução aproximada e obter, através de algoritmo adequado, aproximações
melhores, até que um determinado índice de precisão seja atingido. Este esquema tem
vantagens, mas é inconveniente. Quando se armazenam apenas os elementos não nulos da
matriz de coeficientes, o armazenamento é compacto e é feito por linhas e colunas, isto é,
os elementos ficam distribuídos aleatoriamente, dificultando a sua localização para se obter
a matriz transposta de A ou a inversa da transposta. Apresenta-se também excessivo tempo
de computação, quando houver necessidade de soluções repetidas, pois cada solução
envolve a repetição com novos dados de todo o processo iterativo.
Em suma, os métodos iterativos têm desvantagens em relação aos métodos
diretos[Monticelli,99] , como também a inversão implícita da matriz A.
c) Transformar o sistema linear dado, através de operações aritméticas elementares, até que
seja obtido um sistema equivalente, mas de solução imediata, ou pelo menos mais simples
(métodos diretos). O método direto principal é o método de Gauss que apresenta mais
vantagens sobre os outros, pois é mais preciso e a esparsidade da matriz A é menos
afetada.
Considera-se o seguinte algoritmo:
16
1.2 Algoritmo – Método de Gauss
Seja A = (aij ) nxn
Notações:
− índices inferiores indicam ordem de linha e coluna do elemento respectivamente;
− índices superiores representam o estado atual do elemento, isto é, o número de
transformações sofridas por este elemento no decorrer de um processo qualquer;
− índice superior 0 (zero) indica que o elemento se encontra com o seu valor inicial.
1) Fazendo k ← 1 (iniciação do contador de linha);
2) Soma-se a linha k de A, multiplicada por − a (jkk −1) às linhas j = k+1, ..., n, ou seja,
calculam-se:
a (jkk ) = a (jkk −1) − a (jkk −1) .a ki( k ) :
j ,i = k+1,..., n ;
b (j k ) = b (j k −1) − a (jhk −1) .b (j k )
Isto resulta na eliminação de x k das (n-k ) equações subseqüentes à k-ésima.
3) Fazendo k ← k + 1 e indo para 2.
4) Substituindo os valores dos xi (i = k+1, ..., n), já calculados na k-ésima equação e
resolvendo para x k . Isto equivale a se calcular: x k =
n
1 (k )
(bk − ∑ a ki( k ) xi )
a kk
i = k +1
5) Fazendo k ← k − 1 . Se k= 0, FIM. Senão, → 4
Este algoritmo de Gauss, normalmente apresentado na literatura de ágebra linear
[Leon ,98], exige que a matriz a ser triangularizada seja armazenada integralmente na
17
memória do computador. O armazenamento linha por linha, no caso de matrizes esparsas,
dificulta a localização dos elementos de uma dada coluna.
Pode-se verificar que o número de operações aritméticas envolvido no método de
Gauss é menor que o número de operações necessárias à inversão de uma matriz. Em se
tratando de uma matriz esparsa, as vantagens são bem mais significativas.
1.2 Exemplo
Aplicação do algoritmo a um sistema de quatro incógnitas.
2 x1
− 2 x
1
− 4 x1
− 2 x2
− 4 x3
+ 4 x2
− 6 x2
+ 6 x3
− 8 x3
=1
− 6 x4
=0
− 8 x4
+ 18 x 4
=2
=1
Pode-se escrever o sistema na forma matricial A x = b.
Como as operações de eliminação envolvem apenas os coeficientes aij ( i ,j = 1,2,3,4),
e os termos independentes bi ( i = 1,2,3,4), para fim de exposição, será suficiente trabalhar
com a matriz aumentada [ A / b] .
2 − 2 − 4 0 1
− 2 4
0 − 6 0
[ A / b]=
, Aplicando o algoritmo tem-se:
6 − 8 2
−4 0
0 − 6 − 8 18 1
18
2 − 2 − 4 0 1
0 2 − 4 − 6 1
T (2); T (3) →
42
0 − 4 − 2 − 8 4 32
0 − 6 − 8 18 1
2 − 2 − 4 0 1
− 2 4
0 − 6 0
T21 (1); T31 (2) →
− 4 0
6 − 8 2
0 − 6 − 8 18 1
0
2 − 2 − 4
0 2
−4 −6
0 0 − 10 − 20
0
0 0 − 20
x1
x
Cálculo de x= 2 :
x3
x4
x4 =
x2 =
1
1
T43 ( −2) →
6
4
sendo
1
b4
=−
5
a44
x3 =
2 x1
0
1
2 − 2 − 4
0 2 − 4 − 6
1
0 0 − 10 − 20 6
0
40 − 8
0 0
− 2 x2
− 4 x3
2 x2
− 4 x3
− 6 x4
=1
− 10 x3
− 20 x4
40 x4
=6
= −8
1
(b3 −
a33
4
1
1
(b2 − a2i xi ) = a22
2
i =3
∑
4
∑a
3i xi
3
5
)= − +
i =4
x1 =
=1
2
1
=−
5
5
4
1
2
(b1 − a1i xi ) = a11
5
i=2
∑
1.4 Decomposição L U e L D U ∗ .[ Monticelli,83]
A pode ser reduzida à forma triangular superior:
a1n
a11
0
a2 n
a3n → A(1) = 0
...
...
0
ann
a11 a12
a
21 a22
A = a31 a32
... ...
an1 an 2
a13
a23
a33
...
an 3
...
...
...
...
...
a11
0
= 0
...
0
a12
a13
(1)
a22
0
(1)
a23
( 2)
a33
...
...
...
...
0
an( 23)
...
→A
( 2)
A(n ) = U
...
...
a12
(1)
a22
(1)
a32
...
an(12)
a13
(1)
a23
(1)
a33
...
an(13)
a11
a1n
(1)
a2 n
0
( n)
( 2)
a3n .... → A = 0
...
...
( 2)
0
ann
...
...
...
...
...
a1n
a2(1n)
a3(1n)
...
(1)
ann
a12
a13
...
(1)
a22
...
0
(1)
a23
( 2)
a33
...
...
...
0
0
...
...
a1n
a2(1n)
a3( 2n)
...
(n)
ann
19
Em geral, se a matriz A n x n poder ser colocada na forma triangular superior U sem
trocar as linhas, então A pode ser fatorada em L U, onde L é triangular inferior e tem todos os
elementos da diagonal iguais a um (1). O elemento lij de L abaixo da diagonal principal vai
ser um múltiplo da i-ésima que foi subtraído da j-ésima linha durante o método de Gauss.
Pode-se verificar que L U = A.
Para se ver como esta fatoração funciona, analisa-se o processo em termos de matrizes
elementares. Isto é equivalente a multiplicar a matriz A à esquerda por n matrizes
elementares: E1, E2 ,..., En. Então, En ...E3 E2 E1 A = U e como as matrizes elementares são
invertíveis, tem-se: A = ( E1−1.E2−1.E3−1. ... .En−1 ).U , logo A = L U.
Dada a fatoração L U de uma matriz A, é possível prosseguir e fatorar U em um
produto D U ∗ , onde D é a diagonal e U ∗ é triangular superior, com todos os elementos da
diagonal iguais a 1.
U11 0
0 u
22
D U∗ = 0
0
... ...
0
0
0
...
0
...
u33 ...
...
...
0
...
0
0
0
...
unn
u12
1 u
11
0 1
0 0
... ...
0 0
u13
u11
u 23
u 22
...
...
1
...
...
...
0
...
u1n
u11
u2n
u22
u3 n
u33
...
1
∗
Tem-se A = L D U .
Tem-se, portanto, que:
1) A matriz U é ∆ - superior e contém os elementos quando se escalona a matriz A (método
de Gauss).
20
2) A matriz L é ∆ - inferior e representa as etapas de eliminação. Seus elementos lij , i<j são
os simétricos dos multiplicadores λ que aparecem nas operações Tij (λ ) ; os elementos da
diagonal de L são todos iguais a 1.
3) Além disso, a solução de A x = b é equivalente à solução dos dois sistemas L y = b e U
x = y, sendo o primeiro por substituição progressiva e o segundo por substituição
regressiva.
4) Assim, o método de eliminação de Gauss é equivalente a fatorar A no produto L U, onde
L é ∆ - inferior com diagonal unitária e U é ∆ -superior.
∗
∗
5) Introduzindo a matriz diagonal D, pode-se escrever A = L U = L D U , onde U tem
também diagonal unitária; esta decomposição é única. [Leon 98, cap. 6].
Pode-se mostrar que o número de operações aritméticas envolvido na solução de A =
L U, L y = b , U x = y é o mesmo que o envolvido no método de eliminação de Gauss.
∗
Na prática, o método L U (ou L D U ) oferece algumas vantagens.
Por exemplo, se ao resolver A x = b para vários valores de b, A permanece inalterada,
o método L U é o mais conveniente.
Caso de matriz simétrica:
t
Supõe-se que na equação linear A x = b, A é matriz simétrica ( isto é A= A ).
∗
t
Se A = L D U = A = U
∗t
D Lt , resulta U ∗ = Lt ,de modo que A= L D Lt .
Por exemplo,
2 − 2
4
A = 2 − 8 − 4
− 2 − 4 3
1
1
=
2
− 1
2
0
1
1
3
1
0 4 0 0
0 0 − 9 0 0
0 0 3
0
1
1
2
1
0
1
−
2
1
3
1
21
Seja U uma matriz tiangular superior com todos os elementos da diagonal principal
diferentes de zero, e L, uma matriz triangular inferior com os elementos da diagonal principal
iguais a um.
Decompondo A em L U e substituindo-os no sistema A x = b, tem-se:
L U x = b , e fazendo Ux = y, tem-se Ly = b.
Como L é uma matriz tringular inferior, resolve-se, diretamente, calculando y por
substituição de cima para baixo. Uma vez determinado o valor de y, e usando a matriz
triangular superior U, resolve-se a outra equação U x = y por substituição de baixo para cima.
Sejam as matrizes:
1 0
l
21 1
l31 l32
l
l
L= 41 42
.
.
.
.
.
.
l n1 l n 2
0
0
0
0
1
0
l 42
1
.
.
.
.
.
.
l n3
ln 4
0 0 ... 0
u11 u12
0 u
0 0 ... 0
22
0
0
0 0 ... 0
0
0
0 0 ... 0
,U =
.
.
. . . .
.
. . . .
.
.
.
. . . .
0
. . . 1
0
u13
u14
u 23
u 24
u33
u34
0
u 44
.
.
.
.
.
.
0
0
. . . u1n
. . . u 2 n
. . . u 3n
. . . u4n
. . . .
. . . .
. . . .
. . . u nn
Para resolver um sistema A x = b, faz-se as seguintes considerações:
y1
y
2
.
A = L U; L y = b; U x = y toma-se y =
.
.
y n
j −1
e calcula-se y usando a fórmula y j = b j − ∑ l jk y k ,
k =1
22
x1
x
2
.
e tomando-se x = , calcula-se x,
.
.
xn
usando a fórmula
xj=
n
1
(b j −
u jk xk )
u jj
k = j +1
∑
1.5 Exemplo
2 − 2 − 4 0 x1
1
− 2 4
0
0 − 6 x2
=
− 4 0
2
6 − 8 x3
0 − 6 − 8 18 x4
1
2 − 2
2 −2 −4 0
0 2
− 2 4
0 − 6 T21 (1)
A=
→
0 − 4
− 4 0
6 − 8 T31 (2)
0 − 6
0 − 6 − 8 18
−4
0
− 4 − 6 T32 (2)
→
− 2 − 8 T42 (3)
− 8 18
0
0
2 − 2 − 4
2 − 2 − 4
0 2
− 4 − 6
−4 −6
0 2
T (−2) →
=U
0 0 − 10 − 20
0 0 − 16 − 20 43
0
0
40
0 0
0 0 − 20
0
1
−1 1
Então : L =
− 2 − 2
0 −3
L y = b, tem-se:
0 0
0 0
1 0
2 1
0
1
−1 1
− 2 − 2
0 −3
0 0
0 0
1 0
2 1
1
y1
y
2 = 0
y3
2
1
y4
1
1
logo y =
6
− 8
23
0
2 − 2 − 4
0 2 − 4 − 6
U x = y , tem-se:
0 0 − 10 − 20
0
40
0 0
x1
1
x
2 = 1
6
x3
− 8
x4
2
− 5
1
−
logo, x = 2
− 1
5
1
−
5
∗
Como a matriz A pode também ser fatorada em L D U , onde L e U
∗
são matrizes
triangulares inferior e superior, respectivamente, com elementos da diagonal principal
unitários e D é uma matriz diagonal, tem-se também:
D U ∗ = Uentão L D U
∗
∗
x = b , DU x = y e L y = b.
1.6 Variante do método de Gauss apresentada por W. F. Tinney e seus
colaboradores [Tinney,67;85].
Para exemplificar, considera-se um sistema A x = b de quatro equações e quatro
∗
variáveis. Decompondo A em L D U tem-se:
1 0
l
21 1
l31 l32
l 41 l 42
0
0
1
l 43
0
0
0
1
d11
0
0
0
0
0
d 22
0
0
0
d 33
0
0
0
0
d 44
1 u12
0 1
0 0
0 0
u13
u 23
1
0
u14
u 24
u34
1
x1 b1
x b
2 = 2
x3 b3
x 4 b4
∗
Fazendo y = U x, resolve-se o sistema L D y = b, usando a substituição progressiva
por colunas (de cima para baixo) e determinando o vetor y.
Coluna 1
(1)
1
y
b
= 1
d11
y 2(1) = b2 − l 21b1
Coluna 2
y
( 2)
2
y 2(1)
=
d 22
y 3( 2 ) = y 3(1) − l 32 y 2(1)
Coluna 3
y
( 3)
3
y 3( 2)
=
d 33
y 4(3) = y 4( 2) − l 43 y 3( 2)
Coluna 4
y
( 4)
4
y 4(3) ( 4) y 4( 3)
=
y4 =
d 44
d 44
24
y 3(1) = b3 − l 31b1
y 4( 2 ) = y 4(1) − l 43 y 2(1)
y 4(1) = b4 − l 41b1
y1(1)
( 2)
y
A solução do sistema L D y = b é y= 2( 3)
y3
( 4)
y 4
De fato, tem-se:
1) d 11 y1 = b1
y1 =
b1
= y1(1)
d11
2) l 21 d11 y1 + d 22 y 2 = b2 ; d 22 y 2 = b2 l 21 d 11 y1 ; y 2 =
b2 − l 21 d11 y1 b2 − l 21 y1(1) y 2(1)
=
=
= y 2( 2)
d 22
d 22
d 22
3) l31 d11 y1 + l32 d 22 y 2 + d 33 y 3 = b 3 , d 33 y 3 = b3 − l31b1 − l32 y 2(1) = y 3(1) − l32 y 2(1)
y3 =
(1)
y 32
− l 32 y 2(1) y 3( 2)
=
= y 3( 3)
d 33
d 33
4) l 41 d11 y1 + l 42 d 22 y 2 + l 43 d 33 y 3 + d 44 y 4 = b4
d 44 y 4 = b4 − l 41l 32 d11 y1 − l 42 d 22 y 2 − l 43 d 33 y 3
d 44 y 4 = b4 − l 41b1 − l 42 y 2
(1)
− l 42 y 3
(2 )
d 44 y 4 = y 4(1) − l 42 y 2(1) − l 42 y 3( 2 )
d 44 y 4 = y 4( 2 ) − l 42 y 3( 2)
d 44 y 4 = y 4(3)
(3 )
y
(4 )
y4 = 4 = y4
d 44
Agora, resolve-se o sistema U ∗ x = y , usando a substituição por linhas.
25
1 u12
0 1
0 0
0 0
u13 u14 x1
y1
y
u23 u24 x2
= 2
y3
1 u34 x3
0
1 x4
y4
y1
y
2
y3
y4
linha 4 : x4 = y 4
linha 3 : x3 = y 3 − u 34 x 4
linha2: x2 = y2 − u12 x3 − u24 x4
linha1: x1 = y1 − u12 x2 − u13 x3 − u14 x4
1.7 Exemplo
2 x1
− 2 x
1
Seja o sistema
4
x
−
1
− 2 x2
− 4 x3
+ 4 x2
− 6 x2
=1
− 6 x4
=0
+ 6 x3
− 8 x4
=2
− 8 x3
+ 18 x4
=1
O sistema pode ser colocado na forma de uma equação matricial A x = b.
∗
Decompondo A em L D U :
0
1
−1 1
− 2 − 2
0 −3
0 0
0 0
1 0
2 1
2
0
0
0
0
2 0
0
0 − 10 0
0
0
40
0
0
1 − 1 − 2 0 x1 1
0 1 − 2 − 3 x 0
2 =
0 0
1
2 x3 2
0
1 x4 1
0 0
Calcula-se primeiro y, resolvendo o sistema L D y = b(por colunas de cima para baixo)
0
1
−1 1
− 2 − 2
0 −3
0 0
0 0
1 0
2 1
2
0
0
0
0
2 0
0
0 − 10 0
0
0
40
0
0
y1 1
y 0
2 =
y3 2
y 4 1
coluna 1 ) y1(1) =
1
; y 2 (1) = 0 + 1.1 = 1 ;
2
coluna 2) y 2 (2 ) =
1
; y 3 ( 2 ) = 4 + 2 .1 = 6 ; y 4 ( 2 ) = 1 + 3 .1 = 4
2
coluna 3) y3 (3) =
6
3
( 3)
= − ; y 4 = 4 − 2.6 = 8
5
− 10
y3
(1)
= 2 + 2 .1 = 4 ; y 4
(1)
= 1 + 0.1 = 1
26
coluna 4 ) y 4 ( 4) =
8
1
=−
5
− 40
∗
Resolve-se o sistema U x = y para encontrar o x, que é a solução final:
1 − 1 − 2 0
0 1 − 2 − 3
0 0
1
2
0
1
0 0
x4 = −
x1 =
1
5
1
x1 12
x
2 = 2
x3 3
− 5
x4 1
−
5
;
tem-se:
3
1
1
x3 = − − 2(− ) = −
5
5
5
;
x2 =
1
1
1
1
+ 2(− ) + 3(− ) = −
2
5
5
2
1
1
1
2
+ (− ) + 2(− ) = − .
2
2
2
5
Então a solução é:
2
x = −
5
1
−
2
1
−
5
1
−
5
t
1.8 Grafos
Neste trabalho introduz-se um procedimento para se estudar o efeito da esparsidade na
eliminação de Gauss, que habilita a escolher-se um esquema de eliminação ótimo e/ou
determinar se a eliminação é uma abordagem prática. Isto é feito pela associação de um grafo
linear (GA) a uma matriz A. Admite-se que a matriz A envolvida é suficientemente bem
condicionada de modo que a eliminação de Gauss seja sempre possível. Uma rede de energia
elétrica pode ser usada como sendo seu próprio grafo não direcional [Szwarcfiter ,83].
No sistema A x = b associado a uma rede de energia elétrica, a estrutura de A está
intimamente ligada à configuração geométrica da rede, ou seja, à maneira pela qual as
diversas barras estão interconectadas, independendo complemente dos valores dessas
conexões. Desta forma, é suficiente para descrever a estrutura geométrica de uma rede
27
substituir seus componentes por segmentos de reta. Estes segmentos são denominados ramos
e seus terminais nós. Um nó e um ramo são incidentes quando o nó for um dos terminais do
ramo. Os nós podem incidir em um ou mais ramos. Ao conjunto de nós e ramos que
descrevem a estrutura topológica de uma rede dá-se o nome de grafo.
Quando um nó é eliminado, a conexão que existia antes precisa ser preenchida. Isto
pode requerer a adição de novo ramo no grafo reduzido. Os novos ramos correspondem aos
preenchimentos (fill-ins) na matriz.
Cada nó é considerado uma ligação própria (um laço), isto é, tem admitância própria,
representada por aij onde i = j. Se um nó é ligado a outro, tem-se a admitância mútua
representada também por aij , mas com i ≠ j . O nó i é ligado ao nó j e vice-versa. Quando
não há ligações entre os nós, o termo é igual a zero.
Um caminho de um nó v a um nó w é uma cadeia formada por ramos distintos que têm
origem em v e término em w. O caminho pode ser fechado, se v = w, ou aberto, caso
contrário.
Um ciclo é um caminho fechado vi ,..., vk , sendo vi = vk e k ≥ 3.
Um grafo é conexo, quando existe um caminho entre cada par de nós, isto é, a partir de
um nó arbitrário do grafo é possível alcançar todos os demais. Caso contrário, o grafo é
desconexo.
Alterando-se uma matriz, têm-se modificações nos grafos e serão afetadas suas
linhas/colunas. As linhas/colunas afetadas podem ser encontradas partindo do caminho do
grafo da matriz. Quando mais de uma linha/coluna é modificada, as linhas/colunas afetadas
são originadas pela união dos caminhos [Betancourt ,86;88]. Um grafo é chamado de árvore
se ele for conexo e acíclico.
28
1
5
2
3
4
Figura 1 – Grafo de uma rede de energia elétrica com cinco nós
X
0
A = X
0
X
0
X
0
X
0
X
0
X
X
X
X
X
X
X
0
X
X
X
0
X
FIGURA 2 – Matriz associada ao grafo da figura 1
1
2
3
4
5
FIGURA 3 – Árvore associada à matriz da Figura 2 partindo dos nós 1 e 2
Neste capítulo foi apresentado o método de Gauss para solução do sistema A x = b, a
decomposição de A em L D U e uma variante do método de Gauss feita por W. F. Tinney .
Foi introduzido o conceito de grafos associado à matriz A e dado um exemplo tanto do
método de Gauss como o de Tinney.
29
O CAPÍTULO 2
MÉTODO DE COMPENSAÇÃO
O capítulo apresenta a dedução e análise do método de compensação para a solução
eficiente de problemas de rede envolvendo modificação na matriz A do sistema A x = b.
Mostra-se que as modificações podem ser feitas por ramo ou nó orientados. No final, chega-se
a uma equação que se altera a cada modificação, mas calcula-se e armazena-se a inversa de A
somente uma vez [Alsaç, Stott, Tinney, 83].
2.1 Método de Compensação
Este método torna o processo de cálculo da solução de problemas de redes de energia
elétrica mais rápido e eficiente, e é mais apropriado, quando:
a) as modificações não são grandes;
b) as modificações não são permanentes;
c) as equações modificadas não precisam ser resolvidas repetidamente .
Considera-se a seguinte equação matricial ( A+ ∆A) x = b, sendo A, a matriz
admitância da rede de energia elétrica de dimensão n x n, ∆A também de dimensão nxn a
matriz modificação envolvendo um ou mais elementos da rede, x o vetor das tensões e b o
vetor das correntes.
31
Para todas as modificações incidentes simétricas da matriz A, toma-se ∆A = M∆YM t ,
onde:
∆ Y é uma matriz m x m contendo as modificações de A , sendo m ≤ n .
M é a matriz conexão, de dimensão n x m.
Quando a modificação é por ramo orientado, M tem colunas de entradas +1 e –1 nas
posições relevantes. No caso de modificações por nó orientado, as colunas de M têm uma
única entrada +1.
2.2 Modificação por ramo orientado
Quando m ramos são modificados simultaneamente, ∆Y é uma matriz diagonal m x m
de mudança da admitância, e M, como já foi dito, tem m colunas, com as entradas +1 e –1
nas posições relevantes.
Primeiro,considera-se apenas um ramo modificado (FIG. 4).
2
1
5
∆y1
3
4
y11
y
21
A = y31
0
0
y12
y13
0
y 22
0
y 24
0
y33
y 44
y 42
y 43
y 44
y52
0
y54
0
y 25
0
y 45
y55
FIGURA 4 – Rede de energia elétrica com modificação de um ramo (1 - 4)
∆y1
+ 1
0
0
∆A = 0 [∆y1 ][+ 1 0 0 − 1 0] = 0
− ∆y1
− 1
0
0
0 0 − ∆y1 0
0 0
0
0
0 0
0
0
0 0 ∆y1 0
0 0
0
0
32
y11 + ∆y1
y
21
A + ∆A = y31
− ∆y1
0
y12
y13
− ∆y1
y 22
0
y 24
0
y33
y34
y 42
y 43
y 41 + ∆y1
y52
0
y54
0
y 25
0
y 45
y55
Observa-se que a admitância entre os nós i e j muda de ∆y1 , adicionando ∆y1 a y ii e a
y jj e subtraindo-o de y ij e y ji .
Consideram-se dois ramos modificados (FIG. 5):
1
2
5
∆y1
3
4
∆y
y
∆
2
2
FIGURA 5 – Rede de energia elétrica com modificações de dois ramos (1-4) e (3-5)
+ 1 0
0
0
∆y
∆A = 0 + 1 1
0
− 1 0
0 − 1
∆y1
0
∆A = 0
− ∆y1
0
0 + 1 0 0 − 1 0
= M∆YM t
∆y2 0 0 + 1 0 − 1
0
0
− ∆y1
0
0
0
0
∆y 2
0
0
0
∆y1
0 − ∆y 2
0
0
0
− ∆y 2
0
∆y 2
Logo, a matriz com dois ramos modificados será:
33
y11 + ∆y1
y
21
A + ∆A = y31
− ∆y1
0
y12
y13
− ∆y1
y 22
0
y 24
0
y33 + ∆y 2
y34
y 42
y 43
y 44 + ∆y1
y52
− ∆y 2
y54
y 25
− ∆y 2
y 45
y55 + ∆y 2
0
Observa-se que os ∆ y são adicionados a
y11 , y 33 , y 44 ey55 e subtraídos de
y14 , y 41 , y 35, y 53 que, particularmente, são nulos, porque não há ligações entre estes nós.
2.3 Modificação por nó orientado
Esta é uma modificação mais geral, pois podem ser feitas simultaneamente alterações
em todos os elementos da matriz associados a um dado conjunto de nós. A matriz ∆Y tem a
dimensão m x m, onde m é o número de nós envolvidos. A matriz M tem m colunas cada uma
com uma única entrada +1 na posição relevante como se vê no exemplo a seguir.
Considera-se a mesma matriz A vista no exemplo anterior .
Quando dois nós são envolvidos, tem-se:
1
0
∆A = 0
0
0
0
0
∆y1
0.
− ∆y1
1
0
∆y1
0
− ∆y1 1 0 0 0 0
= 0
.
∆y1 0 0 0 1 0
− ∆y1
0
0 0 − ∆y1
0 0
0
0 0
0
0 0
∆y1
0 0
0
0
0
0
0
0
∆A é a mesma matriz encontrada quando a modificação foi feita no ramo 1-4 (FIG. 4).
Quando quatro nós são envolvidos, considerando a matriz A dada no exemplo anterior
e as modificações que agora envolverão os nós 1, 3, 4 e 5, tem-se:
34
1
0
∆A = 0
0
0
0 0 0
∆y1
0 0 0
− ∆y1
0 1 0.
0
1 0 0
0
0 0 1
∆y1
0
∆A = 0
− ∆y1
0
0
0
0
0
0 ∆y 2
0
0
0 − ∆y 2
− ∆y1
∆y1
0
0
0
0
∆y 2
− ∆y 2
− ∆y1
0
0
∆y1
0
0 1
0 0
.
− ∆y 2 0
∆y 2 0
0
0
0
0
0
0
1
0
0
0
0
1
0
1
0
0
0
0
− ∆y 2
0
∆y 2
Mais uma vez, conclui-se que ∆ A continua sendo a mesma anterior, quando se fez a
modificação nos ramos 1-4 e 3-5 (FIG. 5).
Sendo a matriz A uma matriz de admitância de uma rede de energia elétrica, ela será
simétrica e esparsa, principalmente se tratando de um sistema de grande porte. A matriz
modificação ∆A é ainda mais esparsa e poderá também ser decomposta num produto LU,
onde L é uma matriz triangular inferior com elementos unitários na diagonal principal e U é
uma matriz triangular superior:
∆A =
∆y1
0
0
− ∆y1
0
0
0
0
0
0
0 − ∆y1
0
0
0
0
0 ∆y1
0
0
0 1
0 0
0 = 0
0 − 1
0 0
0
1
0
0
0
∗
0
0
1
0
0
0
0
0
1
0
0
0'
0
0
1
∆y1
0
0
0
0
0
0
0
0
0
0 − ∆y1
0
0
0
0
0
0
0
0
∗
0
0
0
0
0
U = D U , onde D é matriz diagonal e U é uma matriz triangular superior com os elementos
∗
da diagonal iguais a um. Então, ∆ A = L D U .
1
0
∆ A = 0
− 1
0
0
1
0
0
0
0
0
1
0
0
0
0
0
1
0
0 ∆y1
0 0
0. 0
0 0
1 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 1
0 0
0.0
0 0
0 0
0
1
0
0
0
0 −1
0 0
1 0
0 1
0 0
0
0
0
0
1
35
O que determina o melhor método a ser usado (nó orientado ou ramo orientado) é o
gasto computacional necessário para o processamento de cada um deles. Este gasto aumenta
com a ordem da matriz modificação ∆Y .
Inversa da matriz modificada
Para resolver a equação (A + ∆A) x = b, deve-se calcular a inversa ( A + ∆A) −1 .
Tem-se:
( A + ∆A) −1 = A−1 ( I + ∆A. A−1 ) −1
e desenvolvendo-se formalmente em série de potências,
( I + ∆A. A−1 ) −1 = I − ∆A. A−1 + (∆A. A−1 ) 2 − (∆A. A−1 )3 + ...
Portanto,
[
]
( A + ∆A) −1 = = A−1 ( I + ∆A. A−1 ) −1 = A−1. I − ∆A. A−1 + (∆A. A−1 ) 2 − ( ∆A. A−1 )3 + ... =
[
]
= A −1 − A −1 ∆A. A −1 − (∆AA −1 ) 2 + (∆A. A −1 ) 3 − ... =
[
]
= A −1 − A −1 . I − (∆A. A −1 ) + (∆A. A −1 ) 2 − (∆AA −1 ) 3 + ... .∆A. A −1 =
= A −1 − A −1 ( I + ∆A. A −1 ) −1 .∆A. A −1 ,
[
]
então, x = A −1 − A −1 ( I + ∆A. A −1 ) −1 ∆A. A −1 .b
Se ∆A = M∆YM t , então, ( A + ∆A) −1 = A−1 − A−1 ( I + M∆YM t A−1 ) −1 M .∆YM t . A−1
Simplificando, mais a identidade acima, toma-se:
C= ( I + ∆YM t A −1 M ) −1 .∆Y , donde:
∆Y = ( I + ∆YM t A −1 M ).C , logo:
36
M∆Y = ( M + M∆YM t A −1 M ).C
M∆Y = ( I + M∆YM t A −1 ).M .C , donde :
MC = ( I + M∆YM t A −1 ) −1 M∆Y e , tem − se :
MCM t = ( I + M∆YM t A −1 ) −1 M∆YM t
Então, tem-se que: ( A + ∆A) −1 = A−1 − A−1.M .C.M t . A−1
Pondo :
X = A −1 .M
Z = M t A −1 M = M t X
Vem: C = ( I + ∆Y .Z ) −1 ∆Y ,
Após a modificação a solução é: x = ( A−1 − A−1M .C.M t . A−1 ).b .
Pode-se
escrever também: x = A−1 ( I − M .CM t . A−1 ).b . Se
A = LU, vem A −1 = U −1 .L−1 ,
então:
x = U −1 L−1 ( I − M .C.M t . A −1 ).b
x = U −1 ( L−1 − L−1 .M .C.M t .U −1 .L−1 ).b
x = U −1 ( I − L−1 .M .C.M t .U −1 ) L−1 .b
A idéia é tornar o processo de cálculo de x rápido e eficiente. A inversa de A será
calculada somente uma vez e armazenada. Verifica-se que as inversas, que se devem calcular,
são altamente esparsas e as outras operações necessárias são somente adição e multiplicação
de matrizes.
37
São mostradas abaixo as estruturas da admitância nodal com cinco nós, antes e depois
de se fazer a modificação, isto é, antes de ligar o nó 4 ao nó 5. Pode-se verificar como ficaram
também as estruturas das matrizes admitâncias antes e depois da modificação.
4
5
2
3
1
FIGURA 6 – Estrutura da matriz admitância nodal: rede com cinco nós
5
4
2
3
1
FIGURA 7 – Estrutura da matriz admitância nodal: rede de cinco nós, com a conexão
dos nós 4 e 5
38
1
X
X
X
2
X
X
X
3
4
X
X
X
X
X
X
X
5
1
2
1
2
X3
4
X 5
X
X
X
X
X
Estrutura da matriz admitância: A
X
3
4
X
X
X
X
X
X
X
X
5
1
2
X3
X4
X 5
Estrutura da matriz admitância: A + ∆A
2.4 Exemplo
Considerando-se as estruturas das figuras 6 e 7 dá-se o exemplo seguinte, atribuindo
valores numéricos;
0 x1 0
6 −3 −3 0
− 3 6
0 − 3 0 x2 6
− 3 0
6 − 1 − 2 x3 = 0
0 x4 − 2
0 − 3 −1 5
0
0 −2 0
3 x5 − 3
Considera-se uma rede de energia elétrica de 5 nós, portanto, a matriz A é 5 x 5 e as
matrizes x e b são 5 x 1 cada uma. Pode-se também, neste caso, ter cada matriz x e b até 5
x 5, isto é, se A é de ordem n x n, as matrizes x e b podem ser n x m, onde m ≤ n . É
interessante observar que se têm várias combinações possíveis, principalmente se n é um
número grande.
0
6 −3 −3 0
− 3 6
0 − 3 0
A= − 3 0
6 − 1 − 2 e decompondo A em L U, tem-se:
0
0 − 3 −1 5
0
0 −2 0
3
39
1
1
− 2
1
−
L= 2
0
0
0
0
0
1
0
0
1
0
1
3
2
−
3
−
0
1
2
1
−
2
−
6 − 3 − 3
5
3
0 2 − 2
4
U = 0 0
0
0 0
0
0 0
1
−
1
2
0
0
0
0
1
0
−3 0
− 2 − 2
2 − 1
3
0
2
e
0
e
1
1
2
2
−1
L = 3
2
3
2
3
1
6
0
−1
U = 0
0
0
1
9
2
9
0
0
0
0
1
0
0
1
0
1
3
5
6
7
12
1
6
1
12
1
4
0
0
0
0
1
2
3
4
1
3
5
12
1
4
1
2
0
1
1
2
0
0
0
0
1
4
9
7
18
1
2
1
3
2
3
Considera-se, por exemplo, ∆Y = [1] fazendo uma ligação do nó 4 ao nó 5, e
0
0
como foi visto acima a matriz de conexão é M = 0 e, como C = [I + ∆Y Z ]−1 ∆Y , tem-se:
1
− 1
2
23 19
27 27
3
19 89
7
27 108 12
2
7
3
t −1
Z = M A M = [0 0 0 1 − 1]
3
12
4
5
11 1
18
2
9
7
1
4
18
2
9
3
2
C = ( [1] + [1]. ) −1 [1] =
3
5
5
4
9 9
11 7 0
18 18 0
2
1 1
0 =
2 2 3
2 1 1
3 3 − 1
1
2
3 3
40
0
0
0
MCM t =
0
0
0 0
0
0 0
0
0
0
0
3 .
−
5
3
5
0 0
0
3
0 0
5
3
0 0 −
5
E:
x = [U −1 ( I − L−1MCM −1U −1 ) L−1 ] b, sendo que x é o vetor das tensões e b o vetor das correntes.
Resolvendo-a por partes, tem-se:
L−1MCM t =
1
1
2
2
3
2
3
2
3
0
0
0
1
0
0
1
0
1
3
5
6
7
12
1
2
3
4
1
1
2
0
0
0
0
0 0
0
0
0
1
0 0
0
0 0
0
0 0
0
3
0 0
5
3
0 0 −
5
0 0
0 0
0 0
3 =
− 0
5
3
0
5
0 0
0
0 0
0
0
0
0
3
−
5
3
10
0 0
0
3
0 0
5
3
0 0 −
10
L−1MCM tU −1 =
0
0
0
0
0
0 0
0
0 0
0
0 0
0
3
0 0
5
3
0 0 −
10
1
0 6
0 0
0
3 0
−
5
3 0
10
0
1
0
0
I − L−1 M C M t U −1 =
0
0
(
)
U −1 I − L−1MCM tU −1 =
1
9
2
9
0
1
6
1
12
1
4
0
0
0
0
0
0 0
0
1 0
0
0 1
1
3
5
12
1
4
1
2
0
7
0 0
10
3
0 0
20
4
9 0
7
0
18
1 0
=
2 0
1
3 0
2
3
0
0
0
1
5
9
10
0 0
0
0 0
0
0 0
0
3
0 0
10
3
0 0 −
20
0
0
0
1
−
5
1
10
41
1
6
0
0
0
0
1
9
2
9
1
6
0
0
0
0
0
1
6
1
12
1
4
4
9
7
18
1
2
1
3
2
3
1
5
5
12
1
4
1
2
0
0
0
0
0
1
9
2
9
1
6
1
12
1
4
3
10
7
20
1
4
2
5
1
10
0
0
0
0
0
1
0
0
0
0
0 0
0
1 0
0
0 1
0
7
0 0
10
3
0 0
10
0
0
0
1 =
5
9
10
7
15
13
30
1
2
2
5
3
5
U −1 ( I − L−1MCM tU −1 ) L−1 = ( A + ∆A) −1
1
6
0
0
0
0
1
9
2
9
0
1
6
1
12
1
4
0
0
0
0
3
10
7
20
1
4
2
5
1
10
7
1
15 1
13
2
30 2
1
3
2 2
2
3
5 2
3
3
5
0
0
0
1
0
0
1
0
1
3
5
6
7
12
1
2
3
4
1
1
2
38 31 2
0
3
45 45
0 31 143 7
45 180 12
0 2
7
3
=
3
12
4
0 8
17
1
15 30
2
1 7
13
1
2
15 30
Agora, pode-se calcular as tensões, usando a equação:
x = {U −1 ( I − L−1MCM tU −1 ) L−1} b
x = ( A + ∆A ) −1 b
38 31 2
45 45 3
31 143 7
45 180 12
2
7
3
x =
3
12
4
8
17
1
2
15 30
13
1
7
15 30
2
8
15
17
30
1
2
3
5
2
5
7
15
5
13 0 2
30 6 7
1
0 = 3
2 1
2 − 2
1
5 − 3
0
3
5
8
15
17
30
1
2
3
5
2
5
7
15
13
30
1
2
2
5
3
5
42
5
2
7
3
As tensões serão, portanto: x1 = , x2 = , x3 = 1, x4 = 1 e x5 = 0
Para cada caso específico de modificação da rede, o processo pode ser dividido em
duas fases:
Fase preparatória: cálculo da inversa de A, cálculo de L e U e suas inversas e cálculo das
matrizes C e Z;
Fase solução: encontrar o vetor x através da equação: x = U −1 ( I − L−1MCM tU −1 ) L−1 b .
Para outras modificações no sistema não será necessária a inversão da nova matriz da
rede, a matriz inicial A será invertida e armazenada uma única vez como também as inversas
de L e U e as eventuais modificações serão calculadas conforme as fases acima. Como se
pode observar, efetua-se simplesmente a multiplicação de matrizes a cada nova modificação.
Neste capítulo foram consideradas modificações da matriz A, modificações estas feitas
por nós ou ramos orientados. A partir da equação matricial ( A + ∆A) x = b, chegou-se à
fórmula x = [ A −1 − A −1 ( I + ∆A. A −1 ) −1 ∆AA −1 ] b através de deduções matemáticas. A idéia
principal é tornar o cálculo de x rápido e eficiente. Isto é possível, pois basta armazenar
A −1 uma única vez e usar este fato a cada modificação efetuada. Verificou-se que pode-se
armazenar as inversas de L e U que também permitem o cálculo de x da equação matricial.
CAPÍTULO 3
MÉTODOS DE REFATORIZAÇÃO PARCIAL
Quando o número de alterações da matriz A no sistema A x = b não é pequeno, as
modificações são permanentes e não há necessidade de se resolver o sistema repetidamente,
os métodos de refatorização parcial são mais eficientes.
A refatorização parcial é aplicável a qualquer matriz A não singular, mas neste
capítulo considera-se que ela seja a matriz dos coeficientes de uma rede nodal.
Vê-se neste capítulo como se faz a decomposição L U por blocos para calcular a
refatorização da submatriz que foi modificada.
3.1 Método de refatorização parcial simples
Este método refatora somente a submatriz que contém os elementos modificados da
matriz. Os elementos da matriz fora da submatriz não são alterados[Brandwajn,86].
O problema deste método é que o tamanho da submatriz, e consequentemente o
trabalho de refatorização são influenciados pela posição dos elementos modificados na matriz
A. Se os elementos modificados estiverem próximos do topo da matriz, pouca ou nenhuma
vantagem é alcançada neste método.
y11
y
21
0
Seja A =
0
y51
0
y12
0
0
y15
y 22
0
y 24
0
0
y33
y34
0
y 42
y 43
y 44
0
0
0
0
y55
y62
y63
0
y65
1
0
y26
y36
0
y56
y66
2
4
∆y
5
6
FIGURA 8 – Rede referente à matriz A da rede nodal modificada
3
44
y 44
Seja B = 0
0
0
y55
y56
0
y56 a submatriz de A que contém os
y66
elementos que serão modificados, e eles se encontram na parte inferior, o que é vantajoso,
pois pode-se refatorar somente a submatriz que neste caso é de ordem menor que a matriz A.
1
∆y 0 − ∆y
∆B = 0 [∆y ] [1 0 − 1] = 0
0
0 e tem-se portanto:
− 1
− ∆y 0 ∆y
y 44 + ∆y
B + ∆B = 0
− ∆y
0
y55
y65
− ∆y
y56
y66 + ∆y
3.2 Decomposição L U em blocos
Para se fatorar apenas a submatriz B, faz-se a decomposição da matriz A em blocos:
A
A = 11
A21
A12
A22
onde:
A é n x n, invertível, A11 é r x r, invertível, A12 é r x (n-r) e A21 é (n-r) x (n-r).
Prova-se que:
Ir
A = LU = __
A21 A11−1
Demonstração
Tem-se:
|
__
|
0
__
I r
A11
__
0
|
__
|
A12
__ , onde C = A22 − A21 A11−1 A12
C
45
| 0 A11 | A12 A11
Ir
__
| __ __ | __ = __
:
A21 A11−1 | I r 0 | C A21
|
__
|
____________ = A
A21 A11−1 A12 + C
A12
−1
donde: C = A22 − A21 A11 A12
A11
A21
Supõe-se que: A + ∆A =
A12
, ou seja, modifica-se apenas o bloco A22
B22
| 0 A11 | A12
Ir
Então A + ∆A = __
| __ __ | __
A21 A11−1 | I r 0 | D
−1
Onde: D = B22 − A21 A11 A12
•
isto é, A + ∆A = L U
•
•
com L = L
A11 | A12
U = __ | __
0 | D
•
e
De ( A + ∆A) y = b vem L U • y = b . Pondo z = U • y , vem L z = b.
Assim, para resolver ( A + ∆A) y = b , acha-se A= L U ( em bloco ) e resolve Z1
, onde
Z 2
se L z = b para achar z. Achando-se z , resolve-se U • y = z =
Z1 é r x1 e Z 2 é (n-r) x 1.
A11
U = __
0
•
|
__
|
A12
__ e
D
−1
D = B22 − A21 A11
A21
Y1
Se y = ___ onde Y1 é r x 1 e, Y2 é (n-r) x 1. Então, U • y = z dá:
Y2
A11Y1 + A12Y2 = z1
. Após achar , Y2 tem-se
DY2 = z 2
−1
Y1 = A11
( z1 − A12Y2 ) ,
46
Em D Y2 = Z 2 , faz-se
D = LoU o ,
U oY2 = w eLo w = Z 2
calculado w pode − se achar Y2 , resolvendo U oY2 = w
3.3 Exemplo
Seja a resolução do sistema ( A + ∆A) y = b onde foram feitas modificações apenas
no bloco A22 .
0
6 −3 −3 0
− 3 6
0 − 3 0
A11
Seja A = − 3 0
6 − 1 − 2 =
A21
0
0 − 3 −1 5
0
0 −2 0
3
A12
, onde:
A22
0
6 − 3 − 3
0
5 0
0 − 3 − 1
A11 = − 3 6
e A22 =
0 , A12 = − 3 0 , A21 =
0 0 − 2
0 3
− 3 0
6
− 1 − 2
Tem-se:
1 1
3 6
1 1
A11−1 =
6 4
1 1
6 12
1
2
6
−
1
−1
; A21 A11 = 3
1
12
−
1
2
4
5
6
1
−
6
−
1
−
2 e A A−1 A = 3 1
1 21 11 12 1 1
−
2
Considera-se que a modificação foi a ligação do nó 4 ao nó 5, onde ∆y = 1 , então
0
6 −3 −3 0
− 3 6
0 − 3 0
6 − − 1
( A + ∆A) = − 3 0
6 − 1 − 2 , com B22 =
−1 4
0
3
1
6
1
−
−
−
0
0 − 2 − 1 4
decompondo A + ∆A em L• U • (em blocos) pode-se escrever a equação:
L• U • y = b , mas como L• = L tem-se LU • y = b e fazendo U • y = z ,vem L z = b.
47
Cálculo de z, na equação L z = b;
1
0
0
2
−
3
− 1
2
0
1
0
5
−
6
1
−
6
0
0
1
1
−
2
1
−
2
0 0
z1 0
0 0
z2
6
0 0
z3 = 0
1 0
z
−2
4
0 1 z5 − 3
0
6
então z = 0
3
− 2
Como U • y = z , passa-se ao cálculo de y , onde :
A11
U = __
0
•
|
__
|
A12
__
D
D = B22 − A21 A11−1 A12
0
6 −3 −3 | 0
− 3 6 − 3 | − 3 0 y1 0
y 6
2
−
−
−
−
3
6
|
1
2
U•y =
y3 = 0
__ __ __ | __ __ y 3
4
0
0
0 | 3 − 2
y5 − 2
0
0 | − 2 3
0
D Y2 = Z 2 , fatorando D em Lo U o , pode-se resolver a equação LoU o y = Z 2 e U o y = w
3 − 2
5
Uo =
0 3
1
Lo = 2
− 3
0
1
1
Z2 = 2
− 3
Lo w = Z 2
0 w 3
1 =
,
1 w2 − 2
3 − 2 y 4 3
0 5 y = 0
5
Uoy = w
w1 = 3
w2= 0
y4 = 1
y5 = 0
Y1 = A11−1 ( Z1 − A12Y2 )
1
y1 3
y = 1
2 6
y3 1
6
1
6
1
4
1
12
5
3
1
5
3
6 0 0
0
7
1
1
( 6 − − 3 0 ) = ,
0
12
3
1 0 − 1 − 2
1
4
7
3
portanto, a solução será: y1 = , y2 = , y3 = 1, y4 = 1 e y5 = 0
48
3.4 Método de refatorização parcial simples com arranjo especial
Este método supera o problema do método anterior, forçando os elementos que
precisam ser modificados para a parte inferior da matriz. O arranjo especial aumenta a
eficiência da refatorização; porém, ele pode degradar a esparsidade e aumentar o trabalho
computacional para operações progressivas e regressivas. Além disso, este método restringese a aplicações onde as posições dos elementos da matriz modificada podem ser previstas
antecipadamente [Brandwajn ,86] .
3.5 Observações do que ocorre quando se decompõe A em L D U ∗
∗
Na transformação de A em L D U , a matriz triangular inferior L e a matriz triangular
superior U
∗
não mantêm o mesmo grau de esparsidade da matriz A. Aparecem em L
elementos não nulos lij ( i < j ) onde aij correspondente de A é nulo. Estes elementos lij que
surgem em L são os “fill-ins”, preenchimentos que correspondem a novas ligações ao se
eliminar alguns nós.
Para se conseguir minimizar o número de “fill-ins”, existem técnicas que podem ser
usadas em conjunto com a fatoração de A, procurando preservar a estrutura esparsa de L.
As matrizes de incidência simétricas podem ter uma estrutura corretamente descrita
por grafos não orientados, se aij ≠ 0 implica que a ji também é diferente de zero.
Sempre que houver um segmento orientado de i para j haverá, necessariamente, outro
orientado de j para i, fato que pode ser descrito por um segmento não orientado.
49
O capítulo três mostra que o método de refatorização simples é mais conveniente, pois
o método de refatorização com arranjo especial degrada a esparsidade da matriz A. Mostra-se
que para a solução do sistema deve-se decompor a matriz A em L D U por blocos, também,
como se decompõe por blocos uma matriz e, finaliza-se com um exemplo aplicando o
método.
CAPÍTULO 4
MÉTODO DE BENNETT
O capítulo mostra um método engenhoso, devido a Bennett [Bennett,65] , utilizado
∗
para atualizar os fatores L D U de uma matriz modificada. Chega-se à atualização dos
fatores, sem alterar a matriz A do sistema A x = b. Nele também pode-se ver o algoritmo que
∗
permite estas alterações a partir das inversas das matrizes L .D e U .
4.1 Método de Bennett
Quando adiciona-se ∆A à matriz A, pelo método de composição, conclui-se que:
( A + ∆A) −1 = A−1 − A−1 ( I + ∆A A−1 ) −1 ∆A A−1 . Com a substituição de ∆A por XCY t , chega-se a
outra identidade, efetuando as operações abaixo:
( A + XCY t ) −1 =
A −1 − A −1 ( I + XCY t A −1 ) −1 XCY t A −1 =
A −1 − A −1 [ ( XC ) −1 ( I + XCY t A −1 )
[
A −1 − A −1 ( XC ) −1 + Y t A −1
]
−1
]−1 Y t A −1 =
Y t A −1 =
A −1 − A −1 [C −1 X −1 + Y t A −1 ] −1 Y t A −1 =
A −1 − A −1 [ (C −1 + Y t A −1 X ) X −1
[
A −1 − A −1 X C −1 + Y t A −1 X
]
−1
]−1 Y −1 A −1 =
Y t A −1
51
Pode-se ver como fica a inversa de uma matriz modificada quando acrescenta-se à
matriz o termo XCY t . Não é necessário armazenar a inversa da matriz a cada modificação,
basta apenas armazenar a inversa de A uma única vez; com isto, o número de operações
diminui, o que é fundamental para a computação.
Supõe-se que A seja decomposta no produto L D U. Se A for modificada pela adição
de XCY t , é desejável também ser capaz de modificar L D U, de modo a se obter a
•
•
•
decomposição L D U de A + XCY t .
Para se evitar a inversa completa da matriz modificada, toda vez que há uma alteração,
trabalha-se com os fatores triangulares.
•
•
•
Para iniciar o processo do cálculo de L D U da matriz modificada, introduze-se
as anotações:
Li → representa a matriz obtida da matriz unidade I, pela substituição da i-ésima coluna de I
pela i-ésima coluna de L.
U i → representa a matriz obtida da matriz unidade I, pela substituição da i-ésima linha de I
pela i-ésima linha de U.
K i ∗ → é a i-ésima linha de K
K ∗i → é a i-ésima coluna de K
Também são válidas as propriedades:
n
A) L = I n + ∑ ( Li − I n )
i =1
De fato, tem-se:
52
1 0 0
l
21 1 0
1=
L = l31 0 1
... ... ...
ln1 0 0
... 0
1 0 0 ...
0 1 0 ...
... 0
... 0 , L2 = 0 l32 1 ...
... ...
... ... ... ...
0 ln 2 0 ...
... 1
0
1
0
0
0 , ..., Ln = 0
...
...
0
1
0
1
0
...
0
0 ... 0
0 ... 0
1 ... 0
... ... ...
0 0 1
Logo:
0 0
l
n
21 0
i
( L − I n ) = l31 l32
i =1
... ...
ln1 ln 2
0
0
∑
0
...
ln 3
... 0
... 0
... 0 e, portanto: L = I n +
... ...
... 0
∑ (L − I )
n
i
n
i =1
B) ( Li ) −1 = 2 I n − Li ,
Se A é n x n e B é a matriz obtida de A pela aplicação da transformação elementar T,
isto é , B = T(A), então B = E A, onde E é a matriz T( I n ).
1
0
...
0
j
Tem-se: L = .
0
0
...
0
0 ...
0
0 ... 0
1 ...
0
0 ... 0
... ...
...
... ... ...
0 ...
1
0 ... 0
0 ... l j +1. j
0 ... l j +2, j
1 ... 0
... ...
...
... ... ...
0 ...
ln , j
0 ... 0
0 ... 0
0
0
...
0
0
0
...
1
Portanto:
L j = Enj (ln , j )...E j +1, j (l j +1, j ) I n , onde Eij (lij ) é a matriz elementar obtida de
linha i a linha j multiplicada por lij .
I n somando à
53
1
0
...
0
j −1
Donde ( L ) = E j +1, j (−l j +1. j )...Enj (−lnj ) I n
0
0
...
0
0 ...
0
0 ... 0
1 ...
0
0 ... 0
... ...
...
... ... ...
0 ...
1
0 ... 0
0 ... − l j +1, j
0 ... − l j +2, J
1 ... 0
... ...
...
... ... ...
0 ...
− ln , j
0 ... 0
0 ... 0
Exemplo: para j = 1
2 I n − L1 =
2
0
0
...
0
0
2
0
...
0
0 0
... 0 1
... 0 − l21 1 0
... 0 = − l31 0 1
... ... ... ... ...
... 0 − ln1 0 0
0 ... 0 1 0 0
0 ... 0 l21 1 0
2 ... 0 − l31 0 1
... ... ... ... ... ...
0 ... 2 Ln1 0 0
... 0
... 0
... 0 = ( L1 ) −1
... ...
... 1
Valem também as propriedades:
n
C) U = I n + ∑ (U i − I n )
i =1
D) (U i ) −1 = 2 I n − U i
∗
Seja A uma matriz de ordem n fatorada em um produto L D U .
a11
a
21
A = a31
...
an1
d11
0
D= 0
...
0
a12
a13
a22
a23
a32
a33
...
...
an 2
an 3
0
0
d 22
0
0
...
d 33
...
0
0
... a1n
... a2 n
... a3n
... ...
... ann
,
1 0
l
21 1
L = l31 l32
... ...
ln1 ln 2
0
1 u12
0 1
0 0
0 0 e U ∗ = 0 0
... ...
... ...
0 0
0 d 55
0
0
...
0
...
1
...
...
...
ln3 ...
u13
u23
1
...
0
0
0
0 ,
...
lnn
... u1n
... u2 n
... u3n
... ...
... 1
0
0
...
0
= 2I n − Lj
0
0
...
1
54
Fazendo: d11 = a11 ,
1
l
21
l31
1
obtém-se L∗1 =
A∗1 = .
a11
.
.
l n1
e
U 1∗ =
1
A1∗ = [1 u12
a11
u13 ... u1n ] .
a
F
Além disso, se A é particionada em submatrizes tem-se A = 11
G H
a11
Então, ( L ) A(U ) = 0
∗1 −1
1 −1
a11
A =L
0
1
0
e tem-se ainda:
1
H−
GF
a11
0
U ∗1
1
H−
GF
a11
Referindo-se A como A (1) , escreve − se A ( 2 ) = H −
1
GF .
a11
Então:
1 0
l
21 1
(1)
A = l31 l32
... ...
ln1 ln 2
0
0
1
...
ln 3
... 0
... 0
a
... 0 11
0
... ...
... 1
1 u12
0 1
0
0 0
A( 2)
... ...
0 0
As operações acima repetidas em A
( 2)
u13
u 23
1
...
0
... u1n
... u 2 n
... u3n
... ...
... 1
fornecem L(∗22) ,U 2( 2∗) e d 22 .
∗
Com repetições adicionais, completa-se o processo e determina-se L D U .
55
4.2 Prova do processo
Tem-se:
1 a11
L
0
1
− −
l21
= .
.
.
l
n1
1
− − −
l 21
0 ∗1
U = .
A( 2 )
.
.
l
n1
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0
− −
a11 | a12
− − | − −
I n−1 0 |
0
− − −
a11
|
− − − −
|
I n−1 0
...
−−
A
( 2)
a11
− −
a1n a21
− − = .
.
.
.
a21a12
l21a12 ... l21a1n a
11
Onde K = ... ... ... = ...
ln1a12 ... ln1a1n an1a12
a11
a21
.
G = . e F = [a12 . . . a1n ] .
.
a2 n
A( 2 ) = H −
0 1
|
u12 ...
− − − − − − − − − −
A( 2 ) 0
|
I n−1
| a12
...
| −−
−−
|
|
|
K + A( 2 )
|
|
a1n
− −
a21a1n
a11
1
...
... =
GF pois
an1a1n a11
...
a11
...
Portanto H = A( 2) +
1
GF . Assim ,
a11
a11 |
a11 | 0
0
1 −1
∗
1 −1
−−
( L ) A(U ) = − − |
= − − | − −
0 | H − 1 GF 0 | A( 2)
a11
1
GF e
a11
u1n
− − =
56
4.3 Exemplo
A= A
(1)
4 − 6
2
∗
= − 4 − 7 10 . Decompondo A no produto L D U tem-se :
2
7 − 9
1 0 0
1 2 − 3
2 0 0
∗
L = − 2 1 0 , U = 0 1 − 2 e D = 0 1 0
1 3 1
0 0 1
0 0 3
Toma-se:
d11 = 2 ,
U1• =
L1•1
2 1
1
1
A•1 = − 4 = − 2 e
=
a11
2
2 1
1
1
A1• = [2 4 − 6] = [1 2 − 3]
a11
2
2
A=
G
F
− 4
− 16 24
, F = [4 − 6] e G = ; log o GF =
H
− 12
2
8
− 7 10
− 7 10 − 8 12 1 − 2
1
H =
,H −
GF =
−
=
= A( 2 )
a11
7 − 9
7 − 9 4 − 5 3 − 3
( 2)
d 22 = 1 = a11
H−
F = [− 2]
G = [3]
1
GF = −3 + 6 = 3 = A(3)
( 2)
a11
1
L(•21 ) =
3
U1(•2) = [1 − 2]
então:
H = [− 3] e GF = [− 6]
57
2 0 0
1 2 − 3
1 0 0
∗
L = − 2 1 0 , U = 0 1 − 2 eD = 0 1 0
0 0 3
0 0 1
1 3 1
No procedimento com a matriz modificada A + XCY t , a primeira linha e a primeira
•
•
•
•
coluna levam imediatamente a L•1 , d11 e U 1• , que são a primeira coluna de L , o primeiro
•
•
termo da matriz D e a primeira linha de U , considerando-se que A + XCY t = L• D • U • .
d11•
0
.
B ( 2)
Fazendo-se A + X C Y t = B , tem-se: B (1) =
0
Como é possível obter L D U de A conhecendo-se
•
•
•
A(1) , A( 2 ) ,..., A( n ) , o mesmo
t
pode ser feito para L , D eU da matriz A + X C Y = B.
Exemplo:
Para i = 1
0
1 β
, (U • )1 =
I n−1
0 I n−1
1
( L• )1 =
α
( )
( L• ) −1 B (1) U •
−1
d •
= 11
0
0
A
( 2)
+X
( 2)
C
( 2)
(Y
( 2) t
)
t
Se XC = D, YC = E , então:
( DY t )11 =
m
∑
D1k Ykt1 =
k =1
m
∑
Y1k (Ckt j ) X 1 j =
j ,k =1
m
m
∑ ∑
Y1k
k =1
m
∑E
j =1
t
1 j X j1
X 1 j C jk =
j =1
m
m
∑∑ Y
1k
X 1 j C jk =
j =1 k =1
(
= X CY
t
)
11
;
d11• = a11 + Y1• C t X 1t• e como a11 = d11 tem − se que d11• = d11 + Y1• C t X 1t•
Por conveniência computacional, incluem-se os vetores:
58
p (i ) = (C (i ) ) t ( X 1(•i ) ) t ou ( p i ) t = X 1(•i )C (i )
q (i ) = C (i ) (Y1(•i ) )t
ou (q (i ) ) t = Y1(•i ) (C (i ) ) t
4.4 Algoritmo
Define − se C 1 = C , X 1 = X , Y 1 = Y , p 0 = 0, q 0 = 0
i =1
→
( p i ) t = X 1i• C i ;
p i → p i −1
d ii = Dii + Y1i• p i ;
d ii → Dii
→→→→→→→→→→ i < n
; i = n, fim
↓
↓
q i → q i −1
(q i ) t = Y1i• (C i ) t ;
j =1
→→→→→→→→→→
X ij•+1 = − Li + j ,i X 1i• + X ij +1,• ;
X ij+•1 → X ij +1,•
Y ji•+1 = − X i ,i + j Y1i• + Y ji+1,• ;
Y ji•+1 → Y ji+1,•
L•i + j ,i = Li + j ,i + (
U i•,i + j = U i ,i + j (
1
) X ij+•1q i ;
d ii
L•i + j ,i → Li + j ,i
1 i +1 i
)Y j • p ;
d ii
U i•,i + j → U i ,i + j
j+i → j
j < n − i +1
j = n − i +1
↑←←←←←←
↓
C i+ j = C i − (
1 i i t
)q ( p ) ; C i +1 → C i
d ii
↓
i +1→ i
retorna
o processo,
59
Se A e C forem simétricas e X=Y, o algoritmo se simplifica e o número de operações é
dividido em duas partes. Isto acontece nos problemas que envolvem as redes de energia
elétrica, pois as matrizes admitâncias são simétricas, X é a matriz de conexão M e
Y t = M t , isto é , Y = M = X .
4.5 Exemplo
Sejam as matrizes
4 − 6
2
A = − 4 − 7 10 ,
2
7 − 9
X
(1)
1 2
= X = 0 − 3 ,
8 6
1 0 0
L = − 2 1 0 ,
1 3 1
C
(1)
2 3
=C =
,
1 1
Cálculo de X ( 2) e Y ( 2 )
X (j•i +1) = − Li + j ,i X 1(•i ) + X (ji+)1,•
Y j(•i +1) = −U i ,i + j Y1•( i ) + Y j(+i1) ,•
X 1(•2) = 2[1 2]+ [0 − 3] = [2 1]
X 2(•2) = −1[1 2] + [8 6] = [7 4]
Y1•( 2) = −2[8 − 5] + [17 − 10] = [1 0]
Y2(•2) = 3[8 − 5] + [− 26 16] = [− 2 1]
1 0
Y ( 2) =
− 2 1
2 1
X ( 2) =
7 4
•
Cálculo de: L•1
e U 1••
2 0 0
D = 0 1 0 ,
0 0 3
Y
(1)
1 2 − 3
U = 0 1 − 2
0 0 1
−5
8
= Y = 17 − 10
− 26 16
60
L•21 = −2 +
L•31 = 1 +
1
1
[2 1]. = − 13
9
9
3
1
1
[7 4]. = 28
9
3 9
U 12• = 2 +
4
1
[1 0] = 22
9
5 9
U 13• = −3 +
4
1
[− 2 1]. = − 10
9
3
5
1
13
L••1 = −
9
28
9
22
U 1•• = 1
9
−
10
3
Cálculo de C ( 2 ) , p ( 2) , q ( 2) , d 22 , X (3) e Y (3)
C
( 2)
14
2 3 1 1
=
− [4 5] = 9
3
1 1 9 3
−
9
14
( p ( 2 ) ) t = [2 1] 9
3
−
9
(q
14
) = [1 0] 9
22
9
( 2) t
22
9 = 25
6
− 9
9
3
−
9 = 14
6
− 9
9
X (3) = −3[2 1] + [7 4] = [1 1]
Y (3) = 2[1 0] + [− 2 1] = [0 1]
22
9
6
−
9
38
9
3
−
9
p ( 2)
q
( 2)
25
=9
38
9
14
= 9
3
−
9
61
Cálculo de L••2 e U 2••
L•32
14
113
9
= 3 + [1 3] 9 =
3
34
− 34
9
•
U 23
25
9
15
= −2 + [0 1] 9 = −
38
34
17
9
L••2
0
= 1
113
34
15
U 2•• = 0 1 −
17
Cálculo de C (3) , p (3) , q (3) e d 33
C ( 3)
14
= 9
3
−
9
22
14
9
9 − 9 25
3
6
− 34 − 9
9
9
14
( p ) = [1 1] 34
3
−
34
24
34 = 11
10
− 34
34
( 3) t
14
(q (3) ) t = [0 1] 34
24
34
3
34 = 24
10 34
−
34
−
14
38 34
=
9 − 3
34
24
34
10
−
34
11
14 (3) 34
,p =
14
34
34
24
10 (3) 34
− ,q =
10
34
−
34
11
58
d 33 = 3 + [0 1] 34 =
10
− 17
34
•
Tem-se: L
1
0 0 9 0
D • U • = − 13 1 0 0 34
9
9
28 113 1 0 0
34
9
Pode-se verificar que A + X C Y t = L• D • U •
1
0
0 0
58 0
7
22
9
1
0
10
3
15
−
17
1
−
62
Nste capítulo, mostra-se que John M. Bennett, substituindo ∆A por X C Y t na
matriz ( A + ∆A) −1 , chegou à identidade ( A + X C Y t ) −1 = A −1 X [C −1 + Y t A −1 X ] −1 Y −1 A −1 .
Mostra-se também que não é necessário armazenar a inversa da matriz A a cada
modificação, é suficiente armazenar a inversa de A uma única vez, o que é fundamental
para a computação. Com o desenvolvimento de um algoritmo, chegou-se à
decomposição da matriz A + X C Y t em L• D •U • a partir de L D U da matriz A. É
possível evitar a inversa da matriz modificada A + X C Y t usando L D U da matriz A,
permitindo resolver o sistema matricial a cada modificação e, também, o método foi
exemplificado.
CAPÍTULO 5
MÉTODOS DE REFATORIZAÇÃO PARCIAL I E II
Os métodos de refatorização parcial I e II foram impulsionados pela observação de
que a alteração em uma linha/coluna de uma matriz afeta somente um subgrupo de
linhas/colunas subseqüentes. Neste capítulo considera-se também a vantagem que existe em
se ordenar a matriz A, antes de fatorá-la em L U, para minimizar os "fill-ins"[ Brandwajn,86].
5.1 Métodos de refatorização parcial I e II
Alguns melhoramentos tornam os métodos de refatorização parcial I e II mais
eficientes:
1) A matriz original é organizada sem levar em consideração as modificações futuras.
Nenhum nó é forçado para fora da organização direcionada pela esparsidade da matriz.
2) Somente os fatores alterados na matriz são atualizados, não importando a posição dos
elementos na mesma.
Se a matriz A é definida positiva e simétrica, a matriz U é transposta da matriz L, mas
não é simétrica; é idêntica em relação à estrutura esparsa.
A refatorização parcial pode ser aplicada a qualquer matriz A não singular, mas neste
capítulo considera-se que ela seja a matriz dos coeficientes de uma rede nodal.
64
5.2 Ordenação dos nós por grau mínimo
Se A é a matriz admitância de uma rede, a ordenação dos nós tem o propósito de
minimizar o aparecimento de elementos não nulos, que surgem nas matrizes L e U em
posições onde na matriz original A, existem elementos nulos: são os "fill-ins" [Tinney,85].
O processo de ordenação por grau mínimo consiste em enumerar os nós da rede em
ordem crescente de seus graus. O grau de um nó é dado pelo número de ligações que ele
possui. Isto é o mesmo que enumerar as linhas/colunas da matriz, de acordo com o número de
elementos não nulos fora da diagonal.
Após ordenar a matriz usando o método de grau mínimo, ela será organizada em
níveis em forma de L invertido. Cada um desses níveis será referido com linha/coluna
[Brandwajn,86],fig 9.
Elementos de U (linhas horizontais)
↓
Elementos de L
(linhas verticais) →
b
↑ Parte da matriz que será trabalhada
Linha/coluna que está sendo alterada
FIGURA 9 – Esquema de fatorização linha superior /coluna inferior
Ao se modificar a linha/coluna em negrito somente um subgrupo das linhas/colunas
subseqüentes será afetado.
65
Para utilização dos novos métodos, utiliza-se um esquema que permite um cálculo
eficiente do caminho de refatorização. Neste capítulo é usado o esquema linha superior/coluna
inferior.
É interessante observar que a alteração de uma linha/coluna afeta somente um
subgrupo de linhas/colunas subseqüentes. Isto pode ser visto no exemplo dado abaixo para
uma rede de doze nós (FIG. 10).
4
1
9
7
6
10
11
12
3
8
5
2
FIGURA 10 – Rede nodal com 12 nós
Pode-se formar uma tabela e um grafo onde se verifica o caminho, isto é, a seqüência
que cada nó tem na rede (Tabela abaixo).
O grafo do caminho é fornecido pela tabela como pode ser conseguido, também,
através da matriz L. O grafo tem um sentido do nó de menor numeração para o de maior
numeração. Quando o nó não é conectado diretamente ao nó seguinte, usa-se o "fillin"[Tinney,85].
66
•
•
•
•
•
•
•
•
•
•
•
8
•
•
•
9
•
•
•
10
•
•
o •
11
•
• •
•
12
•
• •
• o •
1
2
3
4
5
6
7
• elemento original
o elemento novo − fill − ins
FIGURA 11 – Estrutura esparsa da matriz L
Tabela de caminhos
Nó Pr óximo Nó
1
2
3
4
5
6
4
5
8
9
8
9
Nó Pr óximo Nó
7
8
9
10
11
12
10
11
10
12
12
0
FIGURA 12 – Tabela de caminhos referentes à rede nodal da FIG. 10
6
1
2
4
5
9
8
3
7
10
11
12
FIGURA 13 – Grafo do caminho com 12 nós (árvore)
67
Considerando, como exemplo, a alteração dos elementos na linha/coluna 7, os fatores
da matriz afetados podem ser encontrados diretamente a partir da matriz L (FIG. 11). Na
verdade, os fatores modificados estão localizados nas linhas/colunas 10, 11, 12.
As linhas/colunas afetadas podem, também, ser encontradas a partir do caminho do
grafo da matriz (FIG. 13). O caminho consiste de nós encontrados durante o traçado do
caminho, partindo de um nó até o último nó em ordem crescente.
Quando mais de uma linha/coluna é modificada, as linhas/colunas afetadas são
originadas pela união dos caminhos. Por exemplo: sejam as alterações nas linhas/colunas 6 e
1. Os fatores modificados estarão localizados nas linhas/colunas 4, 9, 10 , 11, 12. As outras
linhas/colunas permanecerão inalteradas.
5.3 Método de refatorização parcial I
O método de refatorização parcial I trabalha todas as linhas/colunas no caminho
composto:
Primeiro passo: identificar todas as linhas/colunas no caminho composto;
Segundo passo: encontrar o caminho composto para essas linhas/colunas, podendo assim
saber quais serão as linhas/colunas subseqüentes que serão alteradas;
Terceiro passo: trabalhar apenas com a submatriz onde se encontram as linhas/colunas que
serão modificadas.
Pode se perceber a vantagem que se consegue neste método em relação ao método de
compensação. Faz-se antes uma ordenação da matriz antes de fatorá-la, para se minimizar ao
máximo os "fill-ins", ao mesmo tempo que se identifica a submatriz que será formada pelas
linhas/colunas que serão modificadas.
68
5.4 Método de refatorização parcial II
O método de refatorização parcial II é uma variante do método de Bennett, no qual o
caminho de fatorização é utilizado para identificar as linhas/colunas que devem ser
atualizadas. É importante também ordenar a matriz A para minimizar os "fill-ins" e, para que
as matrizes L e U tenham o maior número de elementos nulos possíveis.
Ao se atualizarem os fatores L D U da matriz modificada, trabalha-se apenas as
linhas/colunas que serão atualizadas devido à modificação de uma ou mais linha/coluna.
5.5 Exemplo: (método de refatorização parcial I)
Considera-se uma rede de energia elétrica com 11 nós, (FIG. 14), onde a ordenação,
não foi efetuada.
5
4
7
6
3
2
1
8
10
9
11
FIGURA 14 – Rede de energia elétrica com 11 nós onde não foi feita uma ordenação
prévia dos mesmos
69
•
•
•
L = •
•
•
•
• •
o •
• •
• •
•
o •
•
• •
o •
•
•
•
• • •
• elemento original
o elemento novo ( fill − ins )
FIGURA 15 – Estrutura da matriz L conseguida ao se fatorar a matriz A da rede nodal
acima, apresentando 3 fill-ins
1
2
10
3
4
5
6
7
8
9
11
FIGURA 16 – Grafo do caminho sem ordenação (árvore de 11 nós)
70
Com a ordenação de grau mínimo tem-se a ordenação mostrada na FIG. 16. Pode-se
verificar que o número de fill-ins diminuiu de três para um.
1
4
7
3
11
9
8
10
6
5
2
FIGURA 17 – Rede de energia elétrica com 11 nós onde foi feita
a ordenação de grau mínimo
•
•
•
•
•
•
•
•
•
•
•
•
• • •
•
• •
•
•
• •
• • • o •
• elemento original
o elemento novo ( fill − in )
FIGURA 18 – Estrutura da matriz L ao se fatorar a matriz A, após a ordenação
observando o grau mínimo. Neste caso aparece apenas um fill-in
71
6
3
2
1
7
8
5
10
4
9
11
FIGURA 19 – Grafo após se fazer uma ordenação de grau mínimo (árvore com 11 nós)
Pode-se representar por um sistema a rede de energia elétrica considerada acima após
a ordenação. Atribuindo-se valores às admitâncias e às correntes, considerando como
incógnitas as tensões, tem-se:
− x4
4 x1
− 2 x7
− x5
4 x2
− 2 x6
− 2 x10
6 x3
− x1
3 x4
− x2
− 2 x2
− x10
− x10
4 x5
6 x6
− 2 x1
4 x7
4 x8
− x7
− 2 x3
− x3
− x5
− x4
− x11
− x11
− x6
− x7
E, portanto, tem-se a equação matricial A x = b onde:
− x8
− x8
− x8
− x9
− x9
− x10
6 x9
6 x10
− 2 x9
− x11
− x11
=1
=1
=3
=1
=2
=3
=0
=1
− 2 x11 = 2
=1
6 x11 = 0
72
0
0 −1 0
0 −2 0
0
0
0
4
1
0
1
4
0
0 −1 − 2 0
0
0
0
0
0
3
0
6
0
0
0
0
0
0 − 2 − 1
0
3 0
0
0
0
0
0 − 1
−1 0
1
0 −1 0
2
0
4
0
0
0
0 −1 0
A= 0 −2 0
0
0
6
0
0
0 − 1 0 e b = 3
− 2 0
0
0
0
0
0
4
0 − 1 0 − 1
0
0
0
0
0
0
4 − 1 − 1 − 1
0
1
0
0
0
0
0 −1 −1 6
0 − 2
0
2
0
1
0 − 2 0 −1 −1 0 −1 0
6
0
0 −1 −1 0
0 −1 −1 − 2 0
6
0
0
Para exemplificar suponha-se que foi feita uma alteração na linha/coluna de número 7.
Observando as FIG. 17 ou 18 pode-se ver que esta alteração irá modificar apenas as
linhas/colunas 9 e 11. Isto permite trabalhar apenas na submatriz que envolve estas
linhas/colunas. Decompõe-se a matriz A em L U por blocos.
0
4
0
4
0
0
−1 0
A11 =
0 −1
0 −2
− 2 0
0
0
0
0
0
0
A12 =
0
0
− 1
− 1
0
0
−2
0
−1
−1
0
−1
0
0
6
0
0
0
0
0
3
8
0
0 − 2 0
−1 0
0
0 − 1 − 2 0 0
0
0
0
0 0
1
3 0
0
0 0 −1 8
, A11 =
0
4
0
0 0
0
0
0
6
0 0
0
0
0
0
4 0
0
0
0
0 4
3
16
0
0
0
12
37
1
8
0
0
1
6
0
3
16
0
0
0
3
37
4
37
0
0
0
0
0
3
8
0
0
0
0
0
0
3
37
15
74
0
0
1
16
10
37
1
37
1
16
0
0
11
32
0
0
0
0
0
0
3
37
4
37
0
0
0
0
0
0
0
0
0
0
1
4
0
0
− 1
0 0 0 − 1 − 1
6 0 − 2
0 0 0
− 1
, A21 = 0 0 − 2 0 − 1 − 1 0 − 1 e A22 = 0 6 0
0
− 2 0 6
0 0 − 1 − 1 0 0 − 1 − 1
0
− 1
− 1
73
3
− 16
A21. A11−1 = 0
− 5
16
0
−
7
37
1
3
1
−
6
−
0
19
32
1
A21. A11−1. A12 =
4
21
32
−
0
1
16
0
−
0
−
7
16
−
0
11
37
−
0
17
74
11
32
0
−
0
13
32
1
−
4
1
−
4
1
−
4
1
85
−
4
32
2023
7
−
444
12
7
455
−
12
96
1
21
173
32
4
32
1
641 7
, A22 − A21. A11−1. A12 = −
444 12
4
7
121
− 85
32
12
96
−
A x = b e A = L U em blocos então L U x = b. Fazendo U x = y, tem-se L y = b.
Cálculo de y:
3
−
16
0
5
−
16
0
−
7
37
0
0
1
3
1
−
6
I8
1
−
16
−
0
−
0
−
7
16
11
37
0
−
0
17
74
0
11
−
32
0
−
13
32
1
−
4
1
−
4
1
−
4
0
I3
y1 1
y 1
2
y3 3
y 4 1
y5 2
y6 = 3
y 0
7
y8 1
y9 2
y10 1
y11 0
Tem-se:
y1 = 1, y 2 = 1, y3 = 3, y 4 = 1, y5 = 2, y6 = 3, y7 = 0, y8 = 1, y9 =
39
551
3
, y10 =
e y11 =
16
148
2
Cálculo de x:
A11
U x = y, então tem-se :
0
173
32
1
−
4
85
−
32
A12
1
−
4
2023
444
7
−
12
85
−
32
7
−
12
455
96
X 1 Y1
X = Y
2 2
74
1
1
39
3
16
551
1
onde Y1 = e Y2 =
2
148
3
3
2
0
1
Como
A11. X 1 + A12 . X 2 = Y1 e B22 . X 2 = Y2 ,então :
173
32
B22 . X 2 = − 1
4
− 85
32
matriz
1
85
39
−
16
4
32
551
2023
7
, para calcular X 2 decompõe-se a
− .X 2 =
444
12
148
7
455
3
−
2
12
96
−
B22 em LoU o . Pode-se escrever então LoU o X 2 = Y2 e fazendo
U o X 2 = z , Lo z = Y2 .
39
16
551
0
0
1
o
L z = − 0,0462
1
0 z = 148 ,
3
− 0,4913 − 0,1554 1
2
o
U X2 =
2,4375
logo z = 3,8355
3.2935
5,4066 − 0,25 − 2,6563
2,4375
0
4,5447 − 0,7061 X 2 = 3,8355 ,
0
2,2935
0
3,3248
1
x10 = 0,9978 fazendo os arredondamentos tem-se X 2 = 1
1
x11 = 0,9837
x9= 0,9906
A11−1 . X 1 = Y1 − A12 . X 2
A11−1 . A11 . X 1 = A11−1 (Y1 − A12 . X 2 ) ,
X 1 = A11−1 (Y11 − A12 . X 2 )
então:
75
3
8
0
0
1
X1 = 8
0
0
3
16
0
0
0
12
37
1
8
0
0
1
6
0
3
16
0
0
0
3
37
4
37
0
0
0
0
0
3
8
0
0
0
0
0
0
1
37
1
37
0
0
1
16
10
37
15
37
1
16
0
0
11
32
0
0
0
0
0
0
3
37
4
37
0
0
0
0
0
0
0
(
0
0
0
1
4
1 0
1 0
3 0
1 − 0
2 0
3 0
0 − 1
1 − 1
0
0
1
1
0
0
1
− 2 − 1
1
0 − 1
1
.1 ) = Temos,
1
−1 0
1
−1 0
1
1
0 −1
− 1 − 1
1
portanto: x1 = x2 = x3 = ... = x9 = x10 = x11 = 1 , que é a solução do sistema dado.
Mostra-se neste capítulo que ao se alterar uma ou mais linhas/coluna de uma matriz são
afetados somente um subgrupo de linhas/colunas subseqüentes e, também, a vantagem
de ordenar uma matriz antes de fatorá-la em L U, pois proporciona menos
preenchimentos(“fill-ins”) das posições de elementos nulos na matriz A que aparecem
nas matrizes L e U, quando se efetua a decomposição.
CONCLUSÃO
O método de compensação é apropriado às aplicações onde o número de mudanças é
pequeno, as modificações não são permanentes e as equações modificadas não precisam ser
resolvidas repetidamente.
O método da refatorização simples refatora somente a submatriz que contém os
elementos da matriz modificada e neste método, se os elementos modificados estiverem no
topo da matriz, não existe vantagem em usá-lo.
O método de refatorização simples, com arranjo especial, força os elementos para a
parte inferior da matriz, porém, pode alterar seriamente a esparsidade da matriz e restringe-se
às aplicações onde a posição dos elementos da matriz modificada pode ser prevista
antecipadamente.
O método de Bennett atualiza os fatores L D U da matriz modifica e a vantagem é que
não precisa armazenar a inversa da matriz modifica a cada modificação feita, basta apenas
armazenar a inversa uma única vez. Pode-se evitar a inversa completa da matriz modificada e
chega-se aos fatores L D U da matriz modificada através dos fatores L D U da matriz inicial.
Os métodos de refatorização parcial I e II são significativamente mais rápidos que os
métodos anteriores. A eficiência deriva dos conceitos de caminho recentemente desenvolvido.
O método de refatorização parcial I refatora seletivamente um subgrupo de linhas/colunas da
matriz modificada, isto é, somente o subgrupo que sofre modificações ao se alterar uma ou
mais linha/colunas. O método de refatorização parcial II é o de Bennett modificado que utiliza
o caminho para alterar seletivamente fatores matriciais, portanto, os dois últimos métodos são
mais usados.
REFERÊNCIAS BIBLIOGRÁFICA
[Alsaç, 83] O.Alsaç ,B. Stott &W.F. Tinney . “ Sparsity-0riented Compensation Mehods For
Modified Network Solutions”. IEEE Transactions on Power Apparatus and Systems ,vol.Pas102 , N.5, May 1983.
[Bennett, 65] Bennett, John M. "Triangular Factors of Modified Matrices", Numerische
Mathematik. 1965-pp.217-221.
[Betancourt, 86] Betancourt, Ramon and Alvarado, Fernando L. ,“Parallel inversion of sparse
matrices” – IEEE Transactions on Power Systems, vol. PWRS – 1 – Feb-1986 – pp.74-81.
[Betancourt, 88] Betancourt, Ramon, “An efficient heuristic ordering algorithm for partial
matrix refactorization” , IEEE Transactions on Power systems, vol. , August, 1988.
[Brandwajn, 86] Brandwajn,Vladimir & Chan, M. Sherman ."Partial Matrix Refactorization"IEEE Transactions on Power Systems,vol.pWRS-1, February 1986- pp.193-199.
[Dierker, 86] Dierker, Paul F. and Voxman, William L., “Discrete Mathematic”, Harcourt
Brace Jovanovich, International Edition, 1986.
[Duff ,86] Duff, J. F., Erisman, A. M. and Reid, J. K., “Direct Methods for Sparse Matrices”,
Oxford University Press, Oxford, 1986.
[George, 81] George, A. and Liw, J. W. H., “Computer Solution of Large Sparse Positive
Definite Systems”, Prendice-Hall, Englewood Clifs, N. J.,1981.
[Golumbic, 80] Golumbic, Martin Chales, “Algorithmic Graph Teory and Perfect Graphs”,
Academic Press, 1980.
[Gustavson, 72] Gustavson, F. G., “Some Basic Techiniques for Solving Sparse Systems of
Equations in Sparse Matrix and their Applications”, Rose, D. J. and Willowghby, R. A., eds.,
Plenum Press, New York, 1972, pp. 41-52.
[Leon, 98] Leon, Steven J., “Linear Álgebra with Applications”, Prentice-Hall, Englewood
Clifs, N. J., 1998.
[Monticelli, 83] Monticelli,Alcir José. “Fluxo de Carga em Rede de Energia Elétrica”.Ed.
Edgard Blücher Ltda –São Paulo. 1983.
[Monticelli,99] A. B. Alves,E.N. Asada and A. Monticelli. “Critical Evaluation of Direct and
Iterave Methods for Solving Ax = b Systems im Power Flow Calculations and Contigency
Analysis”. IEEE Transactions on Power Systems,Vol.14.No.2, May 1999.
[Morozowski, 81] Morozowski, Marciano Filho, “Matrizes Esparsas em Redes de Potência”,
Técnicas de Operação- Livros Técnicos e Científicos Editora S/A, 1981.
[Sato, 63] Sato, N. & Tinney, W. F., “Techniques for exploiting the sparsity of the network
admit New York,Dec.1963.
78
[Szwarcfiter, 83] Szwarcfiter, Jayme Luiz, “Grafos e algoritmos computacionais”, Editora
Campos Ltda., Segunda Edição, 1983.
[Tinney,67] Tinney, W. F. and Walker, J. W., “Direct Solutions of Sparse Network Equation
by Optimally Ordered Triangular Factorization”, Proceedings of the Institute of Electrical and
Electronics Engineers, New York, 55 (11):18.
[Tinney ,85] Tinney, W. F., Brandwajn, V. and Chan, S. M., “Sparse Vector Methods”, IEEE
Transactions on Power Apparatus and Systems, vol.: PAS-104, n. 2, Feb. 1985.
79
APÊNDICE I
a
0
1
Processo detalhado para mostrar que L1 11
( 2) U = A
0
A
a 22
a
H = 32
...
an 2
... a 2 n
... a3n
... ...
... a nn
a 23
a33
...
an3
a 21
a
G = 31 , F = [a12
...
a n1
a13 ... a1n ]
a 21a12
a 22 − a
11
a31a12
GF
a −
= 32
H−
a11
a11
...
a a
a n 2 − n1 12
a11
a
L1 11
0
a 21a13
a11
a31a13
a33 −
a11
...
a a
a n 3 − n1 13
a11
a 23 −
a 21a1n
a11
a a
... a3n − 31 1n
( 2)
a11 = A
...
...
a n1a1n
... a nn =
a11
... a 2 n −
0
=
A (2 )
1 0 0
l
21 1 0
= l31 0 1
... ... ...
l n1 0 0
...
...
...
...
...
0
0
0
...
1
a11
l 21a11
= l31a11
...
l a
n1 11
a11
0
0
0
0
0
a a
a 22 − 21 12
a11
a a
a32 − 31 12
a11
...
a a
a n 2 − n1 12
a11
0
a a
a 22 − 21 12
a11
a a
a32 − 31 12
a11
...
a a
a n 2 − n1 12
a11
0
a a
a 23 − 21 13
a11
a a
a33 − 31 13
a11
...
a a
a n 3 − n1 13
a11
0
a a
a 23 − 21 13
a11
a a
a33 − 31 13
a11
...
a a
a n 3 − n1 13
a11
0
a a
... a 2 n − 21 1n
a11
a a
... a3n − 31 1n =
a11
...
...
a a
... a nn − n1 1n
a11
...
0
a 21a1n
... a 2 n −
a11
a a
... a3n − 31 1n
a11
...
...
a a
... a nn − n1 1n
a11
...
81
a
L1 11
0
0 1
1 a11
( 2) U = L
A
0
u13 ... u1n
0 ... 0
1 ... 0 = A
... ... ...
0 0 1
1 u12
0 1
0
0 0
A( 2)
... ...
0 0
a11u12
a11
a21a12
l21a11 l21a11u12 + a22 −
a11
a31a12
= l31a11 l31a11u12 − a32 −
a11
...
...
ln1a11 ln1a11u12 − a3n + an1a12
a11
...
a11u13
a a
l21a11u13 + a23 − 21 13
a11
a31a13
l31a11u13 + a33 −
a11
...
a a
ln1a11u13 + an 3 − n1 13
a11
...
...
...
...
Pode-se verificar que:
l21 a11 = a21 ⇒ l21 =
a21
a11
a11u12 = a12 ⇒ u12 =
a12
a11
l31a11 = a31 ⇒ l31 =
a31
a11
a11u13 = a13 ⇒ u13 =
a13
a11
.
. . . . . . . .
a
ln1a11 = an1 ⇒ ln1 = n1
a11
. . . .
. . . . .
a
a11u1n = a1n ⇒ u1n = 1n
a11
entâo, tem − se :
a a
a a
a a
l21a11u12 = 21 12 , l21a11u13 = 21 13 , ... , ln1a11u1n = n1 1n ,
a11
a11
a11
conclue − se :
ai 1
a1 j
li 1 =
, u1 j =
a11
a11
e
li 1a11u1 j =
ai 1a1 j
a11
a21a1n
l21a11u1n + a2 n −
a11
a a
l31a11u1n + a3n − 31 1n
a11
...
a a
ln1a11u1n + ann − n1 1n
a11
a11u1n
82
APÊNDICE II
Implementação computacional do Método de Bennett
Linguagem Delphi, Borland do Brasil Ltda.
unit UFrmPrincipal;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, Buttons, ExtCtrls, RXDice, TimerLst, MPlayer, Menus,
XP_MainMenu;
const size = 15;
type Tmatrix = array[1..size,1..size] of double;
type lc = array[1..2] of integer;
type
TCounter = class(TThread)
private
hour,min,sec,msec: word;
protected
procedure Execute; override; // Método anulado
procedure contar;
public
constructor Create();
end;
type
Tprinti = class(TThread)
private
protected
procedure Execute; override; // Método anulado
procedure contar;
public
constructor Create();
end;
type
TFrmLDU = class(TForm)
Panel1: TPanel;
BitBtn1: TBitBtn;
Edit1: TEdit;
83
Edit2: TEdit;
Edit3: TEdit;
Edit4: TEdit;
Edit5: TEdit;
Edit6: TEdit;
Label1: TLabel;
Label2: TLabel;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
Memo1: TMemo;
BitBtn2: TBitBtn;
Label7: TLabel;
Timer1: TTimer;
ttfXPMainMenu1: ttfXPMainMenu;
Arquivo1: TMenuItem;
Importao1: TMenuItem;
N1: TMenuItem;
Sair1: TMenuItem;
OpenDialog1: TOpenDialog;
procedure BitBtn1Click(Sender: TObject);
procedure BitBtn2Click(Sender: TObject);
procedure Sair1Click(Sender: TObject);
procedure Importao1Click(Sender: TObject);
private
parathread: boolean;
function MAdicao(lina, cola, linb, colb: integer; ma, mb: TMatrix;
var linr, colr: integer; param: integer): TMatrix;
function MMult(lina, cola, linb, colb: integer; ma, mb: TMatrix;
var linr, colr: integer): TMatrix;
function MMultValor(lin, col: integer; m: TMatrix; var linr,
colr: integer; valor: double; param: integer): TMatrix;
function MTransposta(lin, col: integer; var linr, colr: integer;
m: TMatrix): TMatrix;
Function LeDados(matrix: string; dados: string; var linm, colm: integer; var m: tmatrix):
boolean;
function InsereEspaco(valor: integer): string;
{ Private declarations }
public
{ Public declarations }
end;
var
FrmLDU: TFrmLDU;
implementation
{$R *.dfm}
function TFrmLDU.InsereEspaco(valor: integer): string;
84
var i: integer;
s: string;
begin
s:= '';
for i:= 1 to valor do
s:= s + ' ';
result:= s;
end;
Function TFrmLDU.MMult(lina,cola,linb,colb: integer; ma, mb: TMatrix; var linr, colr:
integer): TMatrix;
var mr: TMatrix;
u, t, v: integer;
x: double;
begin
if cola <> linb then
exit;
linr:= lina;
colr:= colb;
for u:= 1 to lina do
for t:= 1 to colb do
begin
x:= 0;
for v:= 1 to cola do
x:= x + (ma[u,v] * mb[v,t]);
mr[u,t]:= x;
end;
result:= mr;
end;
function TFrmLDU.MTransposta(lin, col: integer; var linr,colr: integer; m: TMatrix):
TMatrix;
var mr: TMatrix;
i, j: integer;
begin
linr:= col;
colr:= lin;
for i:= 1 to linr do
for j:= 1 to colr do
mr[i,j]:= m[j,i];
result:= mr;
end;
Function TFrmLDU.MAdicao(lina, cola, linb, colb: integer; ma, mb: TMatrix; var linr, colr:
integer; param: integer): TMatrix;
var mr: TMatrix;
i, j: integer;
begin
if (lina <> linb) or (cola <> colb) then
exit;
85
linr:= lina;
colr:= colb;
for i:= 1 to lina do
for j:= 1 to cola do
if param = 1 then
mr[i,j]:= ma[i,j] + mb[i,j]
else
mr[i,j]:= ma[i,j] - mb[i,j];
result:= mr;
end;
Function TFrmLDU.MMultValor(lin, col: integer; m: TMatrix; var linr, colr: integer; valor:
double; param: integer): TMatrix;
var mr: TMatrix;
i, j: integer;
begin
linr:= lin;
colr:= col;
for i:= 1 to lin do
for j:= 1 to col do
if param = 1 then
mr[i,j]:= valor * m[i,j]
else
if param = 2 then
mr[i,j]:= m[i,j] / valor
else
if m[i,j] <> 0 then
mr[i,j]:= valor / m[i,j]
else
mr[i,j]:= 0;
result:= mr;
end;
Function TFrmLDU.LeDados(matrix: string; dados: string; var linm, colm: integer; var m:
tmatrix): boolean;
var i, x, j, verro: integer;
neg: boolean;
res: string;
begin
result:= false;
neg:= false;
x:= 1; j:= 1;
res:= '';
for i:= 1 to length(dados) do
begin
if dados[i] in ['0'..'9',',','.'] then
begin
res:= res + dados[i];
end
86
else
if (dados[i] = ';') then
begin
if neg then
begin
m[x,j]:= strtofloat('-' + res);
neg:= false;
end
else
m[x,j]:= strtofloat(res);
res:= '';
inc(j);
inc(x);
if x = 2 then
verro:= j
else
if j <> verro then
begin
result:= true;
showmessage('Matrix ' + matrix + ' inválida.');
exit;
end;
j:= 1;
end
else
if dados[i] = '-' then
neg:= true
else
if dados[i] = ' ' then
begin
if res <> '' then
begin
if neg then
begin
m[x,j]:= strtofloat('-' + res);
neg:= false;
end
else
m[x,j]:= strtofloat(res);
res:= '';
inc(j);
end;
end;
end;
if res <> '' then
begin
if neg then
begin
m[x,j]:= strtofloat('-' + res);
neg:= false;
end
87
else
m[x,j]:= strtofloat(res);
res:= '';
end;
linm:= x;
colm:= j;
end;
procedure TFrmLDU.BitBtn1Click(Sender: TObject);
const dist: integer = 7;
var l, d, u, l1, d1, u1, raux, raux1: tmatrix;
c, x, y, p, q, pt, qt: array[1..size] of Tmatrix;
lcl, lcd, lcu, lcl1, lcd1, lcu1: lc;
lcc, lcx, lcy, lcp, lcq, lcpt, lcqt: array[1..size] of lc;
i, j, a, f, g, laux, caux, laux1, caux1: integer;
resp, h: string;
begin
parathread:= false;
tcounter.Create;
if ledados('L',edit1.Text,lcl[1],lcl[2],l) then
exit;
if ledados('D',edit2.Text,lcd[1],lcd[2],d) then
exit;
if ledados('U',edit3.Text,lcu[1],lcu[2],u) then
exit;
if ledados('X',edit4.Text,lcx[1][1],lcx[1][2],x[1]) then
exit;
if ledados('C',edit5.Text,lcc[1][1],lcc[1][2],c[1]) then
exit;
if ledados('Y',edit6.Text,lcy[1][1],lcy[1][2],y[1]) then
exit;
//Cálculo de x2 e c2
i:= 1;
memo1.Lines.Clear;
memo1.Lines.Add('RESULTADO');
repeat
memo1.Lines.Add('------------------');
memo1.Lines.Add('Para I = ' + inttostr(i));
memo1.Lines.Add('------------------');
for j:= 1 to lcx[i][1] - 1 do
begin
raux:= MMultvalor(1,lcx[i][2],x[i],laux,caux,-l[i+j][i],1);
for a:= 1 to lcx[i][2] do
x[i+1][j][a]:= raux[1][a] + x[i][j+1][a];
lcx[i+1][1]:= j;
lcx[i+1][2]:= a - 1;
raux:= MMultvalor(1,lcy[i][2],y[i],laux,caux,-u[i][i+j],1);
for a:= 1 to lcy[i][2] do
y[i+1][j][a]:= raux[1][a] + y[i][j+1][a];
lcy[i+1][1]:= j;
lcy[i+1][2]:= a - 1;
88
end;
memo1.Lines.Add('X' + inttostr(i) + ' = ');
resp:= ' ';
for f:= 1 to lcx[i+1][1] do
begin
for g:= 1 to lcx[i+1][2] do
begin
resp:= resp + h + formatfloat('#,##0.##',(x[i+1][f][g]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(x[i+1][f][g]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
memo1.Lines.Add('Y' + inttostr(i) + ' = ');
resp:= ' ';
for f:= 1 to lcy[i+1][1] do
begin
for g:= 1 to lcy[i+1][2] do
begin
resp:= resp + h + formatfloat('#,##0.##',(y[i+1][f][g]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(y[i+1][f][g]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
// cálculo de l**1 e u*1*
for f:= 1 to lcl[1] do
for g:= 1 to lcl[2] do
if f = g then
begin
l1[f,g]:= 1;
u1[f,g]:= 1;
end;
pt[i]:= MMult(lcx[i][1],lcx[i][2],lcc[i][1],lcc[i][2],x[i],c[i],lcpt[i][1],lcpt[i][2]);
lcpt[i][1]:= 1;
p[i]:= MTransposta(lcpt[i][1],lcpt[i][2],lcp[i][1],lcp[i][2],pt[i]);
raux:= MMult(1,lcy[i][2],lcp[i][1],lcp[i][2],y[i],p[i],laux,caux);
d1[i][i]:= d[i][i] + raux[laux][caux];
raux:= MTransposta(lcc[i][1],lcc[i][2],laux,caux,c[i]);
qt[i]:= MMult(1,lcy[i][2],laux,caux,y[i],raux,lcqt[i][1],lcqt[i][2]);
q[i]:= MTransposta(lcqt[i][1],lcqt[i][2],lcq[i][1],lcq[i][2],qt[i]);
lcl1[1]:= lcl[1]; lcl1[2]:= lcl[2];
89
lcu1[1]:= lcu[1]; lcu1[2]:= lcu[2];
lcd1[1]:= lcd[1]; lcd1[2]:= lcd[2];
for j:= 1 to lcx[i][1] - 1 do
begin
raux:= MMult(j,lcx[i+1][2],lcq[i][1],lcq[i][2],x[i+1],q[i],laux,caux);
l1[i+j][i]:= l[i+j][i] + ((1 / d1[i][i]) * raux[laux][caux]);
raux:= MMult(j,lcy[i+1][2],lcp[i][1],lcp[i][2],y[i+1],p[i],laux,caux);
u1[i][i+j]:= u[i][i+j] + ((1 / d1[i][i]) * raux[laux][caux]);
end;
j:= 1;
memo1.Lines.Add('D*(' + inttostr(i) + ',' + inttostr(i) + ') = ');
resp:= ' ';
for f:= 1 to i do
begin
for g:= 1 to i do
begin
resp:= resp + h + formatfloat('#,##0.##',(d1[f][g]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(d1[f][g]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
memo1.Lines.Add('L*(*,' + inttostr(i) + ') = ');
resp:= ' ';
for f:= 1 to lcl1[2] do
begin
for g:= 1 to i do
begin
resp:= resp + h + formatfloat('#,##0.##',(l1[f][g]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(l1[f][g]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
memo1.Lines.Add('U*(*,' + inttostr(i) + ') = ');
resp:= ' ';
for f:= 1 to i do
begin
for g:= 1 to lcu1[2] do
begin
resp:= resp + h + formatfloat('#,##0.##',(u1[f][g]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(u1[f][g]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
90
end;
memo1.Lines.Add('');
//Cálculo de c i + j
raux:= MMult(lcq[i][1],lcq[i][2],lcpt[i][1],lcpt[i][2],q[i],pt[i],laux,caux);
raux1:= MMultvalor(laux,caux,raux,laux1,caux1,(1 / d1[i][i]),1);
c[i+j]:= MAdicao(lcc[i][1],lcc[i][2],laux1,caux1,c[i],raux1,lcc[i+j][1],lcc[i+j][2],2);
memo1.Lines.Add('C' + inttostr(i) + ' = ');
resp:= ' ';
for f:= 1 to lcc[i+j][1] do
begin
for g:= 1 to lcc[i+j][2] do
begin
resp:= resp + h + formatfloat('#,##0.##',(c[i+j][f][g]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(c[i+j][f][g]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
inc(i);
until i = lcx[1][1];
pt[i]:= MMult(lcx[i][1],lcx[i][2],lcc[i][1],lcc[i][2],x[i],c[i],lcpt[i][1],lcpt[i][2]);
lcpt[i][1]:= 1;
p[i]:= MTransposta(lcpt[i][1],lcpt[i][2],lcp[i][1],lcp[i][2],pt[i]);
raux:= MMult(1,lcy[i][2],lcp[i][1],lcp[i][2],y[i],p[i],laux,caux);
d1[i][i]:= d[i][i] + raux[laux][caux];
// close;
for i:= 1 to lcl1[1] do
for j:= 1 to lcl1[2] do
begin
if i = j then
begin
l1[i][j]:= 1;
u1[i][j]:= 1;
end;
if i > j then
begin
u1[i][j]:= 0;
d1[i][j]:= 0;
end;
if i < j then
begin
l1[i][j]:= 0;
d1[i][j]:= 0;
end;
end;
91
memo1.Lines.Add('-------------------------------');
memo1.Lines.Add('RESULTADO');
memo1.Lines.Add('L* = ');
resp:= ' ';
for i:= 1 to lcl1[1] do
begin
for j:= 1 to lcl1[1] do
begin
resp:= resp + h + formatfloat('#,##0.##',(l1[i][j]));
h:= insereespaco(10 - length(formatfloat('#,##0.##',(l1[i][j]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
memo1.Lines.Add('D* = ');
resp:= ' '; h:= '';
for i:= 1 to lcl1[1] do
begin
for j:= 1 to lcl1[1] do
begin
resp:= resp + h + formatfloat('#,##0.##',d1[i][j]);
h:= insereespaco(10 - length(formatfloat('#,##0.##',(d1[i][j]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
memo1.Lines.Add('');
memo1.Lines.Add('U* = ');
resp:= ' '; h:= '';
for i:= 1 to lcl1[1] do
begin
for j:= 1 to lcl1[1] do
begin
resp:= resp + h + formatfloat('#,##0.##',u1[i][j]);
h:= insereespaco(10 - length(formatfloat('#,##0.##',(u1[i][j]))));
end;
memo1.Lines.Add(resp);
resp:= ' ';
h:= '';
end;
// timer1.Enabled:= false;
parathread:= true;
end;
procedure TFrmLDU.BitBtn2Click(Sender: TObject);
begin
TPRINTI.Create;
92
end;
procedure TCounter.Execute;
begin
contar;
end;
procedure TPRINTI.Execute;
begin
contar;
end;
procedure TCounter.contar;
begin
msec:= 0; min:= 0; sec:= 0; hour:= 0;
while not frmldu.parathread do
begin
inc(msec);
if msec = 1000 then
begin
msec:= 0;
inc(sec);
end;
if sec = 60 then
begin
sec:= 0;
inc(min);
end;
if min = 60 then
begin
min:= 0;
inc(hour);
end;
FrmLDU.Label7.Caption:=
formatdatetime('hh:nn:ss:zzz',encodetime(hour,min,sec,msec));
FrmLDU.Label7.Refresh;
end;
end;
procedure tprinti.contar;
begin
TRY
frmldu.Memo1.Lines.SaveToFile('LPT1');
EXCEPT
END;
end;
constructor TCounter.Create();
begin
inherited Create(True); { Chama o contrutor herdado. Ele irá temporariamente colocar o
thread em estado de espera para depois executá-lo. }
93
FreeOnTerminate := True; // Libera o objeto após terminar.
Priority := TpLower; { Configura sua prioridade na lista de processos do Sistema
operacional. }
Resume; // Inicia o Thread.
end;
constructor tprinti.Create();
begin
inherited Create(True); { Chama o contrutor herdado. Ele irá temporariamente colocar o
thread em estado de espera para depois executá-lo. }
FreeOnTerminate := True; // Libera o objeto após terminar.
Priority := TpLower; { Configura sua prioridade na lista de processos do Sistema
operacional. }
Resume; // Inicia o Thread.
end;
procedure TFrmLDU.Sair1Click(Sender: TObject);
begin
self.Close;
end;
procedure TFrmLDU.Importao1Click(Sender: TObject);
var f: textfile;
s: string;
begin
if opendialog1.Execute then
if opendialog1.FileName <> '' then
begin
assignfile(f,opendialog1.FileName);
reset(f);
Readln(F, S);
Edit1.Text := S;
Readln(F, S);
edit2.Text:= s;
Readln(F, S);
edit3.Text:= s;
Readln(F, S);
edit4.Text:= s;
Readln(F, S);
edit5.Text:= s;
Readln(F, S);
edit6.Text:= s;
end;
end;
end.
94
EXEMPLO:
0
5 −2 −3 0
− 2 5
0 − 3 0
Seja a matriz A = − 3 0
5 − 1 − 1 , A = L D U
0
0 − 3 −1 5
0
0 −1 0
2
1
2
− 5
3
−
L= 5
0
0
0
0
1
0
2
7
5
−
7
−
1
13
20
7
−
20
−
0
2
1 − 5
0 1
U =
0 0
0 0
0 0
0
1
− 1 0
X =0 1
0 0
0 − 1
0
5
0
0
0
0
0
0
D=
0
1
0
13
−
−1
0
33
0
3
5
2
−
7
−
1
0
0
0
21
5
0
0
0
0
0
20
7
0
0
0
33
20
0
0
0
0
5
−
0
7
13
7
−
−
7
20
− 13
1
33
0
1
0
1 0
D=
0 1
0
1
− 1 0
Y =0 1
0 0
0 − 1
0
0
0
0
46
33
95
Aplicação do programa no exemplo da página anterior:
RESULTADO
-----------------Para I = 1
-----------------X1 =
-0,5
0
0,6
1
0
0
0
-1
Y1 =
-0,5
0,6
0
0
0
1
0
-1
D*(1,1) =
6
L*(*,1) =
1
-0,58
-0,5
0
0
U*(*,1) =
1
C1 =
0,83
0
-0,58
0
1
-----------------Para I = 2
-----------------X2 =
0,45
1
-0,35
0
0
-1
-0,5
0
0
96
Y2 =
0,45
-0,35
0
1
0
-1
D*(2,2) =
6
0
0
4,41
L*(*,2) =
1
-0,58
-0,5
0
0
0
1
-0,33
-0,68
0
U*(*,2) =
1
0
-0,58
1
-0,5
-0,33
0
-0,68
0
0
0
-0,68
-0,49
0
0
-0,53
C2 =
0,79
0
0
1
-----------------Para I = 3
-----------------X3 =
-0,06
0,65
0,16
-0,65
Y3 =
-0,06
0,18
0,65
-0,61
D*(3,3) =
6
0
0
0
4,41
0
0
0
4,02
L*(*,3) =
1
-0,58
-0,5
0
0
0
1
-0,33
-0,68
0
0
0
1
-0,49
-0,5
U*(*,3) =
1
0
0
-0,58
1
0
-0,5
-0,33
1
97
C3 =
0,76
-0,09
-0,09
0,75
-----------------Para I = 4
-----------------X4 =
0,14
-0,4
Y4 =
0,15
-0,36
D*(4,4) =
6
0
0
0
0
4,41
0
0
0
0
4,02
0
0
0
0
1,98
L*(*,4) =
1
-0,58
-0,5
0
0
0
1
-0,33
-0,68
0
0
0
1
-0,49
-0,5
0
0
0
1
-0,5
U*(*,4) =
1
0
0
0
-0,58
1
0
0
-0,5
-0,33
1
0
0
-0,68
-0,49
1
0
0
-0,53
-0,49
0
0
0
1
-0,5
0
0
0
0
1
0
0
0
1,98
0
0
0
0
0
1,5
C4 =
0,76
-0,06
-0,06
0,63
------------------------------RESULTADO
L* =
1
0
0
-0,58
1
0
-0,5
-0,33
1
0
-0,68
-0,49
0
0
-0,5
D* =
6
0
0
0
0
0
4,41
0
0
0
0
0
4,02
0
0
98
U* =
1
0
0
0
0
-0,58
1
0
0
0
-0,5
-0,33
1
0
0
0
-0,68
-0,49
1
0
0
0
-0,53
-0,49
1
Download
