Capı́tulo 5

Matrizes e Sistemas lineares

Neste capı́tulo estudaremos alguns métodos para calcular a solução de sistemas de equações
lineares. Apenas nos preocuparemos com sistemas quadrados, isto é, aqueles em que o
número de equações é igual ao número de incógnitas. Supõe-se que as noções básicas de
álgebra matricial, como adição e multiplicação de matrizes, matriz inversa e identidade,
determinante de uma matriz etc., sejam conhecidas do leitor.

5.1 Introdução
Um sistema de equações algébricas de ordem n, que é um conjunto de n equações com n
incógnitas,
8
>
>a11 x1 + a12 x2 + . . . + a1n xn = b1
>
>
<a21 x1 + a22 x2 + . . . + a2n xn = b2
> ..
>
> .
>
:
an1 x1 + an2 x2 + . . . + ann xn = bn
pode ser representado através de uma equação matricial
Ax = b ,
onde 2 3
a11 a12 ··· a1n
6 .. .. .. .. 7 é matrix dos coeficientes,
A=4 . . . . 5
an1 an2 ··· ann
2 3
x1
6 .. 7
x = 4 . 5 é o vetor colunar das incógnitas,
xn
2 3
b1
6 .. 7
b = 4 . 5 é o vetor dos termos independentes.
bn
Em todo o texto, salvo menção em contrário, sempre indicaremos um sistema linear
genérico de ordem n por Ax = b. Para facilidade de notação usaremos indistintamente
2 3
x1
6 .. 7
x = 4 . 5 ou x = (x1 , . . . , xn ) .
xn

60
5.2 Exemplos de Aplicação
1
5.2.1 Provetas
Considere o seguinte problema: quatro tipos de materiais particulados estão distribuı́dos
por quatro provetas, e em cada proveta os materiais são dispostos em camadas, não mis-
turadas, de modo que seja possı́vel medir facilmente o volume de cada material em cada
uma delas. Dado que possamos medir a massa total de cada proveta, e que saibamos a
massa da proveta vazia, queremos calcular a densidade de cada um dos materiais.
Para colocar o problema em termos matemáticos, chamemos os materiais de A, B, C e
D, e suas densidades respectivas de ⇢A , ⇢B , ⇢C e ⇢D . Essas são as incógnitas do problema,
números que queremos descobrir.
Entre os dados disponı́veis para resolvê-lo estão a massa conjunta dos quatro materiais
em cada uma das provetas (numeradas de 1 a 4), que chamaremos de m1 , m2 , m3 e m4 , já
descontada a tara das provetas.

Além disso, temos o volume de cada um dos materiais em cada uma das provetas.
Chamaremos de v1A , v1B , v1C e v1D o volume dos materiais A, B, C e D na Proveta 1,
v2A , v2B , v2C e v2D o volume dos materiais A, B, C e D na Proveta 2, e assim por diante.
Como a densidade é a razão entre massa e volume, a massa do material A na Proveta
1 é v1 ⇥ ⇢A . Estendendo esse raciocı́nio para os demais materiais, obtemos que a massa
total m1 contida na Proveta 1 é

v1A ⇥ ⇢A + v1B ⇥ ⇢B + v1C ⇥ ⇢C + v1D ⇥ ⇢D .

Considerando as quatro provetas, obteremos quatro equações:

8
>
>v1A ⇥ ⇢A + v1B ⇥ ⇢B + v1C ⇥ ⇢C + v1D ⇥ ⇢D = m1
>
>
<v ⇥ ⇢ + v ⇥ ⇢ + v ⇥ ⇢ + v ⇥ ⇢ = m2
2A A 2B B 2C C 2D D
>
> v 3A ⇥ ⇢A + v 3B ⇥ ⇢B + v 3C ⇥ ⇢C + v 3D ⇥ ⇢D = m3
>
>
:v ⇥ ⇢ + v ⇥ ⇢ + v ⇥ ⇢ + v ⇥ ⇢ = m4
4A A 4B B 4C C 4D D

Trata-se de um sistema linear de quatro equações e quatro incógnitas.

1
Extraı́do de Asano & Coli 2009

61
 Uma possı́vel aplicação em geologia seria a seguinte. Uma sonda faz o papel das
provetas, e uma coluna de material é retirada, contendo materiais diferentes dispostos
em camadas (pode ser até uma sonda coletando material congelado). A sonda permitiria
medir a dimensão de cada camada, mas não poderı́amos desmanchar a coluna para medir
a densidade de cada material isoladamente, sob o risco de alterar a compactação.

5.2.2 Resolução do Cı́rculo

Vamos agora concluir o exemplo iniciado no Capı́tulo 2. Nosso problema era o seguinte: da-
das as coordenadas de três pontos quaisquer, (x1 , y1 ), (x2 , y2 ), (x3 , y3 ), resolver a equação
do cı́rculo que passa por estes três pontos,

(x a)2 + (y b)2 = r2 ,

de forma a determinar o centro do cı́rculo (a, b) e seu raio, r. Temos três incógnitas, de
forma que são necessárias três equações para resolver o problema. São elas
8
> 2 2 2
<(x1 a) + (y1 b) = r
(x2 a)2 + (y2 b)2 = r2
>
:
(x3 a)2 + (y3 b)2 = r2 .

Vamos inicialmente manipular a primeira equação. Expandindo os termos quadráticos

obtemos
x21 + a2 2ax1 + y12 + b2 2by1 r2 = 0 .
Definindo
k ⌘ a 2 + b2 r2
obtemos
2x1 a + 2y1 b k = x21 + y12 .
Manipulando as demais equações da mesma forma, obtemos o seguinte sistema de equações
lineares
8
> 2 2
<2x1 a + 2y1 b k = x1 + y1
2x2 a + 2y2 b k = x22 + y22
>
:
2x3 a + 2y3 b k = x23 + y32 .

que, escrito em forma matricial, fica

2 32 3 2 2 3
2x1 2y1 1 a x1 + y12
42x2 2y2 15 4 b 5 = 4x22 + y22 5 . (5.1)
2x3 2y3 1 k x23 + y32
O problema resume-se, agora, em resolver o sistema acima para obter a, b e k.

5.2.3 Calculando as populações do H em uma região H II

5.3 Método de Cramer

Um método para resolver sistemas lineares, talvez já conhecido do leitor, é o método de
Cramer. Nele a solução do sistema Ax = b é dada por
det(Ai )
xi = , i = 1, 2, . . . , n
det(A)

62
onde det(A) é o determinante da matriz A, e Ai é a matriz obtida de A substituindo a
sua i-ésima coluna pelo vetor b dos termos independentes.
O determinante de uma matriz A de ordem n pode ser calculado através do desenvol-
vimento por linhas (regra de Laplace):
n
X
det(A) = ( 1)i+j aij det(Aij )
j=1

onde i é o ı́ndice de uma linha qualquer e Aij é a matriz obtida de A retirando-se a i-ésima
linha e a j-ésima coluna.
Observe que se det(A) 6= 0 então o sistema Ax = b tem uma única solução. Se
det(A) = 0 então podem ocorrer dois casos:

1. o sistema não possui solução (sistema inconsistente);

2. o sistema possui infinitas soluções (sistema indeterminado).

