Upload
tranngoc
View
219
Download
3
Embed Size (px)
Citation preview
Cálculo Automático de Estruturas
MÉTODOS NUMÉRICOS
J P Moitinho de Almeida
E M B Ribeiro Pereira
2006
Nota introdutória
Estes apontamentos foram editados pela primeira vez em Outubro de 1986 para o Curso de Cálculo Automático de Estruturas, organizado pelo Departamento de Engenharia Civil do IST, com o apoio do FSE.
Em 1994 foram reeditados e revistos, com a preciosa colaboração do Engº Luís Castro.
Esta edição de 2006, com algumas gralhas corrigidas, vai ser colocada 'online'.
Índice
Capítulo 1 Técnicas de Armazenamento de Matrizes...........................................................................1
1.1 Matrizes Não Simétricas.................................................................................................11.1.1 Dimensionamento Directo............................................................................................11.1.2 Matrizes Bandeadas......................................................................................................21.1.3. Matrizes Esparsas........................................................................................................3
1.2 Matrizes Simétricas.........................................................................................................41.2.1 Matrizes Não Esparsas.................................................................................................41.2.2. Matrizes Bandeadas....................................................................................................51.2.3 Matrizes Esparsas.........................................................................................................61.2.4 Armazenamento em Perfil...........................................................................................6
1.3 Operação de Matrizes em Computador..........................................................................71.3.1 Somar duas matrizes simétricas cheias : Cij = Aij + Bij............................................81.3.2 Multiplicar uma matriz simétrica bandeada por um vector ci = Aij bj....................91.3.3 Multiplicar uma matriz simétrica em perfil por um vector ci=Aijbj.......................11
Capítulo 2 Métodos de Resolução de Sistemas de Equações..............................................................14
2.1 Métodos de Gauss e GaussJordan...............................................................................14
2.2 Métodos de Factorização de Matrizes.........................................................................20
2.3 Métodos Iterativos de Resolução de Sistemas.............................................................242.3.1 Critérios de Paragem para Métodos Iterativos.........................................................242.3.2 Métodos de Relaxação...............................................................................................25Método de Jacobi................................................................................................................26Método de GaussSeidel.....................................................................................................27Método SOR (Successive OverRelaxation)......................................................................272.3.3 Métodos Projectivos..................................................................................................28Método do Gradiente..........................................................................................................29Método dos Gradientes Conjugados..................................................................................30
2.4 Adaptação dos Métodos de Resolução à Forma de Armazenamento das Matrizes...312.4.1 Matrizes Simétricas....................................................................................................312.4.2 Matrizes Simétricas Bandeadas................................................................................322.4.3 Matrizes Armazenadas em Perfil..............................................................................352.4.4 Algoritmos para os Métodos Iterativos.....................................................................372.4.5 Solução Frontal..........................................................................................................38
Capítulo 3 Integração Numérica...........................................................................................................39
3.1 Fórmulas de NewtonCottes.........................................................................................39
3.2 Fórmula de Gauss..........................................................................................................41
3.3 Resultados Numéricos..................................................................................................43
3.4 Integração em Área e em Volume................................................................................46
3.5 Implementação da Integração Numérica.....................................................................47
Capítulo 4 Valores e Vectores Próprios................................................................................................48
4.1 Métodos de Determinante.............................................................................................48
4.2 Métodos de Potências...................................................................................................50
Capítulo 1
Técnicas de Armazenamento de Matrizes
1.1 Matrizes Não Simétricas
1.1.1 Dimensionamento Directo
Todas as linguagens de programação de alto nível permitem a utilização de variáveis indexadas. Uma matriz é assim referida como uma variável com dois índices.
Dado que a memória do computador trabalha numa só dimensão, é função do compilador (ou do interpretador) transformar variáveis pluridimensionais em cadeias unidimensionais.
Consoante as linguagens, podemos encontrar diferentes formas de efectuar esta transformação.
Para matrizes (variáveis bidimensionais) é possível efectuar uma transformação por linhas ou por colunas. Exemplificando, a matriz
43332313
42322212
41312111
=
aaaa
aaaa
aaaa
A
corresponde ao vector
[ ]433323134232221241312111= aaaaaaaaaaaaA'
em queA ( i , j ) = A' ( k ) com k = ( i 1 ) × ncol + j .
1
1 Técnicas de Armazenamento de Matrizes
Numa arrumação por colunas, a matriz A corresponde ao vector
[ ]434241333231232221131211= aaaaaaaaaaaaA ''
em queA ( i , j ) = A'' ( k ) com k = ( j 1 ) × nlin + i .
Advém da norma da linguagem FORTRAN que a arrumação das matrizes seja feita por colunas. Esta norma permite que uma variável seja referida de mais do que uma forma (com um número diferente de dimensões) em módulos diferentes, ou até no mesmo módulo.
Nestas situações, e para uma ''optimização'' da memória e do tempo de cálculo, pode convir ao programador utilizar o armazenamento de uma matriz na sua forma unidimensional.
Noutras linguagens como o PASCAL e (em geral) o BASIC, não é possível tirar partido destas transformações.
1.1.2 Matrizes Bandeadas
São de especial importância no cálculo científico, pela redução de memória e de tempo de cálculo que, quando convenientemente aproveitadas, podem proporcionar, as matrizes em que se sabe à priori que existe um grande número de elementos nulos.
De entre estas destacamse as matrizes bandeadas, que se caracterizam por serem nulos todos os elementos que estão afastados mais do que um certo número de colunas da diagonal.
À faixa diagonal, de elementos predominantemente não nulos, chamase banda da matriz, e à sua dimensão chamase largura de banda da matriz. Por exemplo, na matriz
665646
65554535
6454443424
5343332313
42322212
312111
000
00
0
0
00
000
=
aaa
aaaa
aaaaa
aaaaa
aaaa
aaa
A
a largura de banda é igual a 5. A forma mais simples de armazenar estas matrizes consiste em guardálas numa outra matriz com o mesmo número de linhas e com o número de colunas igual à largura de banda.
2
1 Técnicas de Armazenamento de Matrizes
À matriz A anteriormente apresentada corresponde
665646
65554535
6454443424
5343332313
42322212
312111
aaa
aaaa
aaaaa
aaaaa
aaaa
aaa
=A'
em que
contrário caso ,
1 =)(
0
+ + < < +,',
lbandlsbanijlsbaniselsbanijiAjiA
, )(
indicando lsban a coluna onde é armazenada a diagonal de A na matriz A'. Se a diagonal é centrada então lsban = (lband + 1) / 2 .
Para efeitos de programação é conveniente que os termos que aparecem nos extremos desta matriz e que não têm correspondência na matriz original sejam postos a zero.
1.1.3. Matrizes Esparsas
Quando a densidade de elementos nulos numa matriz for apreciável, mas a sua disposição irregular, não é rentável utilizar o armazenamento em banda, pois este obrigaria a trabalhar com um grande número de elementos nulos.
Podese então recorrer a outras formas de armazenamento, que proporcionam uma economia sensível de tempo de cálculo e/ou de memória utilizada.
A sua rentabilidade depende do índice de esparsidade da matriz, sendo este definido como o quociente entre o número de elementos nulos e o número de elementos da matriz.
Todos estes métodos têm como base a descrição da posição dos elementos não nulos através de informação armazenada em vectores auxiliares (vectores guias).
A forma mais simples que os mesmos podem tomar consiste em utilizar dois vectores guias para guardar as posições dos elementos não nulos, que são armazenados num terceiro vector.
3
1 Técnicas de Armazenamento de Matrizes
Por exemplo, a matriz
0000
00000
000
0000
5414
63
523222
4111
=
aa
a
aaa
aa
A
com número de elemento nulos = 16 e índice de esparsidade = 16 / ( 4× 6) =2 / 3 , pode ser armazenada através dos seguintes vectores
[ ][ ]
[ ]5414635232224111
51653241
44322211
aaaaaaaaA
C
L
='==
Para referenciar o elemento A ( i , j ) é necessário percorrer os vectores guias L e C até encontrar um índice k tal que L ( k ) = i e C ( k ) = j , sendo então A ( i , j ) = A ' ( k ) . Se não se encontrar um índice que satisfaça a condição referida, então o elemento A ( i , j ) é nulo.
Caso cada elemento da matriz e dos vectores guias ocupe a mesma memória de computador, este processo só será rentável em termos de espaço para índices de esparsidade superiores a 2/3.
A rentabilidade em termos de tempo de cálculo só poderá ser definida em função das operações a efectuar. No entanto, chamase a atenção para o facto de a ''procura'' do índice de um dado elemento ser uma operação longa, repetitiva, e dificilmente optimizável.
1.2 Matrizes Simétricas
Uma matriz simétrica é, como se sabe, uma matriz quadrada (nº de linhas = nº de colunas = dimensão) em que A ( i , j ) = A ( j , i ) . Os processos antes definidos para operar matrizes em geral podem ser adaptados à sua estrutura de forma a reduzirse ainda mais o tempo de cálculo e/ou a memória necessária para as operar.
1.2.1 Matrizes Não Esparsas
Neste caso, é necessário guardar em memória todos os elementos do triângulo superior (ou inferior) e a diagonal da matriz, num total de (dimensão + 1) × dimensão / 2 elementos. O armazenamento é feito num vector com esta dimensão.
4
1 Técnicas de Armazenamento de Matrizes
Exemplificando, a matriz
44433442244114
433332233113
4232222112
41312111
= = =
= =
= =
aaaaaaa
aaaaaa
aaaaa
aaaa
A
será armazenada no vector
[ ]44434241333231222111= aaaaaaaaaaA'
em que:
≤×≥×
jiseijjA
jisejiiAjiA
+/)((' +/)(('
, , ,
)21
)21=)(
Este processo é semelhante ao utilizado para armazenar matrizes não simétricas em vector, mas exige um maior número de operações. É agora necessário efectuar um teste lógico, duas adições e duas multiplicações para encontrar um elemento, enquanto no caso anterior bastava fazer duas adições e uma multiplicação.
A escolha da utilização ou não deste armazenamento está dependente do habitual binómio tempo de cálculo/memória.
1.2.2. Matrizes Bandeadas
Para operar matrizes simétricas bandeadas basta guardar em memória os elementos da banda que pertencem à diagonal e ao triângulo superior. Costuma referirse este conjunto como a semibanda, e o seu número de colunas como a largura da semibanda.
Exemplificando, a matriz
5554455335
544443344224
53433332233113
4232222112
312111
==
==
==
=
=
00
0
0
00
aaaaa
aaaaaa
aaaaaaa
aaaaa
aaa
A
5
1 Técnicas de Armazenamento de Matrizes
com largura de banda = 5 e largura da semibanda = 3, pode ser armazenada em
55
5444
534333
423222
312111
=
a
aa
aaa
aaa
aaa
A'
em que
≥≥
casos restantes nossese
0
)1(
)1(
= )( jilsbanjjijA
ijlsbaniijiA
jiA > ++,' > ++,'
,
1.2.3 Matrizes Esparsas
Um processo de armazenar matrizes simétricas esparsas é proceder de uma forma semelhante à anteriormente indicada, mas guardando unicamente os elementos não nulos do triângulo superior.
Para encontrar o elemento A ( i , j ) deve verificarse se j é maior ou igual a i. Nesse caso, procurase i no vector L e j no vector C. Caso contrário, devese procurar j no vector L e i no vector C.
1.2.4 Armazenamento em Perfil
As matrizes de rigidez de sistemas estruturais têm normalmente uma estrutura em banda mas, se a numeração dos nós não é feita convenientemente, aparecem alguns elementos que ''alargam'' a faixa central e obrigam a definir uma semibanda com um grande número de elementos nulos.
Existe um grande número de métodos que permitem armazenar e operar eficientemente este tipo de matrizes. Deles, o mais vulgar é o armazenamento por altura efectiva de perfil ou ''skyline''.
Este método, ao contrário do anteriormente exposto, exige um vector guia com dimensão igual à dimensão da matriz a guardar. Em contrapartida, pode ser necessário armazenar elementos nulos.
6
1 Técnicas de Armazenamento de Matrizes
Exemplificando, na matriz
66
6555
5444
5333
423222
412111
Simétrica
0
00
00
000
=
a
aa
aa
aa
aaa
aaa
A
o vector guia indica quais as posições ocupadas no vector auxiliar pelos elementos da diagonal da matriz.
[ ]14129531=L
[ ]665544332211
14129531aaaaaa 65545342413221 0
13111087642
aaaaaaaA ='
Se se pretende obter um determinado elemento da matriz, o processo envolve uma busca no vector guia, relativamente complicada mas não tão repetitiva como para o armazenamento esparso. Se se quer ter uma listagem dos elementos que podem não ser nulos numa dada coluna col , basta fazer:
A' ( k ) , com k desde L(col1) +1 até L(col),
excepto no caso da coluna 1 que contém unicamente o elemento a11, necessáriamente guardado
em A'(1).
Esta forma de armazenamento tem a grande vantagem de, no caso (que é o mais vulgar) de se querer resolver um sistema de equações por eliminação de Gauss ou decomposição de Choleski, a matriz ''original'' ter uma estrutura igual à que vai sendo obtida durante o processo de solução.
1.3 Operação de Matrizes em Computador
As técnicas de armazenamento descritas anteriormente não devem ser julgadas unicamente com base nas indicações fornecidas.
Efectivamente, não se deve tentar efectuar uma operação sem procurar ''optimizar'' a maneira como os índices são referidos.
7
1 Técnicas de Armazenamento de Matrizes
Apresentamse seguidamente algumas versões possiveis para efectuar operações matriciais simples. O seu interesse principal é chamar a atenção para alguns dos erros e vícios mais comuns.
1.3.1 Somar duas matrizes simétricas cheias : C ij = A ij + B ij
Algoritmo 1 É mau
Entrada: dimensão dimensão das matrizes,matriz A de dimensão × dimensão,matriz B de dimensão × dimensão.
Saída: matriz C de dimensão × dimensão.
PARA i = 1 ATÉ dimensãoPARA j = 1 ATÉ dimensão
SE i > j ENTÃO l = (i 1) × i / 2 + jOU ENTÃO l = (j 1) × j / 2 + i
C(l) = A(l) + B(l)FIM
Algoritmo 2 Funciona, mas não é eficiente
Entrada: igual ao algoritmo 1.Saída: igual ao algoritmo 1.
PARA i = 1 ATÉ dimensãoPARA j = i ATÉ dimensão
(porque j é sempre ≥ i)l = (j 1) × j / 2 + iC(l) = A(l) + B(l)
FIM
8
1 Técnicas de Armazenamento de Matrizes
Algoritmo 3 Um pouco melhor
Entrada: igual ao algoritmo 1.Saída: igual ao algoritmo 1.
PARA j = 1 ATÉ dimensão(l posição do primeiro elemento da coluna j)l = (j 1) × j / 2 + 1PARA i = 1 ATÉ j
C(l) = A(l) + B(l)l = l + 1
FIM
Algoritmo 4 Provavelmente o mais eficaz, mas também o mais difícil de ''entender''
Entrada: igual ao algoritmo 1.Saída: igual ao algoritmo 1.
nº elementos = (dimensão+1) × dimensão/2PARA l=1 ATÉ nº elementos
C(l) = A(l) + B(l)FIM
1.3.2 Multiplicar uma matriz simétrica bandeada por um vector ci = Aij bj
Algoritmo 1 Certo, curto, mas explicado
Entrada: dimensão dimensão da matriz,lsban dimensão da semibanda da matriz,matriz A de dimensão×lsban, contendo os elementos da semibanda da matriz,vector B de dimensão.
Saída: vector C de dimensão.
9
1 Técnicas de Armazenamento de Matrizes
PARA i = 1 ATÉ dimensãoC(i) = 0PARA j = 1 ATÉ dimensão
SE j ≥ i e j < i + lsban ENTÃOcoluna = j i + 1C(i)=C(i)+A(i, coluna) × B(j)
SE i > j e i < j + lsban ENTÃOcoluna = i j + 1C(i) = C(i) + A(j,coluna) × B(j)
FIM
Algoritmo 2 Sendo mais comprido, não passa a vida a fazer testes. Alem disso só dá importância à semibanda
Entrada: igual ao algoritmo 1.Saída: igual ao algoritmo 1.
ls = lsban 1PARA i=1 ATÉ dimensão
soma = 0.(operar o triângulo inferior índices trocados)PARA j = Máximo (1, i ls) ATÉ i 1
l = i j + 1soma = soma + A(j,l) × B(j)
(operar a diagonal sempre na coluna 1)soma = soma + A(i,1) × B(i)(operar o triângulo superior)PARA j = i + 1 ATÉ Mínimo (dimensão , i + ls)
l = j i + 1soma = soma + A(i,l) × B(j)
C(i) = somaFIM
Este algoritmo pode ainda ser optimizado passando a utilizar duas novas variáveis, por exemplo i1 = i + 1 e i2 = i 1, que são substituidas no cálculo do valor de l dentro dos dois ciclos em j.
10
1 Técnicas de Armazenamento de Matrizes
Outras ''optimizações'' irão provavelmente provocar uma maior ilegibilidade do programa. Notese que tanto no algoritmo 1, como no 2 a variável l é sempre calculada de acordo com a mesma fórmula mas, no segundo exemplo,os testes lógicos são eliminados dividindo os ciclos. Além disso, a variável i indica sempre a linha da matriz em causa e j a sua coluna. É sempre de pensar duas vezes antes de pôr um programa muito rápido, mas incompreensível.
Algoritmo 3 Muito rápido, mas pelo menos à primeira vista é mais difícil
Entrada: igual ao algoritmo 1.Saída: igual ao algoritmo 1.
PARA i=1 ATÉ dimensão(operar a diagonal, inicializando C)C(i) = A(i,1) × B(i)
PARA i = 1 ATÉ dimensão 1(l indica a coluna que está armazenada na coluna 2 da linha i)l = i + 1matfim = dimensão i + 1PARA j = 2 ATÉ Minimo(lsban,matfim)
(operar o triângulo inferior)C(l) = C(l)+A(i,j) × B(i)(operar o triângulo superior)C(i) = C(i)+A(i,j) × B(l)l = l+1
FIM
Este algoritmo não efectua os testes lógicos do algoritmo 2. Em tempo de cálculo tem a vantagem de fazer menos ciclos e menos operações.
11
1 Técnicas de Armazenamento de Matrizes
1.3.3 Multiplicar uma matriz simétrica armazenada em perfil por um vector ci=Aijbj
Algoritmo 1 Certo, explicado, mas lento e repetitivo
Entrada: dimensão dimensão da matriz,vector guia do perfil L de dimensão , (necessita da posição L(0) = 0),vector A contendo os elementos do perfil da matriz,vector B de dimensão.
Saída: vector C de dimensão.
PARA i=1 ATÉ dimensãoC(i) = 0.PARA j = 1 ATÉ dimensão
SE j ≥ i ENTÃOposdiag = L(j)ipos = posdiag j + iSE ipos > L(j 1) ENTÃO
C(i) = C(i) + A(ipos) × B(j)SE j < i ENTÃO
posdiag = L(i)ipos = posdiag i + jSE ipos > L(i1) ENTÃO
C(i) = C(i) + A(ipos) × B(j)FIM
Algoritmo 2 Melhor, mas ainda necessitando de testes
Entrada: igual ao algoritmo 1 (não necessita do elemento L(0) definido).Saída: igual ao algoritmo 1.
12
1 Técnicas de Armazenamento de Matrizes
PARA j = 1 ATÉ dimensãoC(j) = A(L(j)) × B(j)
PARA i = 1 ATÉ dimensão 1PARA j = i + 1 ATÉ dimensão
posdiag = L(j)ipos = posdiag j + iSE ipos > L(j1) ENTÃO
C(i) = C(i) + A(ipos) × B(j)C(j) = C(j) + A(ipos) × B(i)
FIM
Algoritmo 3 Melhor, adaptado e seguindo o tipo de arrumação
Entrada: igual ao algoritmo 1 (não necessita do elemento L(0) definido).Saída: igual ao algoritmo 1.
PARA j = 1 ATÉ dimensãoC(j) = A(L(j)) × B(j)
PARA j=2 ATÉ dimensãonelemj = L(j) L(j1)linhaj = j nelemj + 1ipos = L(j1) + 1PARA i = linhaj ATÉ j 1
C(i) = C(i) + A(ipos) × B(j)C(j) = C(j) + A(ipos) × B(i)ipos = ipos + 1
FIM
A única forma de julgar dois algoritmos alternativos (e supostamente certos) é fazêlos correr em contrarelógio na máquina e situação em que vão ser utilizados. O resultado depende da dimensão da matriz, do computador e do compilador, pelo que não se pode garantir que um algoritmo marginalmente mais rápido seja sempre o melhor.
13
Capítulo 2
Métodos Directos e Iterativos de
Resolução de Sistemas de Equações
2.1 Métodos de Gauss e GaussJordan
Para resolver um sistema de equações pelo método de Gauss ou GaussJordan aplicamse três operações distintas nas equações intervenientes. Estas operações, abaixo enunciadas, têm a propriedade de originarem novos sistemas de equações cujas soluções são as mesmas do sistema inicial.
i. A equação Ei pode ser multiplicada por uma constante λ, diferente de zero, e a equação resultante λEi pode ser usada em vez de Ei.
ii. A equação Ei pode ser multiplicada por uma constante λ qualquer, somada à equação Ej, e a equação resultante Ej+ λ Ei pode ser usada em vez de Ej.
iii. As equações Ei e Ej podem ser trocadas de ordem.
Quando o sistema de equações é escrito numa matriz de n linhas por n + m colunas, em que n é a dimensão do sistema e m o número de termos independentes, as operações acima enunciadas implicam que:
i. se pode multiplicar uma linha da matriz por uma constante não nula qualquer;
14
2 Métodos de Resolução de Sistemas de Equações
ii. se pode adicionar a uma linha uma outra linha multiplicada por uma constante qualquer;
iii. se podem trocar duas linhas da matriz.
Estas propriedades são aplicadas de forma a transformar a matriz do sistema numa matriz triangular superior (método de Gauss) ou diagonal (método de GaussJordan).
Exemplificando o sistema de equações:
=+++=+++=+++−
−=++−
44321
34321
24321
14321
9432
023
12342
22324
EXXXX
EXXXX
EXXXX
EXXXX
pode ser escrito:
−
−−
=
94132
01213
123142
22324
A
A propriedade ii pode ser aplicada às linhas 2, 3 e 4 de forma a anular os elementos da primeira coluna destas linhas:
1 Linhaaa4 Linha4 Linha
1 Linhaaa3 Linha3 Linha
1 Linhaaa2 Linha2 Linha
1141
1131
1121
−=−=−=
resultando em:
−−−
−−
=
1032140
232141250
1142530
22324
A
Para eliminar os elementos da segunda coluna das linhas 3 e 4 fazse:
2 Linhaaa4 Linha4 Linha
2 Linhaaa3 Linha3 Linha
2242
2232
−=−=
15
2 Métodos de Resolução de Sistemas de Equações
−−−−−−
−−
=
3143762300
3236233700
1142530
22324
A
Finalmente para anular o elemento da terceira coluna, quarta linha fazse:
3 Linhaaa4 Linha4 Linha 3343−=
−−−
−−
=
1411128111000
3236233700
1142530
22324
A
Para saber qual o valor que as variáveis tomam basta caminhar de trás para a frente. A quarta equação diz que:
214111111281411128111 44 =×=⇒= XX
a terceira que
( ) 026233237332362337 343 =×+−×−=⇒−=−− XXX
e assim sucessivamente.
A solução do sistema é:
.2;0;1;1 4321 ===−= XXXX
O método de Gauss Jordan difere do apresentado por também se eliminarem os elementos acima da diagonal da coluna em operação, resultando numa matriz diagonal cuja solução é imediata. No problema apresentado faziase:
1 Linhaaa4 Linha4 Linha
1 Linhaaa3 Linha3 Linha
1 Linhaaa2 Linha2 Linha
1141
1131
1121
−=−=−=
depois
2 Linhaaa4 Linha4 Linha
2 Linhaaa3 Linha3 Linha
2 Linhaaa1 Linha1 Linha
2242
2232
2212
−=−=−=
16
2 Métodos de Resolução de Sistemas de Equações
3 Linhaaa4 Linha4 Linha
3 Linhaaa2 Linha2 Linha
3 Linhaaa1 Linha1 Linha
3343
3323
3313
−=−=−=
4 Linhaaa3 Linha3 Linha
4 Linhaaa2 Linha2 Linha
4 Linhaaa1 Linha1 Linha
4434
4424
4414
−=−=−=
Depois de ter feito todas estas contas achase o valor dos vários Xi com:
iiii abX =
O número de operações matemáticas necessárias para aplicar estes dois métodos de resolução de sistemas de equações estão indicadas na Tabela 1.
Mult./Divisões Adiç/SubtracçõesGauss 33 23 nmnn −+ ( ) ( )nmnmn −+−+ 61213 23
GaussJordan 22 23 nmnn −+ ( ) ( )nmnmn −+−+ 2112 23
Tabela 1
Verificase que o método de Gauss implica um número de operações menor do que o método de GaussJordan.
Este pode contudo ter a vantagem de ser mais fácil de programar.
Para aplicar o método de Gauss em ''problemas gerais'' convém aplicar um algoritmo que verifique se os elementos da diagonal não se anulam. Caso isso aconteça deve aplicarse a propriedade iii de forma a passar a estar na diagonal um elemento não nulo. Se isso não fôr possível, então a matriz não tem inversa.
Devido a erros de arredondamento pode acontecer que surja um elemento na diagonal que parece ser não nulo, mas que o deveria ser. Se ele é utilizado para divisor aparecem grandes erros numéricos. Dada a dificuldade em distinguir um zero de um valor pequeno, é preferível seguir uma estratégia que evite este problema.
17
2 Métodos de Resolução de Sistemas de Equações
Em problemas delicados é vantajoso aplicar a propriedade iii de forma a colocar na diagonal o elemento com maior valor absoluto da coluna. Garantese assim, à custa de operações adicionais, uma muito maior estabilidade numérica.
O algoritmo que se apresenta procede desta forma. Se fôr aplicado a um problema que se saiba à priori ser bem comportado, então a fase de ''colocar o maior da coluna na diagonal'' pode ser eliminada.
Algoritmo ''GAUSS''
Entrada: n número de variáveis,m número de termos independentes,matriz A tal queA(i,j) = aij (i=1...n, j=1...n)A(i,j+n) = bij (i=1...n, j=1...m)
Saída: matriz A em que os termos dos bij correspondem aos Xij
(Triangularização)PARA i = 1 ATÉ n 1
''Pôr o maior da coluna i na diagonal''PARA j = i + 1 ATÉ n
factor = A(j,i) / A(i,i)PARA k = i ATÉ n + m
A(j,k) = A(j,k) + factor × A(i,k)SE A(n,n) é ''zero'' ENTÃO Erro ''Determinante nulo''
(Retrosubstituição)PARA j = n + 1 ATÉ n + m
PARA i = n ATÉ 1 (passo 1)soma = A(i,j)PARA k = i + 1 ATÉ n
soma = soma A(i,k) × A(k,j)A(i,j) = soma / A(i,i)
FIM.
Rotina para ''Pôr o maior da coluna i na diagonal''maior = A(i,i)
18
2 Métodos de Resolução de Sistemas de Equações
linha = iPARA k = i + 1 ATÉ n
SE Absoluto(A(k,i)) > maior ENTÃOmaior = Absoluto(A(k,i))linha = k
SE linha ð i ENTÃO ''Troca linhas i e linha '' ??????SE maior é ''zero'' ENTÃO Erro ''Determinante nulo''FIM (Pôr o maior da coluna i na diagonal)
Na maior parte dos problemas de análise estrutural a matriz do sistema é simétrica, e também positiva definida. O facto de a matriz ser definida permite que se possa aplicar o método de Gauss sem fazer escolha do elemento a colocar na diagonal. A simetria da matriz permite economizar nas operações aritméticas envolvidas.
Relembrando as matrizes que foram aparecendo na resolução do exemplo anterior:
−−
−
9
0
12
2
4132
1213
3142
2324
−−−
−−
10
23
11
3214
214125
4253
0
0
0
22324
−−
−−−−
−−
314
323
37623
62337
00
00
11
2
425
23
30
24
−
−
−−
−
1411128111000
323
11
2
623
4
3
3700
2530
324
Verificase que as submatrizes ''por eliminar'' continuam a ser simétricas, o que quer dizer que por exemplo quando se passou da primeira para a segunda submatriz se calculou duas vezes o valor 5/2, que surge tanto para a32 como para a23. A mesma coisa se passa com a42 e a24 e ainda com a43 e a34.
19
2 Métodos de Resolução de Sistemas de Equações
O algoritmo seguinte ''Gauss simétricas'' opera tendo em conta esta propriedade. Além disso ''não se dá ao trabalho'' de operar os elementos do triângulo inferior.
Algoritmo ''GAUSS SIMÉTRICAS''
Entrada: Igual a ''GAUSS''
Saida : Igual a ''GAUSS''
(Triangularização)PARA i = 1 ATÉ n 1
PARA j = i + 1 ATÉ nfactor = A(i,j) / A(i,i)PARA k = j ATÉ n + m
A(j,k) = A(j,k) + factor × A(i,k)
(A retrosubstituição é igual à do algoritmo ''GAUSS'')
FIM
Para este algoritmo estão resumidas na Tabela 2 o número de operações aritméticas:
Mult./DivisõesGausssimétricas ( ) ( ) 31216 23 nmnmn +−++Tabela 2
Mais adiante este algoritmo voltará a ser tratado para ''optimizar'' o processamento de matrizes armazenadas na forma simétrica.
2.2 Métodos de Factorização de Matrizes
Todos os métodos de factorização decompõem a matriz de sistema em duas matrizes, uma triangular inferior e uma triangular superior, de tal forma que o seu produto seja igual à matriz original:
ULA=
Os diferentes métodos distinguemse pela escolha que fazem para os elementos da diagonal de L e/ou de U.
O método de Doolittle obriga L(i,i) = 1
20
2 Métodos de Resolução de Sistemas de Equações
O método de Crout obriga U(i,i) = 1O método de Choleski obriga L(i,i) = U(i,i)
Esta escolha faz com que os restantes elementos possam ser imediatamente determinados, como se pode ver no exemplo seguinte.
=
33
2322
131211
333231
2221
11
333231
232221
131211
00
00
00
U
UU
UUU
LLL
LL
L
AAA
AAA
AAA
×+×+××+×××+××+××
×××
333323321331223212311131
23221321222212211121
131112111111
ULULULULULUL
ULULULULUL
ULULUL
Aplicando por exemplo o método de Crout, calculamse os elementos de L e U igualando os elementos de A com o desenvolvimento do produto matricial LU.
( )
233213313333
12313232
3131
2213212323
12212222
2121
111313
111212
1111
332332133133
32123132
3131
2322132123
22122122
2121
131113
121112
1111
1
1
1
1
1
1
ULULAL
ULAL
AL
LULAU
ULAL
AL
LAU
LAU
AL
LULULA
LULA
LA
ULULA
LULA
LA
ULA
ULA
LA
×−×−=×−=
=×−=
×−=====
⇒⇒⇒⇒⇒⇒⇒⇒⇒
×+×+×=×+×=×=
×+×=×+×=×=
×=×=
×=
Um algoritmo geral de factorização de matrizes pode ser:
Algoritmo ''FACTORIZAÇÃO GERAL''
Entrada: n número de variáveis,matriz A tal que A(i,j) = aij (i=1...n, j=1...n)
Saída: matrizes U e L , factores de A .
21
2 Métodos de Resolução de Sistemas de Equações
PARA i = 1 ATÉ nsoma = A(i,i)PARA k = 1 ATÉ i 1
soma = soma L(i,k) × U(k,i)SE soma é ''zero'' ENTÃO Erro ''Determinante nulo''
SE Doolittle ENTÃOL(i,i) = 1U(i,i) = soma
SE Crout ENTÃOL(i,i) = somaU(i,i) = 1
SE Choleski ENTÃOsomaiiL =),(
U(i,i) = L(i,i)PARA j = i + 1 ATÉ n
soma1 = A(i,j)soma2 = A(j,i)PARA k = 1 ATÉ i 1
soma1 = soma1 L(i,k) × U(k,j)soma2 = soma2 L(j,k) × U(k,i)
U(i,j) = soma1/L(i,i)L(j,i) = soma2/U(i,i)
FIM
Notese que os vários métodos são muito semelhantes, sendo até possível escrever um único algoritmo que cubra os três casos. No entanto do ponto de vista computacional não é vantajoso fazêlo. É preferível escolher um método e tentar adaptar o programa à escolha feita.
No caso de matrizes simétricas e positivas definidas é sempre melhor utilizar o método de Choleski, pois as matrizes factores são a transposta uma da outra, pelo que basta calcular uma delas.
Depois de factorizadas as matrizes, o cálculo da solução do sistema de equações passa ainda por duas retrosubstituições.
Querse conhecer L×U×x = b.
Calculase um vector y tal que L×y = b, ou seja U×x = y.
O primeiro passo é calcular y e o segundo calcular x a partir dele. Uma das maiores vantagens da resolução de sistemas de equações por factorização é que estes cálculos são completamente
22
2 Métodos de Resolução de Sistemas de Equações
independentes da factorização propriamente dita, isto é, desde que existam os factores da matriz do sistema ele pode ser resolvido para quantos termos independentes for necessário.
Algoritmo ''FACTORIZAÇÃO GERAL'' (Continuação)
Entrada: n número de variáveis,matrizes U e L , factores de A .vector B tal que B(i) = bi (i=1...n)
Saída: vector da solução X.
(resolver L×y = b)PARA i = 1 ATÉ n
soma = B(i)PARA k = 1 ATÉ i 1
soma = soma L(i,k) × Y(k)Y(i) = soma/L(i,i)
(resolver U×x = y)PARA i = n ATÉ 1 (passo 1)
soma = Y(i)PARA k = i + 1 ATÉ n
soma = soma U(i,k) × X(k)X(i) = soma/U(i,i)
FIM
Do ponto de vista da eficiência, o método de Choleski ''brilha'' para resolver sistemas de equações simétricos de grandes dimensões. Tem o inconveniente de exigir o cálculo de tantas raízes quadradas como a dimensão do sistema, mas como este é um factor que cresce linearmente, não é relevante para sistemas muito grandes.
O número de operações exigidas é o indicado no Tabela 3.
Mult./Divisões Raízes quadradasCholeski ( ) 643 23 nnn −+ nTabela 3
23
2 Métodos de Resolução de Sistemas de Equações
Apresentamse dois algoritmos para a factorização de Choleski. Eles diferem unicamente na ordem com que as contas são feitas. O primeiro calcula a matriz L e segue o raciocínio do algoritmo ''FACTORIZAÇÃO GERAL''. O segundo calcula a matriz U e por ''trocar as contas'' é, como se verá mais adiante, mais eficaz para ser adaptado a outras formas de armazenamento da matriz.
Algoritmo ''CHOLESKI I''
Entrada: n número de variáveis,matriz A tal que A(i,j) = aij (i=1...n, j=1...n)
Saída: matriz L , factor de A .
PARA i = 1 ATÉ nsoma = A(i,i)PARA k = 1 ATÉ i 1
soma = soma L(i,k)2
somaiiL =),(
PARA j = i + 1 ATÉ nsoma = A(j,i)PARA k = 1 ATÉ i 1
soma = soma L(j,k) × L(i,k)L(j,i) = soma/L(i,i)
FIM
Algoritmo ''CHOLESKI II''
Entrada: n número de variáveis,matriz A tal que A(i,j) = aij (i=1...n, j=1...n)
Saída: matriz U , factor de A .
24
2 Métodos de Resolução de Sistemas de Equações
PARA j = 1 ATÉ nPARA i = 1 ATÉ j 1
soma = A(i,j)PARA k = 1 ATÉ i 1
soma = soma U(k,i) × U(k,j)U(i,j) = soma/U(i,i)
soma = A(j,j)PARA k = 1 ATÉ j 1
soma = soma U(k,j)2
somajjU =),(
FIM
2.3 Métodos Iterativos de Resolução de Sistemas
De entre os métodos iterativos para a resolução de sistemas de equações temos de considerar dois grandes grupos: os métodos de relaxação e os métodos projectivos.
Todos eles resolvem um sistema de equações de cada vez, sendo a principal diferença entre estes métodos a estratégia utilizada para as correcções sucessivas.
Destacamse os seguintes métodos de relaxação : Jacobi, GaussSeidel e SOR. Dos métodos projectivos destacase o método do gradiente e o método dos gradientes conjugados.
2.3.1 Critérios de Paragem para Métodos Iterativos
A convergência de um processo iterativo é atingida quando se tem:
∆<− ixAb
sendo ∆ uma tolerância absoluta dada.
Regra geral, preferese utilizar um critério em que se considere um erro relativo δ.
δ<−
b
xAb i
ou
25
2 Métodos de Resolução de Sistemas de Equações
δ<− −
i
ii
x
xx 1
2.3.2 Métodos de Relaxação
Os métodos de relaxação baseiamse em algoritmos que necessitam de decompor a matriz do sistema A em duas outras M e N tais que:
A = M + N
Assim o sistemaA x = b
tomará a forma
M x + N x = b
ou
x = M1 (b N x)
da qual se pode tomar a forma iterativa:
xi = M1 b M1 N xi 1
Os vários métodos diferem fundamentalmente na escolha da forma como se decompõe a matriz do sistema.
Método de Jacobi
No método de Jacobi é considerada para matriz M a matriz que contém os termos diagonais da matriz A.
=
nna
a
a
a
M
33
22
11
00
00
00
26
2 Métodos de Resolução de Sistemas de Equações
=
0
0
0
0
3231
2321
1312
aa
aa
aa
N
A vantagem evidente é a fácil inversão da matriz M, os elementos de M1 são iguais aos inversos algébricos dos elementos da diagonal de M.
A equação matricial
xi = M1 b M1 N xi 1
pode ser escrita indicialmente como
−= ∑
≠ jkkjkj
jjj xab
ax 1-ii 1
A solução inicial costuma ser tomada igual a zero. Daí vem que
x1 = M1 b
ou indicialmente
jj
jj a
bx =1
Método de GaussSeidel
No método de GaussSeidel é considerada para matriz M o triângulo inferior de A, ou seja:
=
nna
aaa
aa
a
M
333231
2221
11
0
00
=
0
000
00
0
23
1312
a
aa
N
27
2 Métodos de Resolução de Sistemas de Equações
Neste caso M1 não é calculada explicitamente, mas as soluções da equação
M xi = b N xi1
são obtidas por retrosubstituição.
Em notação indicial a fórmula iterativa escrevese:
−−= ∑∑
+=
−
=
n
jkkjk
j
kkjkj
jjj xaxab
ax
1
1-i1
1
ii 1
que é facilmente automatizável.
Método SOR (Successive OverRelaxation)
Este método é uma adaptação do método de GaussSeidel e tenta acelerar a convergência por meio da introdução de um parâmetro correctivo que majore o passo.
O passo é dado por:
1-iiijjj xxx −=∆
ou seja
1-i
1
1-i1
1
ii 1j
n
jkkjk
j
kkjkj
jjj xxaxab
ax −
−−=∆ ∑∑
+=
−
=.
Corrigindo o passo por meio de um coeficiente ω diferente de zero podese escrever a fórmula iterativa como:
i1-iijjj xxx ∆+= ω
logo
( )
−−+−= ∑∑
+=
−
=
n
jkkjk
j
kkjkj
jjjj xaxab
axx
1
1-i1
1
i1-ii 1ωω
Provase que para cada matriz existe um parâmetro de sobrerelaxação óptimo que conduz à melhor convergência possível. Este toma em geral valores entre 1 e 2, mas só é possível conhecê
28
2 Métodos de Resolução de Sistemas de Equações
lo em função dos valores próprios da matriz, e o seu cálculo viria a despropósito, pois seria mais pesado do que resolver o sistema de equações sem sobrerelaxação.
Optase por utilizar um processo iterativo em que se estabelece um valor de ω arbitrário (normalmente 1, mas sempre entre 1 e 1.2), que vai sendo corrigido em função da evolução do processo iterativo.
2.3.3 Métodos Projectivos
Os métodos projectivos baseiamse na minimização do resíduo:
r = b A x
através da minimização da função
h2 = rT A1 r
Sendo a matriz A simétrica e positiva definida, quando r for igual a zero então h2 é mínimo e também igual a zero.
Desenvolvendo a expressão de h2,
( ) ( )bAbxbxAxh
bxAAbxAh1TTT2
1T2
2 −
−
+−=−−=
e considerando para x a forma iterativa,
kk1k dxx α+=+
em que α representa o comprimento do passo e dk a direcção seguida na iteração k, temos:
( ) ( ) ( ) ( ) kT1TkTkkTkkTk22 22 xbbAbxAxxAbddAdh −++−−= −αα
ou seja
( ) ( ) ( ) kT1TkTkkTkkTk22 22 xbbAbxAxrddAdh −++−= −αα
Considerando que h2 varia quadraticamente com α então tem um mínimo local quando ( )0
2
=∂α
∂ h
, tem que se procurar o valor de α para o qual
29
2 Métodos de Resolução de Sistemas de Equações
( ) ( ) ( ) 02 kkTk2
=−= rdAdh α
∂α∂
logo( )
( ) kk
kkk
T
T
dAd
rd=α
Todos os métodos projectivos se baseiam nestas conclusões:
( )( ) kk
kkk
kkk1k
T
T
dAd
rd
dxx
=
+=+
α
α
A diferença entre os diferentes métodos projectivos está na escolha da direcção dk seguida em cada iteração.
Método do Gradiente
Neste método é escolhida em cada iteração a direcção correspondente ao resíduo no ponto x que coincide com o gradiente de h2. Assim, podese formular o seguinte algoritmo que minimiza as operações matriciais envolvidas:
Fazer
r0 = b A x0
Iterar
( )[ ] ( )[ ]kkk1k
kkk1k
kkkkk
kk
=
=
=
=TT
urr
rxx
urrr
rAu
αα
α
−+
+
+
O processo deverá parar quando rk for muito pequeno. Convém ter cuidado porque se podem acumular demasiados erros numéricos (rk só é calculado explicitamente na iteração zero). Caso o processo pareça não convergir, convém reiniciar o cálculo com o valor actual de x.
30
2 Métodos de Resolução de Sistemas de Equações
Método dos Gradientes Conjugados
Este método procura que as direcções escolhidas para as várias iterações sejam conjugadas entre si, isto é procura que se verifique:
( ) ji 0ji T ≠= que semprepAp
Um algoritmo possível é, o apresentado por Jennings [1]:
Tomese
p 0 = r 0 = b A x 0
iterese da seguinte forma
( )[ ] ( )[ ]
( )[ ] ( )[ ]kk1+k1k
kkk1+kk
kkk1k
kkk1k
kkkkk
kk
=
-=
=
=
=
=
TT
TT
prp
upur
urr
pxx
uprp
pAu
ββ
αα
α
+
−+
+
+
+
Teoricamente a solução "exacta" é obtida ao fim de n iterações, mas dependendo do sistema em causa ela pode ser atingida antes disso. Pode também acontecer que devido a erros numéricos sejam necessárias mais de n iterações para se atingir a convergência.
2.4 Adaptação dos Métodos de Resolução à Forma de Armazenamento das Matrizes
Os algoritmos até agora apresentados para a resolução de sistemas de equações supõem que se está a operar com matrizes cheias, não simétricas.
No caso de sistemas estruturais aparecemnos normalmente matrizes simétricas, em banda ou esparsas.
Apresentamse algumas adaptações dos algoritmos considerados para estes tipos de armazenamento.
31
2 Métodos de Resolução de Sistemas de Equações
2.4.1 Matrizes Simétricas
O algoritmo ''Gausssimétricas'' é de fácil adaptação. No entanto, não é conveniente transformar directamente os índices (i ,j ) no valor (k ) correspondente em armazenamento simétrico.
Os elementos da matriz referidos são A(i,i) , A(i,j) , A(j,k) , A(i,k) .
A posição de A(i,i) é calculada através de um apontador lii que é posto a zero antes de iniciar o cálculo, e que é incrementado de i dentro de cada ciclo, ficando a valer 1, 3, 6, 10, ... apontando para os elementos da diagonal.
A(i,j) para j =i + 1 até n são os elementos da linha i à direita da diagonal. A posição de um elemento (lij ) pode ser calculada a partir da posição da diagonal (lii ), a que se vai adicionando sucessivamente o número da coluna anterior.
A(i,k) para k = j até n é calculado de um modo semelhante.
A(j,k) para k = j até n começa não na diagonal i (apontada por lii ), mas na diagonal j (apontada por ljj ), mas em tudo o resto funciona de um modo idêntico.
Para tentar seguir o algoritmo é conveniente seguir a Tabela das posições dos elementos.
Algoritmo ''GAUSS ARMAZENAMENTO SIMÉTRICO''
Entrada: n número de variáveis, m número de termos independentes,vector A com a ij em armazenamento simétrico,matriz B com B(i,j) = b ij (i = 1 ..n , j = 1 ..m )
Saída: matriz B em que os termos dos b ij correspondem aos X ij
32
2 Métodos de Resolução de Sistemas de Equações
lii = 0PARA i = 1 ATÉ n
lii = lii + ilij = lii + iljj = lij + 1PARA j = i + 1 ATÉ n
factor = A(lij) / A(lii)lik = lijljk = ljjPARA k = j ATÉ n
A(ljk) = A(ljk) + A(lik) × factorlik = lik + kljk = ljk + k
PARA k = 1 ATÉ mB(j,k) = B(j,k) + B(i,k) × factor
lij = lij + jljj = ljj + j + 1
FIM
O método de Choleski no algoritmo II também é muito fácil de adaptar usando um raciocínio semelhante ao exposto. Deixase isso como um exercício.
2.4.2 Matrizes Simétricas Bandeadas
Tanto o método de Gauss como o de Choleski são fácilmente utilizáveis para matrizes armazenadas desta forma. Isto porque a matriz resultante da sua aplicação tem a mesma forma que a matriz original, ou seja, se a matriz do sistema tem uma semibanda b , a matriz resultante também a tem.
Basta portanto proceder a uma conveniente alteração dos índices para adaptar o algoritmo ''GaussSimétricas'' a matrizes simétricas bandeadas.
33
2 Métodos de Resolução de Sistemas de Equações
Algoritmo ''GAUSS SIMÉTRICAS BANDEADAS''
Entrada: n número de variáveis,m número de termos independentes,b semi banda da matriz,matriz A(n,b) com a ij armazenada em banda,matriz B com B(i,j) = b ij (i = 1..n , j = 1..m ).
Saída: matriz B em que os termos dos b ij correspondem aos X ij
(Triangularização)PARA i = 1 ATÉ n 1
ofset = 0jreal = iPARA j = 2 ATÉ Minimo (b ,n i + 1)
ofset = ofset + 1jreal = jreal + 1SE jreal = n ENTÃO
factor = A(i,j) / A(i,1)PARA k = 1 ATÉ b ofsetA(jreal,k) = A(jreal,k) + factor × A(i,k+ofset)PARA k = 1 ATÉ m
B(jreal,k) = B(jreal,k) + factor × B(i,k)
(Retrosubstituição)PARA i = n ATÉ 1 (passo 1)
factor = 1 / A(i,1)PARA j = 1 ATÉ m
soma = B(i,j)ipos = iPARA k = 2 ATÉ Minimo (b, n i + 1)
ipos = ipos + 1soma = soma A(i,k) × B(ipos,j)
B(i,j) = soma × factorFIM
34
2 Métodos de Resolução de Sistemas de Equações
Algoritmo ''CHOLESKI SIMÉTRICAS BANDEADAS
Entrada: n número de variáveis,m número de termos independentes,b semi banda da matriz,matriz A(n,b) com a ij armazenada em banda,matriz B com B(i,j) = b ij (i = 1..n , j = 1..m ).
Saída: matriz B em que os termos dos b ij correspondem aos X ij
(Factorização)PARA i = 1 ATÉ n
i1 = i + 1soma = A(i,1)PARA k = Maximo (1, i1 b ) ATÉ i 1
coluna = i1 ksoma = soma A(k,coluna) 2
somaiA =)1,(
PARA j = 2 ATÉ bsoma = A(i,j)PARA k = Maximo (1, i1 b , j b + 1) ATÉ i 1
coluna = i1 ksoma = soma A(k,coluna) × A(k,j)
A(i,j) = soma / A(i,1)
(Substituição)factor = 1 / A (1,1)PARA j = 1 ATÉ m
B (1,j ) = B (1,j ) × factorPARA i = 2 ATÉ n
factor = 1 / A (i ,1)linha_i = Maximo (1, i b + 1)PARA j = 1 ATÉ m
soma = B(i,j)icol = 1PARA k = i 1 ATÉ linha_i (passo 1)
icol = icol + 1soma = soma A(k,icol) × B(k,j)
B(i,j) = soma × factor
35
2 Métodos de Resolução de Sistemas de Equações
(Retrosubstituição)PARA i = n ATÉ 1 (passo 1)
factor = 1 / A (i ,1)PARA j = 1 ATÉ m
soma = B(i,j)ipos = iPARA k = 2 ATÉ Minimo (b , n i + 1)
ipos = ipos + 1soma = soma A(i,k) × B(ipos,j)
B(i,j) = soma × factorFIM
O número de operações necessárias para operar uma matriz simétrica em banda é o indicado na Tabela 4.
Mult./Divisões( ) ( ) ( ) 2322 232 nnmbbbbn ++−−−+
Tabela 4
2.4.3 Matrizes Armazenadas em Perfil
Para se tornar eficiente para matrizes armazenadas em perfil é conveniente trabalhar com uma forma alterada do método de Gauss, que em vez de operar por colunas trabalha por linhas. Soriano [2] chamalhe ''Método de Martin'', e apresenta o seu desenvolvimento.
O algoritmo do método de Choleski é de adaptação imediata.
Algoritmo ''CHOLESKI II PERFIL''
Entrada: n número de variáveis,m número de termos independentes,vector A com a ij armazenada em ''skyline'' segundo o vector guia L contendo L (0) = 0,matriz B com B(i,j) = b ij (i = 1..n , j = 1..m ).
Saída: matriz B em que os termos dos b ij correspondem aos X ij
36
2 Métodos de Resolução de Sistemas de Equações
(Factorização)
)1()1( AA =
PARA j = 2 ATÉ naltura_j = L (j ) L (j 1)pos = L (j 1) + 1linha_j = j altura_j + 1PARA i = linha_j ATÉ j
soma = A(pos)altura_i = L(i) L(i 1)linha_i = i altura_i + 1linha_inicial = Maximo (linha_i , linha_j )pos_i = L(i) i + linha_inicialpos_j = L(j) j + linha_inicialPARA k = linha_inicial ATÉ i 1
soma = soma A(pos_i) × A(pos_j)pos_i = pos_i + 1pos_j = pos_j + 1
SE i = j ENTÃO somaposA =)(
OU ENTÃO A(pos) = soma / A(L(i))pos = pos + 1
(Substituição)factor = 1 / A(1)PARA j = 1 ATÉ m
B(1,j ) = B(1,j ) × factorPARA i = 2 ATÉ n
altura_i = L(i) L(i 1)linha_i = i altura_i + 1factor = 1 / A(L(i))PARA j = 1 ATÉ m
pos = L(i 1) + 1soma = B(i,j)PARA k = linha_i ATÉ i 1
soma = soma A(pos) × B(k,j)pos = pos + 1
B(i,j) = soma × factor
37
2 Métodos de Resolução de Sistemas de Equações
(Retrosubstituição)PARA i = n ATÉ 2 (passo 1)
factor = 1 / A(L(i))PARA j = 1 ATÉ m
B(i,j) = B(i,j) × factoraltura_i = L(i) L(i 1)linha_i = i altura_i +1pos = L(i 1) + 1PARA k = linha_i ATÉ i 1
PARA j = 1 ATÉ mB(k,j) = B(k,j) A(pos) × B(i,j)
pos = pos + 1factor = 1 / A (1)PARA j = 1 ATÉ m
B (1,j ) = B (1,j ) × factorFIM
2.4.4 Algoritmos para os Métodos Iterativos
A forma de armazenamento escolhida no caso de se resolver um sistema de equações por um método iterativo só influi nos produtos matriciais necessários para o método em causa.
Os algoritmos apresentados no Capitulo 1 podem servir de base para estas operações. No entanto convém relembrar (nunca é demais) que não é eficiente seguir explicitamente as fórmulas. É sempre melhor pensar um pouco e adaptar o produto matricial ao algoritmo.
Nos métodos de Jacobi e dos gradientes conjugados é indiferente a ordem pela qual os elementos são calculados (excepto do ponto de vista de erros numéricos) pelo que se pode utilizar qualquer algoritmo eficiente para o produto matricial. No método de Gauss Seidel é necessário calcular sequencialmente os elementos do vector produto pelo que o algoritmo 2 do Cap. 1.3 tem de ser utilizado em vez do 3.
Para diminuir os erros numéricos é costume utilizar variáveis de dupla precisão, pelo menos para acumular resultados intermédios.
38
2 Métodos de Resolução de Sistemas de Equações
2.4.5 Solução Frontal
Esta técnica de resolução de sistemas de equações, específica para problemas estruturais e usando memória externa, tem a vantagem de ser muito eficiente, e de simples programação. Pode dizerse que é o mais simples de todos os métodos que utilizam memória auxiliar.
Resumidamente o método de solução frontal é o método de Gauss operando o sistema à medida que ele vai sendo montado. Quando todos os elementos que convergem num nó foram montados, a equação de equilibrio respectiva pode ser eliminada e gravada em disco ou banda. No fim do processo fazse a retrosubstituição a partir da ultima variável eliminada, lendo as equações de trás para a frente.
39
Capítulo 3
Métodos de Integração NuméricaPara calcular o valor de um integral definido de uma função que não tem primitiva explícita, ou cujos valores não sejam facilmente calculáveis, usamse métodos aproximados que são chamados métodos de integração numérica ou de quadratura (o nome advém de terem sido usados para calcular a área de um circulo como se fosse um quadrado).
Fazse assim a aproximação:
)f(x a dx f(x)n
0=i
b
a ii∑∫ ≅
Duas abordagens distintas podem ser seguidas para resolver este problema:
i. Fixar a posição dos xi e achar quais os valores de ai que conduzem a uma melhor aproximação para um dado tipo de funções,
ii. Achar os valores de xi e ai que conduzem a uma melhor aproximação para um dado tipo de funções.
3.1 Fórmulas de NewtonCottes
A primeira abordagem é geralmente aplicada impondo que os pontos xi sejam igualmente espaçados no intervalo [a,b] e que as fórmulas resultantes integrem polinómios do maior grau possível. Conduz aos conhecidos métodos dos Trapézios e de Simpson, quando n = 1 e 2 respectivamente. Em geral denominamse as fórmulas resultantes como fórmulas de integração de NewtonCottes.
40
3 Métodos de Integração Numérica
A determinação dos coeficientes destas fórmulas é feita integrando os polinómios de Lagrange que passam pelos pontos xi.
A função a integrar é aproximada como
∏∑≠
n
i0,=j
n
0=i )x-(x
)x-(x )f(x
ji
ji
pelo que o integral é igual a:
∏ ∫∑≠
n
i0,=j
b
a
n
0=i
dx )x-(x
)x-(x )f(x
ji
ji
Comparando esta expressão com a anteriormente dada, verificase que os pesos de integração ai
são dados por:
dx )x-(x
)x-(x = a
b
a
n
i0,=j∫ ∏
≠ ji
ji
Obtémse as seguintes quatro fórmulas quando n= 1, 2, 3 e 4:
n=1 Regra do Trapézio
][ )(fh
- )f(x+)f(x h
=dx f(x)b
aξ''
122
3
10∫
n=2 Regra de Simpson
][ )(fh
- )f(x+)f(x+)f(x h
=dx f(x)b
aξ(4)
5
210 904
3∫
n=3 NewtonCottes de terceira ordem
][ )(f3h
- )f(x+)3f(x+)3f(x+)f(x 3h
=dx f(x)b
aξ(4)
5
3210 808∫
n=4 NewtonCottes de quarta ordem
][ )(f8h
- )7f(x+)32f(x+)12f(x+)32f(x+)7f(x 2h
=dx f(x)b
aξ(6)
7
43210 94545∫
em que
b<<a ;ih +a = x; n
a-b =h ξi
41
3 Métodos de Integração Numérica
Da análise destas fórmulas pode concluirse que elas permitem a integração exacta de polinómios de grau n.
O erro associado depende de uma derivada de ordem superior a n e de uma potência da mesma ordem de h, pelo que se pode concluir que desde que a função seja regular e o valor de h pequeno se podem obter boas aproximações.
No caso de se querer integrar a função num intervalo grande, é necessário estabelecer uma fórmula de integração para um n grande. Manifestamente isso não é prático, pelo que se recorre à integração composta.
Dividese o intervalo a integrar em vários sub intervalos e calculase o integral total como a soma dos integrais nos sub intervalos.
Dividindo o intervalo [a,b] em m sub intervalos obtêmse as seguintes fórmulas:
n=1
)(fa)h-(b
- )f(x+)f(+)f( h
=dx f(x)1-m
=1j
b
aµ''
122
2
2
j
∑∫ ba
n=2
)(fa)h-(b
- )f(x+)f(x+)f(+)f( h
=dx f(x)m
=1j
1-m
=1j
b
aµ(4)
4
12j2j 18042
3
∑∑∫ ba
em que
b<<a ;jh +a = x; nm
a-b =h µj
Na prática preferese utilizar uma fórmula de ordem inferior, mais simples, utilizando um maior número de subintervalos.
3.2 Fórmula de Gauss
Quando não se impõe qual a posição dos xi, é possível determinar o valor do integral com maior precisão, calculando menos vezes o valor da função. Dado que existem 2n parâmetros que se vão fixar (os ai e os xi ) deve ser possível (e é) integrar exactamente um polinómio de grau 2n1 (que tem 2n coeficientes).
As melhores fórmulas de integração de polinómios determinamse com base nas propriedades de
42
3 Métodos de Integração Numérica
uma série de polinómios conhecidos como polinómios de Legendre. Existe um polinómio de Legendre de grau n para valores de n maiores do que zero, e ele goza, entre outras, das seguintes propriedades:
1. ≠
∫ j=i se 0 >
ji se 0 = dx (x)L (x)L
1
1- ji
2. Todas as n raizes do polinómio Ln(x) são reais e estão contidas no intervalo ]1;1[.
A função que se pretende calcular é aproximada pelos polinómios de Lagrange que passam pelas raizes xk de Ln(x) . Temse assim:
∏∑≠
≅n
i,=j
n
=i )x-(x
)x-(x )f(x f(x)
1 ji
j
1i
Como o polinómio interpolador é de ordem n1, a aproximação é exacta se f(x) fôr um polinómio de grau n1. O integral que se pretende calcular pode ser substituido por:
)f(x a dx f(x)n
=1i
1
1- ii∑∫ ≅
com
dx )x-(x
)x-(x = a
1
1-
n
i1,=j∫ ∏
≠ ji
ji
Sabese que esta fórmula é exacta para f(x)=Pn1. Será que não é para um polinómio P2n1, que seria a precisão máxima teóricamente atingível?
Dividindo P2n1 por Ln temse:
P2n1 = Qk Ln + Rm em que k < n e m < n.
Qk pode escreverse como uma combinação linear dos polinómios de Legendre de grau inferior ou igual a k:
∑k
0=i
L = Q iik α
então
∑ ∫∫k
0=i
n <k , 0 =dx L L =dx L Q poisniink α
43
3 Métodos de Integração Numérica
logo
∑∫∫ −
n
1=i
1
1)(xR a =dx (x)R =dx P
imim
1
1 12n
mas dado que os xi são os zeros de Ln
)(xR = )(xR + )(xL )(xQ = )(xP imiminiki12n
então
ii12n12n a )(xP =dx P n
=1i
1
1- ∑∫
Estas fórmulas de integração denominamse fórmulas de Gauss, sendo os xi referidos como os pontos de Gauss e os ai como os pesos de Gauss. Os seus valores encontramse extensivamente tabelados, para o caso de a=1 e b=1 (por vezes também se encontram para a=0 e b=1).
Recorrendo a uma mudança de coordenadas, calculase qualquer integral:
22com
2b+a
+a-b
= ) x( ))f(x( =dx f(x)1
1-
b
aξξξξ∫∫ − dab
Apresentamse, na Tabela 5, os pontos e pesos de Gauss para n=1 até 5.
3.3 Resultados Numéricos
Na Tabela 6, 7 e 8 indicamse vários integrais calculados pelos métodos descritos.
Interessa particularmente comparar os resultados com a solução exacta em função do número de cálculos da função necessários. Convem chamar a atenção para que:
• As fórmulas apresentadas não se comportam bem para funções com singularidades. Nesse caso devem utilizarse outros tipos de métodos de integração, que tendo em conta o tipo da singularidade extrapolam o valor do integral.
• As fórmulas apresentadas não servem no caso de domínios infinitos.
44
3 Métodos de Integração Numérica
Tabela 5 Pontos e pesos de Gauss
n xi ai
1 0.00000 2.000002 0.57735 1.00000
0.57735 1.000003 0.77460 0.55555
0.00000 0.888890.77460 0.55555
4 0.86114 0.347850.33998 0.652150.33998 0.652150.86114 0.34785
5 0.90618 0.236930.53847 0.478630.00000 0.568890.53847 0.478630.90618 0.23693
Tabela 6 Valores obtidos utilizando os métodos de integração numérica atrás descritos
dx cos(x) 1
0∫ dx cos(x) 2
0∫ 0
0.5∫ 104 × x10 dx dx x10
1
0
10×∫"exacto" 0.841471 0.909292 0.443892 0.909091
Trapézios (1) 0.770151 0.583853 2.441406 5.000000
Simpson (2) 0.841772 0.915021 0.816981 1.673177
NewtCot (3) 0.841604 0.911807 0.642136 1.315094
NewtCot (4) 0.841471 0.909263 0.478178 0.979309
Gauss (1) 0.877583 1.080605 0.004768 0.009766
Gauss (2) 0.841270 0.905451 0.227312 0.465535
Gauss (4) 0.841471 0.909297 0.442585 0.906414
Gauss (6) 0.841471 0.909298 0.443892 0.909091
Gauss (8) 0.841471 0.909297 0.443892 0.909091
45
3 Métodos de Integração Numérica
Tabela 7 Valores obtidos utilizando os métodos de integração numérica atrás descritos
dx e 2x-1.5
0∫ 1
1.5∫ ex 2
dx ∫0.5
0dx )2)-log((5x 2 ∫
2
0dx )/102)-log((5x 2
"exacto" 0.856188 0.109364 0.584112 0.320873
Trapézios (1) 0.829049 0.118320 0.000000 0.554518
Simpson (2) 0.846133 0.109310 0.191788 0.477803
NewtCot (3) 0.852270 0.109340 0.354173 0.412849
NewtCot (4) 0.856722 0.109364 0.664487 0.288723
Gauss (1) 0.854674 0.104806 0.287682 0.439445
Gauss (2) 0.863338 0.109400 1.589030 0.081094
Gauss (4) 0.856210 0.109364 0.403838 0.392983
Gauss (6) 0.856188 0.109364 0.575948 0.324139
Gauss (8) 0.856188 0.109364 0.541608 0.337874
Tabela 8 Valores obtidos utilizando os métodos de integração numérica atrás descritos
∫ ×3
1dx sin(x/10) 10/x 2
"exacto" 0.142602
Trapézios (1) 5.651953
Simpson (2) 5.080399
NewtCot (3) 3.421461
NewtCot (4) 1.197046
Gauss (1) 4.794621
Gauss (2) 2.339978
Gauss (4) 0.221287
Gauss (6) 0.159812
Gauss (8) 0.142602
46
3 Métodos de Integração Numérica
3.4 Integração em Área e em Volume
A integração em área e em volume é imediata no caso de os limites de integração serem constantes. Seja:
∫ ∫ ∑∑b
a
d
c
n
0=i
m
0=j
)y,f(x b a =dy dx y)f(x, jiji
para as fórmulas de NewtonCottes, e
∫ ∫ ∑∑1
1-
1
1-
n
1=i
m
1=j
)y,f(x b a =dy dx y)f(x, jiji
para as fórmulas de Gauss. A mudança de coordenadas necessária para as fórmulas de Gauss é igual à referida anteriormente. Caso os limites de integração não sejam constantes, pode usarse uma mudança de coordenadas que os ''endireite''. Por exemplo:
∫ ∫1
1-
x+1
1-dxdy y)f(x,
fazendo uma mudança de coordenadas tal que:
22; uv
+ u
+ v=y u =x
então:
∫ ∫∫ ∫1
1-
1
1-
1
1-
x+1
1-dvdu v)J(u, v)f(u, = dxdy y)f(x, ||
em que
vu
vu
∂∂
∂∂
∂∂
∂∂
yy
xx
= v)J(u,
Notese que devido a mudanças de coordenadas deste tipo as fórmulas de integração podem deixar de ser exactas devido à presença do termo | J | .
47
3 Métodos de Integração Numérica
3.5 Implementação da Integração Numérica
A escolha da fórmula de integração a utilizar depende do tipo de operações envolvidas. Não é possível dar receitas exaustivas, mas apresentase um exemplo para o qual a escolha pode ser bem definida. Num elemento tipo viga, de inércia variável, para calcular o efeito de uma determinada solicitação de vão, a melhor forma de integração será:
• Uma fórmula tipo NewtonCottes se as propriedades das secções estiverem definidas em posições determinadas (por exemplo décimos de vão), e se além disso interessar posteriormente conhecer os esforços provocados por essa solicitação nessas mesmas secções.
• Uma fórmula de Gauss se as propriedades forem definidas por uma função contínua, e os valores calculados a meio da viga não forem importantes, ou forem fáceis de calcular.
48
Capítulo 4
Problemas de Valores e Vectores Próprios
O problema tipo de valores e vectores próprios a considerar será:
( ) 0=− xBA λ
Duas abordagens distintas são possíveis para este problema:
• achar λ tal que
€
A −λB =0;
• achar x e λ tais que xBxA λ= .
4.1 Métodos de Determinante
A expansão do determinante de A λ B em potências de λ dá origem a um polinómio de grau igual à dimensão das matrizes A e B. Esta expansão implica um número de operações aritméticas que varia factorialmente com a dimensão do problema. O cálculo das raízes do polinómio obtido só pode ser feita (para matrizes de dimensão maior do que 3) por um processo iterativo.
Isto leva a que a expansão directa do determinante não seja numericamente eficiente, e se prefiram em geral outros métodos que não explicitam a expansão do determinante, e que a seguir se apresentam.
Se as matrizes A e B são triangulares do mesmo tipo (ambas superiores ou ambas inferiores) temse que a expansão do determinante consiste unicamente no produto dos termos da diagonal:
( ) res triangulaforem e se1
BAbaBAn
iiiii∏
=
−=− λλ
49
4 Problemas de Valores e Vectores Próprios
Se as matrizes não são triangulares (o que é o caso geral) pode factorizarse a sua soma:
UDLBAULBA =−=− λλ ou
então
UDLBAULBA =−=− λλ ou
sendo os determinantes dos dois factores determinados como se referiu anteriormente para as matrizes triangulares. Se se aplicou um dos métodos de factorização apresentados no capítulo 2, deve terse ULUL === ou ,1ou ,1 .
Este método permite calcular o determinante para um dado λ, a pesquisa dos valores de λ para os quais ele se anula tem de ser feita por um processo iterativo.
Sejam µk1 e µk duas estimativas do valor próprio de BA λ− , a estimativa seguinte será:
( )1-kk1-kk
kk1k µµ
µµµ
ηµµ −−−−
−−=+ BABA
BA
em que η é uma constante escolhida de forma a acelerar a convergência do processo ( = 1).
É necessário aplicar o factor de convergência cuidadosamente porque para valores de η superiores à unidade a convergência não é monotónica.
Para verificar se não ''se saltou'' nenhum valor próprio pode utilizarse a sequência de Sturm:
se as matrizes A e B são positivas definidas então o número de elementos negativos na diagonal de L (se Uii = 1) ou de U (se Lii = 1) é igual ao número de valores próprios inferiores a λ.
Para achar os vectores próprios resolvese a equação
( ) 0=− xBA λ
tendo em conta que é necessário incluir uma condição de normalização.
Outra alternativa é a partir de um valor de λ aproximado usar o método de potências associado a uma translação de frequência, que se apresenta mais adiante.
50
4 Problemas de Valores e Vectores Próprios
4.2 Métodos de Potências
Estes métodos operam sobre estimativas dos valores próprios do sistema, de forma a terse:
xBxA λ=
Regra geral esta equação é escrita na forma:
( )xABx 1−=λ
O vector estimativa pode ser escrito como uma combinação linear dos vectores próprios do sistema:
j
n
jj vx ∑
==
1
α
Multiplicando ambos os termos desta igualdade por ( )kAB 1− temse:
( ) ( ) j
kn
jj
kvABxAB 1
1
1 −
=
− ∑= α
Sabese que
( ) jjj vvAB λ=−1
então
( ) jkj
n
jj
kvxAB λα∑
=
− =1
1
Pondo o maior dos valores próprios em evidência fica:
( ) ( ) jk
j
n
jj
kkvxAB 1
11
1 λλαλ ∑=
− =
como
11<λλ j para j = 2, n e 111 =λλ então:
( ) 1111 limlim vxAB k
k
k
kαλ
∞→
−
∞→= .
Este limite não existe, mas a razão entre dois termos k e k + 1 da sucessão ( )kAB 1− é igual a λ1
para valores de k grandes. Ou seja, se em cada iteração se escalar (sempre da mesma maneira) o vector x, o factor k
1λ desaparece e no limite o produto háde convergir para v1, sendo o factor de
51
4 Problemas de Valores e Vectores Próprios
escala igual a λ1.
O algoritmo para calcular o primeiro vector próprio da matriz pode ser
Algoritmo ''PRIMEIRO VECTOR PRÓPRIO''
Entrada: Matriz S = B1 AVector X, primeira estimativa de v1
t tolerância admitida no erro de λ1
max número máximo de iteraçõesCritério para escalar o vector próprio
Saída: Vector X, valor calculado de v1
lambda, valor calculado de λ1.
nit = 0''Escalar (X, lambda0)''lambda0 = 0. ( A escala inicial não fazia sentido)ENQUANTO nit < max REPETIR
Y = S X ''Escalar (Y, lambda)"SE (Absoluto((lambdalambda0)/lambda) < t) ENTÃO FIMX = Ylambda0 = lambda
ERRO ''Não convergiu''FIM
Para calcular os vectores próprios seguintes, é necessário utilizar um vector tentativa para o qual o coeficiente α1 correspondente ao primeiro modo seja nulo.
Se assim fôr temse:
jk
n
jjj
kk vxAB ) ( = ) (1=
21-
2 ∑ λλαλ
como α1 =0 não importa que o termo (λ1 / λ2) não tenda para zero.
Sabese que os vectores próprios são ortogonais em relação à matriz B.
jivBv ≠ se 0 = jiT
52
4 Problemas de Valores e Vectores Próprios
A estimativa inicial x pode, como se disse anteriormente, escreverse:
∑n
vx1=
jj = j
α
então
11111=
1 = = vBvvBvxBv TTTj
n
jj αα∑
11
11 T
T
= vBv
xBvα
Pelo que
x
xSvBv
BvvI
vxx
)
=
(=
-=
1
11
11
111
T
T
−
α
é uma estimativa para a qual a componente segundo v1 é nula.
Teoricamente bastaria multiplicar a estimativa inicial do segundo modo pela matriz S1 para fazer desaparecer o efeito do primeiro modo. No entanto, devido a erros numéricos acumulados, esta componente pode voltar a aparecer.
Multiplicando em cada iteração o vector estimativa por S1, conseguese controlar esse erro numérico. Na prática, costuma operarse a matriz B1 A já multiplicada por S1.
Em vez de fazer:
x1 = S1 x , y = S x1
calculase uma única vez D = S S1, e depois iterase
y = D x
Para calcular valores e vectores próprios superiores ao segundo utilizase uma matriz Sm que elimina as componentes segundo os m primeiros vectores próprios.
Para valores de m grandes, os erros numéricos tomam proporções exageradas e o processo iterativo deixa de convergir.
Para calcular os menores valores próprios pode iterarse na forma:
53
4 Problemas de Valores e Vectores Próprios
y = ( A1 B ) x
Como
j
jj
1- = λv
vBA
então
€
( A1 B )k x = 1λm
k αj ( λm λjj =1
n
∑ )k v j
e o processo iterativo converge para o valor pretendido pois ( ) 0→kjm λλ para j < m.
Para calcular valores próprios intermédios é costume utilizar um artifício conhecido por translação de frequência.
Se se escrever o valor próprio kk δµλ += , a equação do sistema fica:
( ) 0=−− xBBA kδµ .
O sistema pode ser redefinido chamando
BBBAA ='=' e µ− .
( ) 0'' =− xBA kδ .
Formalmente o sistema é identico ao anterior, mas os valores próprios δk diferem dos λk da constante µ.
Seja λm o valor próprio mais próximo de µ, então
jk
n
jjmjk
m
k vxBA ) ( 1
= ) '' (1=
1- ∑ δδαδ
pelo que o processo iterativo converge para δm e vm. λm é calculado fazendo:
mm δµλ +=
Notese que é necessário que µ não corresponda a nenhum valor próprio, porque nesse caso 0' =−= BAA µ e não existe A'1.
54