Baixe cap4 e outras Provas em PDF para Engenharia Civil, somente na Docsity!
Cap´ıtulo 4
Resolu¸c˜ao de Sistemas de
Equa¸c˜oes Lineares
4.1 Introdu¸c˜ao
A resolu¸c˜ao de sistemas de equa¸c˜oes lineares ´e um dos problemas num´ericos mais comuns em aplica¸c˜oes cient´ıficas. Tais sistemas surgem, por exemplo, em conex˜ao com a solu¸c˜ao de equa¸c˜oes diferenciais parciais, determina¸c˜ao de caminhos ´otimos em redes (grafos) e interpola¸c˜ao de pontos, dentre outros. Consideraremos aqui, inicialmente, a resolu¸c˜ao de um sistema de equa¸c˜oes lineares de n equa¸c˜oes a n vari´aveis (inc´ognitas),
a 11 x 1 + a 12 x 2 +... + a 1 nxn = b 1 a 21 x 1 + a 22 x 2 +... + a 2 nxn = b 2
.. . an 1 x 1 + an 2 x 2 +... + annxn = bn
ou, escrito na forma matricial, Ax = b (4.1)
onde A ´e uma matriz quadrada, de ordem n, e x e b s˜ao vetores de n elementos,
a 11 a 12... a 1 n a 21 a 22... a 2 n .. .
an 1 an 2... ann
x 1 x 2 .. . xn
b 1 b 2 .. . bn
A matriz A pode apresentar, dependendo do problema de onde o sistema foi derivado, uma certa estrutura e esparsidade. Uma matriz ´e dita estruturada se os seus elementos est˜ao dispostos de uma determinada forma como, por exemplo, ao longo de algumas diagonais e/ou colunas/linhas (figuras 4.1-a) e 4.1-b), como um triˆangulo (a matriz em 4.1-c ´e dita triangular inferior) ou, ainda, sem estrutura qualquer (4.1-d). Al´em disso, as matrizes mostradas na figura 4.1 apresentam alguns elementos nulos. Uma matriz ´e dita esparsa se ela cont´em, aproximadamente, em torno de 90% de elementos nulos; caso contr´ario, ela ´e dita densa. Em conseq¨uˆencia, pode-se dizer que um sistema ´e esparso ou denso, dependendo de como ´e a matriz de coeficientes do sistema. Uma das principais metas a se atingir, na resolu¸c˜ao de um sistema de equa¸c˜oes lineares, ´e obter a sua solu¸c˜ao no menor espa¸co de tempo e, se poss´ıvel, sem alterar a sua estrutura e/ou esparsidade. Por isso, existem certos m´etodos e/ou algoritmos espec´ıficos para se resolver alguns sistemas particulares, conforme veremos a seguir.
(a)
× ×
× × ×
× × ×
× ×
,^ (b)
× ×
× ×
× ×
× × × ×
(c)
×
× ×
× × ×
× × × ×
,^ (d)
× ×
× ×
×
× ×
Figura 4.1: Estruturas t´ıpicas de matrizes: (a)tridiagonal, (b)flecha, (c)triangular inferior, (d) n˜ao-estruturada.
4.2 Resolu¸c˜ao de Sistemas Triangulares de Equa¸c˜oes Line-
ares
Se o sistema (4.1) apresenta sua matriz de coeficientes A na forma triangular – seja ela inferior, como mostrado na figura 4.1-c, ou superior – ent˜ao ´e poss´ıvel resolvˆe-lo de forma imediata, atrav´es de substitui¸c˜ao direta, para matrizes triangulares inferiores, e de retro-substitui¸c˜ao, para matrizes triangulares superiores. Suponha ent˜ao um sistema triangular inferior,
Lx = b (4.2)
onde
L =
l 11 l 21 l 22 l 31 l 32 l 33 .. .
ln 1 ln 2 ln 3... lnn
Nesse caso, as inc´ognitas x 1 , x 2 ,.. ., xn, podem ser facilmente determinadas como
x 1 =
b 1 l 11
x 2 =
b 2 − l 21 x 1 l 22
x 3 =
b 3 − l 31 x 1 − l 32 x 2 l 33
.. .
xn =
bn −
∑n− 1 j=1 lnj^ xj lnn
O processo acima, denominado de substitui¸c˜ao direta, pode ser expresso de forma algor´ıtmica como
4.3 Resolu¸c˜ao de Sistemas de Equa¸c˜oes Lineares por Elimi-
na¸c˜ao Gaussiana
Se o sistema de equa¸c˜oes lineares n˜ao apresenta uma forma simples, tal que se possa determinar as inc´ognitas facilmente, ent˜ao podemos efetuar modifica¸c˜oes no sistema de tal forma que o transformamos em um sistema triangular, preservando a solu¸c˜ao do sistema anterior. Uma vez feitas estas modifica¸c˜oes, a solu¸c˜ao ´e obtida de forma imediata, conforme visto na se¸c˜ao anterior. Um processo desse tipo ´e aquele chamado de elimina¸c˜ao Gaussiana, o qual consiste em se aplicar opera¸c˜oes elementares – somas e multiplica¸c˜oes – `as linhas da matriz de coeficientes e do vetor independente b, de tal forma que a matriz passe a ser triangular superior. Suponha, por exemplo, o sistema (^)
x 1 x 2 x 3
cuja solu¸c˜ao ´e x 1 = x 2 = x 3 = 1. Para transformarmos a matriz A em uma matriz triangular superior, devemos eliminar os elementos abaixo da diagonal principal de A. Para tanto, se multiplicamos a primeira linha por a 21 /a 11 = − 1 /4 e subtra´ımo-la da segunda, temos: (^)
x 1 x 2 x 3
Agora, para eliminar o termo a 31 , multiplicamos a primeira linha por a 31 /a 11 = 1/1 e subtra´ımo-la da terceira: (^)
x 1 x 2 x 3
Note que os elementos do vetor independente b s˜ao modificados tamb´em! A matriz agora ´e praticamente triangular superior; falta eliminar o termo a 32. Para tanto, basta multiplicar a segunda linha por a 32 /a 22 = − 2 / 7 , 5 e subtra´ı-la da terceira, de onde
x 1 x 2 x 3
Agora, podemos utilizar o algoritmo de retro-substitui¸c˜ao para determinar as inc´ognitas:
x 3 = 1
x 2 =
x 1 =
Podemos sumarizar o processo ent˜ao da seguinte forma: para se eliminar os elementos abaixo da diagonal na k-´esima coluna (ou seja, os elementos das linhas k + 1, k + 2,.. ., n na coluna k), usamos o elemento akk – chamado de pivˆo – para calcularmos um multiplicador z = a aikkk para cada i-´esima linha abaixo da linha k. Esse multiplicador ser´a utilizado para multiplicar a k-´esima linha e subtra´ı-la da linha i (incluindo, aqui, os elementos do termo independente b). Uma vez eliminados todos os elementos abaixo da diagonal principal de A, resta-nos uma matriz triangular superior, e, ent˜ao, podemos determinar a solu¸c˜ao x usando o algoritmo da retro-substitui¸c˜ao. O processo de elimina¸c˜ao Gaussiana pode ser descrito de forma algoritmica como
Algoritmo 4.3.1 Elimina¸c˜ao Gaussiana
proc elimina¸c˜ao Gaussiana(input: A, b; output: x) for k = 1, 2 ,... , n − 1 do for i = k + 1, k + 2,... , n do z ← a akkik aik ← 0 for j = k + 1, k + 2,... , n do aij ← aij − zakj endfor bi ← bi − zbk endfor endfor call retro substitui¸c˜ao(A, b, x) endproc
4.3.1 Dificuldades
O processo de elimina¸c˜ao Gaussiana, descrito acima, n˜ao consegue resolver todo e qualquer sistema. Considere, por exemplo, o sistema [ 0 1 1 1
] [
x 1 x 2
]
[
]
o qual tem como solu¸c˜ao x 1 = x 2 = 1. No entanto, se formos aplicar elimina¸c˜ao Gaussiana a esse sistema, ele falhar´a, pois o pivˆo a 11 = 0. E ´´obvio, portanto, que os pivˆos n˜ao podem ser nulos. O sistema (4.4) pode, no entanto, ser modificado, procedendo-se a uma troca de linhas - imediatamente temos um sistema triangular superior. No entanto, ´e poss´ıvel que, ao longo do processo de elimina¸c˜ao Gaussiana, surja um zero na diagonal principal e n˜ao seja poss´ıvel, por qualquer troca de linhas, removˆe-lo. Nesse caso, o sistema n˜ao tem solu¸c˜ao^1 ; o algoritmo para a elimina¸c˜ao Gaussiana deve ser modificado adequadamente para se levar em conta tal possibilidade. O pr´oximo exemplo mostra uma outra dificuldade associada ao m´etodo: [ ε 1 1 1
] [
x 1 x 2
]
[
]
onde 0 < ε 1, cuja solu¸c˜ao correta ´e
x 1 = (^1) −^1 ε ≈ 1 x 2 = 11 −−^2 εε ≈ 1
No entanto, se aplicarmos elimina¸c˜ao Gaussiana ao sistema (4.5), obteremos
x 2 = 2 −ε
− 1 1 −ε−^1 ≈^1 x 1 = (1 − x 2 )ε−^1 ≈ 0
o qual obviamente aproxima bem x 2 , mas o valor de x 1 ´e completamente errado! Isso acontece porque, se ε ´e pequeno o suficiente em um determinado computador, tanto 2 − ε−^1 quanto 1 − ε−^1 ser˜ao calculados como −ε−^1 (devido `a perda de d´ıgitos significativos na subtra¸c˜ao). Desse exemplo, tiramos uma outra li¸c˜ao: o pivˆo deve, sempre, ser escolhido como o maior poss´ıvel, em m´odulo. O processo de escolha de pivˆos, chamado de pivotamento, implica na troca de linhas da matriz de coeficientes (bem como do termo independente b). Computacionalmente, no entanto, n˜ao ´e
(^1) Esta ´e, inclusive, uma maneira de se determinar se o sistema ´e singular, isto ´e, a matriz de coeficientes n˜ao tem inversa.
Note que os fatores apik apk k utilizados para se eliminar os elementos abaixo da diagonal s˜ao armazenados na matriz A, onde se colocariam zeros (conforme utilizado no algoritmo da elimina¸c˜ao Gaussiana sem pivotamento). Isso ´e feito de forma a se poder obter, a partir do algoritmo acima, a fatora¸c˜ao LU da matriz A, conforme veremos na se¸c˜ao a seguir. O exemplo abaixo mostra o funcionamento do algoritmo descrito acima:
Exemplo 4.1 Calcule a solu¸c˜ao do sistema
x 1 x 2 x 3
Solu¸c˜ao: Inicialmente, temos p = (1, 2 , 3) (de acordo com o algoritmo)e s = (6, 8 , 3) (verifique, por inspe¸c˜ao). A cada passo, temos:
k j p i z 1 3 (3, 2 , 1) 2 0 , 3333 3 0 , 6667
A =
(^) b =
k j p i z 2 3 (3, 1 , 2) 3 − 1 , 2308
A =
(^) b =
Uma vez efetuada a elimina¸c˜ao, procede-se ao c´alculo das inc´ognitas:
i = 3 : p 3 = 2, x 3 =
i = 2 : p 2 = 1, x 2 =
i = 1 : p 1 = 3, x 1 =
4.3.2 Elimina¸c˜ao Gaussiana e a Fatora¸c˜ao LU
Conforme visto na se¸c˜ao anterior, o algoritmo de elimina¸c˜ao Gaussiana com pivotamento e es- calonamento produz, de forma impl´ıcita, uma matriz triangular inferior, uma matriz triangular superior e um vetor de permuta¸c˜ao. Como essas matrizes foram obtidas por transforma¸c˜oes sobre a matriz A original, podemos de alguma forma relacion´a-las entre si. Primeiramente, analisemos o vetor de permuta¸c˜ao; seus elementos indicam qual linha foi trocada com outra, i.e., se pj = k, isso significa que a linha j foi trocada com a linha k. Essa permuta¸c˜ao pode ser expressa, tamb´em, atrav´es de uma matriz de permuta¸c˜ao, P , a qual tem como elementos apenas o 0 e o 1. No exemplo mostrado na se¸c˜ao anterior, obtemos p = (3, 1 , 2) ao fim; ou seja, a linha 3 est´a no lugar da linha 1; a linha 1 est´a no lugar da linha 2 e, por fim, a linha 2 est´a na linha 3. A matriz de permuta¸c˜ao correspondente ´e
P =
A rela¸c˜ao existente entre A e as matrizes triangular inferior, L; triangular superior, U ; e a matriz de permuta¸c˜ao P ´e a seguinte: P A = LU
onde L ´e triangular inferior com diagonal unit´aria e os seus elementos abaixo da diagonal principal encontram-se armazenados na matriz A, ao final do algoritmo de elimina¸c˜ao Gaussiana com pivotamento e escalonamento, por´em possivelmente permutados. Usando mais uma vez o exemplo anterior, temos:
^ P A^ =^ LU
onde as matrizes L e U foram permutadas adequadamente, usando a matriz P. Pode-se verificar, por inspe¸c˜ao, que o lado direito da igualdade ´e a matriz A com as suas linhas trocadas conforme expresso por P. A fatora¸c˜ao LU ´e ´util quando, para uma mesma matriz de coeficientes A, temos de resolver m sistemas de equa¸c˜oes lineares Ax(j)^ = b(j), com termos independentes b(1), b(2),.. ., b(m). Basta, ent˜ao, obter a fatora¸c˜ao com o algoritmo de elimina¸c˜ao Gaussiana com pivotamento e escalonamento (sem calcular xi – as ´ultimas trˆes linhas do algoritmo), obtendo L, U e P. Valendo- se da igualdade P Ax = P b e, como P A = LU , podemos escrever L(U x) = b, de onde a solu¸c˜ao de um sistema Ax = b ´e obtida resolvendo-se dois sistemas triangulares:
Ly = P b U x = y
A fatora¸c˜ao LU , bem como a solu¸c˜ao dos sistemas triangulares acima, s˜ao expressas pelos algoritmos 4.3.3 e 4.3.4.
Como n − 1 linhas s˜ao processadas dessa forma, temos um total de n(n − 1) ≈ n^2 opera¸c˜oes para a primeira coluna. Para as demais colunas, note que o mesmo racioc´ınio acima ´e v´alido, mas ´e como se a matriz diminu´ısse de uma linha e uma coluna a cada novo valor de k. Assim, para todos os n − 1 pivˆos a serem calculados, teremos:
n^2 + (n − 1)^2 +... + 3^2 + 2^2 =
n^3 +
n^2 +
n − 1 ≈
n^3 +
n^2
a qual ´e obtida usando
∑n k=1 k
6 n(n^ + 1)(2n^ + 1). Para se corrigir o termo independente b, gasta-se n − 1 opera¸c˜oes, depois n − 2, e assim sucessivamente, de onde
(n − 1) + (n − 2) +... + 1 =
n^2 −
n.
Finalmente, o processo de retro-substitui¸c˜ao custa
1 + 2 + 3 +... + n =
n^2 +
n
opera¸c˜oes. Combinando todas as express˜oes, podemos dizer que, para se resolver m sistemas de equa¸c˜oes lineares Ax(i)^ = b(i), usando a fatora¸c˜ao LU , apresenta um custo computacional de aproximada- mente 1 3 n^3 +
n^2
o que mostra que ´e mais eficiente efetuar a fatora¸c˜ao LU apenas uma vez, e depois resolver os m sistemas lineares, do que se resolvˆessemos cada sistema independentemente, pois o custo, nesse caso, seria da ordem de 13 mn^3.
4.3.4 Resolu¸c˜ao de sistemas com m´ultiplos termos independentes
Existem situa¸c˜oes que requerem a solu¸c˜ao de v´arios sistemas lineares, todos eles com a mesma matriz de coeficientes, por´em com diferentes termos independentes. Como visto na se¸c˜ao 4.3.3, ´e mais vantajoso, nesse caso, realizar-se a fatora¸c˜ao LU de A, apenas uma vez; a solu¸c˜ao de todos os sistemas ´e obtida, simplesmente, calculando-se as solu¸c˜oes dos sistemas triangulares Ly(i)^ = P b(i) e U x(i)^ = y(i), onde o ´ındice (i) identifica um sistema espec´ıfico.
4.3.4.1 C´alculo da inversa de uma matriz
Uma dessas situa¸c˜oes ´e o c´alculo da inversa de uma matriz. Note que tal c´alculo n˜ao ´e realizado com o fim de se resolver um sistema de equa¸c˜oes (utilizando-se a rela¸c˜ao x = A−^1 b; aplica¸c˜oes que envolvam certas decomposi¸c˜oes de matrizes exigem que se escreva um vetor v como XDX−^1 , onde X e D s˜ao matrizes. Seja ent˜ao a matriz A, cuja inversa A−^1 ´e desejada. Como, por defini¸c˜ao, o produto entre uma matriz e a sua inversa ´e a matriz identidade I,
AA−^1 = I (4.6)
podemos escrever o problema de determina¸c˜ao da inversa na forma de um sistema de equa¸c˜oes lineares com m´ultiplos termos independentes (e, conseq¨uentemente, m´ultiplas solu¸c˜oes) como
AX = I (4.7)
onde X ≡ A−^1 , i.e., as colunas (X)i s˜ao as colunas de A−^1.
Dessa forma, obtendo-se a fatora¸c˜ao LU de A, a primeira coluna de A−^1 ´e obtida resolvendo-se os sistemas
Ly(i)^ = P
e U (X) 1 = y(1)^ (4.9)
e a segunda coluna ´e obtida como
Ly(2)^ = P
e U (X) 2 = y(2)^ (4.11)
e as demais colunas s˜ao obtidas similarmente. Note que, computacionalmente, basta usar apenas um vetor y, sendo o mesmo reutilizado a cada novo sistema resolvido.
Exemplo 4.2 Obtenha a inversa da matriz
A =
Solu¸c˜ao: Aplicando-se o algoritmo 4.3.3, obtemos os fatores L, U e P :
L =
U =
P =
Agora, aplica-se o algoritmo 4.3.4 usando-se como termo independente o vetor (1, 0 ,... , 0)T^ , i.e., resolve-se
Ly(1)^ = P
(^) y(1)^ =
y(1)^ =
Assim, A−^1 ´e dada por
A−^1 = X =
e pode-se verificar que
^ AA−^1 =^ I
4.4 Resolu¸c˜ao Iterativa de Sistemas de Equa¸c˜oes Lineares
Em certos casos, n˜ao ´e conveniente se resolver o sistema Ax = b atrav´es de um m´etodo direto como a elimina¸c˜ao Gaussiana. Considere, por exemplo, a matriz A derivada da discretiza¸c˜ao em diferen¸cas-finitas (com estˆencil de 5 pontos) do operador diferencial ∇^2 , cuja estrutura ´e mostrada na figura 4.2; se aplicarmos a fatora¸c˜ao LU sobre A, alguns dos elementos que eram nulos em A passar˜ao a ser diferentes de zero, tanto em L como em U (figura 4.3). Note que A tem 64 elementos n˜ao-nulos, ao passo que L e U apresentam um total de 134 elementos n˜ao-nulos. A
0 2 4 6 8 10 12 14 16
0 2 4 6 8 10 12 14 16 nz = 64
Estrutura da matriz A
Figura 4.2: Estrutura da matriz A derivada da discretiza¸c˜ao em diferen¸cas-finitas do operador diferencial ∇^2.
elimina¸c˜ao Gaussiana est´a, nesse caso, destruindo a estrutura e/ou a esparsidade da matriz, o que n˜ao ´e aconselh´avel, principalmente para matrizes grandes (n > 10000). Por outro lado, mesmo quando a matriz ´e densa – ou seja, a inser¸c˜ao de elementos n˜ao-nulos n˜ao implicar´a em aumento consider´avel do uso da mem´oria – pode n˜ao ser aconselh´avel utilizar um m´etodo direto, se a solu¸c˜ao desejada necessita apenas um n´umero pequeno de d´ıgitos corretos.
0 2 4 6 8 10 12 14 16
0 2 4 6 8 10 12 14 16 nz = 67
Estrutura da matriz L
0 2 4 6 8 10 12 14 16
0 2 4 6 8 10 12 14 16 nz = 67
Estrutura da matriz U
Figura 4.3: Estruturas da matriz L (a esquerda)e U (a direita), resultantes da fatora¸c˜ao LU de A.
Uma outra raz˜ao, que justifica o uso de m´etodos iterativos (ver [13]), ´e o fato de seu custo computacional ser proporcional a n^2 (e, `as vezes, at´e mesmo a n), o que os torna bastante competitivos, se comparados a um m´etodo direto (cujo custo ´e proporcional a n^3 ).
4.4.1 Normas de vetores e de matrizes
Como todo processo iterativo, ´e necess´ario saber quando se alcan¸cou a convergˆencia do processo
- em nosso caso, obteve-se uma estimativa xk que aproxima suficientemente x = A−^1 b. Fazendo uma analogia com o m´etodo da bissec¸c˜ao (ver se¸c˜ao 2.2), onde se detectava a convergˆencia quando o comprimento do intervalo era menor do que uma tolerˆancia pr´e-especificada, aqui vamos tamb´em calcular um comprimento de um vetor (em IRn). Para se calcular esse comprimento, utiliza-se uma norma. Uma norma de um vetor x per- tencente a um espa¸co vetorial V ´e uma fun¸c˜ao || x || : V −→ IR+^ que obedece aos seguintes postulados: || x || > 0 , se x = 0, x ∈ V || λx || = | λ | || x ||, se λ ∈ IR, x ∈ V || x + y || ≤ || x || + || y || se x, y ∈ V (desigualdade triangular) A norma de um vetor ´e o seu ”comprimento” no espa¸co vetorial V ; ´e uma generaliza¸c˜ao da no¸c˜ao de valor absoluto de um n´umero real. Para o espa¸co vetorial IRn, a norma mais conhecida ´e a chamada norma Euclidiana, definida por
|| x || 2 =
( (^) n ∑
i=
x^2 i
onde x = (x 1 , x 2 ,... , xn)T^. Particularmente, em IR^2 , temos || x || 2 =
x^21 + x^22 , que ´e a express˜ao para a distˆancia de um ponto com coordenadas (x 1 , x 2 ) em rela¸c˜ao `a origem do sistema de eixos cartesianos. Existem outras normas que s˜ao bastante usadas em c´alculos num´ericos, como a norma-l∞
|| x ||∞ = n max i=
| xi | (4.13)
e a norma-l 1 ,
|| x || 1 =
∑^ n
i=
| xi | (4.14)
as quais s˜ao bem mais simples e menos onerosas de se calcular do que a norma Euclidiana.
4.4.2 Normas de matrizes
Uma vez especificada uma norma de um vetor, a norma matricial subordinada ´e definida como
|| A || = sup || Au || : u ∈ IRn, || u || = 1 (4.15)
para uma matriz A n × n. Pode-se verificar que
|| Ax || ≤ || A || || x ||, x ∈ IRn
Por exemplo, a norma matricial subordinada da norma vetorial || · ||∞ ´e dada por
|| A ||∞ =
n max i=
∑^ n
j=
| aij | (4.16)
4.4.4 Erros computacionais e condicionamento
Qualquer solu¸c˜ao de um sistema linear deve ser considerada uma solu¸c˜ao aproximada, em virtude de erros de arredondamento e outros. O m´etodo mais natural para determina¸c˜ao da precis˜ao de uma solu¸c˜ao ´e verificar qu˜ao bem esta solu¸c˜ao satisfaz o sistema original, calculando o vetor res´ıduo. Se a solu¸c˜ao aproximada ˜x for uma boa aproxima¸c˜ao, pode-se esperar que cada componente de r = b − Ax˜ seja pequeno, pelo menos em um conceito relativo. H´a sistemas de equa¸c˜oes, contudo, em que o resto n˜ao proporciona uma boa medida da precis˜ao da solu¸c˜ao. S˜ao sistemas nos quais pequenas altera¸c˜oes nos dados de entrada conduzem a mudan¸cas significativas na solu¸c˜ao. Estes s˜ao denominados sistemas inst´aveis ou mal-condicionados.
Exemplo 4.4 A solu¸c˜ao exata do sistema { x 1 + x 2 = 2 1 , 01 x 1 + x 2 = 2, 01
´e x 1 = x 2 = 1. Supondo que, devido a erros, a solu¸c˜ao calculada fosse
x 1 = 0 x 2 = 2 , 005
o vetor res´ıduo neste caso seria RT^ = [− 0 , 005; 0, 005]. Entretanto, o erro em cada resposta, x 1 e x 2 , ´e de aproximadamente uma unidade. Por outro lado, os coeficientes tamb´em podem conter erros. Supondo que algum tipo de erro tenha mudado as equa¸c˜oes acima para
{ x 1 + x 2 = 2 1 , 0001 x 1 + x 2 = 2, 007
at´e mesmo uma solu¸c˜ao bem diferente da anterior, como x 1 = 100 e x 2 = − 98 produziria um res´ıduo bem pequeno, RT^ = [0; − 0 , 003].
Erros deste tipo, ao contr´ario daqueles causados pela acumula¸c˜ao de erros de arredondamento, n˜ao podem ser evitados por uma programa¸c˜ao cuidadosa. Como, ent˜ao, determinar quando um problema ´e mal-condicionado?
Figura 4.4: Os sistemas que descrevem a intersec¸c˜ao das retas s˜ao, da esquerda para a direita: bem-condicionado, mal-condicionado e singular.
Em geral, tem-se a situa¸c˜ao mostrada na figura 4.4, para o caso de duas retas. Quando o sistema ´e ordem maior, no entanto, deve-se recorrer a medidas algebricas para se estimar o mal- condicionamento do sistema. Uma dessas medidas ´e o chamado determinante normalizado da matriz dos coeficientes. Para obtˆe-lo, normaliza-se a matriz de coeficientes A dividindo-se cada
linha de A pela raiz quadrada da soma dos quadrados dos elementos de cada linha,
norm| A | =
a 11 α 1
a 12 α 1 · · ·^
a 1 n a^ α^1 21 α 2
a 22 α 2 · · ·^
a 2 n α 2 .. .
an 1 αn
an 2 αn · · ·^
ann αn
| A |
α 1 α 2 ·... · αn
onde αi =
a^2 i 1 + a^2 i 2 +... + a^2 in. Diz-se, ent˜ao, que uma matriz A ´e mal-condicionada se o n´umero norm| A | for pequeno, comparado com a unidade.
Exemplo 4.5 Seja a matriz
A =
[
]
verifique se ela ´e mal-condicionada. Solu¸c˜ao: Calculando-se α 1 =
1 e α 2 =
2 , 0201 , podemos obter
norm| A | =
α 1 α 2
ou seja, A ´e dita ser mal-condicionada.
4.4.5 M´etodos iterativos
Dado um sistema n˜ao-singular de n equa¸c˜oes lineares Ax = b, um m´etodo iterativo para resolver esse sistema ´e definido por um conjunto de fun¸c˜oes Φk(x 0 , x 1 ,... , xk, A, b), onde x 0 = Φ 0 (A, b) ´e uma estimativa inicial para a solu¸c˜ao x = A−^1 b e x 1 , x 2 ,... s˜ao as aproxima¸c˜oes sucessivas para a solu¸c˜ao,
x 1 = Φ 1 (x 0 , A, b) x 2 = Φ 2 (x 0 , x 1 , A, b) .. . xk = Φk(x 0 , x 1 ,... , xk, A, b)
As fun¸c˜oes Φk nos definem os m´etodos iterativos. Diz-se que um m´etodo ´e estacion´ario se, para um m > 0, Φn n˜ao depende de n para todo n ≥ m, ou seja, Φ = Φm = Φm+1 =... Nesse caso, xn+1 depende de, no m´aximo, m vetores anteriores, xn, xn− 1 ,.. ., xn−m+1. Por exemplo, para m = 2, temos
x 0 = Φ 0 (A, b) (4.19) x 1 = Φ 1 (x 0 , A, b) (4.20) xk = Φ(xk− 2 , xk− 1 , A, b), k = 2, 3 ,... (4.21)
O grau de um m´etodo estacion´ario ´e ˆm (para ˆm ≤ m) se, para n ≥ m − 1, xn+1 depende de xn, xn− 1 ,.. ., xn− mˆ+1 mas n˜ao para k < n − mˆ+. O grau de um m´etodo iterativo definido pelas equa¸c˜oes (4.19)-(4.21) ´e 2. Um m´etodo iterativo ´e dito linear se todas as fun¸c˜oes Φi s˜ao fun¸c˜oes lineares de x 0 , x 1 ,.. ., xn− 1. Assim, um m´etodo iterativo estacion´ario linear de grau 1 pode ser expresso por
xk+1 = Gxk + f (4.22)
onde G ´e uma matriz e f um vetor, escolhidos adequadamente.
Multiplicando r 0 por A−^1 , temos
A−^1 r 0 = A−^1 b − x 0 = x − x 0 = e 0 (4.34)
e, usando essa igualdade, podemos obter uma express˜ao para x:
x = x 0 + A−^1 r 0 (4.35)
Note que a equa¸c˜ao (4.35) envolve A−^1 ; mas, obviamente, n˜ao podemos utiliz´a-la, pois se a calcul´assemos, a solu¸c˜ao do sistema seria imediata! Por outro lado, a equa¸c˜ao (4.34) nos permite escrever Ae 0 = r 0 (4.36)
e, combinando as equa¸c˜oes (4.33), (4.36) e (4.35), podemos descrever o m´etodo do refinamento iterativo como (^)
rk = b − Axk resolve Aek = rk xk+1 = xk + ek
, k = 0, 1 ,... (4.37)
O m´etodo do refinamento iterativo ´e utilizado em conjunto com um m´etodo direto, como a elimina¸c˜ao Gaussiana. Tendo fatorado A no produto LU e obtido uma solu¸c˜ao ˜x para Ax = b (a qual pode n˜ao ser muito boa, devido a erros de arredondamento), fazemos x 0 = ˜x e refinamos essa solu¸c˜ao, usando (4.37), at´e que xk seja suficientemente bom. Note que a fatora¸c˜ao LU pode, agora, ser utilizada para resolver Aek = (LU )ek = rk. Se consideramos que a solu¸c˜ao obtida com a fatora¸c˜ao LU de A n˜ao foi exata, ent˜ao podemos dizer que U −^1 L−^1 = B ≈ A−^1. Usando a equa¸c˜ao (4.35), escrevemos
xk+1 = xk + A−^1 rk = xk + A−^1 b − A−^1 Axk .. B. ≈ A−^1 ... = xk + B(b − Axk ) (4.38)
De onde podemos mostrar que o m´etodo converge para uma solu¸c˜ao: subtraindo x de ambos os lados da equa¸c˜ao (4.38), temos
xk+1 − x = xk − x + B(b − Axk ) .. b. = Ax ... = xk − x + B(Ax − Axk ) = (I − BA)(xk − x)
e, tomando normas de ambos os lados da igualdade acima, vem, pela desigualdade triangular:
|| xk+1 − x || ≤ || I − BA || ||xk − x || ... || xk − x || = || I − BA || || xk− 1 − x || ... ≤ || I − BA ||^2 ||xk− 1 − x ||
≤
≤ || I − BA ||k^ ||x 0 − x ||
o que nos diz que os erros convergem para 0 se || I − BA || < 1. Assim como nos m´etodos de determina¸c˜ao de ra´ızes de fun¸c˜oes, precisamos definir alguns crit´erios de parada do processo de refinamento. O primeiro desses crit´erios ´e a estipula¸c˜ao de um n´umero m´aximo de itera¸c˜oes (kmax); o segundo pode ser baseado na norma do res´ıduo rk, devidamente escalonada por || b || (usando uma norma qualquer, previamente escolhida). Assim, as itera¸c˜oes proceder˜ao enquanto || rk || < ε|| b || (4.39)
n˜ao for satisfeito; ε ´e um n´umero real, escolhido de acordo com a exatid˜ao requerida. Um algoritmo que expressa o m´etodo do refinamento iterativo pode ser escrito como:
Algoritmo 4.4.1 Refinamento Iterativo
proc refinamento iterativo(input: A, L, U , b, x 0 , kmax, ε; output: xk) t ← ε|| b || for k = 0, 1 ,... , kmax do rk ← b − Axk if || rk || < t then break endif resolve (LU )ek = rk, obtendo ek xk+1 ← xk + ek endfor endproc
O exemplo abaixo [10] mostra o comportamento desse m´etodo.
Exemplo 4.6 Seja o sistema
x^ =
cuja solu¸c˜ao ´e o vetor x = (1, 1 , 1 , 1)T^. Utilizando um computador com apenas 6 casas decimais de precis˜ao, obtemos como solu¸c˜ao inicial, atrav´es da elimina¸c˜ao Gaussiana, com pivotamento, o vetor x = (0, 999988 , 1 , 000137 , 0 , 999670 , 1 , 000215)T
Agora, dispondo dos fatores triangulares da fatora¸c˜ao LU , podemos utilizar o algoritmo 4.4.1 e obter: x = (0, 999994 , 1 , 000069 , 0 , 999831 , 1 , 000110)T x = (0, 999996 , 1 , 000046 , 0 , 999891 , 1 , 000070)T x = (0, 999993 , 1 , 000080 , 0 , 999812 , 1 , 000121)T x = (1, 000000 , 1 , 000006 , 0 , 999984 , 1 , 000011)T
4.4.7 M´etodo iterativo de Jacobi
Suponha o sistema (4.1), com n = 3, sem perda de generalidade. Se os elementos da diagonal de A s˜ao todos n˜ao-nulos, ent˜ao pode-se isolar cada vari´avel x 1 , x 2 e x 3 atrav´es de
x 1 = c 12 x 2 + c 13 x 3 + d 1 x 2 = c 21 x 1 + c 23 x 3 + d 2 x 3 = c 31 x 1 + c 32 x 2 + d 3
onde
cij =
− a aijii , i = j 0 , i = j di = bi/aii
Com essa transforma¸c˜ao, o sistema Ax = b foi transformado em um sistema da forma
(I − C)x = d