Por exemplo, no caso de um sistema linear de ordem 2, cada equação representa uma reta.
Resolver o sistema significa determinar a intersecção das duas retas. Se as duas retas
forem coincidentes, então há infinitos pontos de interseção. Se forem paralelas, não há
nenhum ponto de interseção. Neste texto nos preocuparemos com sistemas lineares que
tenham uma única solução.
Uma das propriedades do determinante é que se uma das linhas da matriz for uma
combinação linear de outra (ou outras), então o determinante será zero. Por exemplo, a
matriz abaixo tem determinante zero

1 2
.
3 6

Sistemas que possuem equações que são combinações lineares são ditos degenerados ou
singulares. Como vimos acima, eles podem ser ou inconsistentes ou indeterminados.
Sistemas não singulares (em que o det(A) 6= 0) possuem sempre uma solução. Entre-
tanto, duas questões numéricas podem impedir que a solução seja obtida

1. embora não sejam combinações lineares exatas de outras, algumas equações podem
ser tão próximas a combinações lineares que erros de arredondamento as tornem
linearmente dependentes em algum estágio da solução. Neste caso, o procedimento
numérico irá falhar.

2. Erros de arredondamento cumulativo podem impedir que a solução seja obtida.

Ao longo deste capı́tulo discutiremos formas de lidar com estas duas questões.
A utilização do método de Cramer para resolver sistemas lineares pode ser inviável, pois
o número de operações aritméticas que devem ser efetuadas aumenta consideravelmente
com um pequeno aumento na ordem do sistema.
Para estimar o número de operações necessárias para a regra de Cramer, vamos con-
siderar o caso de um sistema com n = 20. Para resolvê-lo, precisamos calcular 21 deter-
minantes de ordem 20. Mas, para calcular um determinante de ordem 20, usamos a regra
de Laplace, que decompõe o determinante em uma soma envolvendo 20 determinantes
de ordem 19. Se extrapolarmos o processo até chegarmos em determinantes de ordem 2,
teremos que o número de operações aritméticas será da ordem de 21! ⇡ 5 ⇥ 1019 . Para um
sistema de ordem n, temos que o número de operações será da ordem de (n + 1)!.

63
 Em um computador pessoal de 30 Gflops2 estas 1020 operações levariam 3.3 ⇥ 109 s ou
aproximadamente 100 anos! Na prática, a situação é ainda pior, pois estamos considerando
apenas o tempo para efetuar as operações aritméticas, e não o acesso à memória.
No novo super-computador do IAG, que terá uma capacidade teórica de 20 Tflops, esta
conta levaria “apenas” 57 dias. Embora útil para sistemas de ordem menor, o método de
Cramer é impraticável para sistemas maiores, e outros métodos devem ser empregados
neste caso. Outro aspecto negativo do método de Cramer é que como ele necessita de
muitas operações aritméticas, ele potencialmente gerará mais erros de arredondamento.

Exercı́cio 1: use a regra de Cramer para obter uma solução analı́tica para o problema
do cı́rculo (Eq. 5.1).

5.4 Tarefas da álgebra linear computacional

Há muito mais na álgebra linear do que resolver um único sistema de equações lineares.
Abaixo listamos os principais tópicos abordados neste capı́tulo.

- Solução para a equação matricial Ax = b para um vetor colunar desconhecido, x.

- Solução para mais de uma de uma equação matricial, Axj = bj , para um conjunto
de vetores xj , j = 0, 1, . . . , cada um correspondendo a um dado vetor de termos
independentes, bj . Nesta tarefa a simplificação chave é que a matriz A é mantida
constante, enquanto que os termos independentes variam.

- Cálculo da matriz inversa A 1, que obedece à equação matricial AA 1 = I, onde I

é a matriz identidade.

- Cálculo do determinante de uma matriz quadrada A.

- Melhora iterativa da solução de um sistema.

5.5 Sistemas de acordo com as propriedades das matrizes

Tipicamente, podemos ter dois tipos de sistemas lineares, os sistemas cheios e esparsos.
Nos sistemas cheios, todos, ou ao menos a grande maioria, dos elementos da matriz A é
diferente de zero. Nos sistemas esparsos, uma parte importante dos elementos de A é nula.
Um caso importante são sistemas com matrizes tridiagonais, como ilustrado na Figura 5.1.
Sistemas esparsos possuem soluções particulares e mais rápidas que os sistemas cheios.
Vamos inicialmente estudar os métodos de solução para matrizes cheias.

5.6 Método da Eliminação de Gauss

5.6.1 Sobre o método
É o método mais simples para solução de um sistema de equações. O método de Gauss
possui várias caracterı́sticas que o tornam interessante, mesmo que haja métodos mais
eficientes.
Uma caracterı́stica interessante do método é que quando aplicado para resolver um
conjunto de equações lineares, a eliminação de Gauss produz tanto a solução das equações
2
Flops significa número de operações de ponto flutuante por segundo.

64
 Figura 5.1: Exemplos de matrizes esparsas.

(para um ou mais vetores de termos independentes) quanto a inversa da matriz A (esta

última é obtida quando empregamos uma variante do método, chamada de método de
Gauss-Jordan, seção 5.7). Uma de suas caracterı́sticas mais importantes é que o método é
tão estável quanto qualquer outro método direto (direto, aqui, é usado em contraposição
aos métodos iterativos mostrados no fim do capı́tulo), desde que seja empregado o pivo-
tamento (seções 5.6.4 e 5.7.2)
Algumas deficiências do método são

1. se a matriz inversa não for desejada, o método de Gauss é tipicamente 3 vezes mais
lento que a melhor alternativa disponı́vel (decomposição LU, seção 5.10).

2. quanto o empregamos para mais de uma equação matricial (Axj = bj ), todos os

vetores de termos independentes devem ser armazenados na memória e manipulados
simultaneamente.

A deficiência 1) acima pode suscitar questionamentos, afinal, se temos a matriz inversa,

podemos calcular as incógnitas de um sistema Axj = bj através de:
1
xj = A bj .

Isto realmente funciona, mas este procedimento resulta em uma resposta muito suscetı́vel
a erros de arredondamento, e deve ser evitado.

5.6.2 Procedimento
Vamos ilustrar o procedimento do método de eliminação de Gauss com um exemplo sim-
ples. O objetivo consiste em transformar o sistema Ax = b em um sistema triangular
equivalente. Para isso, usarmos a seguinte propriedade da Álgebra Linear.

Propriedade: A solução de um sistema linear não se altera se subtrairmos de uma

equação outra equação do sistema multiplicada por uma constante.

Considere o seguinte sistema de equações:

8
>
<2x + y + z = 7
4x + 4y + 3z = 21
>
:
6x + 7y + 4z = 32

65
Multiplicando a primeira equação por (-2) e somando na segunda, e multiplicando a pri-
meira equação por (-3) e somando na terceira temos
8
>
<2x + y + z = 7
2y + z = 7
>
:
4y + z = 11
Multiplicando a segunda equação por (-2) e somando na terceira temos
8
>
<2x + y + z = 7
2y + z = 7
>
:
z= 3

Da terceira equação temos z = 3 ) z = 3 .

Substituindo na segunda equação temos 2y + 3 = 7 ) y = 2 .
Substituindo na primeira equação temos 2x + 2 + 3 = 7 ) x = 1 .

