CCI-22

Matemática Computacional

Carlos Alberto Alonso Sanches

CCI-22

3) Raízes de Sistemas Lineares

Eliminação de Gauss, Gauss-Jordan, Decomposição LU, Gauss-Jacobi, Gauss-Seidel

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

MMéétodos de resolutodos de resoluççãoão

� Para a resolução de um sistema linear de equações, há dois grupos de métodos:� Métodos diretos: a solução é obtida através da aplicação de um número finito de operações aritméticas� Regra de Cramer� Eliminação de Gauss e de Gauss-Jordan� Decomposição LU

� Métodos iterativos: a solução é obtida através de uma sequência de aproximações sucessivas, até se alcançar uma resposta que satisfaça a precisão exigida� Gauss-Jacobi� Gauss-Seidel

Sistemas de EquaSistemas de Equaçções Linearesões Lineares

� Forma geral:

nnnnnn

nn

nn

bxaxaxa

bxaxaxa

bxaxaxa

=+++

=+++

=+++

...

...

...

2211

22222121

11212111

M

nnnnnn

nn

nn

bxaxaxa

bxaxaxa

bxaxaxa

=+++

=+++

=+++

...

...

...

2211

22222121

11212111

M

onde:aaijij são os coeficientesxxii são as incógnitasbbii são os termos independentesnn é a ordem do sistema

� Forma matricial:

Ax = bAx = b

=

nnnn

n

n

aaa

aaaaaa

A

L

MOMM

L

L

21

22221

11211

=

nnnn

n

n

aaa

aaaaaa

A

L

MOMM

L

L

21

22221

11211

=

nb

bb

bM

2

1

=

nb

bb

bM

2

1

=

nx

xx

xM

2

1

=

nx

xx

xM

2

1

onde:onde:

ExemploExemplo

1542

2514

5542

321

321

321

−=++

=−+

=−+

xxxxxxxxx

1542

2514

5542

321

321

321

−=++

=−+

=−+

xxxxxxxxx

� Forma matricial:

−

=

−

−

1

2

5

.

542

514

542

3

2

1

xxx

� Forma geral:

CCáálculo das forlculo das forçças em uma trelias em uma treliççaa

� Um exemplo:1

23

4

5

6

7

8

9

10

11

12

13

1415

16

171

2

3

4

5

6

7

8

9

10

F1 F2 F3

Fh Fh

=−−−=

=++−=

=°++°−=

∑∑

∑

0

0

045cos45cos

531

541

541

faffaF

faffaF

fffF

y

x

aax 4342143421

� Condições de equilíbrio:� Na junção 2: � Na junção 3:

=+−=

=+−=

∑∑

0

0

31

62

fFF

ffF

y

x

� Idem para demais junções

� Gerará um sistema de ordem 17

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

Regra de Regra de CramerCramer

� A aplicação da regra de Cramer, em um sistema de ordem n, exige o cálculo de quantos determinantes?� n para os numeradores e 1 para o denominador

Tempo de processamentoTempo de processamento

� Número m de multiplicações, no caso de 17 equações:

18 det17 = 18 ( 17m + 17 det16 )= 18 ( 17m + 17 ( 16m + 16 det15 )) = 18 ( 17m + 17 ( 16m + 16 ( 15m + 15 det14 ))) = 18 ( 17m + 17 ( 16m + 16 ( 15m + 15 ( 14m + 14 (.... ( 3m + 3 ( 2m )....)))))

multiplicações

Lembrando:

Tempo de processamentoTempo de processamento

= 18 ( 17m + 17 ( 16m + 16 ( 15m + 15 ( 14m + 14 (.... ( 3m + 3 ( 2m )....)))))

= m ( 2 x 3 x 4 x 5 x .... x 17 x 18 ++ 3 x 4 x 5 x .... x 17 x 18 ++ 4 x 5 x .... x 17 x 18 + + 5 x .... x 17 x 18 + : :+ 16 x 17 x 18 ++ 17 x 18 )

= 18! (1 + (1/2!) + (1/3!) + ... + (1/16!) ) multiplicações

≈ 9,6 x 1015 multiplicações

Tempo de processamentoTempo de processamento

� Quantidade de multiplicações: ≈ 9,6 x 1015

� Utilizando um supercomputador atual:� 1011 multiplicações por segundo� Tempo gasto: 9,6 x 104 s ≈ 1 dia

� Se o sistema fosse de ordem 20, exigiria cerca de 28 anos de processamento nesse mesmo computador!

� Um algoritmo bem mais eficiente é o Método da Eliminação de Gauss, que gasta tempo O(n3)

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

MMéétodo da Eliminatodo da Eliminaçção de Gaussão de Gauss

� Objetivo� Transformação do sistema linear a ser resolvido

em um sistema linear triangular

� Operações válidas:� Troca da ordem das linhas

� Troca da ordem das colunas (com exceção dos termos independentes)

� Multiplicação de uma equação por um número real não nulo

� Substituição de uma equação por uma combinaçãolinear entre ela mesma e outra equação

Sistemas Lineares TriangularesSistemas Lineares Triangulares� Triangular superior:

� Triangular inferior:

=

nn

n

n

n

a

aaaaaaaaa

A

K

MOMMM

K

K

K

000

00

0

333

22322

1131211

=

nn

n

n

n

a

aaaaaaaaa

A

K

MOMMM

K

K

K

000

00

0

333

22322

1131211

=

nnnnn aaaa

aaaaa

a

A

K

MOMMM

K

K

K

321

333231

2221

11

0

00

000

=

nnnnn aaaa

aaaaa

a

A

K

MOMMM

K

K

K

321

333231

2221

11

0

00

000

ResoluResoluçção de um sistema triangularão de um sistema triangular

� Exemplo:

� Passos da resolução:

2

3154

354

3

3

43

=

=⋅−

=−

xx

xx

2

3154

354

3

3

43

=

=⋅−

=−

xx

xx

1

1122

12

2

2

432

−=

−=⋅−+

−=−+

xx

xxx

1

1122

12

2

2

432

−=

−=⋅−+

−=−+

xx

xxx

1

10125)1(43

10543

1

1

4321

=

−=+⋅−−⋅+

−=+−+

xx

xxxx

1

10125)1(43

10543

1

1

4321

=

−=+⋅−−⋅+

−=+−+

xx

xxxx

12

24 ==x 1

2

24 ==x

22

354

12

10543

4

43

432

4321

=

=−

−=−+

−=+−+

xxxxxxxxxx

22

354

12

10543

4

43

432

4321

=

=−

−=−+

−=+−+

xxxxxxxxxx

PassosPassos� Considere a matriz aumentada [Ab]:

[ ]

=

nnnnnn

n

n

baaaa

baaabaaa

Ab

321

222221

111211

MMOMM

K

K

[ ]

=

nnnnnn

n

n

baaaa

baaabaaa

Ab

321

222221

111211

MMOMM

K

K

� Passo 1: anular os coeficientes de x1 nas linhas L2 a Ln� Substituir a linha L2 pela combinação linear:

11

21211212 ,

aamondeLmL =⋅−

Linha L1Linha L2

Linha Ln

� Se a11 = 0, trocar L1 com Lk, onde ak1 ≠ 0� Se Lk não existir, então o sistema não tem solução

� Continuar analogamente para linhas Lj , 2 < j ≤ n� Passo i, 1 < i < n: anular os coeficientes de xi nas linhas Li+1 a Ln

Exemplo 1Exemplo 1

132

3344

532

321

321

321

−=+−

=−+

=−+

xxxxxxxxx

132

3344

532

321

321

321

−=+−

=−+

=−+

xxxxxxxxx

[ ]

−−

−

−

=

1132

3344

5132

Ab

2,11

212112122 ==⋅−=

aamLmLL 2,

11

212112122 ==⋅−=

aamLmLL [ ] [ ]

[ ]7120

513223344

2

2

−−−=

−⋅−−=

LL

1,11

313113133 ==⋅−=

aamLmLL 1,

11

313113133 ==⋅−=

aamLmLL [ ] [ ]

[ ]6260

513211132

3

3

−−=

−⋅−−−=

LL

[ ]

−−

−−−

−

=

6260

7120

5132

Ab[ ]

−−

−−−

−

=

6260

7120

5132

Ab

Exemplo 1Exemplo 1

[ ]

−−

−−−

−

=

6260

7120

5132

Ab[ ]

−−

−−−

−

=

6260

7120

5132

Ab

3,22

323223233 ==⋅−=

aamLmLL 3,

22

323223233 ==⋅−=

aamLmLL

[ ] [ ][ ]15500

712036260

3

3

=

−−−⋅−−−=

LL [ ]

−−−

−

=

15500

7120

5132

Ab[ ]

−−−

−

=

15500

7120

5132

Ab

=⇒=⇒=−+⇒=−⋅+

=⇒−=−−⇒−=−−

=⇒=

1225362532

273272

3155

111321

2232

33

xxxxxxxxxx

xx

Exemplo 2Exemplo 2

3814222

134311027

57524

321

321

321

=++

=−+

=++

xxxxxx

xxx

3814222

134311027

57524

321

321

321

=++

=−+

=++

xxxxxx

xxx

−=

3814222

134311027

575241

][Ab

[ ] [ ][ ]

[ ] ( ) [ ][ ]12101130860

5752411/223814222

1410140020

575241)1/27(134311027

3

13133

2

12122

−−−=

⋅−=⋅−=

−−=

⋅−−=⋅−=

LLmLL

LLmLL [ ] [ ]

[ ][ ] ( ) [ ]

[ ]12101130860

5752411/223814222

1410140020

575241)1/27(134311027

3

13133

2

12122

−−−=

⋅−=⋅−=

−−=

⋅−−=⋅−=

LLmLL

LLmLL

−−−

−−=

12101130860

1410140020

575241

][Ab

−−−

−−=

12101130860

1410140020

575241

][Ab

Nos cálculos a seguir, consideraremos F(10,3,-5,5):

Exemplo 2Exemplo 2

−−−

−−=

12101130860

1410140020

575241

][Ab

−−−

−−=

12101130860

1410140020

575241

][Ab

[ ] ( ) [ ][ ]618006130000

14101400202/8612101130860

3

23233

−−=

−−⋅−−−−−=⋅−=

LLmLL [ ] ( ) [ ]

[ ]618006130000

14101400202/8612101130860

3

23233

−−=

−−⋅−−−−−=⋅−=

LLmLL

−−

−−=

618006130000

1410140020

575241

][Ab

x3 = -61800/(-61300) = 1,01

x2 = [-1410 – (-1400)⋅1,01]/2 = 0,0

x1 = [57 - 52⋅1,01 - 4⋅0,0]/1 = 4,5

No entanto, a solução exata é:� x1 = 1� x2 = 1� x3 = 1

PivoteamentosPivoteamentos parcial e completoparcial e completo

� Pivôs pequenos geram multiplicadores grandes, que aumentam os erros de arredondamento...

� Uma simples alteração no método de Gauss é escolher como pivô o elemento de maior módulo:� em cada coluna (pivoteamento parcial)� dentre todos os elementos possíveis no processo de eliminação (pivoteamento completo): exige um maior esforço computacional

� Voltemos a resolver o exemplo anterior com precisão de 3 casas decimais, mas com pivoteamento parcial:

−

3814222

134311027

575241

−

3814222

134311027

575241

−

3814222

575241

134311027

−

3814222

575241

134311027

Exemplo 2 com Exemplo 2 com pivoteamentopivoteamento parcialparcial

]715.166,870[

]134311027[)27/22(]3814222[

]521,5207,00[

]134311027[)27/1(]575241[

3

13133

2

12122

−−=

−⋅−=⋅−=

−=

−⋅−=⋅−=

LLmLL

LLmLL

]715.166,870[

]134311027[)27/22(]3814222[

]521,5207,00[

]134311027[)27/1(]575241[

3

13133

2

12122

−−=

−⋅−=⋅−=

−=

−⋅−=⋅−=

LLmLL

LLmLL

−−

−

−

715,166,870

521,5207.00

134311027

−−

−

−

715,166,870

521,5207.00

134311027

−

−−

−

521,5207,00

715,166,870

134311027

−

−−

−

521,5207,00

715,166,870

134311027

]1,521,5200[

]715,166,870[)6,87/07,0(]521,5207,00[

3

23233

=

−−⋅−−=⋅−=

LLmLL

−−

−

1,521,5200

715,166,870

134311027

−−

−

1,521,5200

715,166,870

134311027 x3 = 52,1/52,1 = 1x2 = [-71-16,5⋅1]/(-87,6) = 0,999x1 = [134 – (-3)⋅1 – 110⋅0,999]/27 = 1,00

ExercExercííciocio

622

520002,0

21

21

=+

=+

xxxx

622

520002,0

21

21

=+

=+

xxxx

� Considere o sistema linear abaixo:

� Utilizando aritmética de três dígitos decimais, resolva-o através da eliminação de Gauss:� sem pivoteamento parcial� com pivoteamento parcial

� Confira os resultados encontrados

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

MMéétodo de Gausstodo de Gauss--JordanJordan

� Consiste em efetuar operações sobre as equações do sistema com a finalidade de transformá-lo em um sistema diagonal equivalente, isto é, são nulos todos os coeficientes aik, quando i≠k

=

nna

aa

a

A

K

KOMMM

L

K

K

000

000

000

000

33

22

11

=

nna

aa

a

A

K

KOMMM

L

K

K

000

000

000

000

33

22

11

ExemploExemplo

4232

2325

15

321

321

321

=++

=++

=++

xxxxxx

xxx

4232

2325

15

321

321

321

=++

=++

=++

xxxxxx

xxx[ ]

=

423223251151

Ab

[ ] [ ][ ]6040640L

2325511151LmLL

2

12122

,,,

)/(

=

⋅−=⋅−= [ ] [ ][ ]6040640L

2325511151LmLL

2

12122

,,,

)/(

=

⋅−=⋅−=

[ ] [ ][ ]2380220L

2325524232LmLL

3

13133

,,,

)/(

=

⋅−=⋅−= [ ] [ ][ ]2380220L

2325524232LmLL

3

13133

,,,

)/(

=

⋅−=⋅−=

[ ]

=

238022060406401325

Ab,,,

,,,[ ]

=

238022060406401325

Ab,,,

,,,

423211512325

ExemploExemplo

[ ]

=

238022060406401325

Ab,,,

,,,[ ]

=

238022060406401325

Ab,,,

,,,

[ ] [ ][ ]9132609000L

604064064222380220LmLL

3

23233

,,

,,,),/,(,,,

=

⋅−=⋅−=

[ ]

=

913260900060406401325

Ab,,

,,,[ ]

=

913260900060406401325

Ab,,

,,,

[ ] [ ][ ]31310640L

91326090006090406040640LmLL

2

32322

,,

,,),/,(,,,

−=

⋅−=⋅−=

ExemploExemplo

[ ] [ ][ ]5711305L

313106406421325L

1

1

,

,,),/(

=

−⋅−=

[ ]

−=

9132609000313106401325

Ab,,

,,[ ]

−=

9132609000313106401325

Ab,,

,,

[ ] [ ][ ]7812005L

9132609000609035711305L

1

1

,

,,),/(,

−=

⋅−= [ ] [ ][ ]7812005L

9132609000609035711305L

1

1

,

,,),/(,

−=

⋅−=

[ ]

−

−

=

9132609000313106407812005

Ab,,

,,

,A solução é:

� x1 = -2,556� x2 = -0,2854� x3 = 4,783

Outra aplicaOutra aplicaççãoão

� Uma variação do método de Gauss-Jordan pode ser utilizada para se encontrar a inversa de uma matriz A quadrada de ordem n

� Basta transformar a matriz A na correspondente matriz identidade, aplicando essas mesmas operações em uma matriz identidade de ordem n

[A|I] [I|A-1]Gauss-Jordan

Refinamento por resRefinamento por resííduosduos

� Se x(1) for encontrado como solução do sistema Ax = b, então o erro dessa solução é x – x(1)

� Multiplicando o erro por A:� A(x- x(1)) = Ax – Ax(1) = b – b(1) = r(1)

� O resíduo pode ser utilizado para se encontrar uma solução melhorada x(2):� x(2) = x(1) + δ(1), onde δ(1) é um vetor de correção� Ax(2) = b ⇔ A(x(1) + δ(1)) = b ⇔ Aδ(1) = b - Ax(1) = r(1)� δ(1) é solução do sistema Aδ = r(1)

� Esses cálculos permitem um processo de refinamento da solução do sistema Ax = b

resíduo

ExemploExemplo

3106x521x213x081x021880x411x523x88x353749x145x511x88x524

416x011x39x03x78

4321

4321

4321

4321

,,,,,

,,,,,

,,,,,

,,,,,

−=+−−

−=+−−

−=−+−

=+++

3106x521x213x081x021880x411x523x88x353749x145x511x88x524

416x011x39x03x78

4321

4321

4321

4321

,,,,,

,,,,,

,,,,,

,,,,,

−=+−−

−=+−−

−=−+−

=+++

� Vamos refinar o sistema abaixo:

� Através do método de Gauss, podemos encontrar a solução abaixo:

T1 001970981970x ],,,,[)( −=T1 001970981970x ],,,,[)( −=

� Cálculo do resíduo:

−

=−=

5940594021400420

Axbr 11

,

,

,

,

)()(

−

=−=

5940594021400420

Axbr 11

,

,

,

,

)()(

Não estábom...

ExemploExemplo� Cálculo do vetor de correção δ(1):

−

=

δ

δ

δ

δ

−−

−−

−−

5940594021400420

5212130810214115238835314551188524

011390378

4

3

2

1

,

,

,

,

,,,,

,,,,

,,,,

,,,,

−

=

δ

δ

δ

δ

−−

−−

−−

5940594021400420

5212130810214115238835314551188524

011390378

4

3

2

1

,

,

,

,

,,,,

,,,,

,,,,

,,,,

� Solução:

−=δ

00000029400195002950

1

,

,

,

,

)(

−=δ

00000029400195002950

1

,

,

,

,

)(

� Solução melhorada:

−=δ+=

00001999900000200001

xx 112

,

,

,

,

)()()(

−=δ+=

00001999900000200001

xx 112

,

,

,

,

)()()(

ExemploExemplo� Novo resíduo:

−

−

=−=

0130024001100090

Axbr 22

,

,

,

,

)()(

−

−

=−=

0130024001100090

Axbr 22

,

,

,

,

)()(

Melhor queo anterior

� Cálculo do novo vetor de correção:

−

−

−

=δ

00000000700002000020

2

,

,

,

,)(

−

−

−

=δ

00000000700002000020

2

,

,

,

,)(

� Outra solução melhorada:

−=

00001000010000200001

x 3

,

,

,

,

)(

−=

00001000010000200001

x 3

,

,

,

,

)(

� Novo resíduo:

=

0000

r 3 )(

=

0000

r 3 )(

Melhor aproximaMelhor aproximaççãoão� Dado um sistema Ax = b, sejam y e z duas aproximações da solução exata x. Como saber qual delas é a melhor?

� A estratégia mais lógica parece ser comparar os respectivos resíduos: o menor seria da melhor solução

� Infelizmente, isso nem sempre é verdade...� Exemplo:

=++

=++

=++

640x250x210x150520x240x160x120840x120x360x240

321

321

321

,,,,

,,,,

,,,,

=++

=++

=++

640x250x210x150520x240x160x120840x120x360x240

321

321

321

,,,,

,,,,

,,,,

� Conclusão: nem sempre a aproximação de menor resíduo é a melhor ou a mais exata

� Se a busca por resíduos menores não garante melhores soluções, como saber se o processo de refinamento por resíduos funciona?

−

−=

11425

y

−

−=

11425

y

−

=

043

z

−

=

043

z

=

080000000

ry,

,

,

=

080000000

ry,

,

,

=

250240120

rz,

,

,

=

250240120

rz,

,

,

−

=

143

x

−

=

143

x

Condicionamento de problemasCondicionamento de problemas� Um problema é dito mal condicionado se pequenas alterações nos dados de entrada ocasionam grandes erros no resultado final

� Exemplo:

=+

=+

0600y4210x48101190y8730x9920

,,,

,,,

=+

=+

0600y4210x48101190y8730x9920

,,,

,,,Solução: x=1 e y=-1

� Suponha que os valores desse sistema sejam obtidos experimentalmente, e por isso os termos independentes possam variar de ±0,001:

=+

=+

0600y4210x48101200y8730x9920

,,,

,,,

=+

=+

0600y4210x48101200y8730x9920

,,,

,,,Valor perturbado

Solução: x=0,815 e y=-0,789

Erro na entrada: (|0,119-0,120|/|0,119|) ≈ 0,8%

Erro no resultado: (|1,0-0,815|/|1,0|) ≈ 18,5%

Outro exemploOutro exemplo

� Considere os seguintes sistemas:

Solução: x=2 e y=3

=+

=+

50016y5014x5111y3x

,,,

=+

=+

50016y5014x5111y3x

,,,

=+

=+

50316y5014x5111y3x

,,,

=+

=+

50316y5014x5111y3x

,,,

Solução: x=10,28 e y=0,24

(a)(b)

(a)(c)

x

y

(a)

(c)

(b)

MMéétricas de condicionamentotricas de condicionamento

� Há métricas para o condicionamento de sistemas lineares, baseadas em normas de vetores e matrizes (vide Cláudio & Marins)

� No entanto, esses cálculos são difíceis...� Também é possível identificar o condicionamento de um sistema linear apenas com o uso dos refinamentos:� Se os resíduos r(1), r(2), ..., r(n) são pequenos, mas as correções δ(1), δ(2), ..., δ(n) são grandes, então o sistema é mal condicionado

� Para sistemas bem condicionados, bastam no máximo dois refinamentos (ou seja, δ(2) é muito pequeno)

� Ao longo desse processo, os resíduos e as correções devem ser calculados com precisão dupla

ExemploExemplo

� Resolução de Aδ(1) = r(1):

−=−+

=−+

=++

678900x47941x89652x6951204731x32531x958900x472510647000x62314x62351x47592

321

321

321

,,,,

,,,,

,,,,

−=−+

=−+

=++

678900x47941x89652x6951204731x32531x958900x472510647000x62314x62351x47592

321

321

321

,,,,

,,,,

,,,,

� Considere o sistema abaixo em F(10,5,-98,100):

−

−=

2441900717284061

x1

,

,

,)(

−

−=

2441900717284061

x1

,

,

,)(

� Primeiro refinamento em F(10,10,-98,100):

−

−

−

=

−

−

−

=−=

000076696000005537700001218010

678823304004735537710648218010

678900473106480

Axbr 11

,

,

,

,

,

,

,

,

,)()(

−

−

−

=

−

−

−

=−=

000076696000005537700001218010

678823304004735537710648218010

678900473106480

Axbr 11

,

,

,

,

,

,

,

,

,)()(

−

−

=δ

0000577650000025110000000422820

1

,

,

,)(

−

−

=δ

0000577650000025110000000422820

1

,

,

,)(

� Solução melhorada x(2) = x(1) + δ(1):

−

−=

2441900717284051

x 2

,

,

,)(

−

−=

2441900717284051

x 2

,

,

,)(

Resíduos pequenos

Correções pequenas

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

Uma outra forma de ver...Uma outra forma de ver...

� Consideremos o sistema de 3 equações Ax = b:

)(0

333231

232221

131211

Aaaaaaaaaa

A =

=)(0

333231

232221

131211

Aaaaaaaaaa

A =

=

=

3

2

1

bbb

b

=

3

2

1

bbb

b

=

3

2

1

xxx

x

=

3

2

1

xxx

x

� Após a primeira fase da eliminação de Gauss:

−

−==

=

10m01m001

MondeAMaa0aa0aaa

A

31

21000

133

132

123

122

113

112

111

1 )()()(

)()(

)()(

)()()(

)( ,.

−

−==

=

10m01m001

MondeAMaa0aa0aaa

A

31

21000

133

132

123

122

113

112

111

1 )()()(

)()(

)()(

)()()(

)( ,.

� Após a segunda fase da eliminação de Gauss:

−

==

=

1m0010001

MondeAMa00aa0aaa

A

32

111

233

223

222

213

212

211

2 )()()(

)(

)()(

)()()(

)( ,.

−

==

=

1m0010001

MondeAMa00aa0aaa

A

32

111

233

223

222

213

212

211

2 )()()(

)(

)()(

)()()(

)( ,.

Uma outra forma de ver...Uma outra forma de ver...� Resumindo:

� A = A(0)

� A(1) = M(0).A(0) = M(0).A� A(2) = M(1).A(1) = M(1).M(0).A� A = (M(1).M(0))-1.A(2)

� A = (M(0))-1.(M(1))-1.A(2)

� É fácil comprovar que:

=−−

1

01

001

)()(

3231

21

1)1(1)0(

mmmMM

=−−

1

01

001

)()(

3231

21

1)1(1)0(

mmmMM

� Portanto:

ULa00aa0aaa

1mm01m001

A233

223

222

213

212

211

3231

21 .)(

)()(

)()()(

=

= ULa00aa0aaa

1mm01m001

A233

223

222

213

212

211

3231

21 .)(

)()(

)()()(

=

=

DecomposiDecomposiçção LUão LU� A comprovação anterior pode ser generalizada em um teorema

⋅

==

nn

n333

n22322

n1131211

3n2n1n

3231

21

u000

uu00uuu0uuuu

1mmm001mm001m0001

ULA

K

MOMMM

K

K

K

K

OMMM

K

K

K

.

⋅

==

nn

n333

n22322

n1131211

3n2n1n

3231

21

u000

uu00uuu0uuuu

1mmm001mm001m0001

ULA

K

MOMMM

K

K

K

K

OMMM

K

K

K

.

� Dada uma matriz quadrada de ordem n, seja Ak a matriz constituída das primeiras k linhas e colunas de A. Suponha que det(Ak) ≠ 0, 0 ≤ k ≤ n-1. Então:� Existe uma única matriz triangular inferior L=(mij), com mii = 1, 1 ≤ i ≤ n. Os demais são os multiplicadores da Eliminação de Gauss

� Existe uma única matriz triangular superior U=(uij), tais que L.U = A

� det(A) = u11.u12. ... .unn

DecomposiDecomposiçção LUão LU� Portanto, dados o sistema linear Ax = b e a decomposição (ou fatoração) LU da matriz A, temos:� Ax = b ⇔ (L.U)x = b

� Chamando Ux = y, o sistema original passa a ser Ly = b, ou seja, surgem dois sistemas triangulares

� Por outro lado, é fácil verificar que y = L-1.b é o vetor b acumulando as operações da Eliminação de Gauss

� Por exemplo, no caso de um sistema com 3 equações:� Como L = (M(0))-1.(M(1))-1, então L-1 = M(1).M(0)

� Portanto, y = M(1).M(0).b� Vantagem da decomposição A = L.U: uma vez calculadas as matrizes L e U, resolvemos mais rapidamente qualquer sistema com a matriz A. Isso é útil, por exemplo, no refinamento por resíduos

ExemploExemplo

3x2x3x42x2xx1x4x2x3

321

321

321

=++

=++

=++

=

234211423

A

− 31031032310423

//

//

−4134323131423

/

///

=

11340131001

L/

/

−

=

40032310423

U //

− 3103134323131423

///

///

multiplicadores

ExemploExemplo

3x2x3x42x2xx1x4x2x3

321

321

321

=−+

=++

=++

=

11340131001

L/

/

−

=

40032310423

U //

3yyy342yy31

1y

321

21

1

=++

=+

=

/

/

=

0351

y /

−

=

053

x0x4

35x32x311x4x2x3

3

32

321

=−

=+

=++

///

Ax = b LUx = b

Ly = b

Ux = y

Outra aplicaOutra aplicaççãoão

� A decomposição LU também é útil no cálculo da matriz inversa

� Resolver o sistema AX = B, onde A, X e B são matrizes de ordem n, é o mesmo que resolver n sistemas Ax = b, onde x e b são vetores de tamanho n

� A inversa A-1 da matriz A pode ser encontrada através da resolução do sistema AX = I, onde I é a matriz identidade

� Nesse caso, basta realizar uma única vez a decomposição LU da matriz A, e depois utilizá-la na resolução de n sistemas

DecomposiDecomposiçção LU com pivoteamentoão LU com pivoteamento� É possível incorporar as estratégias de pivoteamento parcial ou completo à decomposição LU

� Uma matriz quadrada P de ordem n é uma matriz de permutação se for obtida da correspondente matriz identidade através de permutações em suas linhas ou colunas

� As eventuais permutações de linhas ou colunas na matriz A(k), obtida em um passo intermediário da Eliminação de Gauss, podem ser realizadas através da multiplicação por uma matriz de permutação

� Exemplo:

=

=

413562951

562951413

001100010

AP k .. )(

=

562951413

A k)(

Exemplo com pivoteamento parcialExemplo com pivoteamento parcial

2x3x43x2x2x9xx4x3

31

321

321

−=−

=++

=+−

−

−

=

304221143

A 0)(

−

−

=

413443411241304

A 1

//

//)(

=

001010100

P 0)(

−

==

143221304

APA 000 )()()(' .

=

010100001

P 1)(

−

−

==

411241413443304

APA 111

//

//. )()()('

Exemplo com pivoteamento parcialExemplo com pivoteamento parcial

−

−

−

=

8352141413443304

A2

///

//)(

−

=

121410143001

L//

/

−

−

=

8350041340304

U/

/

� L.U = A’ = P.A, onde P = P(1).P(0):

−=

−

−

==

221143304

304221143

010001100

APA ..'

� A’x = b’ ⇔ PAx = Pb ⇔ LUx = Pb

Exemplo com pivoteamento parcialExemplo com pivoteamento parcial

−

=

4352212

y/

/

−=

211

x

2x3x43x2x2x9xx4x3

31

321

321

−=−

=++

=+−

−

=

121410143001

L//

/

−

−

=

8350041340304

U/

/

−

=

− 239

010001100

yyy

121410143001

3

2

1

..

//

/

=

−

−

435221

2

xxx

8350041340304

3

2

1

/

/.

/

/

Ax = b LUx = Pb

Ly = Pb

Ux = y

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

MMéétodos iterativostodos iterativos

� Como foi inicialmente comentado, os métodos iterativos para resolução de sistemas lineares consistem em encontrar uma sequência de aproximações sucessivas

� Dada uma estimativa inicial x(0), calcula-se a sequência x(1), x(2), x(3) ..., até que determinado critério de parada seja satisfeito

� O sistema Ax = b é transformado em x(k) = Cx(k-1) + g, k>0, onde C é uma matriz e g um vetor

� Possíveis critérios de parada:� máximo erro absoluto ou relativo� número de iterações

� Principais métodos: Gauss-Jacobi e Gauss-Seidel

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

MMéétodo de Gausstodo de Gauss--JacobiJacobi

� Considere o sistema linear em sua forma inicial:

nnnn22n11n

2nn2222121

1nn1212111

bxaxaxa

bxaxaxa

bxaxaxa

=+++

=+++

=+++

...

...

...

MMOMM

nnnn22n11n

2nn2222121

1nn1212111

bxaxaxa

bxaxaxa

bxaxaxa

=+++

=+++

=+++

...

...

...

MMOMM

� Isolando a i-ésima incógnita na i-ésima equação:� x1 = (b1 - a12 x2 - ... - a1n x2)/a11� x2 = (b2 - a21 x1 - ... - a2n xn)/a22...� xn = (bn - an1 x1 - ... - an,n-1 xn-1)/ann

MMéétodo de Gausstodo de Gauss--JacobiJacobi

� Dessa forma, para x(k) = Cx(k-1) + g:

−−

−−

−−

=

0aaaa

aa0aaaaaa0

C

nn2nnn1n

22n22221

11n11112

L

MOMM

L

L

//

//

//

−−

−−

−−

=

0aaaa

aa0aaaaaa0

C

nn2nnn1n

22n22221

11n11112

L

MOMM

L

L

//

//

//

=

nnn

222

111

ab

abab

g

/

/

/

M

� Exemplos de critérios de parada:� Erro absoluto: d(k) = maxi |x(k) – x(k-1)| < ε� Erro relativo: dr(k) = d(k)/(maxi |x(k)|) < ε

ExemploExemplo

−=

606170

x 0

,

,

,)(ε = 0,05

−−

−−

−−

=

010351510511011020

C//

//

//

=

10658107

g/

/

/

6x10x3x28xx5x7xx2x10

321

321

321

=++

=++

=++

−=+=

940861960

gCxx 01

,

,

,)()(

|x1(1) – x1(0)| = 0,26

|x2(1) – x2(0)| = 0,26

|x3(1) – x3(0)| = 0,34

dr(1) = 0,34/(max |xi(1)|) = 0,1828 > ε

ExemploExemplo

−−

−−

−−

=

010351510511011020

C//

//

//

=

10658107

g/

/

/

−=

940861960

x 1

,

,

,)(

−=+=

96609819780

gCxx 12

,

,

,)()( dr(2) = 0,12/1,98 = 0,0606 > ε

dr(3) = 0,0324/1,9888 = 0,0163 < ε

−=+=

998409888199940

gCxx 23

,

,

,)()(

CritCritéério das linhasrio das linhas

� Em um método iterativo, a convergência para a solução exata não é garantida: é preciso que o sistema satisfaça alguns requisitos

� Há uma condição suficiente para a convergência do Método de Gauss-Jacobi, conhecido como o critério das linhas :

n21iparaaa ii

n

ij1j

ij ,...,,, =<∑≠=

ExemplosExemplos

6x10x3x28xx5x7xx2x10

321

321

321

=++

=++

=++

� Considere o exemplo anterior:2+1 < 101+1 < 5

2+3 < 10

Garantia de convergência

3x3x3xx

21

21

−=−

=+

� Considere o exemplo abaixo:

1 = 11 < 3

Não há garantia de convergência

� No entanto, o método de Gauss-Jacobi converge neste sistema para a solução exata x1 = x2 = 3/2. Verifique!

� Isso mostra que o critério das linhas é suficiente, mas não necessário

DemonstraDemonstraççãoão� Sejam:

� x* = [x1, x2, ..., xn]T: a solução exata de Ax = b� x(k) = [x(k)1, x(k)2, ..., x(k)n]T: a k-ésima aproximação de x*� e(k) = x(k) – x*: o erro na k-ésima aproximação

� Queremos garantir que limk→∞ e(k)i = 0, 1≤i≤n� No Método de Gauss-Jacobi, pode-se constatar que:

� e(k+1)1 = -(a12e(k)2 + a13e(k)3 + ... + a1ne(k)n)/a11� e(k+1)2 = -(a21e(k)1 + a23e(k)3 + ... + a2ne(k)n)/a22� e(k+1)n = -(an1e(k)1 + an2e(k)2 + ... + an(n-1)e(k)n-1)/ann

� Sejam:� E(k) = máx1≤i≤n{|e(k)i|}� αi = (|ai1| + ... + |ai(i-1)| + |ai(i+1)| + ... + |ain|)/|aii|, 1≤i≤nQuando o critério das linhas é satisfeito, αi < 1

DemonstraDemonstraçção (continuaão (continuaçção)ão)

� Quando k→∞, x(k)→x* é equivalente a E(k)→0� Demonstraremos que E(k+1) ≤ α.E(k), onde α = máx1≤i≤n{αi}� Para 1 ≤ i ≤ n:

� e(k+1)i = -(ai1e(k)1 + ... + ai(i-1)e(k)i-1 + ai(i+1)e(k)i+1 + ... + aine(k)n)/aii� |e(k+1)i| ≤ (|ai1|.|e(k)1| + ... + |ai(i-1)|.|e(k)i-1| + |ai(i+1)|.|e(k)i+1| + ... + |ain|.|e(k)n|)/|aii|

� |e(k+1)i| ≤ (|ai1| + ... + |ai(i-1)| + |ai(i+1)| + ... + |ain|).E(k)/|aii|� |e(k+1)i| ≤ αi.E(k)

� Portanto, E(k+1) ≤ α.E(k)

� Consequentemente, E(k+1)/E(k) ≤ αComo α<1, então E(k)→0 quando k→∞: há convergência!

Mais um exemploMais um exemplo

6x8x63x2x2x52xx3x

32

321

321

−=+

=++

−=++

� Considere o sistema a seguir:3+1 > 15+2 > 2

6 < 8

Não há garantia de convergência

� No entanto, uma permutação entre as duas primeiras linhas garante a convergência:

6x8x62xx3x3x2x2x5

32

321

321

−=+

−=++

=++ 2+2 < 51+1 < 3

6 < 8

Garantia de convergência

� Quando o critério das linhas não for satisfeito, convém tentar uma permutação de linhas e/ou colunas

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

MMéétodo de Gausstodo de Gauss--SeidelSeidel

� Analogamente ao Método de Gauss-Jacobi, calcula-se x(k) = Cx(k-1) + g:

−−

−−

−−

=

0aaaa

aa0aaaaaa0

C

nn2nnn1n

22n22221

11n11112

L

MOMM

L

L

//

//

//

−−

−−

−−

=

0aaaa

aa0aaaaaa0

C

nn2nnn1n

22n22221

11n11112

L

MOMM

L

L

//

//

//

=

nnn

222

111

ab

abab

g

/

/

/

M

� No entanto, utiliza-se no cálculo de :� valores calculados na mesma iteração:� valores da iteração anterior:

)1( +kjx

)1(

1

)1(

1 ,..., +−

+ kj

k xx)()(

1,...,kn

kj xx +

ExemploExemplo

=

000

x 0)(ε = 0,05

0x6x3x36xx4x35xxx5

321

321

321

=++

=++

=++

� Processo iterativo:

)()()(

)()()(

)()()(

,,

,,,

,,

1k2

1k1

1k3

k3

1k1

1k2

k3

k2

1k1

x50x50xx250x75051x

x20x201x

+++

++

+

−−=

−−=

−−=

ExemploExemplo

� Primeira iteração (k=0):

)()()(

)()()(

)()()(

,,

,,,

,,

1k2

1k1

1k3

k3

1k1

1k2

k3

k2

1k1

x50x50xx250x75051x

x20x201x

+++

++

+

−−=

−−=

−−=

875075050150x7500175051x

1001x

13

12

11

,,.,.,

,.,,

)(

)(

)(

−=−−=

=−−=

=−−=

|x1(1) – x1(0)| = 1

|x2(1) – x2(0)| = 0,75

|x3(1) – x3(0)| = 0,875

dr(1) = 1/(max |xi(1)|) = 1 > ε

−

=

87507501

x 1

,

,)(

ExemploExemplo

� Segunda iteração (k=1):

)()()(

)()()(

)()()(

,,

,,,

,,

1k2

1k1

1k3

k3

1k1

1k2

k3

k2

1k1

x50x50xx250x75051x

x20x201x

+++

++

+

−−=

−−=

−−=

9875095050025150x9508750250025175051x

0251875020750201x

23

22

21

,,.,,.,

,),.(,,.,,

,),.(,,.,

)(

)(

)(

−=−−=

=−−−=

=−−−=

|x1(2) – x1(1)| = 0,025

|x2(2) – x2(1)| = 0,20

|x3(2) – x3(1)| = 0,1125

dr(2) = 0,2/(max |xi(2)|) = 0,1951 > ε

−

=

987509500251

x 2

,

,

,)(

ExemploExemplo

� Terceira iteração (k=2):

)()()(

)()()(

)()()(

,,

,,,

,,

1k2

1k1

1k3

k3

1k1

1k2

k3

k2

1k1

x50x500xx250x75051x

x20x201x

+++

++

+

−−=

−−=

−−=

9993099120500075150x99120987502500075175051x

007519875020950201x

33

32

31

,,.,,.,

,),.(,,.,,

,),.(,,.,

)(

)(

)(

−=−−=

=−−−=

=−−−=

|x1(3) – x1(2)| = 0,0175

|x2(3) – x2(2)| = 0,0412

|x3(3) – x3(2)| = 0,0118

dr(3) = 0,0412/(max |xi(3)|) = 0,0409 < ε

−

=

999309912000751

x 3

,

,

,)(

InterpretaInterpretaçção geomão geoméétricatrica

� No caso de um sistema linear de ordem 2, épossível visualizar a convergência do método:

x1

x2

x*

� Os pontos (x1(k+1), x2(k)) satisfazem a primeira equação, enquanto os pontos (x1(k+1), x2(k+1)) satisfazem a segunda

x1

x2

x*� No mesmo sistema, a convergência pode não ocorrer...

CritCritéério de Sassenfeldrio de Sassenfeld� Sejam os seguintes valores:

� Se β < 1, então o Método de Gauss-Seidel gera uma sequência convergente, qualquer que seja x(0)� Quanto menor for β, mais rápida será a convergência

∑=

⋅=βn

2jj1

111 a

a1

+β⋅⋅=β ∑∑

+=

−

=

n

1ijij

1i

1jjij

iii aa

a1

, para 1 < i ≤ n

β = max {βj}, 1 ≤ j ≤ n

ExemploExemplo

010x4x80x21x4001x20xx20x1087x30x60x3x60

40x20x20xx2

4321

4321

4321

4321

,,,,

,,,,

,,,,

,,,

−=+++

=++−−

−=−−+

=+−+

( )

( )

( )

( )

170

2736035808044021704041

35802044020701011

4403060706031

702020121

4

3

2

1

<=β

=⋅+⋅+⋅⋅=β

=+⋅+⋅⋅=β

=++⋅⋅=β

=++⋅=β

,

,,,,,,,

,,,,,,

,,,,,

,,,

DemonstraDemonstraççãoão� Sejam:

� x* = [x1, x2, ..., xn]T: a solução exata de Ax = b� x(k) = [x(k)1, x(k)2, ..., x(k)n]T: a k-ésima aproximação de x*� e(k) = x(k) – x*: o erro na k-ésima aproximação

� Queremos garantir que limk→∞ e(k)i = 0, 1≤i≤n� No Método de Gauss-Seidel, pode-se constatar que:

� e(k+1)1 = -(a12e(k)2 + a13e(k)3 + ... + a1ne(k)n)/a11� e(k+1)2 = -(a21e(k+1)1 + a23e(k)3 + ... + a2ne(k)n)/a22� e(k+1)n = -(an1e(k+1)1 + an2e(k+1)2 + ... + an(n-1)e(k+1)n-1)/ann

� Sejam:� E(k) = máx1≤i≤n{|e(k)i|}� β1 = (|a12| + |a13| + ... + |a1n|)/|a11|� βi = (β1.|ai1| + ... + βi-1.|ai(i-1)| + |ai(i+1)| + ... + |ain|)/|aii|, 1<i≤n

DemonstraDemonstraçção (continuaão (continuaçção)ão)� Quando k→∞, x(k)→x* é equivalente a E(k)→0� Demonstraremos por indução em i (1≤i≤n) que E(k+1) ≤ β.E(k), onde β = máx1≤i≤n{βi}

� Base (i=1): � |e(k+1)1| ≤ (|a12|.|e(k)2| + |a13|.|e(k)3| + ... + |a1n|.|e(k)n|)/|a11|� |e(k+1)1| ≤ [(|a12| + |a13| + ... + |a1n|)/|a11|]. máx1≤i≤n{|e(k)i|}� |e(k+1)1| ≤ β1.E(k) ≤ β.E(k)

� Hipótese: |e(k+1)i-1| ≤ βi-1.E(k) ≤ β.E(k)

� Passo:� |e(k+1)i| ≤ (|ai1|.|e(k+1)1| + ... + |ai(i-1)|.|e(k+1)i-1| + |ai(i+1)|.|e(k)i+1| + ... + |ain|.|e(k)n|)/|aii|

� |e(k+1)i| ≤ (|ai1|.β1 + ... + |ai(i-1)|.βi-1 + |ai(i+1)| + ... + |ain|).E(k)/|aii|� |e(k+1)i| ≤ βi.E(k) ≤ β.E(k)

� Portanto, E(k+1)/E(k) ≤ β� Como β<1, então E(k)→0 quando k→∞: há convergência!

ExemplosExemplos

3x3x3xx

21

21

−=−

=+

� Considere o sistema abaixo, anteriormente visto:

� No entanto, o Método de Gauss-Seidel converge neste sistema para a solução exata x1 = x2 = 3/2. Verifique!

� Isso mostra que o critério de Sassenfeld, como o critério das linhas, é suficiente, mas não necessário

( )1

31311111

2

1

=β

==β

==β

//.

/

18x2x623xx10

21

21

=−

=+

� Considere outro sistema:

( )130

30210610101

2

1

<=β

==β

==β

,

,/,.

,/

� Neste caso, o critério de Sassenfeld garante a convergência, mas o critério das linhas, não...

Não há garantia de convergência

CCICCI--2222

� Introdução

� Métodos diretos� Regra de Cramer

� Eliminação de Gauss

� Gauss-Jordan

� Decomposição LU

� Métodos iterativos� Gauss-Jacobi

� Gauss-Seidel

� Considerações finais

RelaRelaçção entre os critão entre os critéériosrios� Se um sistema satisfaz o critério das linhas, então também satisfará o critério de Sassenfeld

� Demonstração:� Seja α = máx1≤i≤n{αi} < 1, onde αi = (|ai1| + ... + |ai(i-1)| + |ai(i+1)| + ... + |ain|)/|aii|

� Vamos provar por indução em i que βi ≤ αi < 1, 1≤i≤n� Base (i=1):

� β1 = (|a12| + |a13| + ... + |a1n|)/|a11| = α1 < 1

� Hipótese: βi ≤ αi < 1, 1≤i<n� Passo (1≤i≤n):

� βi = (β1.|ai1| + ... + βi-1.|ai(i-1)| + |ai(i+1)| + ... + |ain|)/|aii|� βi ≤ (|ai1| + ... + |ai(i-1)| + |ai(i+1)| + ... + |ain|)/|aii| = αi < 1

� Portanto, α < 1 ⇒ β < 1

� A volta nem sempre é verdadeira...

ConsideraConsideraçções finaisões finais

� Tanto o critério das linhas como o critério de Sassenfeld são condições suficientes para a convergência, mas não necessárias

� Em sistemas esparsos (com grande número de coeficientes nulos), o Método da Eliminação de Gauss não é apropriado, pois não preserva esta vantajosa qualidade. Nesses casos, convém utilizar métodos iterativos

� Os métodos iterativos são menos suscetíveis ao acúmulo de erros de arredondamento

MMéétodos diretos todos diretos versusversus iterativositerativos

� Convergência� Diretos: não faz sentido considerar essa questão, pois calculam a solução exata

� Iterativos: ocorre sob determinadas condições� Esparsidade da matriz de coeficientes

� Diretos: alteram a estrutura da matriz� Iterativos: utilizam sempre a matriz inicial

� Erros de arredondamento� Diretos: ocorrem a cada etapa e acumulam-se� Iterativos: somente os erros da última etapa afetam a solução