Vemos que o método de Gauss é uma forma sistemática de triangularizar um sistema

linear. A solução é obtida em dois passos:
1. Eliminação (forward elimination): triangularização propriamente dita.
2. Substituição (back substitution): obtenção da solução final (vetor x).
Se usarmos a notação matricial, estamos resolvendo a equação
2 32 3 2 3
2 1 1 x 7
44 4 35 4y 5 = 4215
6 7 4 z 32
transformando-a em 2 3 2 3 2 3
2 1 1 x 7
4 2 1 5 4 y = 75.
5 4
1 z 3
| {z }
matriz
triangular superior
Entretanto, podemos trabalhar somente com os números sem escrever as equações.
Para tanto é conveniente escrever a chamada matriz aumentada
2 3
2 1 1 7
44 4 3 215 .
6 7 4 32
Como antes multiplicamos a primeira equação por 2 e subtraı́mos da segunda; multiplica-
mos a primeira equação por 3 e subtraı́mos da terceira:
2 3
2 1 1 7
40 2 1 7 5 .
0 4 1 11
Então multiplicamos a segunda equação por 2 e subtraı́mos da terceira
2 3
2 1 1 7
40 2 1 75.
0 0 1 3

66
5.6.3 Estimativa do número de operações realizadas
Vamos estimar o número de operações realizadas na obtenção da solução x. Estimaremos
separadamente o número de operações feitas durante a eliminação e a substituição.

1) Processo de eliminação
Para estimar o número de operações realizadas durante a triangulação da matriz,
calcularemos quantas adições e multiplicações são necessárias em cada etapa do processo.
Por exemplo, para eliminarmos a primeira coluna, temos (n 1) linhas onde para cada
uma delas são calculadas n + 1 multiplicações e n adições.

eliminação da: multiplicações adições

1a coluna (n 1)(n + 1) (n 1)n
2a coluna (n 2)n (n 2)(n 1)
.. .. ..
. . .
(n 1)a coluna (1)(3) (2)(1)
Pn 1 Pn 1
Total i=1 i(i + 2) i=1 (i + 1)i

O total de multiplicações é
n
X1 n
X1 n
X1
i(i + 2) = i2 + 2 i.
i=1 i=1 i=1

Avaliando cada uma das somatórias

n
X1 n
X
2 n(n + 1)(n + 2) n3 n2 n
i = i2 n2 = n2 = +
6 3 2 6
i=1 i=1
n
X1 n
X n(n + 1) n2 n
i= i n= n= ,
2 2 2
i=1 i=1

que implica
n
X1 n
X1 ✓ ◆
2 n3 n2 n n2 n n3 n2 5n
i +2 i= + +2 = + .
3 2 6 2 2 3 2 6
i=1 i=1

O número total de adições pode ser obtido de forma análoga:

n
X1 n3 n
(i + 1)i = .
3 3
i=1

Obtemos, assim, que o número total de operações de ponto flutuante para o processo
de eliminação é
2n3 n2 7n
Nelim = + .
3 2 6
Para um valor de n suficientemente grande, temos que os termos n3 dominam nas ex-
pressões acima, de forma que o total de operações na eliminação será O(2n3 /3).

2) Processo de substituição
Vamos agora estimar quantas operações de ponto flutuante são feitas durante o cálculo
da solução final a partir da matriz triangularizada (back substitution).

67
 passo multiplicações adições
linha n 1 0
linha n 1 2 1
.. ..
. .
linha 1 n n 1
Pn Pn 1
Total 1 i 1 i

Obtemos que o número de operações para esta fase

Nsubst = n2 .

Chegamos, assim, ao o número total de operações necessárias para resolver um sistema

de ordem n pelo método de Gauss

2n3 3n2 7n
NGauss = + .
3 2 6
Concluı́mos que para valores altos de n o processo de eliminação necessita de um número
muito maior de operações que a substituição e que, neste caso, o total de operações é

2n3
NGauss ⇡ .
3
Por exemplo, um sistema matricial de 20 ⇥ 20 implica em aproximadamente 2 · 203 /3 ⇡
5 · 103 flop. Com um PC de 30 Gflops o problema será resolvido em

5 · 103 flop 7
t= ⇡ 2 · 10 s!
30 · 109 flops
Esta estimativa é muito otimista, pois consideramos que cada operação de ponto flu-
tuante é efetuada em um ciclo da CPU. Isto é válido para adições, mas não para mul-
tiplicações, que tipicamente requerem da ordem de dez ciclos de CPU. Além disso não
consideramos fatores como a perda de eficiência devido ao acesso à memória. De qualquer
maneira, vemos que o método de Gauss é imensamente mais eficiente que o método de
Cramer.

5.6.4 Pivotamento parcial

Seja o sistema 2 32 3 2 3
10 7 0 x1 7
4 3 2.099 65 4x2 5 = 43.9015 (5.2)
5 1 5 x3 6
cuja solução é x = [0, 1, 1].
Vamos considerar que nosso operador de ponto flutuante tenha apenas 5 algarismos
significativos, e vamos resolver o sistema acima pelo método de Gauss.
Multiplicando a 1a equação por 0.3 e somando na 2a ; multiplicando a 1a equação por
-0.5 e somando na 3a , obtemos
2 3
10 7 0 7
40 0.001 6 6.0015 . (5.3)
0 2.5 5 2.5

68
 Multiplicando a 2a equação por 2.5/ 0.001 = 2500 e somando na 3a equação
2 3
10 7 0 7
40 0.001 6 6.001 5 .
0 0 15005 15004

Note que, devido à restrição de 5 algarismos significativos, tivemos que truncar as seguintes
operações

6.001 ⇥ 2500 = 15002.5

15002.5 + 2.5 = 15004.5 = 15004 .

Ao efetuarmos a substituição obteremos

15004
x03 = = 0.99993
15005
6.001 6 ⇥ 0.99993 6.001 5.99958
x02 = = = 1.5
0.001 0.001
7 + 7 ⇥ ( 1.5)
x01 = = 7 10.510 = 0.35
10

Comparando este vetor x0 = (0.99993, 1.5, 0.35) obtido com o vetor x = (1, 1, 0)
solução, vemos o quão grande foi o erro gerado pela restrição de 5 algarismos significativos!
O que causou este problema? O primeiro elemento da linha que está sendo usada para
eliminar os termos das demais é chamado de pivô. Na primeira etapa da eliminação acima
(Eq. 5.3), o pivô, ( 0.001), tornou-se muito pequeno em relação aos outros coeficientes,
resultando num enorme multiplicador (2500) que fez aparecerem erros de arredondamento.
Estes erros por sua vez são ampliados na fase de substituição, onde apareceram subtrações
de números muito próximos divididas por números muito pequenos, o que amplifica enor-
memente o erro (por exemplo, veja o cálculo de x02 , acima).
Uma solução simples e eficiente para este problema é empregar pivotamento parcial no
método de Gauss, que consiste em trocar linhas de forma que tenhamos sempre o maior
valor absoluto possı́vel para o pivô. Isto garantirá multiplicadores . 1 em módulo.
No exemplo acima, empregamos o pivotamento parcial já na primeira etapa
2 3 2 3
10 7 0 7 10 7 0 7
pivotamento
40 0.001 6 6.0015 ! 40 2.5 5 2.5 5 .
parcial
0 2.5 5 2.5 0 0.001 6 6.001

O multiplicador será ( 0.001)/( 2.5) = 0.0004. Multiplicando a 2a equação por este valor
e somando na 3a equação, obtemos a matriz estendida
2 3
10 7 0 7
4 0 2.5 5 2.5 5 ,
0 0 6.002 6.002

que resulta na solução exata x0 = (1, 1, 0).

Uma regra importante a ser seguida: o pivotamento parcial sempre deve ser empregado
no método de Gauss!

69
5.6.5 Solução simultânea de várias equações matriciais
Vimos na Seção 5.4 uma tarefa corriqueira da álgebra linear é resolver um conjunto de
equações matriciais, Axj = bj , j = 1, . . . , m. Neste conjunto as equações matriciais com-
partilham a matriz A e possuem cada uma um dado dado vetor de termos independentes,
bj . Neste caso, em vez de fazer a mesma eliminação m vezes, podemos “guardar” a
sequência de operações aplicadas na triangulação da matriz A para depois aplicar em bj ,
j = 1, . . . , m.
Por exemplo, seja a matriz 2 3
2 6 2
A = 41 3 45
3 6 9
e o vetor de pivotamento que contém o número da linha que foi pivotada,
2 3

p=4 5.

Com o pivotamento da 3a linha, temos

2 3 2 3
3 6 9 3
41 3 5
4 ; 4 5.
2 6 2

Fazendo 1a linha ⇥ ( 1/3) + 2a linha e 1a linha ⇥ ( 2/3) + 3a linha,

2 3 2 3
3 6 9 3
6 1/3 1 77 4 5
4 5 ; .
2/3 2 8

Nesta operação, preservamos os multiplicadores que foram utilizados para eliminar os

primeiros coeficientes das linhas que não eram o pivô. Para concluir a triangularização da
matriz, novamente utilizamos o pivotamento da 3a linha:
2 3 2 3
3 6 9 3
6 1/3 2 8 7 4 5
4 5; 3 .
2/3 1 7

Fazendo 2a linha ⇥ ( 1/2) + 3a linha,

2 3 2 3
3 6 9 3
6 1/3 87
4 2 5 ; 435 .
2/3 1/2 3

Para ilustrar como utilizar os multiplicadores armazenados e o vedor de pivotamento,

vamos encontrar a solução para o vetor b = [4, 7, 39].

- 1o passo: trocar a linha 1 com a linha p(1) = 3 ! b = [39, 7, 4];

- 2o passo: multiplicar a 1a linha por 1/3 e somar na 2a linha; multiplicar a 1a

linha por 2/3 e somar na 3a linha, ! b = [39, 20, 22];

70
 - 3o passo: trocar a linha 2 com a linha p(2) = 3 ! b = [39, 22, 20];

- 4o passo: Multiplicar a 2a linha por 1/2 e somar na 3a linha, b = [39, 22, 9].

Assim, o vetor x com a solução do sistema será dado por

2 32 3 2 3
3 6 9 x1 39
40 2 85 4x2 5 = 4 225 ,
0 0 3 x3 9
que pode ser facilmente resolvido por substituição

x3 = 9/( 3) = 3
x2 = [ 22 + 8 ⇥ 3]/2 = 1
x1 = [39 6⇥1 9 ⇥ 3]/3 = 2 .

.
O procedimento ilustrado pode ser repetido para um número arbitrário de vetores
bj . Uma sugestão para uma implementação eficiente do método de Gauss é fazer uma
subrotina para a eliminação, que retorna a matriz triangularizada com os coeficientes de
eliminação, segundo procedimento acima, e outra para a substituição.

Abaixo delineamos um possı́vel algoritmo para implementar a eliminação de Gauss

computacionalmente, mantendo os multiplicadores para uso posterior:
8
>
> de i = 1 até n 1 faça
>
>
>
> determine o ı́ndice do pivoteamento l i
>
>
>
>
> troque as linhas i e l, das colunas i até n
>
>
>
>
>
< registre
8 o i-ésimo pivotamento: p(i) = l
>
> de j = i + 1 até n faça
>
> >
>
>
> < calcule o multiplicador para a linha j
>
> laço
>
> >
>
> >
> guarde-o no lugar do elementos eliminado
>
> >
:
>
> adicione múltiplos da linha l à linha j
>
>
:
fim do laço

5.6.6 Cálculo do determinante de uma matriz A

Pelas propriedades do determinante, o determinante não se altera se somarmos um múltiplo
de uma linha da matriz à outra, ou seja, se efetuarmos uma combinação linear entre as
linhas. Assim, a eliminação de Gauss para se obter uma matriz triangular superior não
afeta o valor do determinante
2 3 2 0 3
a11 a12 · · · a1n a11 a012 · · · a01n
6 a21 a22 · · · a2n 7 6 0 a0 a02n 7
6 7 6 22 · · · 7
det 6 . .. .. 7 = det 6 .. .. .. 7
4 .. . . 5 4 . . . 5
an1 an2 · · · ann 0 0 ··· a0nn
| {z } | {z }
A U

Mas o determinante de uma matriz triangular é o produto dos elementos da matriz.

Portanto,
det A = det U = a011 a022 . . . ha0nn .

71
 Atenção: cada operação de pivotamento troca o sinal do determinante! Assim, se para
implementar a eliminação de Gauss foram realizados np pivotamentos, o determinante será

det A = ( 1)np a011 a022 . . . aa0nn .

Exemplo: calcule o determinante da matriz abaixo usando o método de Gauss:

2 3
2 1 1
A = 4 4 4 35 .
6 7 4
Usando a eliminação de Gauss sem pivotamento, obtemos
2 3
2 1 1
U= 04 2 1 5 ! det A = det U = 2 ⇥ 2 ⇥ ( 1) = 4.
0 0 1

Exercı́cio 2: Calcule, usando o método de Gauss com pivotamento parcial a solução

do sistema 2 32 3 2 3
4 3 2 2 x1 5
62 1 1 27 6 x 2 7 6 87
6 76 7 = 6 7
42 2 2 45 4x3 5 435
6 1 1 4 x4 1

5.7 Método de Gauss-Jordan

Este método é uma variante do método de Gauss, onde são eliminados todos os elementos
acima e abaixo do pivô. O resultado da eliminação de Gauss-Jordan é uma matriz diagonal.
Para ilustrar método, vamos aplicá-lo à solução de um sistema Ax = b e à obtenção
simultânea da matriz inversa de A.
Considere a equação matricial

A · [x1 _ x2 _ Y ] = [b1 _ b2 _ I] . (5.4)

onde A e Y são matrizes quadradas, xi e bi são vetores colunares, I é a matriz identidade,

o operador (·) significa um produto de matrizes e o operador _ significa o aumento de
matriz, ou seja, a remoção dos parênteses das matrizes para fazer uma matriz mais larga.
É fácil perceber, da equação acima, que os x1 e x2 são simplesmente a solução das equações
matriciais

A · x 1 = b1 ,
A · x 2 = b2 ,

e que a matriz Y é a inversa de A, ou seja

A·Y =I.

Para simplificar, mas sem perda de generalidade, vamos podemos escrever explicita-
mente a Eq. (5.4) usando matrizes de ordem 3
2 3 82 3 2 3 2 39 82 3 2 3 2 39
a11 a12 a13 < x11 x12 y11 y12 y13 = < b11 b12 1 0 0 =
4a21 a22 a23 5· 4x21 5 _ 4x22 5 _ 4y21 y22 y23 5 = 4b21 5 _ 4b22 5 _ 40 1 05
: ; : ;
a31 a32 a33 x31 x32 y31 y32 y33 b31 b32 0 0 1

72
 Usando a notação de matrizes aumentada, o sistema acima fica
2 3
a11 a12 a13 b11 b12 1 0 0
4a21 a22 a23 b21 b22 0 1 05
a31 a32 a33 b31 b32 0 0 1

5.7.1 Exemplo de aplicação: inversão de matrizes

A · [x _ Y ] = [b _ I] . (5.5)
2 3
a11 a12 a13 b1 1 0 0
4a21 a22 a23 b2 0 1 05
a31 a32 a33 b3 0 0 1
2 3
2 1 1 7 1 0 0
44 4 3 21 0 1 05
6 7 4 32 0 0 1
2 3
2 1 1 7 1 0 0
40 2 1 7 2 1 05
0 4 1 21 3 0 1
A partir daqui, elimina-se os elementos superiores e inferiores ao pivô:
2 3
2 0 1/2 7/2 2 1/2 0
40 2 1 7 2 1 05
0 0 1 3 1 2 1
2 3
2 0 0 2 5/2 3/2 1/2
40 2 0 4 1 1 1 5
0 0 1 3 1 2 1
Finalmente, normaliza-se a matriz, de forma que à esquerda ficamos com uma matriz
identidade. Obtém-se, assim, o vetor solução do problema e a matriz inversa de A.
2 3 2 3
1 0 0 1 5/4 3/4 1/4
40 1 05 2 4 1/2 3/2 1/25 .
0 0 1 3 1 2 1
|{z} | {z }
⌘X ⌘A 1

5.7.2 Pivotamento total

Vimos acima (seção 5.6.4) um exemplo em que o pivotamento parcial foi usado para
evitar erros de arredondamento que podem aparecer quando o multiplicador fica muito
pequeno. Partindo do princı́pio que os erros de arredondamento são tão menores quanto
maiores forem os multiplicadores, em geral obtém-se melhores resultados empregando-se o
pivotamento total, em que se pivotam colunas, além das linhas, de forma a sempre manter
o maior termo de uma data linha como pivô.
O pivotalmento total pode ser empregado devido à seguinte propriedade da Álgebra
Linear:

Propriedade: a solução de um sistema linear não é alterada quando trocamos de

lugar duas colunas i e j de A, desde que se troquem as duas linhas correspondentes nos
vetores x e na matriz Y (Eq. 5.5).

73
 Desta forma, vemos que as operações de troca de colunas embaralham o(s) vetor(es)
das incógnitas e a matriz inversa, de forma que para empregar o pivotamento total devemos
manter um registro das trocas de colunas efetuadas para podermos colocar a solução final
na ordem correta. Para exemplificar, vamos resolver novamente o sistema da Eq. (5.2).
Após a eliminação da primeira coluna, obtemos
2 3
10 7 0 7
40 0.001 6 6.0015 .
0 2.5 5 2.5

Vamos agora trocar de lugar as colunas 2 e 3, de forma que o coeficiente 6 será o pivô
2 3
10 0 7 7
40 6 0.001 6.0015 .
0 5 2.5 2.5

Eliminando o segundo termo da terceira linha, obtemos

2 3
10 0 7 7
40 6 0.001 6.001 5 ,
0 0 2.50083 2.50083

que pode ser facilmente resolvido para obtermos xp = (0, 1, 1). Como fizemos uma troca
de colunas, a solução final é obtida trocando-se de lugar as linhas 2 e 3 de xp , ou seja,
x = (0, 1, 1).

5.8 Refinamento da Solução

Seja o sistema
Ax = b . (5.6)
Resolvendo por Gauss ou Gauss-Jordan, obtemos a solução x(0) . Sabemos que erros de
arredondamento podem ocorrer quando se resolve um sistema linear pelo método de eli-
minação, podendo compromer o resultado obtido. Mesmo utilizando pivoteamento total,
não se pode assegurar que a solução obtida seja exata.
Inicialmente, notemos que é trivial verificarmos se a solução de um sistema está correta,
para isso basta multiplicarmos a matriz A pela solução obtida, x(0) , e o resultado deve ser
b. Numericamente, esta verificação deve ser feita impondo-se um critério de convergência
do tipo:
(0)
bi bi
< ✏, i = 1, . . . , n ,
bi

onde b(0) é o vetor obtido do produto Ax(0) .

O que fazemos quando o resultado obtido x(0) não passa pelo critério de convergência?
Uma possibilidade é fazermos o refinamento da solução, como delineado a seguir. Vamos
chamar de erro a diferença entre o valor verdadeiro, x, e o valor obtido, e(0) = x x(0) .
Substituindo no sistema (5.6), temos

A(x(0) + e(0) ) = b (5.7)

Ae(0) = b Ax(0) ⌘ r(0) . (5.8)

74
 Resolvendo o sistema (5.8) determinamos e(0) , a partir do qual podemos fazer uma nova
estimativa para a solução: x(1) = x(0) + e(0) . Caso x(1) obedeça ao critério de convergência
estipulado, teremos encontrado nossa solução. Caso contrário, podemos refinar novamente
a solução obtendo uma estimativa para o erro de x(1) resolvendo o sistema

Ae(1) = b Ax(1) ⌘ r(1) . (5.9)

Este processo pode ser executado quantas vezes desejarmos. É fundamental que as
operações envolvidas nos cálculos dos resı́duos sejam feitas em dupla precisão.

Importante: Como o refinamento envolve a resolução de vários sistemas que com-

partilham a mesma matriz A, podemos empregar o procedimento descrito na seção 5.6.5
para tornar o processo mais eficiente.

5.9 Sistemas mal-condicionados

Estes sistemas também são conhecidos pelo termo mal-condicionados (“ill conditioned”,
em inglês). Vejamos um exemplo.
(
x+ y = 1
99x + 100y = 99.5
A solução deste sistema é única e exata: x = 0.5, y = 0.5. Agora considere o sistema
(
x+ y = 1
,
99.4x + 99.9y = 99.2

cuja solução única e exata é x = 1.4, y = 0.4, muito diferente da do sistema anterior,
apesar dos coeficientes serem parecidos.
Graficando as retas no plano (x, y) vemos porque isto acontece: as retas corresponden-
tes às equações são quase paralelas, ou seja, as equações são quase linearmente dependen-
tes.
Uma maneira de se “medir” o condicionamento de uma traz é calcular seu determi-
nante. Entretanto, o valor do determinante depende da norma (módulo) de cada um dos
seus vetores. Assim, devemos normalizar os vetores para depois calcular o determinante.
Isto produzirá um determinante com valor (em módulo) no intervalo [0, 1]. Se o módulo
do determinante for próximo de zero, então o sistema é mal-condicionado.
Por exemplo, vamos considerar o primeiro exemplo acima. Normalizando os vetores
da matriz 
1 1
,
99 100
obtemos  p p
1/ 2 1/ 2
.
99/140.716 100/140.716
O módulo do determinante da matriz normalizada é
1 100 1 99 3
p p ⇡ 5 ⇥ 10 ,
2 140.716 2 140.716

o que demonstra, quantitativamente, que a matriz original é mal-condicionada.

75
 Há outras medidas do condicionamento de uma matriz, assim como há fórmulas que
relacionam o erro cometido no método de Gauss ou Gauss-Jordan com essas medidas e
o número de algarismos significativos utilizado. Isto, porém, está além do escopo des-
tas notas. Veja, por exemplo, a seção sobre singular value decomposition no Numerical
Recipes.

5.10 Decomposição LU
Suponhamos que se possa escrever a matriz A como o produto de duas matrizes

A=L·U, (5.10)

onde L é uma matriz triangular inferior (tem elementos somente na diagonal e abaixo) e
U é uma matriz diagonal superior (com elementos somente na diagonal e acima). Para o
caso de uma matriz 4 ⇥ 4, Eq. (5.10) ficaria
2 3 2 32 3
a11 a12 a13 a14 l11 0 0 0 u11 u12 u13 u14
6a21 a22 a23 a24 7 6l21 l22 0 07 6 7
6 7 6 7 6 0 u22 u23 u24 7 .
4a31 a32 a33 a34 5 = 4l31 l32 l33 0 5 4 0 0 u33 u34 5
a41 a42 a43 a44 l41 l42 l43 l44 0 0 0 u44
| {z }| {z }
L U
triangular triangular
inferior superior
Pode-se usar esta decomposição para resolver o conjunto de equações lineares

Ax = (LU )x = L(U x) = b .

Inicialmente resolvemos para o vetor y tal que

Ly = b (5.11)

e depois resolvemos
Ux = y , (5.12)
para obter a solução final.
Qual a vantagem em quebrarmos um sistema em dois sistemas sucessivos? A vantagem
é que a solução de um sistema triangular é trivial. Dessa forma, o sistema (5.11) é resolvido
por substituição para frente:
b1
y1 = (5.13)
l11
2 3
i
X
1 4
yi = bi lij yi 5 , i = 2, 3, . . . , n , (5.14)
lii
j=1

e o sistema (5.12) por substituição para trás:

yn
xn = (5.15)
unn
2 3
n
X
1 4
xi = yi uij xj 5 , i = n 1, n 2, . . . , 1 . (5.16)
uii
j=i+1

76
5.10.1 Efetuando a Decomposição LU
Como podemos achar L e U dado A? Abaixo vamos delinear um algoritimo bastante
utilizado, que pode ser estudado em mais detalhes no Numerical Recipes. Vamos escrever
explicitamente o componente i.j da Eq. (5.10). Este componente é sempre uma soma que
começa co,
li1 u1j + · · · = aij .
O número de termos da soma depende se i < j, i > j, ou i = j. De fato, temos os três
casos acima

i < j : li1 u1j + li2 u2j + · · · + lii uij = aij (5.17)

i = j : li1 u1j + li2 u2j + · · · + lii ujj = aij (5.18)
i > j : li1 u1j + li2 u2j + · · · + lij ujj = aij (5.19)

As Eqs. (5.17) — (5.19) perfazem n2 equações para n2 + n incógnitas (note que a diagonal
está representada duas vezes). Trata-se, assim, de um sistema indeterminado. Para resol-
ver este sistema, deve-se, assim, especificar arbitrariamente valores para n incógnitas. Um
procedimento muito usado para resolver a decomposição é o Algoritmo de Crout, que re-
solve de forma trivial as equações acima para todos os l’s e u’s simplesmente rearranjando
as equações em determinada ordem. O algoritmo é como se segue:

- Faça lii , i = 1, · · · , n (de forma a reduzir o número de incógnitas para n2 );

- Para cada j = 1, · · · , n, faça os dois procedimentos seguintes:

– Primeiramente, para i = 1, · · · , j, use Eqs. (5.17) e (5.18), e a condição acima,

para determinar os uij , ou seja
i 1
X
uij = aij lik ukj . (5.20)
k=1

Quando i = 1 a soma anterior é tomada como zero.

– Em segundo lugar, para i = j +1, · · · , n use (5.17) para achar os lij , da seguinte
maneira !
j 1
X
1
lij = aij lik ukj . (5.21)
ujj
k=1

Certifique-se de executar ambos os procedimentos antes de passar para o próximo j!

A chave para compreender o procedimento acima é que os l’s e u’s que ocorrem no lado
direito das equações (5.20) e (5.21) já estão sempre determinados no momento em que são
necessários (por isso que o método de Crout é basicamente uma ordem em que as equações
devem ser resolvidas). Vemos, também, que cada aij é usado apenas uma vez durante o
processo. Isso significa que os l’s e u’s podem ser armazenados nos lugares que os termos
a’s ocupavam na memória do computador. Ou seja, o método de Crout substitui a matriz
original A por uma matriz combinada dos elementos de L e U :
2 3
u11 u12 u13 u14
6 l21 u22 u23 u24 7
6 7
4 l31 l32 u33 u34 5 .
l41 l42 l43 u44

77
 O pivotamento, tal como no caso dos métodos de Gauss e Gauss-Jordan, é essencial
para a estabilidade do método de Crout. Quando se emprega o pivotamento parcial, na
realidade não se decompõe a matriz A na sua forma LU , mas sim uma permutação das
linhas de A. Para ver como efetuar o pivotamento no método de Crout, consulte o capı́tulo
2 do Numerical Recipes.
Qual a vantagem da Decomposição LU sobre o método de Gauss? Como listado na
seção 5.12, o número de operações necessárias para efetuar a decomposição é da ordem
de 1/3 n3 , exatamente o mesmo número de passos necessários para fazer a eliminação de
Gauss. Na literatura, frequentemente cita-se uma vantagem da Decomposição LU que é o
fato de que uma vez tendo-se L e U é trivial obter a solução para um número arbitrário de
vetores de termos independentes (ou seja, resolve-se facilmente um conjunto de sistema de
equações lineares). Entretanto, o mesmo procedimento pode ser feito de forma igualmente
eficiente à partir do procedimento delineado na seção 5.6.5.
Conclusão: o método de Gauss e o método da Decomposição LU são igualmente
eficientes quando se trata de resolver um sistema de equações lineares, ou um conjunto de
sistemas de equações lineares.

5.10.2 Um caso especial: decomposição LU de matrizes tridiagonais

Um caso particular em que a decomposição LU oferece uma solução eficiente é no caso de
matrizes tridiagonais. Suponha que A seja uma matriz na forma
2 3
b1 c 1
6 a 2 b2 c2 7
6 7
6 .. .. .. 7
A=6 . . . 7,
6 7
4 a n 1 bn 1 c n 1 5
an bn

ou seja, apenas os elementos da diagonal principal e das diagonais imediatamente acima

e abaixo são não nulos. Neste caso, a decomposição LU então tem uma forma simples
2 3
1
6 l2 1 7
6 7
6 .. .. 7
L=6 6 . . 7,
7
6 .. .. 7
4 . . 5
ln 1
2 3
u1 v1
6 u2 v2 7
6 7
6 .. .. 7
U =6
6 . . 7.
7
6 .. 7
4 . vn 15
un

É fácil demonstrar que a determinação dos l’s, u’s e v’s é feita através da seguintes
relações de recorrência:
8
>
< u 1 = b1
lj = aj /uj 1
>
:
uj = bj lj cj 1 , j = 2, 3, . . . , n

78
Note que nas relações acima está implı́cito que vi = ci .
Tendo determinado as matrizes L e U através das relações acima, o procedimento para
determinar o vetor solução x do sistema Ax = d (note que aqui usamos d para evitar
confusão com os b’s da matriz tridiagonal) é simples. Inicialmente calcula-se a solução do
sistema Ly = d através da substituição para frente:
(
y1 = d 1
yi = di li yi 1 , i = 2, . . . , n

Calcula-se, então, o vetor solução x resolvendo-se o sistema U x = y por substituição para

trás: (
xn = uynn
x
xk = yk ck uk+1 k
, k = n 1, . . . , 1
O procedimento descrito acima é muito eficiente do ponto de vista computacional e
pode ser implementado com facilidade em duas subrotinas, uma para o cálculo da decom-
posição e outra para a solução do sistema. Note que o fato de que a decomposição LU de
uma matriz tridiagonal também é tridiagonal simplifica muito as substituições para frente
e para trás. Veremos no Cap. 6 que esta solução para um sistema tridiagonal será muito
útil para calcular os coeficientes da interpolação por spline cúbica.

5.11 Forma alternativa para o cálculo da matriz inversa

Denotemos a matriz inversa de A por B, tal que:

AB = I ,

onde I é a matriz identidade. Como usar a decomposição LU ou o método descrito na

seção (5.6.5) para encontrar B? Isso é feito simplesmente escrevendo a equação matricial
acima para cada uma das colunas de B, ou seja,
2 3 2 3
b11 1
6b21 7 607
A6 7 6 7
4b31 5 = 405 ,
b41 0
2 3 2 3
b12 0
6b22 7 617
A6 7 6 7
4b32 5 = 405 ,
b42 0
e assim por diante. Ou seja, o cálculo da matriz inversa reduz-se a resolver um conjunto
de n sistemas lineares em que os vetores de termos independentes são as diferentes colunas
da matriz identidade.

5.12 Comparando Gauss, Gauss-Jordan, e Decomposição

LU
Concluı́mos esta parte do capı́tulo comparando a eficiência dos três métodos diretos estu-
dados. Temos no quadro abaixo a comparação do número de operações empregadas em
cada método para as diferentes tarefas da álgebra linear.

79
 método operações
Solução de Gauss 1/3 n3
Sistemas Gauss-Jordan 1/2 n3
Lineares LU 1/3 n3
Inversão Gauss 5/6 n3
de Gauss-Jordan n3
Matriz LU 5/6 n3
1 3 1 2
m Gauss 3 n + 2 mn
1 3 2
lados Gauss-Jordan 2 n + mn
1 3 1 2
direitos LU 3 n + 2 mn

5.13 Métodos Iterativos

Os métodos que vimos até agora (Gauss, Gauss-Jordan, Decomposição LU) são conhecidos
como métodos diretos, pois a solução é obtida através da manipulação direta das equações
do sistema. Tais métodos podem se tornar ineficientes quando o número de equações fica
muito grande (n & 100), pois o número de operações de ponto-flutuante é O(n3 ). Mais
detalhes em Blum(1972), p.131.
Nos métodos ditos iterativos (também chamados de métodos indiretos), arbitra-se um
vetor inicial x(0) para a solução e calcula-se uma nova estimativa da solução, x(1) como
função de x(0) e assim sucessivamente, ou seja,

xk+1 = g(xk ) ,

onde k é a k-ésima iteração e g representa uma função qualquer. O processo é repetido

até obter a precisão desejada, que se traduz em uma diferença muito pequena entre xk+1
e xk .
Nota: não confunda métodos iterativos com a melhora iterativa da solução, apresentada
na seção 5.8.

5.13.1 Método de Jacobi

Seja um sistema linear de ordem n
8
>
> a11 x1 + a12 x2 + . . . + a1n xn = b1
>
>
<a21 x1 + a22 x2 + . . . + a2n xn = b2
> .. ..
>
> . .
>
:
an1 x1 + an2 x2 + . . . + ann xn = bn

Podemos reescrevê-lo na seguinte forma

8 1
>
>x1 = a11 (b1 a12 x2 a13 x3 . . . a1n xn )
>
>
<x2 = 1 (b2 a21 x1 a23 x3 . . . a2n xn )
a22
> ..
>
> .
>
: 1
xn = ann (bn an1 x1 an2 x2 . . . an,n 1 xn 1)

No método de Jacobi, escolhemos arbitrariamente um vetor inicial x(0) e substituı́mos

no lado direito das equações acima obtendo um novo vetor x(1) . Repetindo-se o pro-
cesso k vezes, vemos que a k-ésima estimativa da solução é obtida da seguinte relação de
recorrência:

80
 8 (k+1) 1 (k) (k)
>
> x1 = a11 (b1 a12 x2 ... a1n xn )
>
>
<x(k+1) = 1 (k) (k)
2 a22 (b2 a21 x1 ... a2n xn )
.
>..
>
>
>
: (k+1) (k) (k)
1
xn = ann (bn an1 x1 ... an,n 1 xn 1 )

ou, de forma mais compacta,

0 1
n
X
(k+1) 1 BBb i (k) C
xi = aij xj C A.
aii @
j=1
j6=i

Dizemos que o processo iterativo converge se, para a sequência de aproximações gerada,
dado ✏ > 0, existir um j tal que para todo k > j e i = 1, 2, . . . , n, |xki x̄i  ✏, onde x̄ é
a solução do sistema. Como na prática não a conhecemos, torna-se necessário um critério
de parada para o processo iterativo. Um possı́vel critério é impor que a variação relativa
entre duas aproximações consecutivas seja menor que ✏. Dado x(k+1) e x(k) , tal condição
é escrita como ( (k+1) )
(k)
xi xi
max (k)
, i = 1, . . . , n  ✏ . (5.22)
xi

Exemplo: considere o seguinte sistema de equações

8
>
<4x1 + 2x2 + x3 = 11
x1 + 2x2 =3 ,
>
:
2x1 + x2 + 4x3 = 16

cuja solução é x = (1, 2, 3). Rescrevendo as equações como

8 x3
11 1
>
< x1 = 4 2 x2 4
x2 = 32 + 12 x1 ,
>
: 1 1
x3 = 4 2 x1 4 x2

temos que as relações de recorrência, pelo método de Jacobi, são

(k)
(k+1) 11 1 (k) x3
x1 = x ,
4 2 2 4
(k+1) 3 1 (k)
x2 = + x , (5.23)
2 2 1
(k+1) 1 (k) 1 (k)
x3 =4 x x .
2 1 4 2
Começando com um vetor arbitrário x(0) = [1, 1, 1] obtemos
(1) 11 1 1
x1 = ·1 ·1=2
4 2 4
(1) 3 1
x2 = + ·1=2
2 2
(1) 1 1 13
x3 = 4 .1 ·1=
2 4 4

81
Substituindo x(1) do lado direito do sistema (5.26) obtemos

(2) 11 1 1 13 15
x1 = ·2 · =
4 2 4 4 16
(2) 3 1 5
x2 = + ·2=
2 2 2
(2) 1 1 5
x3 = 4 .2 ·2=
2 4 2
Na tabela abaixo listamos os resultados para as 5 primeiras iterações. Vemos que a
sequência converge, e atinge uma precisão de aproximadamente 5% em 5 iterações.
(k) (k) (k) (k) (k 1) (k)
k x1 x2 x3 max{|(xi xi )/xi |, i = 1, . . . , n}
0 1 1 1 -
1 2 2 13/4 9/13
2 15/16 5/2 5/2 3/10
3 7/8 63/62 93/32 17/63
4 133/128 31/16 393/128 7/131
5 519/512 517/256 767/256 21/517

5.13.2 Convergência do Método de Jacobi

(REVISAR)
Vamos escrever a matriz A como

A=L+D+U

onde L é a “lower triangular matrix” (sem diagonal); U “upper triangular matrix” (sem
diagonal); e D a matriz diagonal.
Desta forma,
Ax = (L + D + U )x = b
Dx = (L + U )x + b
1
x=D [ (L + U )x + b]
x = Jx + c
onde J = D 1 (L + U) e c = D 1 b.
Aplicando o método iterativo teremos
2 a12 a13 a1n 3
0 a11 a11 ... a11
6 a21 0 a23
... a2n 7
6 a22 a22 a22 7
6 . .. 7
x(k+1) = Jx + c, onde J = 6
(k) . . 7
6 . 0 7
6 .. .. 7
4 . 0 . 5
an1 an2 an,n 1
ann ann ... ann 0

Partindo de x(0) e fazendo sucessivamente a iteração temos

x(k) = J k x(0) + [1 + J + J 2 + . . . + J k
|{z}
1
]c (5.24)
elevado
ak
Para que convirja, requer que
lim J k = [0]
k!1

82
 O que implica que limk!1 [1 + J + J 2 + . . . + J k 1 ] = (1 J) 1 . Assim quando (5.24) é
satisfeita, x = limk!1 x(k) existe e x = 0 + (1 J) 1 c, isto é, (1 J)x = c ou x = Jx + c.
Mas a condição (5.24) é válida se e somente se todos os auto valores da matriz J forem
em módulo < 1.
Seja ⇢s = max | 1 |, | 2 |, . . . , | n | onde | i | são os autovalores da matrix J. ⇢s é também
chamado de raio espectral (“spectral radius”).
Então para atingir precisão p após k iterações devemos ter

p ln 10
⇢ks ⇡ 10 p
! k⇡
ln ⇢s

Assim se ⇢s estiver próximo de 1 a convergência será muito lenta. Existem métodos

de aceleração. Ver Quinney e NR seção 19.5.
Determinar os auto valores da matriz J requerirá outro algoritmo, em geral. Na
prática, muitas vezes é mais fácil testar numericamente a convergência.

Critério das linhas

Uma condição mais simples de convergência , porém apenas suficiente, é que o sistema
possua diagonal principal estritamente dominante, ou seja,
X
|aii | > |aij |, i = 1, . . . , n (5.25)
j=1
j6=i

que é chamado de critério das linhas. Note que por este critério, os elementos da diagonal
principal nunca podem ser nulos.

Exercı́cio: mostre que a matriz do sistema

2 32 3 2 3
4 2 1 x1 11
4 1 2 05 4x2 5 = 4 3 5
2 1 4 x3 16

satisfaz o critério das linhas.

Por ser um critério apenas suficiente, um sistema que não satisfaz o critério das linhas
pode convergir. Além disso, alterando a ordem das linhas ou colunas pode-se tornar um
sistema convergente em divergente e vice-versa.

5.13.3 Método de Gauss-Seidel

O método de Gauss-Seidel é muito semelhante ao método de Jacobi, mas em geral apre-
senta uma convergência mais rápida. Neste método, aproveita-se os valores já calculados
(k+1)
em uma iteração (ex:, x1 ) para a estimativa dos termos seguintes.
As relações de recorrência tomam a seguinte forma
8 (k+1) (k) (k) (k)
>
>
> x1 = a111 (b1 a12 x2 a13 x3 . . . a1n xn )
>
> (k+1) (k+1) (k) (k)
>
< x2
> = a122 (b2 a21 x1 a23 x3 . . . a2n xn )
(k+1) (k+1) (k+1) (k)
x3 = a122 (b3 a31 x1 a32 x2 . . . a2n xn )
>
>
>...
>
>
>
>
: (k+1) 1 (k+1) (k+1)
xn = ann (bn an1 x1 . . . an,n 1 xn 1 )

83
ou, de forma mais compacta,
0 1
i 1
X n
X
(k+1) 1 @ (k+1) (k)
xi = bi aij xj aij xj A .
aii
j=1 j=i+1

Exemplo: vamos considerar novamente o sistema

8
>
<4x1 + 2x2 + x3 = 11
x1 + 2x2 =3 .
>
:
2x1 + x2 + 4x3 = 16
As relações de recorrência, pelo método de Gauss-Seidel, são
(k+1) 11 1 (k) 1 (k)
x1 = x x ,
4 2 2 4 3
(k+1) 3 1 (k+1)
x2 = + x1 , (5.26)
2 2
(k+1) 1 (k+1) 1 (k+1)
x3 =4 x x .
2 1 4 2
Começando novamente com o vetor x(0) = [1, 1, 1] obtemos sucessivamente
0 1 0 1 0 1 0 1
2 29/32 1033/1024 1.0087
x(1) = @ 5/2 A , x(2) = @ 125/64 A , x(3) = @ 4095/2048 A ⇡ @1.9995A .
19/8 783/256 24541/8192 2.9957
Note que, neste exemplo, a taxa de convergência é muito maior.

5.13.4 Convergência do Método de Gauss-Seidel

O critério das linhas também pode ser aplicado ao método de Gauss-Seidel, mas, como no
método de Jacobi, trata-se apenas de uma condição suficiente.
Para o método de Gauss-Seidel existe um outro critério, menos restritivo que o critério
das linhas, chamado critério de Sassenfeld. Seja
M = max i,
1in

onde os i são definidos por

|a12 | + |a13 | + · · · + |a1n |
1 = ,
|a11 |
Pi 1 Pn
j=1 j |aij | + j=i+1 |aij |
i = .
|aii |
A condição M < 1 é suficiente para que as aproximações sucessivas pelo método de Gauss-
Seidel convirjam.
Muito importante: A convergência (ou não) dos métodos de Jacobi ou Gauss-Seidel
independe do vetor inicial escolhido.

Exercı́cio: use o método de Gauss-Seidel para resolver o sistema abaixo. Verifique se

a matriz satisfaz o critério das linhas e o critério da Sassenfeld.
2 32 3 2 3
10 2 2 1 x1 3
6 2 5 1 1 7 6x2 7 6 5 7
6 76 7 = 6 7
4 1 1/2 6 1 5 4x3 5 4 95
1 1 0 20 x4 17

84