Notas de Aula M¶etodos Num¶ericos para Equa»c~oes ...rodney/notas_de_aula/metodos_numericos_EDPs... · encial parcial, procedemos µa discretiza»c~ao desta. Usando diferen»cas

Notas de Aula

Metodos Numericos para EquacoesDiferenciais Parciais Elıpticas

Rodney Josue Biezuner 1

Departamento de MatematicaInstituto de Ciencias Exatas (ICEx)

Universidade Federal de Minas Gerais (UFMG)

Notas de aula do curso Topicos em Analise: Metodos Numericos para EDPs Elıpticas do Programa

de Pos-Graduacao em Matematica, ministrado durante o primeiro semestre do ano de 2007.

15 de junho de 2007

1E-mail: [email protected]; homepage: http://www.mat.ufmg.br/∼rodney.

Sumario

1 Metodo de Diferencas Finitas 31.1 O Caso Unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.1.1 Series de Taylor e Diferencas Finitas em Uma Dimensao . . . . . . . . . . . . . . . . . 31.1.2 Discretizacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.1.3 Resolucao Numerica do Problema de Autovalor Unidimensional . . . . . . . . . . . . . 5

1.2 O Caso Bidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2.1 A Formula dos Cinco Pontos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2.2 Existencia e Unicidade da Solucao Discreta – Autovalores do Problema Bidimensional 101.2.3 Princıpio do Maximo Discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.2.4 Convergencia da Solucao Discreta para a Solucao Classica . . . . . . . . . . . . . . . . 15

1.3 Discretizacoes de Ordem Superior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.3.1 Caso Unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.3.2 Caso Bidimensional: A Formula dos Nove Pontos Compacta . . . . . . . . . . . . . . 20

1.4 Diferencas Finitas em Coordenadas Polares . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251.5 Domınios Arbitrarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 291.6 Exercıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2 Existencia e Unicidade de Solucoes Discretas 332.1 Normas Matriciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.2 Matrizes Diagonalmente Dominantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.3 Teorema dos Discos de Gershgorin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.4 Propriedade FC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 402.5 Matrizes Irredutıveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.6 Invertibilidade de Matrizes de Discretizacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

2.6.1 Esquemas de Diferencas Finitas para o Intervalo e para o Retangulo . . . . . . . . . . 482.6.2 Esquema de Coordenadas Polares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 482.6.3 Esquema de Shortley-Weller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3 Metodos Iterativos para a Resolucao de Sistemas Lineares 503.1 Metodos Iterativos Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.1.1 Metodo de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.1.2 Metodo de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.1.3 Metodo SOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.1.4 Comparacao da Velocidade de Convergencia dos Tres Metodos . . . . . . . . . . . . . 533.1.5 Metodo de Jacobi Amortecido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3.2 Analise de Convergencia dos Metodos Iterativos Lineares . . . . . . . . . . . . . . . . . . . . . 553.2.1 Convergencia dos Metodos Iterativos Lineares . . . . . . . . . . . . . . . . . . . . . . . 563.2.2 Velocidade de Convergencia dos Metodos Iterativos Lineares . . . . . . . . . . . . . . 583.2.3 Convergencia para Matrizes Simetricas Positivas Definidas . . . . . . . . . . . . . . . . 60

1


3.3 Convergencia dos Metodos Iterativos Lineares para as Matrizes de Discretizacao . . . . . . . 613.3.1 Convergencia do Metodo de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.3.2 Convergencia do Metodo de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . 633.3.3 Convergencia do Metodo SOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.3.4 Convergencia do Metodo de Jacobi Amortecido . . . . . . . . . . . . . . . . . . . . . . 733.3.5 Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

3.4 Metodo do Gradiente Conjugado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 743.4.1 Metodos de Descida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 753.4.2 Metodo da Descida Mais Acentuada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 773.4.3 Metodo do Gradiente Conjugado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

3.5 Convergencia do Metodo do Gradiente Conjugado . . . . . . . . . . . . . . . . . . . . . . . . 82

4 Metodos Multigrid 854.1 A Malha de Multigrid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.2 Frequencias Altas e Baixas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.3 Suavizacao do Erro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

4.3.1 Metodo de Jacobi Amortecido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.4 O Ciclo de Duas Malhas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 914.5 O Ciclo Multigrid: Ciclos V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

5 Metodo dos Volumes Finitos 945.1 Leis de Conservacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

5.1.1 Lei de Conservacao Unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 945.1.2 Lei de Conservacao em Varias Dimensoes . . . . . . . . . . . . . . . . . . . . . . . . . 965.1.3 Relacoes Constitutivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

5.2 O Caso Unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.3 O Caso Bidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.4 Linearizacao do Termo Fonte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5.4.1 Termo Fonte do Tipo f (u) = Au + B com A < 0 . . . . . . . . . . . . . . . . . . . . . 1065.4.2 Termo Fonte do Tipo f (u) = Au + B com A > 0 . . . . . . . . . . . . . . . . . . . . . 1065.4.3 Termo Fonte do Tipo f (u) com f ′ (u) < 0 . . . . . . . . . . . . . . . . . . . . . . . . . 106

Capıtulo 1

Metodo de Diferencas Finitas

1.1 O Caso Unidimensional

Nesta secao, desenvolveremos um metodo numerico de diferencas finitas para resolver o problema de Dirichletpara a equacao de Poisson em uma dimensao

−u′′ = f (x) em [0, L] ,u (0) = a, u (L) = b.

1.1.1 Series de Taylor e Diferencas Finitas em Uma Dimensao

Seja ∆x > 0. Considere as seguintes expansoes de Taylor de uma funcao u em torno de um ponto x0,respectivamente a direita e a esquerda de x0:

u(x0 + ∆x) = u(x0) + u′(x0)∆x +12!

u′′(x0)∆x2 +13!

u′′′(x0)∆x3 + . . . , (1.1)

u(x0 −∆x) = u(x0)− u′(x0)∆x +12!

u′′(x0)∆x2 − 13!

u′′′(x0)∆x3 + . . . (1.2)

Daı,

u′(x0) =u(x0 + ∆x)− u(x0)

∆x− 1

2!u′′(x0)∆x− 1

3!u′′′(x0)∆x2 − . . . ,

u′(x0) =u(x0)− u(x0 −∆x)

∆x+

12!

u′′(x0)∆x− 13!

u′′′(x0)∆x2 + . . .

Isso fornece duas aproximacoes possıveis para a primeira derivada u′(x0) de u em x0:

u′(x0) ≈ u(x0 + ∆x)− u(x0)∆x

, (1.3)

u′(x0) ≈ u(x0)− u(x0 −∆x)∆x

. (1.4)

A primeira e chamada uma diferenca progressiva e a segunda e uma diferenca regressiva. Pela Formulade Taylor com Resto, o erro destas aproximacoes e dado por

ε = ±12u′′(ξ)∆x = O(∆x),

onde x0 6 ξ 6 x0 + ∆x no primeiro caso, e x0 −∆x 6 ξ 6 x0 no segundo caso.

3


Por outro lado, se subtrairmos (1.2) de (1.1), obtemos

u′(x0) =u(x0 + ∆x)− u(x0 −∆x)

2∆x− 1

3!u′′′(x0)∆x2 − 1

5!u(5)(x0)∆x4 − . . .

o que da uma outra aproximacao possıvel para a primeira derivada u′(x0) de u em x0:

u′(x0) ≈ u(x0 + ∆x)− u(x0 −∆x)2∆x

(1.5)

com erroε = −1

6u′′′(ξ)∆x2 = O(∆x2),

para algum x0−∆x 6 ξ 6 x0 +∆x. Esta aproximacao por diferenca finita e chamada diferenca centrada.Ela e uma melhor aproximacao que as aproximacoes laterais (progressiva e regressiva).

Se, ao inves, adicionarmos (1.1) e (1.2), obtemos

u′′(x0) =u(x0 + ∆x) + u(x0 −∆x)− 2u(x0)

∆x2− 2

4!u(4)(x0)∆x2 − 2

5!u(6)(x0)∆x4 − . . .

o que fornece uma aproximacao para a derivada segunda u′′(x0) de u em x0:

u′′(x0) ≈ u(x0 + ∆x) + u(x0 −∆x)− 2u(x0)∆x2

(1.6)

com erroε = − 1

12u(4)(ξ)∆x2 = O(∆x2),

onde x0 − ∆x 6 ξ 6 x0 + ∆x. Esta aproximacao e tambem chamada uma diferenca centrada para aderivada segunda.

1.1.2 Discretizacao

Dividimos o intervalo [0, L] em n subintervalos de comprimento ∆x = L/n atraves de n−1 pontos interioresuniformemente espacados:

x0 = 0, x1 = ∆x, x2 = 2∆x, . . . , xn−1 = (n− 1)∆x, xn = n∆x = L,

de modo que [0, L] = [x0, x1] ∪ [x1, x2] ∪ . . . ∪ [xn−1, xn]. Introduzimos a notacao:

ui = u(xi),fi = f (xi) .

Esta e uma discretizacao uniforme do intervalo [0, L]. Uma vez discretizado o domınio da equacao difer-encial parcial, procedemos a discretizacao desta. Usando diferencas centradas para cada ponto interior xi,1 6 i 6 n− 1, temos

−ui−1 + 2ui − ui+1

∆x2= fi. (1.7)

Para os pontos de fronteira, a condicao de Dirichlet implica

u0 = a e un = b. (1.8)

Portanto, para encontrar a solucao discretizada temos que resolver o sistema linear com n − 1 equacoes an− 1 incognitas:

∆x−2 (2u1 − u2) = f1 + a∆x−2

∆x−2 (−u1 + 2u2 − u3) = f2

...∆x−2 (−un−3 + 2un−2 − un−1) = fn−2

∆x−2 (−un−2 + 2un−1) = fn−1 + b∆x−2

,


ou seja,

1∆x2

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 2

u1

u2

...

...un−2

un−1

=

f1 + a∆x−2

f2

...

...fn−2

fn−1 + b∆x−2

.

Esta e uma matriz tridiagonal simetrica, esparsa. Alem disso, como veremos na proxima subsecao, ela epositiva definida (isto e, seus autovalores sao positivos) e portanto possui uma inversa, o que garante aexistencia e unicidade da solucao. Dada sua simplicidade, ela pode ser resolvida por eliminacao gaussianaou sua inversa pode ser efetivamente calculada. Por exemplo, para n = 4, 5, 6 temos

2 −1 0−1 2 −1

0 −1 2

−1

=

1 12

13

0 1 23

0 0 1

12 0 00 2

3 00 0 3

4

1 0 012 1 013

23 1

=

14

3 2 12 4 21 2 3

,

2 −1 0 0−1 2 −1 0

0 −1 2 −10 0 −1 2

−1

=

1 12

13

14

0 1 23

24

0 0 1 34

0 0 0 1

12 0 0 00 2

3 0 00 0 3

4 00 0 0 4

5

1 0 0 012 1 0 013

23 1 0

14

24

34 1

=

15

4 3 2 13 6 4 22 4 6 31 2 3 4

2 −1 0 0 0−1 2 −1 0 0

0 −1 2 −1 00 0 −1 2 −10 0 0 −1 2

−1

=

1 12

13

14

15

0 1 23

24

25

0 0 1 34

35

0 0 0 1 45

0 0 0 0 1

12 0 0 0 00 2

3 0 0 00 0 3

4 0 00 0 0 4

5 00 0 0 0 5

6

1 0 0 0 012 1 0 0 013

23 1 0 0

14

12

34 1 0

15

25

35

45 1

=16

5 4 3 2 14 8 6 4 23 6 9 6 32 4 6 8 41 2 3 4 5

.

A forma da inversa no caso geral pode ser facilmente adivinhada.

1.1.3 Resolucao Numerica do Problema de Autovalor Unidimensional

Os autovalores de Dirichlet do laplaciano em [0, L] devem ser aproximados pelos autovalores da matriz(n− 1)× (n− 1)

A =1

∆x2

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 2

quando n →∞ e correspondentemente ∆x → 0.Lembrando que as autofuncoes de Dirichlet do laplaciano no intervalo [0, L] sao as funcoes

Uj (x) = senjπx

L,


este fato sugere que os autovetores uj da matriz A sao os vetores de coordenadas

Uj (x1) , Uj (x2) , . . . , Uj (xn−2) , Uj (xn−1) = Uj (∆x) , Uj (2∆x) , . . . , Uj ((n− 2) ∆x) , Uj ((n− 1) ∆x) ,

ou seja, como ∆x = L/n, os vetores12

sinθ

2= cos θ

uj =(

senjπ

n, sen

2jπ

n, . . . , sen

(n− 2) jπ

n, sen

(n− 1) jπ

n

).

Usando identidades trigonometricas, vamos verificar que isso de fato acontece:

1.1 Lema. Os n− 1 autovalores da matriz A sao

λj =2

∆x2

(1− cos

jπ

n

)=

4∆x2

sen2 jπ

2n, j = 1, . . . , n− 1, (1.9)

e os autovetores correspondentes sao

uj =(

senjπ

n, sen

2jπ

n, . . . , sen

(n− 2) jπ

n, sen

(n− 1) jπ

n

), j = 1, . . . , n− 1. (1.10)

Prova. Temos

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 2

senjπ

n

sen2jπ

n...

sen(n− 2) jπ

n

sen(n− 1) jπ

n

=

2 senjπ

n− sen

2jπ

n

− senjπ

n+ 2 sen

2jπ

n− sen

3jπ

n...

− sen(n− 3) jπ

n+ 2 sen

(n− 2) jπ

n− sen

(n− 1) jπ

n

− sen(n− 2) jπ

n+ 2 sen

(n− 1) jπ

n

= 2(

1− cosjπ

n

)

senjπ

n

sen2jπ

n...

sen(n− 2) jπ

n

sen(n− 1) jπ

n

,

pois

2 senjπ

n− sen

2jπ

n= 2 sen

jπ

n− 2 sen

jπ

ncos

jπ

n= 2

(1− cos

jπ

n

)sen

jπ

n,


− sen(n− k − 1) jπ

n+ 2 sen

(n− k) jπ

n− sen

(n− k + 1) jπ

n

= − sen[(n− k) jπ

n− jπ

n

]+ 2 sen

(n− k) jπ

n− sen

[(n− k) jπ

n+

jπ

n

]

= − sen(n− k) jπ

ncos

jπ

n+ cos

(n− k) jπ

nsen

jπ

n+ 2 sen

(n− k) jπ

n

− sen(n− k) jπ

ncos

jπ

n− cos

(n− k) jπ

nsen

jπ

n

= 2(

1− cosjπ

n

)sen

(n− k) jπ

n,

e

− sen(n− 2) jπ

n+ 2 sen

(n− 1) jπ

n

= − sen[(n− 1) jπ

n− jπ

n

]+ 2 sen

(n− 1) jπ

n

= − sen(n− 1) jπ

ncos

jπ

n+ cos

(n− 1) jπ

nsen

jπ

n+ 2 sen

(n− 1) jπ

n

= − sen(n− 1) jπ

ncos

jπ

n− sen

(n− 1) jπ

ncos

jπ

n+ 2 sen

(n− 1) jπ

n

= 2(

1− cosjπ

n

)sen

(n− 1) jπ

n,

onde na penultima identidade usamos o fato que

cos(n− 1) jπ

nsen

jπ

n= − sen

(n− 1) jπ

ncos

jπ

n

porque

0 = sen jπ = sen[(n− 1) jπ

n+

jπ

n

]= sen

(n− 1) jπ

ncos

jπ

n+ cos

(n− 1) jπ

nsen

jπ

n.

¥Os autovalores de A sao positivos, portanto A e uma matriz positiva definida. Observe que, fixado j, se n earbitrariamente grande entao

cosjπ

n≈ 1− j2π2

2n2,

pois o desenvolvimento em serie de Taylor da funcao cosseno em torno da origem e

cosx = 1− 12x2 + O

(x3

);

tomando x = jπ/n para n suficientemente grande e desprezando os termos de terceira ordem, obtemos aaproximacao acima. Daı,

2∆x2

(1− cos

jπ

n

)=

2n2

L2

(1− cos

jπ

n

)≈ 2n2

L2

(1−

[1− j2π2

2n2

])=

j2π2

L2,

de forma que os menores autovalores da matriz A sao uma boa aproximacao para os menores autovalores deDirichlet do laplaciano no intervalo [0, L]. Ja o maior autovalor da matriz A e

λn−1 =2

∆x2

(1− cos

(n− 1)π

n

)=

2n2

L2

(1− cos

(n− 1)π

n

)≈ 4n2

L2,


que nao e uma boa aproximacao para um autovalor do laplaciano. Vemos que se aumentarmos o numero depontos de discretizacao (malha mais refinada) obteremos melhores aproximacoes e uma quantidade maior deautovalores proximos aos autovalores do laplaciano. Para comparar, veja a tabela a seguir para os autovaloresdo laplaciano no intervalo [0, π]; na primeira coluna temos os autovalores exatos do laplaciano, enquanto que

na demais colunas os autovalores da matriz A, λj =2n2

π2

(1− cos

jπ

n

), com a linha superior indicando o

numero n de subintervalos na malha

n = 11 n = 21 n = 31 n = 51 n = 101 n = 10011 0.993 221 21 0.998 136 38 0.999 144 44 0.999 683 82 0.999 919 37 0.999 999 184 3.892 419 95 3.970 248 82 3.986 325 21 3.994 943 16 3.998 710 15 3.999 986 879 8.462 720 39 8.849 945 24 8.930 889 79 8.974 415 97 8.993 471 18 8.999 933 5116 14.333 863 96 15.528 221 28 15.782 100 25 15.919 213 41 15.979 370 36 15.999 789 8725 21.030 205 54 23.855 895 28 24.469 653 89 24.802 991 47 24.949 649 29 24.999 486 9936 28.009 247 34 33.646 940 78 34.904 404 68 35.592 050 94 35.895 629 79 35.998 936 2249 34.705 588 92 44.682 641 99 46.979 277 93 48.245 465 23 48.806 722 35 48.998 029 2364 40.576 732 50 56.716 479 58 60.570 369 11 62.715 235 6 63.670 436 30 63.996 637 9781 45.147 032 93 69.479 637 52 75.538 215 24 78.946 473 26 80.472 391 97 80.994 614 71100 48.046 231 68 82.687 007 94 91.729 225 95 96.877 607 56 99.196 334 56 99.991 792 02

1.2 O Caso Bidimensional

Nesta secao, desenvolveremos um metodo numerico de diferencas finitas para resolver o problema de Dirichletpara a equacao de Poisson no retangulo (0, a)× (0, b)

−∆u = f (x, y) em (0, a)× (0, b) ,u = 0 sobre ∂ ((0, a)× (0, b)) ,

e para o problema de autovalor de Dirichlet para o laplaciano no retangulo −∆u = λu em (0, a)× (0, b) ,

u = 0 sobre ∂ ((0, a)× (0, b)) .

1.2.1 A Formula dos Cinco Pontos

Vamos estabelecer alguma notacao. Denote

Ω = (0, a)× (0, b) =(x, y) ∈ R2 : 0 < x < a, 0 < y < b

.

Ao discretizar Ω atraves dos pontos

(xi, yj) = (i∆x, j∆y) , 0 6 i 6 n, 0 6 j 6 m

onde∆x =

a

n, ∆y =

b

m,

substituımos o domınio Ω pela malha (ou gride) uniforme

Ωd = (x, y) ∈ Ω : x = i∆x, y = j∆y, 1 6 i 6 n− 1, 1 6 j 6 m− 1 .

Sua fronteira discretizada e o conjunto

∂Ωd = (x, y) ∈ ∂Ω : x = i∆x, y = j∆y, 0 6 i 6 n, 0 6 j 6 m ,


de forma queΩd =

(x, y) ∈ Ω : x = i∆x, y = j∆y, 0 6 i 6 n, 0 6 j 6 m

.

A equacao de Poisson−uxx − uyy = f (x, y)

pode ser agora discretizada. Denotamos

ui,j = u (xi, yj) ,

fi,j = f (xi, yj) .

Aproximamos cada derivada parcial de segunda ordem pela sua diferenca centrada, obtendo

−uxx ≈ −ui−1,j + 2ui,j − ui+1,j

∆x2,

−uyy ≈ −ui,j−1 + 2ui,j − ui,j+1

∆y2.

Portanto, a equacao de Poisson discretizada toma a forma

−ui−1,j + 2ui,j − ui+1,j

∆x2+−ui,j−1 + 2ui,j − ui,j+1

∆y2= fi,j . (1.11)

Como a funcao u e calculada em cinco pontos, esta equacao e chamada a formula dos cinco pontos.Para cada ponto interior da malha obtemos uma equacao, logo temos um sistema linear de (n− 1) (m− 1)

equacoes com o mesmo numero de incognitas. Diferente do caso unidimensional, no entanto, nao existe umamaneira natural de ordenar os pontos da malha, logo nao podemos obter imediatamente uma representacaomatricial para o problema discretizado. Precisamos antes escolher uma ordenacao para os pontos da malha,e como existem varias ordenacoes possıveis, existem varias matrizes associadas.

Talvez a mais simples ordenacao e a ordem lexicografica induzida de Z2. Nesta ordem, os pontos damalha sao percorridos linha por linha, da esquerda para a direita, de baixo para cima:

u1,1, u2,1, . . . , un−1,1, u1,2, u2,2, . . . , un−1,2, . . . . . . , u1,m−1, u2,m−1, . . . , un−1,m−1.

Neste caso, a matriz associada ao sistema linear e uma matriz (n− 1) (m− 1) × (n− 1) (m− 1) que podeser escrita como uma matriz de (m− 1)× (m− 1) blocos de dimensao (n− 1)× (n− 1) na forma

A =

B − 1∆y2

I

− 1∆y2

I B − 1∆y2

I

− 1∆y2

I. . . . . .

. . . . . . − 1∆y2

I

− 1∆y2

I B − 1∆y2

I

− 1∆y2

I B

(m−1)×(m−1)


onde I e a matriz identidade (n− 1)× (n− 1) e B e a matriz (n− 1)× (n− 1) dada por

2(

1∆x2

+1

∆y2

)− 1

∆x2

− 1∆x2

2(

1∆x2

+1

∆y2

)− 1

∆x2

− 1∆x2

. . . . . .

. . . . . . − 1∆x2

− 1∆x2

2(

1∆x2

+1

∆y2

)− 1

∆x2

− 1∆x2

2(

1∆x2

+1

∆y2

)

(n−1)×(n−1)

Observe que

aii = 2(

1∆x2

+1

∆y2

)

para todo 1 6 i 6 (n− 1) (m− 1), enquanto que

aij = − 1∆y2

se o ponto j e vizinho a esquerda ou a direita do ponto i e

aij = − 1∆x2

se o ponto j e vizinho acima ou abaixo do ponto i. Por exemplo, no caso especial ∆x = ∆y, se n = 4 e m = 6(ou seja 3× 5 = 15 pontos internos na malha e uma matriz 15× 15), temos

A =1

∆x2

4 −1 0 −1 0 0 0 0 0 0 0 0 0 0 0−1 4 −1 0 −1 0 0 0 0 0 0 0 0 0 0

0 −1 4 0 0 −1 0 0 0 0 0 0 0 0 0−1 0 0 4 −1 0 −1 0 0 0 0 0 0 0 0

0 −1 0 −1 4 −1 0 −1 0 0 0 0 0 0 00 0 −1 0 −1 4 0 0 −1 0 0 0 0 0 00 0 0 −1 0 0 4 −1 0 −1 0 0 0 0 00 0 0 0 −1 0 −1 4 −1 0 −1 0 0 0 00 0 0 0 0 −1 0 −1 4 0 0 −1 0 0 00 0 0 0 0 0 −1 0 0 4 −1 0 −1 0 00 0 0 0 0 0 0 −1 0 −1 4 −1 0 −1 00 0 0 0 0 0 0 0 −1 0 −1 4 0 0 −10 0 0 0 0 0 0 0 0 −1 0 0 4 −1 00 0 0 0 0 0 0 0 0 0 −1 0 −1 4 −10 0 0 0 0 0 0 0 0 0 0 −1 0 −1 4

Observe que a matriz A e uma matriz simetrica, pentadiagonal e esparsa.

1.2.2 Existencia e Unicidade da Solucao Discreta – Autovalores do ProblemaBidimensional

Denotaremos por ud a funcao u|Ωd, isto e, ud e a discretizacao da funcao u no domınio discretizado Ωd.

Vamos definir o operador laplaciano discreto obtido a partir da formula dos cinco pontos por

−∆dud = −(

ui−1,j − 2ui,j + ui+1,j

∆x2+

ui,j−1 − 2ui,j + ui,j+1

∆y2

). (1.12)


de modo que a discretizacao do problema −∆u = f em Ω,

u = 0 sobre ∂Ω,

e o problema −∆dud = fd em Ωd,ud = 0 sobre ∂Ωd.

(1.13)

Para estabelecer a existencia e unicidade da solucao discreta, provaremos que a matriz de discretizacao A,que e uma matriz simetrica, e tambem uma matriz positiva definida, pois isso implica em particular que Ae invertıvel.

Lembrando que as autofuncoes de Dirichlet do laplaciano no retangulo [0, a]× [0, b] sao as funcoes

Ukl (x, y) = senkπx

asen

lπy

b,

este fato sugere que os autovetores ukl da matriz A na ordem lexicografica sao os vetores de coordenadas

Ukl (x1, y1) , Ukl (x2, y1) , . . . , Ukl (xn−1, y1) ,

Ukl (x1, y2) , Ukl (x2, y2) , . . . , Ukl (xn−1, y2) ,

...Ukl (x1, ym−1) , Ukl (x2, ym−1) , . . . , Ukl (xn−1, ym−1)

= Ukl (∆x, ∆y) , Ukl (2∆x, ∆y) , . . . , Ukl ((n− 1)∆x, ∆y) ,

Ukl (∆x, 2∆y) , Ukl (2∆x, 2∆y) , . . . , Ukl ((n− 1)∆x, 2∆y) ,

...Ukl (∆x, (m− 1)∆y) , Ukl (2∆x, (m− 1)∆y) , . . . , Ukl ((n− 1)∆x, (m− 1)∆y) ,

ou seja, como ∆x = a/n e ∆y = b/m, os vetores

ukl =(

senkπ

nsen

lπ

m, sen

2kπ

nsen

lπ

m, . . . , sen

(n− 1) kπ

nsen

lπ

m,

senkπ

nsen

2lπ

m, sen

2kπ

nsen

2lπ

m, . . . , sen

(n− 1) kπ

nsen

2lπ

m,

. . . ,

senkπ

nsen

(m− 1) lπ

m, sen

2kπ

nsen

(m− 1) lπ

m, . . . , sen

(n− 1) kπ

nsen

(m− 1) lπ

m

).

1.2 Lema. Os (n− 1)× (m− 1) autovalores da matriz A sao

λkl = 2[

1∆x2

(1− cos

kπ

n

)+

1∆y2

(1− cos

lπ

m

)]= 4

(1

∆x2sen2 kπ

2n+

1∆y2

sen2 lπ

2m

), (1.14)

k = 1, . . . , n− 1, l = 1, . . . ,m− 1, e os autovetores correspondentes sao

ukl =(

senkπ

nsen

lπ

m, sen

2kπ

nsen

lπ

m, . . . , sen

(n− 1) kπ

nsen

lπ

m,

senkπ

nsen

2lπ

m, sen

2kπ

nsen

2lπ

m, . . . , sen

(n− 1) kπ

nsen

2lπ

m, (1.15)

. . . ,

senkπ

nsen

(m− 1) lπ

m, sen

2kπ

nsen

(m− 1) lπ

m, . . . , sen

(n− 1) kπ

nsen

(m− 1) lπ

m

),

k = 1, . . . , n− 1, l = 1, . . . , m− 1.


Prova. Embora a demonstracao deste lema possa ser feita de maneira analoga a do Lema 1.1, usandoidentidades trigonometricas, daremos uma demonstracao diferente. Lembrando que as autofuncoes e osautovalores de Dirichlet do laplaciano no retangulo sao facilmente obtidos atraves do metodo de separacaode variaveis, encontraremos os autovalores da matriz A usando um metodo de separacao de variaveis discretopara achar os autovalores do laplaciano discreto

−(

ui−1,j − 2ui,j + ui+1,j

∆x2+

ui,j−1 − 2ui,j + ui,j+1

∆y2

)= λui,j . (1.16)

Em particular, este metodo nao depende da maneira como os pontos da malha sao ordenados (nao dependeda matriz A usada para representar o laplaciano discreto). Como no metodo de separacao de variaveiscontınuo, assumimos que as solucoes da equacao discreta acima sao produtos da forma

ui,j = F (i)G (j) , (1.17)

onde F e G sao funcoes de uma variavel inteira. Substituindo esta expressao na equacao de Helmholtzdiscreta, obtemos

F (i− 1)G (j)− 2F (i)G (j) + F (i + 1) G (j)∆x2

+F (i) G (j − 1)− 2F (i) G (j) + F (i)G (j + 1)

∆y2= −λF (i)G (j) .

Dividindo esta equacao por F (i)G (j), segue que

F (i− 1)− 2F (i) + F (i + 1)∆x2F (i)

+G (j − 1)− 2G (j) + G (j + 1)

∆y2G (j)= −λ.

Separando as variaveis, concluımos que cada um dos quocientes acima e independente de i ou de j, isto e,eles sao constantes:

F (i− 1)− 2F (i) + F (i + 1)F (i)

= A, (1.18)

G (j − 1)− 2G (j) + G (j + 1)G (j)

= B, (1.19)

onde as constantes α, β estao relacionadas pela identidade

A

∆x2+

B

∆y2= −λ. (1.20)

Estas equacoes podem ser escritas como formulas de recorrencia (analogas as equacoes diferenciais ordinariasobtidas no metodo de separacao de variaveis contınuo)

F (i + 1)− (A + 2) F (i) + F (i− 1) = 0,

G (j − 1)− (B + 2) G (j) + G (j + 1) = 0.

Para resolve-las, e mais conveniente trabalhar com as constantes

2α = A + 2, 2β = B + 2.

Desta forma, as equacoes para F e G tornam-se

F (i− 1)− 2αF (i) + F (i + 1) = 0, (1.21)G (j − 1)− 2βG (j) + G (j + 1) = 0. (1.22)


Observe que

λ = 2(

1− α

∆x2+

1− β

∆y2

). (1.23)

Vamos resolver a equacao para F , ja que a equacao para G e completamente analoga. Substituindo em(1.21) uma solucao da forma

F (i) = zi (1.24)

obtemoszi−1 − 2αzi + zi+1 = 0,

donde, dividindo por zi−1 extraımos a equacao quadratica (analoga a equacao indicial)

z2 − 2αz + 1 = 0. (1.25)

As duas raızes saoz± = α±

√α2 − 1,

com z+ + z− = 2α e z+z− = 1. Portanto, a solucao geral para a equacao (1.21) e

F (i) = c1zi+ + c2z

i−

para algumas constantes c1, c2. Para determinarmos estas constantes e tambem α, aplicamos as condicoesde fronteira, que implicam

F (0) = F (n) = 0.

A primeira destas por sua vez implica que c1 = −c2, logo

F (i) = c(zi+ − zi

−). (1.26)

Como a equacao para F e homogenea, a constante c e arbitraria. Aplicando a segunda, segue que

zn+ = zn

−,

ou, como z+z− = 1,z2n+ = 1

Consequentemente, z+ e uma 2n-esima raiz complexa de 1:

z+ = eijπ/n (1.27)

para algum inteiro 1 6 k 6 2n − 1, onde i =√−1. Como z− = 1/z+, podemos restringir 0 6 k 6 n − 1 e

(1.26) produz todas as solucoes nao-triviais F de (1.21).Portanto,

α =z+ + z−

2=

eiπk/n + e−iπk/n

2= cos

kπ

n, 0 6 k 6 n− 1,

e, escolhendo c = 1/2,

Fk (i) = eiπki/n − e−iπki/n = senikπ

n.

Analogamente,

β = coslπ

m, 0 6 l 6 m− 1,

eGl (j) = sen

jlπ

m.


Segue que os autovalores sao

λkl = 2[

1∆x2

(1− cos

kπ

n

)+

1∆y2

(1− cos

lπ

m

)]

e as coordenadas das autofuncoes associadas sao dadas por

(ukl)i,j = Fk (i) Gl (j) = senikπ

nsen

jlπ

m.

¥

1.3 Teorema. (Existencia e Unicidade da Solucao Discreta) Seja Ω = (0, a) × (0, b). Entao o problemadiscretizado −∆dud = fd em Ωd,

ud = 0 sobre ∂Ωd,

possui uma unica solucao.

Prova. Pelo lema anterior, os autovalores da matriz simetrica A sao positivos, logo ela e uma matrizinvertıvel. ¥

1.2.3 Princıpio do Maximo Discreto

Para obter uma estimativa a priori para a equacao de Poisson discretizada, e com isso provar a convergenciada solucao discreta para a solucao classica, usaremos um princıpio do maximo discreto que enunciaremos eprovaremos nesta subsecao.

1.4 Lema. (Propriedade do Valor Medio) Se ∆dud = 0, entao para pontos interiores vale

ui,j =∆x2 (ui,j−1 + ui,j+1) + ∆y2 (ui−1,j + ui+1,j)

2 (∆x2 + ∆y2).

Em particular, se ∆x = ∆y, entao para pontos interiores vale

ui,j =ui,j−1 + ui,j+1 + ui−1,j + ui+1,j

4.

1.5 Teorema. (Princıpio do Maximo Discreto) Se ∆dud > 0, o maximo de ud em Ωd e atingido na fronteira∂Ωd; se o maximo de ud e atingido no interior, entao ud e constante.

Se ∆dud 6 0, o mınimo de ud em Ωd e atingido na fronteira ∂Ωd; se o mınimo de ud e atingido nointerior, entao ud e constante.

Prova. Primeiro provaremos para ∆x = ∆y, para ilustrar a analogia com o caso contınuo. ∆dud > 0 implica

ui,j 6 ui,j−1 + ui,j+1 + ui−1,j + ui+1,j

4.

Logo, um ponto interior e um maximo local, isto e,

ui,j > ui,j−1, ui,j+1, ui−1,j , ui+1,j

(ou seja, e um maximo em relacao aos seus quatro vizinhos), somente se cada um dos seus quatro vizinhosassume este mesmo valor maximo, e a desigualdade torna-se uma identidade. Aplicando este argumento atodos os pontos da malha, concluımos que ou nao existe um maximo interior, e portanto o maximo e atingido


na fronteira, ou existe um maximo interior e todos os pontos da malha assumem o mesmo valor, isto e, ud econstante.

No caso geral ∆x 6= ∆y, se ∆dud > 0 temos(

1∆x2

+1

∆y2

)ui,j 6 1

2

(ui,j−1 + ui,j+1

∆y2+

ui−1,j + ui+1,j

∆x2

).

Se ui,j e um maximo local, segue que(

1∆x2

+1

∆y2

)ui,j 6 1

2

(ui,j + ui,j

∆y2+

ui,j + ui,j

∆x2

)=

12

(1

∆x2+

1∆y2

)ui,j ,

logo nenhum dos seus quatro vizinhos pode assumir um valor menor que ui,j , isto e, cada um dos quatrovizinhos assume o mesmo valor maximo e o argumento prossegue como no caso anterior. O caso ∆dud 6 0e provado considerando-se −ud. ¥

1.2.4 Convergencia da Solucao Discreta para a Solucao Classica

Por simplicidade, trabalharemos no quadrado unitario, isto e, Ω = (0, 1) × (0, 1). Consideraremos a normado maximo discreta para funcoes vd definidas no domınio discretizado Ωd:

‖vd‖∞ = max06i6n06j6m

|vi,j | .

Em primeiro lugar, obtemos uma estimativa a priori discreta (que tambem pode ser visto como um resultadode regularidade discreto) para solucoes da equacao de Poisson discreta com condicao de Dirichlet homogenea:

1.6 Lema. (Estimativa a Priori) Seja Ω = (0, 1)2. Seja ud uma solucao de −∆dud = fd em Ωd,

ud = 0 sobre ∂Ωd.

Entao‖ud‖∞ 6 1

8‖∆dud‖∞ . (1.28)

Prova. Considere a funcao

w (x, y) =14

[(x− 1

2

)2

+(

y − 12

)2]

e sua versao discretizada wd definida por

wi,j =14

[(xi − 1

2

)2

+(

yj − 12

)2]

. (1.29)

Entaow > 0 e ∆w = 1,

e tambemwd > 0 e ∆dwd = 1, (1.30)


pois

∆dwd =wi−1,j − 2wi,j + wi+1,j

∆x2+

wi,j−1 − 2wi,j + wi,j+1

∆y2

=14

[(xi−1 − 1

2

)2 +(yj − 1

2

)2 − 2(xi − 1

2

)2 − 2(yj − 1

2

)2 +(xi+1 − 1

2

)2 +(yj − 1

2

)2

∆x2

+

(xi − 1

2

)2 +(yj−1 − 1

2

)2 − 2(xi − 1

2

)2 − 2(yj − 1

2

)2 +(xi − 1

2

)2 +(yj+1 − 1

2

)2

∆y2

]

=14

[(xi−1 − 1

2

)2 − 2(xi − 1

2

)2 +(xi+1 − 1

2

)2

∆x2+

(yj−1 − 1

2

)2 − 2(yj − 1

2

)2 +(yj+1 − 1

2

)2

∆y2

]

=14

[(xi −∆x− 1

2

)2 − 2(xi − 1

2

)2 +(xi + ∆x− 1

2

)2

∆x2+

(yj −∆y − 1

2

)2 − 2(yj − 1

2

)2 +(yj + ∆y − 1

2

)2

∆y2

]

=14

[(x2

i + ∆x2 + 14 − 2xi∆x− xi + ∆x

)− 2(x2

i − xi + 14

)+

(x2

i + ∆x2 + 14 + 2xi∆x− xi −∆x

)

∆x2

+

(y2

j + ∆y2 + 14 − 2yj∆y − yj + ∆y

)− 2(y2

j − yj + 14

)+

(y2

j + ∆y2 + 14 + 2yj∆y − yj −∆y

)

∆y2

]

=14

(2∆x2

∆x2+

2∆y2

∆y2

)= 1.

Considere agora a funcaoud − ‖∆dud‖∞ wd. (1.31)

Temos entao

∆d (ud − ‖∆dud‖∞ wd) = ∆dud − ‖∆dud‖∞∆dwd

= ∆dud − ‖∆dud‖∞6 0.

Segue do Princıpio do Maximo Discreto que a funcao ud − ‖∆dud‖∞ wd assume o seu mınimo na fronteira.Este ultimo e igual a −‖∆dud‖∞max∂Ωd

wd. Por sua vez, o maximo de wd na fronteira e menor ou igual aomaximo de w em ∂Ω, dado por

max06x61

14

(x− 1

2

)2

= max06x61

14

(y − 1

2

)2

=18.

Portanto, concluımos que

ui,j > ui,j − ‖∆dud‖∞ wi,j > −18‖∆dud‖∞ (1.32)

para todos i, j. Analogamente,∆d (ud + ‖∆dud‖∞ wd) > 0

e a funcao ud + ‖∆dud‖∞ wd assume o seu maximo na fronteira, igual a ‖∆dud‖∞max∂Ωdwd 6 1

8a, donde

ui,j 6 ui,j − ‖∆dud‖∞ wi,j 6 18‖∆dud‖∞ (1.33)

para todos i, j. Reunindo as duas desigualdades, segue que

|ui,j | 6 18‖∆dud‖∞

para todos i, j, o que conclui a demonstracao. ¥


1.7 Teorema. Seja Ω = (0, 1)2. Sejam u ∈ C4(Ω

)uma solucao classica para o problema de Dirichlet

−∆u = f em Ω,u = 0 sobre ∂Ω,

e vd uma solucao do correspondente problema discretizado −∆dvd = fd em Ωd,

vd = 0 sobre ∂Ωd.

Entao existe uma constante C > 0 independente de u tal que

‖ud − vd‖∞ 6 C∥∥D4u

∥∥L∞(Ω)

(∆x2 + ∆y2

). (1.34)

Prova. A hipotese f ∈ C2,α(Ω

)garante que u ∈ C4

(Ω

). Lembre-se que

∥∥D4u∥∥

L∞(Ω)= sup

(x,y)∈Ωp+q=4

∣∣∣∣∂4u

∂xp∂yq(x, y)

∣∣∣∣ .

Pela Formula de Taylor,

∂2u

∂x2(xi, yj) =

u(xi −∆x, yj)− 2u(xi, yj) + u(xi + ∆x, yj)∆x2

− 24!

∂4u

∂x4(xi, yj)∆x2 − 2

5!∂6u

∂x6(xi, yj)∆x4 − . . .

=ui−1,j − 2ui,j + ui+1,j

∆x2− 2

4!∂4u

∂x4(xi, yj)∆x2 − 2

5!∂6u

∂x6(xi, yj)∆x4 − . . . ,

∂2u

∂y2(xi, yj) =

u(xi, yj −∆y)− 2u(xi, yj) + u(xi, yj + ∆y)∆y2

− 24!

∂4u

∂y4(xi, yj)∆y2 − 2

5!∂6u

∂y6(xi, yj)∆y4 − . . .

=ui,j−1 − 2ui,j + ui,j+1

∆y2− 2

4!∂4u

∂y4(xi, yj)∆y2 − 2

5!∂6u

∂y6(xi, yj)∆y4 − . . . ,

donde

∆u (xi, yj) = (∆dud)ij −13!

(∂4u

∂x4(xi, yj)∆x2 +

∂4u

∂y4(xi, yj)∆y2

)+ O

(∆x4, ∆y4

). (1.35)

Como−∆u (xi, yj) = f (xi, yj) ,

temos que

− (∆dud)i,j = (fd)i,j −13!

(∂4u

∂x4(xi, yj)∆x2 +

∂4u

∂y4(xi, yj)∆y2

)+ O

(∆x4, ∆y4

). (1.36)

Subtraindo desta equacao a equacao− (∆dvd)i,j = (fd)i,j ,

obtemos

− (∆dud −∆dvd)i,j = − 13!

(∂4u

∂x4(xi, yj)∆x2 +

∂4u

∂y4(xi, yj)∆y2

)+ O

(∆x4, ∆y4

),

o que implica

‖∆d (ud − vd)‖∞ 6 13!

∥∥D4u∥∥

L∞(Ω)

(∆x2 + ∆y2

)+ O

(∆x4, ∆y4

)

6 C∥∥D4u

∥∥L∞(Ω)

(∆x2 + ∆y2

).

Usando a estimativa a priori do lema anterior, obtemos finalmente o resultado desejado. ¥


Definicao. Dizemos que as solucoes do problema discretizado −∆dvd = fd em Ωd,vd = 0 sobre ∂Ωd,

convergem para a solucao exata u do problema de Poisson −∆u = f em Ω,u = 0 sobre ∂Ω,

com relacao a norma ‖·‖ se‖ud − vd‖ → 0

quando ∆x, ∆y → 0. Dizemos que a convergencia e de ordem k (ou que o esquema de diferencasfinitas e convergente de ordem k) se

‖ud − vd‖ = O(∆xk,∆yk

).

O Teorema 1.7 diz que o esquema de diferencas finitas da formula de cinco pontos e um esquema convergentena norma do sup de ordem 2, se u ∈ C4

(Ω

). Maior regularidade da solucao u nao causa melhor convergencia

no metodo. Na verdade, a ordem de convergencia da formula de cinco pontos ainda e 2 mesmo sob hipotesesmais fracas sobre a regularidade de u: basta assumir u ∈ C3,1

(Ω

), ao inves de u ∈ C4

(Ω

). No entanto,

regularidade menor que esta em u afeta negativamente a ordem de convergencia da formula de cinco pontos.Em geral, pode-se provar que se u ∈ Ck,α

(Ω

), 2 6 k 6 4, entao existe uma constante C = C (k, α) tal que

‖ud − vd‖∞ 6 C(∆xk+α−2 + ∆yk+α−2

) ‖u‖Ck,α(Ω) . (1.37)

Para uma demonstracao destes resultados, veja [Hackbusch], pags. 60-61. Se quisermos uma melhor ordemde convergencia para as solucoes discretizadas, e necessario considerar outras forma de discretizar o laplacianoatraves de diferencas finitas. Isto sera feito na proxima secao.

1.3 Discretizacoes de Ordem Superior

Para obter esquemas de diferencas finitas com melhor ordem de convergencia, em geral e necessario acres-centar mais pontos na formula. O metodo dos coeficientes indeterminados e um metodo simples paraconstruir estes esquemas.

1.3.1 Caso Unidimensional

Vamos obter um esquema de diferencas finitas convergente de ordem 4 para o caso unidimensional. Oesquema envolvendo tres pontos, que obtivemos no inıcio do capıtulo atraves da aproximacao da derivadasegunda em um ponto por uma diferenca finita centrada (que envolve o ponto e seus dois vizinhos, a esquerdae a direita), e convergente de ordem 2 (isso que pode ser provado de maneira semelhante a como fizemos paraa formula de cinco pontos). Para obter um esquema com uma maior ordem de convergencia, acrescentamosmais dois pontos a formula de diferencas finitas do esquema, que denotaremos por δui:

δui = c1ui−2 + c2ui−1 + c3ui + c4ui+1 + c5ui+2. (1.38)

Cada termo tem sua expansao em serie de Taylor:

u(xi − 2∆x) = u(xi)− 2u′(xi)∆x +42!

u′′(xi)∆x2 − 83!

u′′′(xi)∆x3 +164!

u(4)(xi)∆x4 − 325!

u(5)(xi)∆x5 + O(∆x6

),

u(xi −∆x) = u(xi)− u′(xi)∆x +12!

u′′(xi)∆x2 − 13!

u′′′(xi)∆x3 +14!

u(4)(xi)∆x4 − 15!

u(5)(xi)∆x5 + O(∆x6

),

u(xi + ∆x) = u(xi) + u′(xi)∆x +12!

u′′(xi)∆x2 +13!

u′′′(xi)∆x3 +14!

u(4)(xi)∆x4 +15!

u(5)(xi)∆x5 + O(∆x6

),

u(xi + 2∆x) = u(xi) + 2u′(xi)∆x +42!

u′′(xi)∆x2 +83!

u′′′(xi)∆x3 +164!

u(4)(xi)∆x4 +325!

u(5)(xi)∆x5 + O(∆x6

).


Substituindo estas expressoes na formula acima, obtemos:

δui = (c1 + c2 + c3 + c4 + c5)u (xi)+ ∆x (−2c1 − c2 + c4 + 2c5)u′(xi)

+ ∆x2

(2c1 +

12c2 +

12c4 + 2c5

)u′′(xi)

+ ∆x3

(−4

3c1 − 1

6c2 +

16c4 +

43c5

)u′′′(xi)

+ ∆x4

(23c1 +

124

c2 +124

c4 +23c5

)u(4)(xi)

+ ∆x5

(− 4

15c1 − 1

120c2 +

1120

c4 +415

c5

)u(5)(xi)

+ O(∆x6

).

Como procuramos um esquema de diferencas finitas com ordem de convergencia maior que 2, queremos obteruma solucao nao-nula para o sistema

c1 + c2 + c3 + c4 + c5 = 0−2c1 − c2 + c4 + 2c5 = 0

2c1 +12c2 +

12c4 + 2c5 =

1∆x2

−43c1 − 1

6c2 +

16c4 +

43c5 = 0

23c1 +

124

c2 +124

c4 +23c5 = 0

;

isso implicaria em princıpio em um esquema com ordem de convergencia pelo menos igual a 3:

δui = u′′(xi) + O(∆x3

).

Como a matriz

1 1 1 1 1−2 −1 0 1 2

212

012

2

−43

−16

016

43

23

124

0124

23

tem determinante igual a 1, ela e invertıvel e o sistema possui a solucao unica

c1 = − 112

1∆x2

,

c2 =43

1∆x2

,

c3 = −52

1∆x2

c4 =43

1∆x2

,

c5 = − 112

1∆x2

.


Incidentalmente, esta solucao tambem implica

− 415

c1 − 1120

c2 +1

120c4 +

415

c5 = 0

o que permite obter um esquema com ordem de convergencia igual a 4:

δui = u′′(xi) + O(∆x4

),

aproximando a derivada segunda u′′ pela diferenca finita

u′′ =− 1

12ui−2 +

43ui−1 − 5

2ui +

43ui+1 − 1

12ui+2

∆x2

ou−u′′ =

ui−2 − 16ui−1 + 30ui − 16ui+1 + ui+2

12∆x2. (1.39)

1.3.2 Caso Bidimensional: A Formula dos Nove Pontos Compacta

Um esquema de ordem 4 para a equacao de Poisson em duas dimensoes e a formula de nove pontos compacta.Se buscassemos uma formula de nove pontos simplesmente a partir da formula de cinco pontos unidi-

mensional obtida na subsecao precedente (como obtivemos a formula de cinco pontos bidimensional a partirda formula de tres pontos unidimensional), escreverıamos

−∆dud =ui−2,j − 16ui−1,j + 30ui,j − 16ui+1,j + ui+2,j

12∆x2+

ui,j−2 − 16ui,j−1 + 30ui,j − 16ui,j+1 + ui,j+2

12∆y2,

(1.40)que pode ser resumida na forma

−∆dud =

− 112∆y2

− 1612∆y2

− 112∆x2

− 1612∆x2

30(

112∆x2

+1

12∆y2

)− 16

12∆x2− 1

12∆x2

− 1612∆y2

− 112∆y2

.

Embora este esquema seja de fato de ordem 4, ele apresenta dificuldades para pontos interiores adjacentes afronteira do retangulo (por exemplo, se considerarmos o ponto (x1, y1), os pontos (x−1, y1) e (x1, y−1) estaofora do retangulo). Uma possibilidade para resolver este problema seria aplicar a formula dos cinco pontosnos pontos interiores adjacentes a fronteira e aplicar a formula dos nove pontos apenas nos pontos interioresmais distantes da fronteira. No entanto, como a formula de cinco pontos e de segunda ordem, a convergenciadeste metodo misto nao deve ser de ordem 4.

Vamos tentar encontrar uma formula de nove pontos compacta, em que os nove pontos estao dispostosem tres linhas e tres colunas, de modo que nao ha problemas em usa-la nos pontos interiores adjacentes afronteira. Aplicando o metodo dos coeficientes indeterminados, buscamos nove coeficientes para a diferencafinita

−∆dud = c1ui−1,j−1 + c2ui,j−1 + c3ui+1,j−1

+ c4ui−1,j + c5ui,j + c6ui+1,j (1.41)+ c7ui−1,j+1 + c8ui,j+1 + c9ui+1,j+1.


Observe a distribuicao dos nove pontos. Alem dos cinco usuais, foram acrescentados os quatro pontos queocupam as posicoes diagonais. Para os quatro pontos vizinhos horizontais ou verticais do ponto central, aformula de Taylor produz

u(xi −∆x, yj) = u(xi, yj)− ∂u

∂x(xi, yj)∆x +

12!

∂2u

∂x2(xi, yj)∆x2 − 1

3!∂3u

∂x3(xi, yj)∆x3 +

14!

∂4u

∂x4(xi, yj)∆x4

− 15!

∂5u

∂x5(xi, yj)∆x5 + O

(∆x6

)

u(xi + ∆x, yj) = u(xi, yj) +∂u

∂x(xi, yj)∆x +

12!

∂2u

∂x2(xi, yj)∆x2 +

13!

∂3u

∂x3(xi, yj)∆x3 +

14!

∂4u

∂x4(xi, yj)∆x4

+15!

∂5u

∂x5(xi, yj)∆x5 + O

(∆x6

)

u(xi, yj −∆y) = u(xi, yj)− ∂u

∂y(xi, yj)∆y +

12!

∂2u

∂y2(xi, yj)∆y2 − 1

3!∂3u

∂y3(xi, yj)∆y3 +

14!

∂4u

∂y4(xi, yj)∆y4

− 15!

∂5u

∂x5(xi, yj)∆x5 + O

(∆x6

)

u(xi, yj + ∆y) = u(xi, yj) +∂u

∂y(xi, yj)∆y +

12!

∂2u

∂y2(xi, yj)∆y2 +

13!

∂3u

∂y3(xi, yj)∆y3 +

14!

∂4u

∂y4(xi, yj)∆y4

+15!

∂5u

∂x5(xi, yj)∆x5 + O

(∆x6, ∆y6

)

enquanto que para os quatro pontos diagonais temos

u(xi + ∆x, yj + ∆y)

= u(xi, yj) +[∂u

∂x(xi, yj)∆x +

∂u

∂y(xi, yj)∆y

]+

12!

[∂2u

∂x2(xi, yj)∆x2 + 2

∂2u

∂x∂y(xi, yj)∆x∆y +

∂2u

∂y2(xi, yj)∆y2

]

+13!

[∂3u

∂x3(xi, yj)∆x3 + 3

∂3u

∂x2∂y(xi, yj)∆x2∆y + 3

∂3u

∂x∂y2(xi, yj)∆x∆y2 +

∂3u

∂y3(xi, yj)∆y3

]

+14!

[∂4u

∂x4(xi, yj)∆x4 + 4

∂4u

∂x3∂y(xi, yj)∆x3∆y + 6

∂4u

∂x∂y3(xi, yj)∆x2∆y2 + 4

∂3u


∂4u

∂y4(xi, yj)∆y4

]

+15!

[∂5u

∂x5(xi, yj)∆x5 + 5

∂5u

∂x4∂y(xi, yj)∆x4∆y + 10

∂5u

∂x3∂y2(xi, yj)∆x3∆y2 + 10

∂5u

∂x∂y4(xi, yj)∆x2∆y3

+5∂5u


∂5u

∂y5(xi, yj)∆y5

]+ O

(∆x6, ∆y6

),

u(xi −∆x, yj −∆y)

= u(xi, yj)−[∂u

∂x(xi, yj)∆x +

∂u

∂y(xi, yj)∆y

]+

12!

[∂2u

∂x2(xi, yj)∆x2 + 2

∂2u


∂2u

∂y2(xi, yj)∆y2

]

− 13!

[∂3u

∂x3(xi, yj)∆x3 + 3

∂3u

∂x2∂y(xi, yj)∆x2∆y + 3

∂3u


∂3u

∂y3(xi, yj)∆y3

]

+14!

[∂4u

∂x4(xi, yj)∆x4 + 4

∂4u

∂x3∂y(xi, yj)∆x3∆y + 6

∂4u

∂x∂y3(xi, yj)∆x2∆y2 + 4

∂3u


∂4u

∂y4(xi, yj)∆y4

]

− 15!

[∂5u

∂x5(xi, yj)∆x5 + 5

∂5u

∂x4∂y(xi, yj)∆x4∆y + 10

∂5u

∂x3∂y2(xi, yj)∆x3∆y2 + 10

∂5u


+5∂5u


∂5u

∂y5(xi, yj)∆y5

]+ O

(∆x6

)


u(xi + ∆x, yj −∆y)

= u(xi, yj) +[∂u

∂x(xi, yj)∆x− ∂u

∂y(xi, yj)∆y

]

+12!

[∂2u

∂x2(xi, yj)∆x2 − 2

∂2u


∂2u

∂y2(xi, yj)∆y2

]

+13!

[∂3u

∂x3(xi, yj)∆x3 − 3

∂3u

∂x2∂y(xi, yj)∆x2∆y + 3

∂3u

∂x∂y2(xi, yj)∆x∆y2 − ∂3u

∂y3(xi, yj)∆y3

]

+14!

[∂4u

∂x4(xi, yj)∆x4 − 4

∂4u

∂x3∂y(xi, yj)∆x3∆y + 6

∂4u

∂x∂y3(xi, yj)∆x2∆y2 − 4

∂3u


∂4u

∂y4(xi, yj)∆y4

]

+15!

[∂5u

∂x5(xi, yj)∆x5 − 5

∂5u

∂x4∂y(xi, yj)∆x4∆y + 10

∂5u

∂x3∂y2(xi, yj)∆x3∆y2 − 10

∂5u


+5∂5u

∂x∂y4(xi, yj)∆x∆y4 − ∂5u

∂y5(xi, yj)∆y5

]+ O

(∆x6, ∆y6

),

u(xi −∆x, yj + ∆y)

= u(xi, yj) +[−∂u

∂x(xi, yj)∆x +

∂u

∂y(xi, yj)∆y

]

+12!

[∂2u

∂x2(xi, yj)∆x2 − 2

∂2u


∂2u

∂y2(xi, yj)∆y2

]

+13!

[−∂3u

∂x3(xi, yj)∆x3 + 3

∂3u

∂x2∂y(xi, yj)∆x2∆y − 3

∂3u


∂3u

∂y3(xi, yj)∆y3

]

+14!

[∂4u

∂x4(xi, yj)∆x4 − 4

∂4u

∂x3∂y(xi, yj)∆x3∆y + 6

∂4u

∂x∂y3(xi, yj)∆x2∆y2 − 4

∂3u


∂4u

∂y4(xi, yj)∆y4

]

+15!

[−∂5u

∂x5(xi, yj)∆x5 + 5

∂5u

∂x4∂y(xi, yj)∆x4∆y − 10

∂5u

∂x3∂y2(xi, yj)∆x3∆y2 + 10

∂5u


−5∂5u


∂5u

∂y5(xi, yj)∆y5

]+ O

(∆x6, ∆y6

).


Substituindo estas expressoes na formula acima, obtemos:

−∆dud = (c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9)u (xi, yj)

+ ∆x (−c1 + c3 − c4 + c6 − c7 + c9)∂u

∂x(xi, yj)

+ ∆y (−c1 − c2 − c3 + c7 + c8 + c9)∂u

∂y(xi, yj)

+ ∆x2

(12c1 +

12c3 +

12c4 +

12c6 +

12c7 +

12c9

)∂2u

∂x2(xi, yj)

+ ∆x∆y (c1 − c3 − c7 + c9)∂2u

∂x∂y(xi, yj)

+ ∆y2

(12c1 +

12c2 +

12c3 +

12c7 +

12c8 +

12c9

)∂2u

∂y2(xi, yj)

+ ∆x3

(−1

6c1 +

16c3 − 1

6c4 +

16c6 − 1

6c7 +

16c9

)∂3u

∂x3(xi, yj)

+ ∆x2∆y

(−1

2c1 − 1

2c3 +

12c7 +

12c9

)∂3u

∂x2∂y(xi, yj)

+ ∆x∆y2

(−1

2c1 +

12c3 − 1

2c7 +

12c9

)∂3u

∂x∂y2(xi, yj)

+ ∆y3

(−1

6c1 − 1

6c2 − 1

6c3 +

16c7 +

16c8 +

16c9

)∂3u

∂y3(xi, yj)

+ ∆x4

(124

c1 +124

c3 +124

c4 +124

c6 +124

c7 +124

c9

)∂4u

∂x4(xi, yj)

+ ∆x3∆y

(16c1 − 1

6c3 − 1

6c7 +

16c9

)∂4u

∂x3∂y(xi, yj)

+ ∆x2∆y2

(14c1 +

14c3 +

14c7 +

14c9

)∂4u

∂x2∂y2(xi, yj)

+ ∆x∆y3

(16c1 − 1

6c3 − 1

6c7 +

16c9

)∂4u

∂x∂y3(xi, yj)

+ ∆y4

(124

c1 +124

c2 +124

c3 +124

c7 +124

c8 +124

c9

)∂4u

∂y4(xi, yj)

+ ∆x5

(− 1

120c1 +

1120

c3 − 1120

c4 +1

120c6 − 1

120c7 +

1120

c9

)∂5u

∂x5(xi, yj)

+ ∆x4∆y

(− 1

24c1 − 1

24c3 +

124

c7 +124

c9

)∂5u

∂x4∂y(xi, yj)

+ ∆x3∆y2

(− 1

12c1 +

112

c3 +112

c7 +112

c9

)∂5u

∂x3∂y2(xi, yj)

+ ∆x2∆y3

(− 1

12c1 − 1

12c3 − 1

12c7 +

112

c9

)∂5u

∂x2∂y3(xi, yj)

+ ∆x∆y4

(− 1

24c1 +

124

c3 − 124

c7 +124

c9

)∂5u

∂x∂y4(xi, yj)

+ ∆y5

(− 1

120c1 − 1

120c2 − 1

120c3 +

1120

c7 +1

120c8 +

1120

c9

)∂5u

∂y5(xi, yj)


Para obter um esquema com ordem de convergencia pelo menos igual a 3, precisarıamos obter uma solucaonao-nula para o sistema

c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9 = 0−c1 + c3 − c4 + c6 − c7 + c9 = 0−c1 − c2 − c3 + c7 + c8 + c9 = 0

c1 + c3 + c4 + c6 + c7 + c9 =1

∆x2

c1 − c3 − c7 + c9 = 0

c1 + c2 + c3 + c7 + c8 + c9 =1

∆y2

−c1 + c3 − c4 + c6 − c7 + c9 = 0−c1 − c3 + c7 + c9 = 0−c1 + c3 − c7 + c9 = 0−c1 − c2 − c3 + c7 + c8 + c9 = 0c1 + c3 + c4 + c6 + c7 + c9 = 0c1 − c3 − c7 + c9 = 0c1 + c3 + c7 + c9 = 0c1 − c3 − c7 + c9 = 0c1 + c2 + c3 + c7 + c8 + c9 = 0

Infelizmente este sistema nao tem solucao pois ele e inconsistente: a sexta e a ultima equacao sao incom-patıveis, assim como a quarta e a decima primeira. Portanto, nao existe uma formula de nove pontoscompacta tal que

−∆dud = −∆u + O(∆x3, ∆y3

).

No entanto, em 1975 o matematico e logico Rosser introduziu a seguinte formula de nove pontos compactano caso especial ∆x = ∆y (em [Rosser1]; veja tambem [Rosser2])

∆dud =ui−1,j−1 + 4ui,,j−1 + ui+1,j−1 + 4ui−1,j − 20ui,j + 4ui+1,j + ui−1,j+1 + 4ui,j+1 + ui+1,j+1

6∆x2, (1.42)

que pode ser resumida na forma

−∆dud =1

6∆x2

−1 −4 −1−4 20 −4−1 −4 −1

, (1.43)

a qual produz um esquema convergente de quarta ordem se a solucao u ∈ C6(Ω

)(ou mesmo se u ∈ C5,1

(Ω

)apenas) dependendo de como a funcao f e discretizada. Para entender como isso ocorre, observe que seu ∈ C8

(Ω

)a formula de Taylor produz

−∆dud = −∆u− ∆x2

12∆2u− ∆x4

360

[∂4

∂x4+ 4

∂4

∂x2∂y2+

∂4

∂y4

]∆u + O

(∆x6

)(1.44)

= −∆u +∆x2

12∆f +

∆x4

360

[∂4

∂x4+ 4

∂4

∂x2∂y2+

∂4

∂y4

]f + O

(∆x6

). (1.45)

O ponto crucial aqui e que o erro e expresso em termos de −∆u e, consequentemente, por f . Ainda enecessario escolher uma discretizacao especial para f :

fd =fi,,j−1 + fi−1,j + 8fi,j + fi+1,j + fi,j+1

12(1.46)

ou

fd =112

11 8 1

1

. (1.47)


Usando a formula de Taylor para f , obtemos que esta discretizacao especial para f satisfaz

fd = f +∆x2

12∆f + O

(∆x4

). (1.48)

Somando esta estimativa com (1.45), e usando −∆dud = fd,−∆u = f , obtemos

−∆dud = −∆u + O(∆x4

)

Para este esquema, pode-se provar (veja [Hackbusch], pag. 64) que existe uma constante C > 0 tal que

‖ud − vd‖∞ 6 C∆x4 ‖u‖C6(Ω) ou ‖ud − vd‖∞ 6 C∆x4 ‖u‖C5,1(Ω) (1.49)

O esquema de Rosser tambem satisfaz o princıpio do maximo. Concluindo, vemos que uma maior regularidadeda solucao permite obter metodos de diferencas finitas com maior ordem de convergencia, embora esta naoseja uma tarefa simples.

1.4 Diferencas Finitas em Coordenadas Polares

Consideraremos nesta secao diferencas finitas em coordenadas polares para domınios com simetria radial.Consideraremos em detalhes os casos do disco e do anel. O primeiro caso inclui a origem no domınio dadefinicao, onde o laplaciano apresenta uma singularidade quando escrito em coordenadas polares, singulari-dade esta que nao existe no problema original, e esta particularidade deve ser tratada com cuidado para naoatrapalhar a ordem de convergencia do esquema obtido.

Considere a equacao de Poisson em coordenadas polares no disco Ω = [0, R)× [0, 2π) :

urr +1rur +

1r2

uθθ = f (r, θ) se 0 6 r < R e 0 < θ < 2π,

u (R, θ) = 0 se 0 6 θ 6 2π.

A solucao exata deste problema deve satisfazer a condicao de continuidade

u (r, 0) = u (r, 2π) para todo 0 6 r 6 R.

Embora esta condicao nao seja uma condicao de fronteira e aparece apenas por causa do sistema de coor-denadas utilizado, ela acaba funcionando como uma condicao de fronteira em muitos metodos numericos (emesmo analıticos), pois nao deixa de ser uma condicao na fronteira do retangulo (0, R)× (0, 2π).

∆r

∆θ

Discretizamos o disco atraves de uma malha polar

Ωd = (ri, θj) ∈ Ω : ri = i∆r, θj = j∆θ, 0 6 i 6 n− 1, 0 6 j 6 m

onde∆r =

R

n, ∆θ =

2π

m.


Sua fronteira discretizada e o conjunto

∂Ωd = (rn, θj) ∈ ∂Ω : rn = n∆r = R, θj = j∆θ, 0 6 j 6 m .

Discretizamos a equacao de Poisson da seguinte forma. Denotamos os valores das discretizacoes ud e fd

em pontos da malha por

ui,j = u (ri, θj) ,

fi,j = f (ri, θj) ,

entendendo que ui,j e fi,j devem satisfazer

u0,0 = u0,j e f0,0 = f0,j (1.50)

para todo 0 6 j 6 m, ja que existe apenas um ponto associado com i = 0 (a origem, correspondente a r = 0).Alem disso, pela condicao de continuidade, devemos ter tambem

ui,0 = ui,2π e fi,0 = fi,2π (1.51)

para todo 0 6 i 6 n. Usando uma diferenca centrada usual para derivadas segundas, o terceiro termo dolaplaciano em coordenadas polares pode ser aproximado para pontos interiores do disco por

(1r2

uθθ

)(ri, θj) ≈ 1

r2i

ui,j−1 − 2ui,j − ui,j+1

∆θ2. (1.52)

Para aproximar os primeiros dois termos, escrevemos

urr +1rur =

1r

(rur)r .

Se (ri, θj) e um ponto interior do disco diferente da origem (isto e, i 6= 0), podemos usar diferencas centradaspara a derivada primeira, tanto na primeira quanto na segunda aproximacoes a seguir, obtendo

1r

(rur)r (ri, θj) ≈ 1ri

(rur) (ri + ∆r/2, θj)− (rur) (ri −∆r/2, θj)2∆r/2

≈ 1ri

ri+1/2u (ri + ∆r, θj)− u (ri, θj)

∆r− ri−1/2

u (ri, θj)− u (ri −∆r, θj)∆r

∆r

=1ri

ri+1/2 (ui+1,j − ui,j)− ri−1/2 (ui,j − ui−1,j)∆r2

. (1.53)

Portanto, a discretizacao da equacao de Poisson no disco para pontos interiores do disco diferentes da origeme

−[

1ri

ri+1/2 (ui+1,j − ui,j)− ri−1/2 (ui,j − ui−1,j)∆r2

+1r2i

ui,j−1 − 2ui,j − ui,j+1

∆θ2

]= fi,j (1.54)

para 1 6 i 6 n − 1 e 1 6 j 6 m − 1. Se j = 0, usando a condicao de continuidade que identifica o ponto(i, 0) com o ponto (i, n), substituımos ui,j−1 por ui,n−1e escrevemos

−[

1ri

ri+1/2 (ui+1,0 − ui,0)− ri−1/2 (ui,0 − ui−1,0)∆r2

+1r2i

ui,n−1 − 2ui,0 − ui,1

∆θ2

]= fi,0 (1.55)

para 1 6 i 6 n − 1. Como este esquema de diferencas finitas foi obtido atraves de diferencas centradas,ele deve ser de segunda ordem. No entanto, devemos ter cuidado ao discretizar a equacao de Poisson na


origem para preservar esta ordem de convergencia. Para isso, multiplicamos a equacao de Poisson por r eintegramos o resultado sobre um pequeno disco Dε centrado na origem de raio ε:

∫ 2π

0

∫ ε

0

fr drdθ =∫ 2π

0

∫ ε

0

r

[1r

(rur)r +1r2

uθθ

]drdθ

=∫ 2π

0

∫ ε

0

(rur)r drdθ +∫ ε

0

1r

∫ 2π

0

uθθ drdθ

=∫ 2π

0

[rur]ε0 dθ +

∫ ε

0

1r

[uθ]2π0 drdθ

= ε

∫ 2π

0

ur (ε, θ) dθ,

onde assumimos u ∈ C2 (Ω) de modo que

uθ (r, 0) = uθ (r, 2π)

para todo 0 6 r < R. Escolhendo ε = ∆r/2, discretizamos a equacao integral

∆r

2

∫ 2π

0

ur (∆r/2, θ) dθ =∫ 2π

0

∫ ∆r/2

0

fr drdθ

aproximando a derivada primeira ur (∆r/2, θ) = (ur)i+1/2,j por diferencas centradas e f por f (0) (pois ∆r

e suposto pequeno), de modo que

ur (∆r/2, θj) ≈ u1,j − u0,j

∆r,

∫ 2π

0

∫ ∆r/2

0

fr drdθ ≈ f (0)∫ 2π

0

∫ ∆r/2

0

r drdθ = 2πf (0)r2

2

∣∣∣∣∆r/2

0

=π

4f (0)∆r2,

e assim∆r

2

m−1∑

j=0

u1,j − u0,j

∆r∆θ =

π

4f (0)∆r2,

donde, como u0 := u0,j independe de j, segue que o valor de u na origem sera dado por

m∆θ

2u0 =

∆θ

2

m−1∑

j=0

u1,j − π

4f (0)∆r2,

ou, usando m∆θ = 2π,4u0

∆r2− 2∆θ

π∆r2

m−1∑

j=0

u1,j = f0. (1.56)

Para escrever essas diferencas finitas em forma matricial

Au = f ,

escolhemos ordenar os pontos da malha discretizada no retangulo polar (ri, θj) : 1 6 i 6 n− 1, 0 6 j 6 mpela ordem lexicografica em (θ, r) e colocando a origem antes de todos estes pontos:.

u = (u0, u1,0, u1,1, . . . , u1,m−1, u2,0, u2,1, . . . , u2,m−1, . . . . . . , un−1,0, un−1,1, . . . , un−1,m−1) . (1.57)


Observe que existem (n− 1)×m + 1 incognitas. Nesta ordenacao, segue que A tem a forma em blocos

A =

α0 ba B1 −β1I

−α2I B2 −β2I. . .

−α3I B3 −β3I. . . . . . . . .

−αn−2I Bn−2 −βn−2I−αn−1I Bn−1

, (1.58)

ondeα0 =

4∆r2

,

a =

−α1

...−α1

m×1

,

αi =1

∆r2

ri−1/2

ri, i = 1, . . . , n− 1,

βi =1

∆r2

ri+1/2

ri, i = 1, . . . , n− 2,

b =[ −β0 . . . −β0

]1×m

,

β0 =2π

∆θ

∆r2,

I = Im,

Bi =

γi −δi 0 −δi

−δi γi −δi

−δi γi −δi

. . . . . . . . .−δi γi −δi

−δi −δi γi

m×m

,

onde

γi =1ri

ri+1/2 + ri−1/2

∆r2+

2r2i

1∆θ2

,

δi =1r2i

1∆θ2

.


A matriz A em geral nao e simetrica. Por exemplo, no caso n = 4 e m = 5 ((n− 1)×m + 1 = 16) temos

α −β0 −β0 −β0 −β0 −β0 0 0 0 0 0 0 0 0 0 0−α1 γ1 −δ1 0 0 −δ1 −β1 0 0 0 0 0 0 0 0 0−α1 −δ1 γ1 −δ1 0 0 0 −β1 0 0 0 0 0 0 0 0−α1 0 −δ1 γ1 −δ1 0 0 0 −β1 0 0 0 0 0 0 0−α1 0 0 −δ1 γ1 −δ1 0 0 0 −β1 0 0 0 0 0 0−α1 −δ1 0 0 −δ1 γ1 0 0 0 0 −β1 0 0 0 0 0

0 −α2 0 0 0 0 γ2 −δ2 0 0 −δ2 −β2 0 0 0 00 0 −α2 0 0 0 −δ2 γ2 −δ2 0 0 0 −β2 0 0 00 0 0 −α2 0 0 0 −δ2 γ2 −δ2 0 0 0 −β2 0 00 0 0 0 −α2 0 0 0 −δ2 γ2 −δ2 0 0 0 −β2 00 0 0 0 0 −α2 −δ2 0 0 −δ2 γ2 0 0 0 0 −β2

0 0 0 0 0 0 −α3 0 0 0 0 γ3 −δ3 0 0 −δ3

0 0 0 0 0 0 0 −α3 0 0 0 −δ3 γ3 −δ3 0 00 0 0 0 0 0 0 0 −α3 0 0 0 −δ3 γ3 −δ3 00 0 0 0 0 0 0 0 0 −α3 0 0 0 −δ3 γ3 −δ3

0 0 0 0 0 0 0 0 0 0 −α3 −δ3 0 0 −δ3 γ3

A primeira linha e a primeira coluna sao diferentes porque os pontos (0, j), j = 0, . . . ,m, sao realmente umunico ponto e este ponto e vizinho a todos os pontos (1, j), j = 0, . . . , m.

A matriz de discretizacao A no caso do anel sera um pouco mais simples, ja que ela sera igual a matrizde discretizacao no caso do disco menos a primeira linha e a primeira coluna.

1.5 Domınios Arbitrarios

Queremos agora discutir a resolucao numerica da equacao de Poisson atraves de diferencas finitas em umdomınio arbitrario.

Seja Ω ⊂ R2 um domınio arbitrario. Se sobrepusermos uma malha uniforme

M = (i∆x, j∆y) ∈ Ω : i ∈ Z e j ∈ Zsobre Ω, obtemos um domınio discretizado definido por

Ωd = (x, y) ∈ Ω : x/∆x ∈ Z e y/∆y ∈ Z . (1.59)

Esta e exatamente a maneira como discretizamos o retangulo. No entanto, o conjunto discretizado dospontos de fronteira ∂Ωd de um domınio arbitrario deve ser tratado de maneira diferente do retangulo, ja quea malha uniforme M em geral nao vai se sobrepor a fronteira de Ω, podendo nao possuir nenhum ponto emcomum com a fronteira ou um numero muito pequeno de pontos em poucas regioes da fronteira.

Uma maneira de tratar este problema e a seguinte. Para determinar se o ponto (xi, yj) ∈ Ωd e adjacentea “fronteira esquerda” de Ω, por exemplo, e ao mesmo tempo encontrar o seu vizinho a esquerda na fronteirase for o caso, basta verificar se o segmento

[xi −∆x, yj ] = (xi − t∆x, yj) : t ∈ [0, 1]esta inteiramente contido em Ω ou nao. Se nao estiver, entao (xi, yj) e um ponto interior adjacente a fronteirae existe um numero tW ∈ (0, 1) tal que

(xi − tW ∆x, yj) ∈ ∂Ω e (xi − t∆x, yj) ∈ Ω para todo t ∈ [0, tW ). (1.60)

Este sera o vizinho a esquerda de (xi, yj) na fronteira discretizada ∂Ωd do domınio. Analogamente, ospontos vizinhos na fronteira discretizada a direita, abaixo e acima de pontos adjacentes a fronteira podemser encontrados; eles satisfazem, respectivamente,

(xi + tE∆x, yj) ∈ ∂Ω e (xi + t∆x, yj) ∈ Ω para todo t ∈ [0, tE). (1.61)


(xi, yj − tS∆y) ∈ ∂Ω e (xi, yj − t∆y) ∈ Ω para todo t ∈ [0, tS). (1.62)

(xi, yj + tN∆y) ∈ ∂Ω e (xi, yj + t∆y) ∈ Ω para todo t ∈ [0, tN ). (1.63)

(os subındices W,E, S, N correspondem aos quatro pontos cardeais oeste, leste, sul, norte em ingles). Defin-imos

∂Ωd = (x, y) ∈ ∂Ω : (x, y) satisfaz (1.60), (1.61), (1.62) ou (1.63) (1.64)

Dependendo da geometria de Ω e concebıvel que um ponto seja simultaneamente adjacente as “quatrofronteiras” de Ω, isto e, que ele tenha os seus quatro vizinhos em ∂Ωd. Alem disso, embora os pontosinteriores da malha estejam distribuıdos uniformemente, esta discretizacao da fronteira do domınio permiteque as vezes dois pontos da malha da fronteira estejam bem proximos um do outro em alguma regiao dafronteira e relativamente distantes em outras (isso ocorre mesmo em domınio regulares como um disco).

Para discretizar a equacao de Poisson nesta malha, observe que pela formula de Taylor temos, para pontosx− < x < x+,

u′′ (x) =2

x+ − x−

(u (x+)− u (x)

x+ − x− u (x)− u (x−)

x− x−

)+ r, (1.65)

onde

|r| 6 13

(x+ − x)2 + (x− x−)2

x+ − x−‖u‖C3([x−,x+]) 6 1

3max (x+ − x, x− x−) ‖u‖C3([x−,x+]) . (1.66)

De fato,

u(x−) = u(x)− u′(x) (x− x−) +12u′′(x) (x− x−)2 − 1

3!u′′′(ξ−) (x− x−)3 ,

u(x+) = u(x) + u′(x) (x+ − x) +12u′′(x) (x+ − x)2 +

13!

u′′′(ξ+) (x+ − x)3 ,

para alguns ξ− ∈ [x−, x] , ξ+ ∈ [x, x+], de modo que

−u (x)− u (x−)x− x−

= −u′(x) +12u′′(x) (x− x−)− 1

6u′′′(ξ−) (x− x−)2 ,

u (x+)− u (x)x+ − x

= u′(x) +12u′′(x) (x+ − x) +

16u′′′(ξ+) (x+ − x)2 ,

donde, somando as duas expressoes,

u (x+)− u (x)x+ − x

− u (x)− u (x−)x− x−

=12u′′(x) (x+ − x−) +

16

[u′′′(ξ+) (x+ − x)2 − u′′′(ξ−) (x− x−)2

].

Assim, podemos aproximar

u′′ (x) ≈ 2x+ − x−

(u (x+)− u (x)

x+ − x− u (x)− u (x−)

x− x−

)

Se x− = x−∆x e x+ = x + ∆x, obtemos a formula de diferencas centradas usual para a derivada segunda.Para aproximar o laplaciano atraves de uma formula de cinco pontos, usamos os quatro pontos vizinhos

(xi − tW ∆x, yj) , (xi + tE∆x, yj) , (xi, yj − tS∆y) , (xi, yj + tN∆y) , com t∗ ∈ (0, 1]


definindo o esquema de diferencas finitas de Shortley-Weller :

∆dud =2

(xi + tE∆x)− (xi − tW ∆x)

(u (xi + tE∆x, yj)− u (xi, yj)

(xi + tE∆x)− xi− u (xi, yj)− u (xi − tW ∆x, yj)

xi − (xi − tW ∆x)

)

+2

(yj + tN∆y)− (yj − tS∆y)

(u (xi, yj + tN∆y)− u (xi, yj)

(yj + tN∆y)− yj− u (xi, yj)− u (xi, yj − tS∆y)

yj − (yj − tS∆y)

)

=2

(tE + tW )∆x

(ui+tE∆x,j − ui,j

tE∆x− ui,j − ui−tW ∆x,j

tW ∆x

)

+2

(tN + tS)∆y

(ui,j+tN∆y − ui,j

tN∆y− ui,j − ui,j−tS∆y

tS∆y

)

ou

−∆dud =2

∆x2

[− 1

tE (tE + tW )ui+tE∆x,j +

1tEtW

ui,j − 1tW (tE + tW )

ui−tW ∆x,j

](1.67)

+2

∆y2

[− 1

tS (tN + tS)ui,j−tS∆y +

1tN tS

ui,j − 1tN (tN + tS)

ui,j+tN∆y

].

Se (xi, yj) e um ponto interior distante da fronteira (isto e, nao adjacente a fronteira), entao t∗ = 1 e para esteponto vale a formula dos cinco pontos usual. Observe que a matriz obtida pelo esquema de Shortley-Wellernao e simetrica, em geral.

Embora a ordem de aproximacao do laplaciano para pontos proximos a fronteira e apenas 1, o esquemade Shortley-Weller e convergente de segunda ordem. No proximo capıtulo, provaremos que o problemadiscretizado possui solucao unica.

1.6 Exercıcios

1. Implemente os metodos discutidos neste capıtulo computacionalmente, verifique a precisao comparandocom a solucao exata e tambem a velocidade de convergencia.

2. Discretize o problema de Poisson com valor de fronteira de Dirichlet a seguir, usando a formula decinco pontos. −∆u = f (x, y) em (0, a)× (0, b) ,

u = g (x, y) sobre ∂ ((0, a)× (0, b)) ,

Implemente alguns exemplos deste problema computacionalmente e compare os resultados obtidos comas solucoes exatas.

3. Prove que a formula dos nove pontos compacta satisfaz o princıpio do maximo discreto.

4. Prove resultados equivalentes ao Lema 1.5 e ao Teorema 1.6 para a formula dos nove pontos compacta.

5. Investigue a ordem de convergencia do esquema de diferencas finitas misto: formula dos nove pontos nospontos interiores distantes da fronteira e formula dos cinco pontos para pontos adjacentes a fronteira.

6. Encontre um esquema de diferencas finitas de segunda ordem para a equacao de laplace tridimensionalem um paralelepıpedo reto. Escolha uma ordenacao apropriada dos pontos da malha e descreva amatriz de discretizacao obtida. Implemente o metodo no computador.

7. Mostre que o esquema de diferencas finitas em coordenadas polares introduzido neste capıtulo satisfazo princıpio do maximo discreto desde que o valor de u0 seja dado pela formula (1.56).


8. Mostre que se ∆d denota o esquema de diferencas finitas em coordenadas polares introduzido nestecapıtulo e Ω e o disco unitario, entao vale a estimativa a priori: se ud e uma solucao de

−∆dud = fd em Ωd,ud = 0 sobre ∂Ωd,

entao‖ud‖∞ 6 1

4‖∆dud‖∞ (1.68)

desde que o valor de u0 seja dado pela formula (1.56). Conclua que este esquema tem ordem deconvergencia 2.

9. Encontre os autovalores da matriz de discretizacao do esquema de diferencas finitas em coordenadaspolares e compare com os autovalores de Dirichlet do laplaciano no disco.

10. Discretize o problema de Poisson com valor de fronteira de Dirichlet para o anel:−∆u = f (r, θ) se R1 < r < R2 e 0 < θ < 2π,u (R1, θ) = g1 (θ)u (R2, θ) = g2 (θ) se 0 6 θ 6 2π.

Implemente alguns exemplos deste problema computacionalmente e compare os resultados obtidos comas solucoes exatas.

11. Mostre que tomando o “quadrado” da formula de tres pontos para o laplaciano unidimensional (es-quema de diferencas centradas para a derivada segunda) obtemos a seguinte formula de cinco pontospara o operador biharmonico unidimensional (esquema de diferencas centradas para a derivada quarta):

δ4ui =ui−2 − 4ui−1 + 6ui − 4ui+1 + ui+2

∆x4(1.69)

Usando a formula de Taylor, obtenha o expoente p tal que

δ4ui = u(4) (xi) + O (∆xp) .

12. O esquema de diferencas finitas mais simples para o operador biharmonico ∆2 em duas dimensoes e aseguinte formula de 13 pontos (para o caso ∆x = ∆y):

∆2u =1

∆x4

12 −8 2

1 −8 20 −8 12 −8 2

1

. (1.70)

Mostre que esta formula pode ser obtida a partir do “quadrado” da formula de cinco pontos parao laplaciano. Como a equacao biharmonica nao satisfaz o princıpio do maximo, a demonstracao daordem de convergencia deste esquema necessita de argumentos diferentes dos usados neste capıtulopara o laplaciano. Na realidade, dependendo de como as duas condicoes de fronteira sao discretizadas,a ordem de convergencia deste metodo pode ser O

(∆x3/2

)ou O

(∆x2

). Veja [Hackbusch], pag. 103 e

pags. 105-109, para detalhes e referencias.

Capıtulo 2

Existencia e Unicidade de SolucoesDiscretas

Determinar a existencia e unicidade de solucoes discretas para as matrizes de discretizacao obtidas viaesquemas de diferencas finitas atraves do calculo de seus autovalores como fizemos no capıtulo anterior paradiferencas centradas em uma dimensao e para a formula de cinco pontos e inviavel em geral (tente calcularos autovalores da matriz de discretizacao para a formula dos nove pontos, para o esquema em coordenadaspolares e para o esquema de Shortley-Weller). Neste capıtulo, desenvolveremos metodos mais gerais e maisfaceis de aplicar.

2.1 Normas Matriciais

Uma norma matricial no espaco vetorial Mn (C) das matrizes complexas n× n e uma norma vetorial quesatisfaz a propriedade submultiplicativa

‖AB‖ 6 ‖A‖ ‖B‖ (2.1)

para todas as matrizes A,B ∈ Mn (C). Algumas das normas mais importantes em Mn (C) sao as seguintes:

1. Norma l1

‖A‖1 =n∑

i,j=1

|aij | . (2.2)

De fato,

‖AB‖1 =n∑

i,j=1

∣∣∣∣∣n∑

k=1

aikbkj

∣∣∣∣∣ 6n∑

i,j,k=1

|aikbkj | 6n∑

i,j,k,l=1

|aikblj | =n∑

i,j=1

|aik|n∑

k,l=1

|blj | = ‖A‖1 ‖B‖1 .

2. Norma l2

‖A‖2 =

n∑

i,j=1

|aij |2

1/2

. (2.3)

Com efeito,

‖AB‖22 =n∑

i,j=1

∣∣∣∣∣n∑

k=1

aikbkj

∣∣∣∣∣

2

6n∑

i,j=1

(n∑

k=1

|aik|2)(

n∑

l=1

|blj |2)

=

n∑

i,k=1

|aik|2

n∑

j,l=1

|blj |2 = ‖A‖22 ‖B‖22 .

33


A norma l2 tambem e chamada norma euclidiana e, mais raramente e somente para matrizes, normade Schur, norma de Frobenius ou norma de Hilbert-Schmidt.

3. Norma l∞ modificada

A norma l∞‖A‖∞ = max

16i,j6n|aij | .

e uma norma vetorial no espaco das matrizes complexas, mas nao e uma norma matricial, pois se

A =[

1 11 1

],

entao

A2 =[

2 22 2

]

e portanto ∥∥A2∥∥∞ = 2 > 1 = ‖A‖∞ ‖A‖∞ .

Mas um multiplo escalar desta norma vetorial e uma norma matricial:

‖A‖n∞ = n max16i,j6n

|aij | . (2.4)

Com efeito,

‖AB‖n∞ = n max16i,j6n

∣∣∣∣∣n∑

k=1

aikbkj

∣∣∣∣∣ 6 n max16i,j6n

n∑

k=1

|aikbkj | 6 n max16i,j6n

n∑

k=1

‖A‖∞ ‖B‖∞

= n ‖A‖∞ n ‖B‖∞ = ‖AB‖n∞ .

4. Norma induzida

Dada uma norma vetorial |·| em Cn, ela induz uma norma matricial atraves da definicao

‖A‖ = max|x|=1

|Ax| = maxx6=0

|Ax||x| . (2.5)

De fato,

‖AB‖ = maxx 6=0

|ABx||x| = max

x6=0

( |ABx||Bx|

|Bx||x|

)6 max

x6=0

|ABx||Bx| max

x 6=0

|Bx||x| 6 max

y 6=0

|Ay||y| max

x 6=0

|Bx||x| = ‖A‖ ‖B‖ .

Esta norma tambem e chamada norma do operador. Ela satisfaz a propriedade muitas vezes util

|Ax| 6 ‖A‖ |x| (2.6)

para todo vetor x ∈ Cn.

5. Norma do maximo das somas das linhas

‖A‖L = max16i6n

n∑

j=1

|aij | . (2.7)

Esta norma e induzida pela norma vetorial l∞. De fato, se x = (x1, . . . , xn), temos

|Ax|∞ = max16i6n

∣∣∣∣∣∣

n∑

j=1

aijxj

∣∣∣∣∣∣6 max

16i6n

n∑

j=1

|aijxj | 6 max16i6n

n∑

j=1

|aij | |x|∞ = ‖A‖L |x|∞ ,


de modo quemax|x|=1

|Ax|∞ 6 ‖A‖L .

Supondo que a k-esima linha de A e nao-nula, definimos o vetor y = (y1, . . . , yn) ∈ Cn por

yi =

akj

|akj | se aij 6= 0,

1 se aij = 0.,

o que implica |y|∞ = 1, akjyj = |akj | e

max|x|∞=1

|Ax|∞ > |Ay|∞ = max16i6n

∣∣∣∣∣∣

n∑

j=1

aijyj

∣∣∣∣∣∣>

∣∣∣∣∣∣

n∑

j=1

akjyj

∣∣∣∣∣∣=

n∑

j=1

|akj | .

Isso vale para todo k, logo

max|x|∞=1

|Ax|∞ > max16k6n

n∑

j=1

|aij | = ‖A‖L .

6. Norma do maximo das somas das colunas

‖A‖C = max16j6n

n∑

i=1

|aij | . (2.8)

Esta norma e induzida pela norma vetorial l1. De fato, escrevendo A em termos de suas colunas

A = [A1 . . . An]

segue que‖A‖C = max

16j6n|Aj |1 .

Se x = (x1, . . . , xn), segue que

|Ax|1 = |x1A1 + . . . + xnAn|1 6n∑

i=1

|xiAi|1 =n∑

i=1

|xi| |Ai|1 6n∑

i=1

|xi| max16j6n

|Aj |1

= ‖A‖C

n∑

i=1

|xi| = ‖A‖C |x|1 ,

dondemax|x|1=1

|Ax|1 6 ‖A‖C .

Agora, se escolhermos y = ej , temos que |y|1 = 1 e

|Ay|1 = |Aj |1para todo k, logo

max|x|1=1

|Ax|1 > |Ay|1 = max16j6n

|Aj |1 = ‖A‖C .

7. p-normas

Este e o nome geral para as normas induzidas pela norma vetorial lp. O caso especial da norma induzidapela norma vetorial l2 (a norma vetorial euclidiana) e tambem chamada a norma espectral e satisfaz

‖|A|‖2 =√

λmax = max√

λ : λ e um autovalor de A∗A

.


De fato, A∗A e uma matriz hermitiana e possui autovalores nao-negativos, pois se A∗Ay = λy, entao

λ |y|22 = 〈y, λy〉2 = 〈y,A∗Ay〉2 = 〈Ay,Ay〉2 = |Ay|22e, alem disso, pela caracterizacao variacional dos autovalores de uma matriz hermitiana temos

λmax = maxx 6=0

〈A∗Ax, x〉2|x|22

= maxx 6=0

|Ax|22|x|22

.

Observe que a 2-norma e diferente da norma matricial l2. Note tambem que se A e uma matrizhermitiana, entao A∗A = A2 e ‖|A|‖2 e portanto o modulo do maior autovalor de A, isto e, a normaespectral de A e o raio espectral de A, definido como sendo o maior valor absoluto dos autovaloresde A:

ρ (A) = maxi=1,...,n

|λi| ,

8. Norma induzida por uma matriz invertıvel

Se ‖·‖ e uma norma matricial qualquer e se S e uma matriz invertıvel, entao

‖A‖S =∥∥S−1AS

∥∥ (2.9)

define uma norma matricial. Com efeito,

‖AB‖S =∥∥S−1ABS

∥∥ =∥∥S−1ASS−1BS

∥∥ 6∥∥S−1AS

∥∥∥∥S−1BS∥∥ = ‖A‖S ‖B‖S .

Lembramos que todas as normas em um espaco vetorial sao equivalentes, e isso vale em particular paranormas matriciais.

2.2 Matrizes Diagonalmente Dominantes

Definicao. Dizemos que uma matriz An×n e diagonalmente dominante se

|aii| >n∑

j=1j 6=i

|aij | para todo i = 1, . . . , n

e estritamente diagonalmente dominante se

|aii| >n∑

j=1j 6=i

|aij | para todo i = 1, . . . , n.

2.1 Proposicao. Se A e uma matriz estritamente diagonalmente dominante, entao A e invertıvel.

Prova. Uma matriz A e invertıvel se existe alguma norma matricial ‖·‖ tal que ‖I −A‖ < 1. De fato, seesta condicao e satisfeita, entao a inversa e dada explicitamente pela serie

A−1 =∞∑

k=0

(I −A)k. (2.10)

A condicao ‖I −A‖ < 1 garante a convergencia desta serie, pois a serie geometrica∑∞

k=0 rk tem raio deconvergencia 1; como para todo N temos

A

N∑

k=0

(I −A)k = [I − (I −A)]N∑

k=0

(I −A)k =N∑

k=0

(I −A)k −N+1∑

k=1

(I −A)k = I − (I −A)N+1,


tomando o limite quando N →∞, concluımos (2.10).Para provar a proposicao, denote por D a matriz diagonal cujas entradas diagonais sao as entradas

diagonais de A. Uma matriz estritamente diagonalmente dominante possui, por definicao, entradas diagonaisnao-nulas, logo D e uma matriz invertıvel. A matriz D−1A tem apenas 1’s na diagonal principal e semostramos que D−1A e invertıvel, isto implicara que A e invertıvel. Para provar isso, considere a matrizI −D−1A. Temos (

I −D−1A)ij

=

0 se i = j,−aij/aii se i 6= j.

Usemos a norma do maximo das somas das linhas. Para cada 1 6 i 6 n temosn∑

j=1

∣∣∣(I −D−1A

)ij

∣∣∣ =n∑

j=1j 6=i

∣∣∣∣aij

aii

∣∣∣∣ =1|aii|

n∑

j=1j 6=i

|aij | < 1,

logo∥∥I −D−1A

∥∥ < 1 e o resultado segue. ¥As vezes, exigir dominancia diagonal estrita em todas as linhas e pedir demais. Para certas matrizes,

dominancia diagonal junto com dominancia diagonal estrita em apenas uma linha e suficiente para garantira sua invertibilidade. As matrizes de discretizacao obtidas no capıtulo anterior satisfazem esta condicao(nas linhas correspondentes a pontos adjacentes a fronteira), e nenhuma delas e estritamente diagonalmentedominante. Por outro lado, esta condicao nao e suficiente para estabelecer a invertibilidade de uma matrizem geral, como o exemplo

4 2 10 1 10 1 1

demonstra. Precisamos de desenvolver varias ideias e ferramentas teoricas antes de provar a invertibilidadedas matrizes de discretizacao do capıtulo anterior.

2.3 Teorema dos Discos de Gershgorin

A primeira ferramenta teorica e o importante Teorema dos Discos de Gershgorin. Ele decorre da seguinteobservacao: se A e uma matriz complexa n × n, podemos sempre escrever A = D + B, onde D = diag(a11, . . . , ann) e a matriz diagonal formada pela diagonal principal de A e B consiste dos elementos restantesde A, possuindo uma diagonal principal nula. Se definirmos Aε = D + εB, entao A0 = D e A1 = A. Osautovalores de D sao a11, . . . , ann, enquanto que os autovalores de Aε devem estar localizados em vizinhancasdos pontos a11, . . . , ann, desde que ε seja suficientemente pequeno. O mesmo deve valer para os autovaloresda matriz A: eles devem estar contidos em discos centrados nos elementos a11, . . . , ann da diagonal principalse os discos sao suficientemente grandes. O Teorema de Gershgorin da uma estimativa precisa e simples decalcular para os raios destes discos em funcao das entradas restantes da matriz A. Denote o disco complexofechado de centro em a e raio R por

DR (a) = z ∈ C : |z − a| 6 R .

2.2 Teorema. (Teorema dos Discos de Gershgorin) Se A ∈ Mn (C) e

Ri (A) =n∑

j=1j 6=i

|aij | (2.11)

denota a soma dos valores absolutos dos elementos da linha i de A excetuando o elemento da diagonalprincipal, entao todos os autovalores de A estao contidos na uniao dos n discos de Gershgorin

G (A) =n⋃

i=1

DRi(A) (aii) . (2.12)


Alem disso, se uma uniao de k destes discos forma uma regiao que e disjunta dos n−k discos restantes,entao existem exatamente k autovalores de A nesta regiao.

Prova. Seja λ um autovalor de A e x = (x1, . . . , xn) 6= 0 um autovetor associado. Seja k um ındice tal que

|xk| > |xj | para j = 1, . . . , n,

isto e, xk e a coordenada de x de maior valor absoluto. Denotando por (Ax)k a k-esima coordenada do vetorAx = λx, temos

λxk = (Ax)k =n∑

j=1

akjxj

que e equivalente a

xk (λ− akk) =n∑

j=1j 6=k

akjxj .

Daı,

|xk| |λ− akk| 6n∑

j=1j 6=k

|akjxj | =n∑

j=1j 6=k

|akj | |xj | 6 |xk|n∑

j=1j 6=k

|akj | = |xk|Rk (A) ,

ou seja,|λ− akk| 6 Rk (A) .

Isso prova o resultado principal do Teorema de Gershgorin (como nao sabemos qual k e apropriado paracada autovalor λ, e um mesmo k pode servir para varios autovalores λ, tudo o que podemos afirmar e queos autovalores estao na uniao dos discos).

Para provar a segunda afirmacao, escreva A = D + B, onde D = diag (a11, . . . , ann) e defina

At = D + tB

para 0 6 t 6 1. Note queRi (At) = Ri (tB) = tRi (A) .

Para simplificar a notacao, assuma que a uniao dos primeiros k discos de Gershgorin

Gk (A) =k⋃

i=1

DRi(A) (aii)

satisfaz Gk (A) ∩ [G (A) \Gk (A)] = ∅. Temos

DRi(At) (aii) = z ∈ C : |z − aii| 6 Ri (At) = z ∈ C : |z − aii| 6 tRi (A) ⊂ DRi(A) (aii) ,

logoGk (At) ⊂ Gk (A)

eGk (A) ∩ [G (At) \Gk (At)] = ∅

para 0 6 t 6 1. Porque os autovalores sao funcoes contınuas das entradas de uma matriz, o caminho

λi (t) = λi (At)

e um caminho contınuo que liga λi (A0) = λi (D) = aii a λi (A1) = λi (A). Como λi (At) ∈ Gk (At) ⊂ Gk (A),concluımos que para cada 0 6 t 6 1 existem k autovalores de At em Gk (A); em particular, fazendo t = 1,


obtemos que Gk (A) possui pelo menos k autovalores de A. Da mesma forma, nao pode haver mais quek autovalores de A em Gk (A), pois os n − k autovalores restantes de A0 = D comecam fora do conjuntoGk (A) e seguem caminhos contınuos que permanecem fora de Gk (A). ¥A uniao G (A) dos discos de Gershgorin e conhecida como a regiao de Gershgorin. Observe que enquantonao podemos em geral afirmar com certeza que cada disco de Gershgorin possui um autovalor, a segundaafirmacao do teorema permite-nos fazer tal conclusao desde que os discos de Gershgorin sejam dois a doisdisjuntos.

O Teorema dos Discos de Gershgorin permite entender o resultado da Proposicao 2.1: se uma matriz A eestritamente diagonalmente dominante, entao os discos de Gershgorin DRi(A) (aii) nao interceptam a origem,logo 0 nao pode ser um autovalor para a matriz A, o que implica que A e invertıvel. Alem disso, se todosos elementos da diagonal principal de A sao reais e positivos, entao os autovalores de A estao localizados nosemiplano direito de C, de modo que se A e tambem simetrica, concluımos que todos os autovalores de Asao positivos.

A aplicacao mais obvia do Teorema dos Discos de Gershgorin e na estimativa dos autovalores de umamatriz, o que e importante se vamos usar os autovalores de matrizes de discretizacao para aproximar osautovalores do laplaciano:

Aplicacao 1. Pelo Teorema dos Discos de Gershgorin, os autovalores da matriz de discretizacao do lapla-ciano no intervalo (0, π) discretizado com n + 1 pontos (esquema de diferencas finitas centradas paraa derivada segunda unidimensional)

A =n2

π2

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 2

estao todos localizados no intervalo (A e simetrica, logo seus autovalores sao todos reais) centrado emx = 2n2/π2 de raio 2n2/π2, ou seja, no intervalo

[0, 4n2/π2

]. Em particular o maior autovalor de A

nao pode exceder 4n2/π2. Como os autovalores do laplaciano neste intervalo sao da forma λj = j2,para termos esperanca em aproximar o autovalor λj por autovalores da matriz A precisamos quej2 6 4n2/π2, isto e, precisamos discretizar o intervalo (0, π) com

n > π

2j

pontos. Isso da uma estimativa bastante grosseira do quao refinada a nossa malha precisa ser paraaproximar os autovalores do laplaciano. Na pratica, vimos que apenas os primeiros autovalores deA aproximam bem os primeiros autovalores do laplaciano e portanto precisamos de uma malha comum numero muito maior de pontos. Observe que uma estimativa semelhante vale para a matriz dediscretizacao M fornecida pela formula de cinco pontos no quadrado (0, π)2 quando tomamos ∆x =∆y = π/n: como os autovalores de M estao localizados no intervalo de centro em x = 4n2/π2 de raio4n2/π2, isto e, em

[0, 8n2/π2

], precisamos de

n > π

2√

2

√i2 + j2

pontos no eixos horizontal e vertical para aproximar o autovalor i2 + j2. Por outro lado, no casobidimensional isso implica em uma matriz de discretizacao da ordem de i2 + j2. ¤

Usos mais refinados do Teorema de Gershgorin permitem obter conhecimento mais preciso sobre ondeos autovalores da matriz se encontram e correspondentemente melhores estimativas para o raio espectral


de uma matriz. Por exemplo, como A e At possuem os mesmos autovalores, existe um teorema dos discosde Gershgorin equivalente para as colunas de uma matriz. Em particular, todos os autovalores de A estaolocalizados na intersecao destas duas regioes: G (A)∩G (At). Isso implica a seguinte estimativa simples parao raio espectral de uma matriz complexa:

2.3 Corolario. Se A ∈ Mn (C), entao

ρ (A) 6 min

max

i=1,...,n

n∑

j=1

|aij | , maxj=1,...,n

n∑

i=1

|aij | = min (‖A‖L , ‖A‖C) .

Prova. O ponto no i-esimo disco de Gershgorin que e mais distante da origem tem modulo

|aii|+ Ri (A) =n∑

j=1

|aij |

e um resultado semelhante vale para as colunas de A. ¥O resultado do Corolario 2.3 nao e surpreendente em vista do raio espectral de uma matriz ser menor quequalquer norma matricial (veja o proximo capıtulo). Um resultado melhor pode ser obtido uma vez quese observa que A e S−1AS tambem possuem os mesmos autovalores, qualquer que seja a matriz invertıvelS. Em particular, quando S = D = diag (p1, . . . , pn) e uma matriz diagonal com todos os seus elementospositivos, isto e, pi > 0 para todo i, aplicando o Teorema de Gershgorin a matriz

D−1AD =(

pj

piaij

)

e a sua transposta, obtemos o seguinte resultado que permite obter uma estimativa arbitrariamente boa dosautovalores de A:

2.4 Corolario. Se A ∈ Mn (C) e p1, . . . , pn > 0, entao todos os autovalores de A estao contidos em

G(D−1AD

) ∩G(DAtD−1

)=

n⋃

i=1

z ∈ C : |z − aii| 6 1pi

n∑

j=1j 6=i

pj |aij |

(2.13)

∩n⋃

i=1

z ∈ C : |z − aii| 6 pj

n∑

i=1i 6=j

1pi|aij |

.

Em particular,

ρ (A) 6 minp1,...,pn>0

max

i=1,...,n

1pi

n∑

j=1

pj |aij | , maxj=1,...,n

pj

n∑

i=1

1pi|aij |

. (2.14)

2.4 Propriedade FC

Na nossa busca por propriedades para matrizes diagonalmente dominantes que garantirao a sua invertibil-idade, uma observacao fundamental e a de que se A e uma matriz diagonalmente dominante, entao 0 naopode ser um ponto interior de nenhum disco de Gershgorin. De fato, se λ e um autovalor de A interior aalgum disco de Gershgorin entao devemos ter desigualdade estrita

|λ− aii| < Ri (A) =n∑

j=1j 6=i

|aij |


para algum i. Se 0 e um autovalor de A interior a algum disco de Gershgorin, entao

|aii| <n∑

j=1j 6=i

|aij |

para algum i e A nao pode ser diagonalmente dominante na linha i.Uma condicao equivalente para que um autovalor λ de A nao seja um ponto interior de nenhum disco de

Gershgorin e que

|λ− aii| > Ri (A) =n∑

j=1j 6=i

|aij | para todo i = 1, . . . , n.

Tais pontos λ na regiao de Gershgorin G (A) (nao necessariamente autovalores de A) constituem precisa-mente a fronteira ∂G (A) da regiao de Gershgorin. Chamaremos a fronteira de um disco de Gershgorinz ∈ C : |z − aii| = Ri (A) um cırculo de Gershgorin.

2.5 Lema. Seja A ∈ Mn (C) e λ um autovalor de A que nao e um ponto interior de nenhum disco deGershgorin. Seja x = (x1, . . . , xn) 6= 0 um autovetor associado a λ e k um ındice tal que

|xk| > |xj | para j = 1, . . . , n.

Se i e qualquer ındice tal que|xi| = |xk|

entao o i-esimo cırculo de Gershgorin passa por λ. Se, alem disso,

aij 6= 0,

entao|xj | = |xk|

e o j-esimo cırculo de Gershgorin tambem passa por λ.

Prova. Como na demonstracao do Teorema de Gershgorin, temos

|xi| |λ− aii| 6n∑

j=1j 6=k

|aijxj | =n∑

j=1j 6=k

|aij | |xj | 6 |xk|n∑

j=1j 6=k

|aij | = |xk|Ri (A) (2.15)

para todo ındice i. Logo, se |xi| = |xk|, temos

|λ− aii| 6 Ri (A) .

Como por hipotese|λ− aii| > Ri (A)

para todo ındice i, segue que|λ− aii| = Ri (A) .

Em geral, |xi| = |xk| implica que as desigualdades em (2.15) sao identidades; em particular,

n∑

j=1j 6=k

|aij | |xj | = |xi|n∑

j=1j 6=k

|aij |


donden∑

j=1j 6=k

|aij | (|xi| − |xj |) = 0.

Esta e uma soma de termos nao-negativos, pois |xi| > |xj |, logo se aij 6= 0 necessariamente devemos ter|xj | = |xi| = |xk|. ¥

Este lema tecnico tem as seguintes consequencias uteis:

2.6 Teorema. Seja A ∈ Mn (C) uma matriz cujas entradas sao todas nao-nulas e seja λ um autovalor deA que nao e um ponto interior de nenhum disco de Gershgorin. Entao todo cırculo de Gershgorinde A passa por λ (isto e, λ esta na intersecao de todos os cırculos de Gershgorin de A) e se x =(x1, . . . , xn) 6= 0 e um autovetor associado a λ entao

|xi| = |xj | para todos i, j = 1, . . . , n.

Prova. Decorre diretamente do lema anterior. ¥

2.7 Corolario. Se A ∈ Mn (C) e uma matriz cujas entradas sao todas nao-nulas e diagonalmente dominante

tal que |aii| >n∑

j=1j 6=i

|aij | para pelo menos alguma linha i, entao A e invertıvel.

Prova. Pois, como A e diagonalmente dominante, se 0 e um autovalor de A entao 0 nao pode ser um pontointerior de nenhum disco de Gershgorin. Por outro lado, pelo teorema anterior, segue que todo cırculo deGershgorin passa por 0. Entretanto, o i-esimo cırculo de Gershgorin centrado em aii e com raio Ri < |aii|nao pode passar por 0. Concluımos que 0 nao e um autovalor de A, logo A e invertıvel. ¥

Na verdade, usando com maior cuidado a informacao dada pelo Lema 2.5 podemos obter resultados aindamelhores:

Definicao. Dizemos que uma matriz A = (aij) ∈ Mn (C) satisfaz a propriedade FC se para todo par deinteiros distintos i, j existe uma sequencia de inteiros distintos i1 = i, i2, i3, . . . , im−1, im = j, com1 6 m 6 n, tais que todas as entradas matriciais

ai1i2 , ai2i3 , . . . , aim−1im

sao nao-nulas.

Por exemplo, a matriz diagonalmente dominante nao-invertıvel

4 2 10 1 10 1 1

,

ja vista anteriormente, nao satisfaz a propriedade FC porque o par 2, 1 nao admite tal sequencia (a unicasequencia possıvel e a23, a31). Ja qualquer par de inteiros distintos i, j tal que aij 6= 0 admite a sequenciatrivial nao-nula aij , de modo que uma matriz cujas entradas nao-diagonais sao todas nao-nulas satisfaz apropriedade FC. O significado da abreviatura “FC”, ou “fortemente conexo”, ficara claro mais adiante.

2.8 Teorema. Seja A ∈ Mn (C) uma matriz que satisfaz a propriedade FC e seja λ um autovalor de A quenao e um ponto interior de nenhum disco de Gershgorin. Entao todo cırculo de Gershgorin de A passapor λ (isto e, λ esta na intersecao de todos os cırculos de Gershgorin de A) e se x = (x1, . . . , xn) 6= 0e um autovetor associado a λ entao

|xi| = |xj | para todos i, j = 1, . . . , n.


Prova. Seja x = (x1, . . . , xn) 6= 0 um autovetor associado a λ e i um ındice tal que

|xi| > |xk| para k = 1, . . . , n.

Pelo Lema 2.5,|λ− aii| = Ri (A) .

Seja j 6= i qualquer outro ındice e i1 = i, i2, i3, . . . , im−1, im = j, com 1 6 m 6 n, ındices tais que todas asentradas matriciais

aii2 , ai2i3 , . . . , aim−1j 6= 0.

Como aii2 6= 0, segue da segunda afirmativa do Lema 2.5 que |xi2 | = |xi|. Mas entao ai2i3 6= 0 e portanto|xi3 | = |xi2 | = |xi|. Prosseguindo desta forma, concluımos que

|xi| = |xi2 | = . . .∣∣xim−1

∣∣ = |xj | .

Em particular, segue novamente do Lema 2.5 que o j-esimo cırculo de Gershgorin passa por λ. Como j earbitrario, isso prova o teorema. ¥

2.9 Corolario. Se A ∈ Mn (C) e uma matriz que satisfaz a propriedade FC e diagonalmente dominante tal

que |aii| >n∑

j=1j 6=i

|aij | para pelo menos alguma linha i, entao A e invertıvel.

Prova. Segue do teorema anterior da mesma forma que o Corolario 2.7 segue do Teorema 2.6. ¥Vamos tentar entender melhor o significado da propriedade FC. Note que ela se refere apenas a localizacao

dos elementos nao-nulos de A fora da diagonal principal – os elementos da diagonal principal e os valoresespecıficos dos elementos fora da diagonal principal sao irrelevantes. Isso motiva as seguintes definicoes:

Definicao. Dada uma matriz A = (aij) ∈ Mn (C) definimos o modulo da matriz A como sendo a matriz

|A| = (|aij |)

cujos elementos sao os modulos dos elementos da matriz A e a matriz indicadora de A como sendoa matriz

M (A) = (µij) ,

onde

µij =

1 se aij 6= 0,0 se aij = 0.

O conceito de uma sequencia de entradas nao-nulas da matriz A que aparece na definicao da propriedadeFC pode ser visualizado em termos de caminhos em um grafo associado a A:

Definicao. Dada uma matriz A ∈ Mn (C), o grafo direcionado de A e o grafo direcionado Γ (A) com nnodos P1, . . . , Pn tais que existe um arco direcionado em Γ (A) de Pi a Pj se e somente se aij 6= 0.

Um caminho direcionado γ em um grafo Γ e uma sequencia de arcos Pi1Pi2 , Pi2Pi3 , . . . em Γ. Ocomprimento de um caminho direcionado e o numero de arcos sucessivos no caminho direcionado. Umciclo e um caminho direcionado que comeca e termina no mesmo no.

Dizemos que um grafo direcionado e fortemente conexo se entre qualquer par de nodos distintosPi, Pj ∈ Γ existir um caminho direcionado de comprimento finito que comeca em Pi e termina em Pj .

Observe que quando Γ e um grafo direcionado com n nodos, se existe um caminho direcionado entre doisnodos de Γ, entao sempre existe um caminho direcionado entre estes dois nodos de comprimento menor queou igual a n− 1.


2.10 Teorema. A ∈ Mn (C) satisfaz a propriedade FC se e somente se Γ (A) e fortemente conexo.

Verificar a propriedade FC a partir do grafo direcionado de A pode ser impraticavel se o tamanho damatriz for muito grande. Existe um metodo computacional mais explıcito para faze-lo:

2.11 Teorema. Sejam A ∈ Mn (C) e Pi, Pj nodos de Γ (A). Existe um caminho direcionado de compri-mento m em Γ (A) de Pi para Pj se e somente se

(|A|m)ij 6= 0

ou, equivalentemente, se e somente se[M (A)m]ij 6= 0.

Prova. Provaremos o teorema por inducao. Para m = 1 a afirmativa e trivial. Para m = 2, temos

(|A|2

)ij

=n∑

k=1

(|A|)ik (|A|)kj =n∑

k=1

|aik| |akj | ,

de modo que(|A|2

)ij6= 0 se e somente se aik, akj sao ambos nao-nulos para algum ındice k. Mas isso e

equivalente a dizer que existe um caminho direcionado de comprimento 2 em Γ (A) de Pi para Pj .Em geral, supondo a afirmativa provada para m, temos

(|A|m+1

)ij

=n∑

k=1

(|A|m)ik (|A|)kj =n∑

k=1

(|A|m)ik |akj | 6= 0

se e somente se (|A|m)ik , akj sao ambos nao-nulos para algum ındice k. Por hipotese de inducao, isso eequivalente a existir um caminho direcionado de comprimento m em Γ (A) de Pi para Pk e um caminhodirecionado de comprimento 1 em Γ (A) de Pk para Pj , isto e, um caminho direcionado de comprimentom + 1 em Γ (A) de Pi para Pj . O mesmo argumento vale para M (A). ¥

Definicao. Seja A = (aij) ∈ Mn (C). Dizemos que A > 0 se aij > 0 para todos 1 6 i, j 6 n e que A > 0 seaij > 0 para todos 1 6 i, j 6 n.

2.12 Corolario. Seja A ∈ Mn (C). Existe um caminho direcionado de comprimento m em Γ (A) de cadanodo Pi para cada nodo Pj se e somente se

|A|m > 0

ou, equivalentemente, se e somente seM (A)m

> 0.

2.13 Corolario. Seja A ∈ Mn (C). A satisfaz a propriedade FC se e somente se

(I + |A|)n−1> 0

ou, equivalentemente, se e somente se

[I + M (A)]n−1> 0.

Prova. Temos

(I + |A|)n−1 = I + (n− 1) |A|+(

n− 12

)|A|2 + . . . +

(n− 1n− 3

)|A|n−1 + |A|n−1

> 0


se e somente se para cada par de ındices i, j com i 6= j pelo menos um dos termos |A| , |A|2 , . . . , |A|n−1

tem uma entrada positiva em (i, j). Pelo Teorema 2.11, isso ocorre se e somente se existe algum caminhodirecionado em Γ (A) de Pi para Pj com comprimento 6 n−1. Isto e equivalente a A satisfazer a propriedadeFC. O mesmo argumento vale para M (A). ¥Em geral, a maneira como uma matriz foi obtida (como as nossas matrizes de discretizacao; veja a ultimasecao do capıtulo) torna clara se elas sao matrizes que satisfazem a propriedade FC ou nao. Se issonao e possıvel, e pretende-se verificar a propriedade FC atraves do Corolario 2.13, e preferıvel calcular[I + M (A)]n−1, ja que M (A) e uma matriz composta apenas de 0’s e 1’s.

2.5 Matrizes Irredutıveis

Lembre-se que uma matriz de permutacao P e uma matriz quadrada cujas entradas sao todas 0 ou 1 e,alem disso, em cada linha e em cada coluna de P existe exatamente um 1. Em particular, P e uma matrizortogonal, de modo que P−1 = P t, isto e, a inversa de P tambem e uma matriz de permutacao. Um casoespecial de uma matriz de permutacao e uma matriz de transposicao, que e uma matriz de permutacao Tigual a matriz identidade exceto em duas posicoes, isto e, para algum par de ındices fixado k, l temos

Tij =

δij se (i, j) 6= (k, l) , (l, k) , (k, k) ou (l, l) ,1 e (i, j) = (k, l) ou se (i, j) = (l, k) ,0 se (i, j) = (k, k) ou se (i, j) = (l, l) .

Matrizes de transposicao sao simetricas. O efeito de multiplicar uma matriz A por uma matriz de transposicaoa esquerda e trocar a posicao de duas linhas da matriz A (no caso acima, as linhas k e l), enquanto que amultiplicacao de A por uma matriz de transposicao a direita muda a posicao de duas colunas de A (no casoacima, as colunas k e l).

TA =

1 0 0 00 0 1 00 1 0 00 0 0 1

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44

=

a11 a12 a13 a14

a31 a32 a33 a34

a21 a22 a23 a24

a41 a42 a43 a44

,

AT =

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44

1 0 0 00 0 1 00 1 0 00 0 0 1

=

a11 a13 a12 a14

a21 a23 a22 a24

a31 a33 a32 a34

a41 a43 a42 a44

.

Pode-se provar que toda matriz de permutacao P e o produto de matrizes de transposicao P = T1 . . . Tm;em particular, P t = Tm . . . T1. A matriz

P tAP = Tm . . . T1AT1 . . . Tm

e portanto obtida atraves da permutacao de linhas e colunas de A, de modo que nenhum novo elemento ecriado ou algum elemento existente de A destruıdo.

Definicao. Dizemos que uma matriz A ∈ Mn (C) e redutıvel se existe alguma matriz de permutacao P ealgum inteiro 1 6 m 6 n− 1 tal que

P tAP =[

B C0 D

]

onde B e uma matriz m×m, D e uma matriz (n−m)× (n−m), C e uma matriz m× (n−m) e 0 ea matriz nula (n−m)×m. Caso contrario, dizemos que A e irredutıvel.


Da definicao vemos que se |A| > 0, entao A e irredutıvel, e para que A seja redutıvel, ela precisa ter pelomenos n− 1 zeros (caso m = 1). A motivacao para este nome e a seguinte. Suponha que queiramos resolvero sistema Ax = b e que A seja redutıvel. Entao, se escrevermos

A = P tAP =[

B C0 D

],

teremos Ax = PAP tx = b ou AP tx = P tb; denotando x = P tx e b = P tb, resolver o sistema Ax = b e entaoequivalente a resolver o sistema

Ax = b.

Escrevendo

x =[

yz

], b =

[b1

b2

]

onde y, b1 ∈ Cm e z, b2 ∈ Cn−m, este sistema e por sua vez equivalente ao sistema

By + Cz = b1

Dz = b2

Se resolvermos primeiro Dz = b2 e utilizarmos o valor de z encontrado na primeira equacao resolvendoBy = b1 − Cz, teremos reduzido o problema original a dois problemas menores, mais faceis de resolver.

2.14 Teorema. Uma matriz A ∈ Mn (C) e irredutıvel se e somente se

(I + |A|)n−1> 0

ou, equivalentemente, se e somente se

[I + M (A)]n−1> 0.

Prova. Para provar o resultado, mostraremos que A e redutıvel se e somente se (I + |A|)n−1 possui pelomenos uma entrada nula.

Assuma primeiramente que A e redutıvel, de modo que para alguma matriz de permutacao P tenhamos

A = P

[B C0 D

]P t =: PAP t.

Observe que|A| = ∣∣PAP t

∣∣ = P∣∣A∣∣ P t,

ja que o efeito de P e apenas trocar linhas e colunas. Alem disso, note que

Ak

=[

Bk Ck

0 Dk

]

para alguma matriz Ck. Logo, como

(I + |A|)n−1 =(I + P

∣∣A∣∣ P t

)n−1= P

(I +

∣∣A∣∣)n−1

P t

= P

[I + (n− 1) |A|+

(n− 1

2

)|A|2 + . . . +

(n− 1n− 3

)|A|n−1 + |A|n−1

]P t

e todos os termos dentro dos colchetes sao matrizes que tem um bloco (n−m)×m nulo no canto esquerdoinferior, segue que (I + |A|)n−1 e redutıvel, logo possui entradas nulas e nao pode ser positiva.


Reciprocamente, suponha que (I + |A|)n−1 possui pelo menos uma entrada nula. Como

(I + |A|)n−1 = I +n−1∑m=1

(n− 1

m

)|A|m ,

(I + |A|)n−1 nao possui entradas diagonais nulas, logo podemos assumir que para algum par i 6= j temos[(I + |A|)n−1

]ij

= 0, o que implica [|A|m]ij = 0 para todo 1 6 m 6 n− 1. Pelo Teorema 2.11 (e observacao

imediatamente posterior a definicao de grafo direcionado), nao existe um caminho direcionado em Γ (A) decomprimento finito entre Pi e Pj . Defina os conjuntos de nodos

S1 := Pk : Pk = Pj ou existe um caminho direcionado em Γ (A) entre Pk e Pj ,

S2 = [ nodos de Γ (A)] \S1.

Por definicao destes conjuntos, nao pode existir nenhum caminho de algum nodo de S2 para algum nodo deS1, logo [|A|m]lk = 0 se Pl ∈ S2 e Pk ∈ S1. E ambos os conjuntos sao nao-vazios, pois Pj ∈ S1 e Pi ∈ S2.Renomeando os nodos de modo que

S1 =

P1, . . . , Pm

,

S2 =

Pm+1, . . . , Pn

,

segue que existe uma matriz de permutacao P tal que

P tAP =[

B C0 D

].

De fato, P e justamente a matriz de permutacao que troca as colunas de tal forma que as variaveis anteriorescorrespondentes aos nodos P1, . . . , Pm no sistema Ax = b sao as novas m primeiras variaveis do sistema linearAx = b; como nao existe nenhum caminho direcionado entre nenhum dos nodos Pm+1, . . . , Pn e qualquer umdos nodos P1, . . . , Pm, temos aij = 0 para m + 1 6 i 6 n e 1 6 j 6 m pelo Teorema 2.11. ¥

2.15 Corolario. Uma matriz A ∈ Mn (C) e irredutıvel se e somente se ela satisfaz a propriedade FC.

2.16 Proposicao. Se A e uma matriz irredutıvel, diagonalmente dominante tal que |aii| >n∑

j=1j 6=i

|aij | para

pelo menos alguma linha i, entao A e invertıvel.

Alem disso, se A e hermitiana e todos os elementos da diagonal principal de A sao positivos, entaotodos os autovalores de A sao positivos.

Prova. O resultado segue do Teorema 2.14, do Corolario 2.9 e do Teorema dos Discos de Gershgorin (vejacomentarios apos o Teorema 2.2). ¥

2.6 Invertibilidade de Matrizes de Discretizacao

Os resultados obtidos nas secoes anteriores fornecem uma demonstracao alternativa de que as matrizesde discretizacao do capıtulo anterior (tanto no caso unidimensional, quanto no caso bidimensional) saoinvertıveis, sem a necessidade de se calcular os seus autovalores.


2.6.1 Esquemas de Diferencas Finitas para o Intervalo e para o Retangulo

E facil ver que todas as matrizes de discretizacao obtidas no capıtulo anterior para o intervalo e para oretangulo (isto e, os esquemas unidimensionais de tres pontos e cinco pontos, e os esquemas bidimensionaisde cinco e nove pontos, compacto ou nao-compacto) sao matrizes diagonalmente dominantes com dominanciadiagonal estrita nas linhas correspondentes a pontos adjacentes a fronteira. Alem disso, elas sao matrizesirredutıveis porque elas satisfazem a propriedade FC. De fato, cada ındice i da matriz corresponde a umponto interior Pi da malha e aij 6= 0 sempre que Pi e Pj sao pontos vizinhos naqueles esquemas. Entao,dados dois pontos distintos Pi, Pj e facil encontrar uma sequencia de ındices i1 = i, i2, i3, . . . , im−1, im = j,com 1 6 m 6 n, tais que todas as entradas matriciais

ai1i2 , ai2i3 , . . . , aim−1im

sao nao-nulas: no caso unidimensional, basta percorrer a malha diretamente de Pi ate Pj (andando a partirde Pi sempre para a direita ou sempre para a esquerda, conforme o caso, ate encontrar Pj), e no casobidimensional basta usar qualquer caminho interior de Pi ate Pj (pode-se usar a ordem lexicografica parapercorrer a malha, ou a ordem lexicografica inversa, dependendo das posicoes relativas de Pi e Pj ; no entanto,estes caminhos sao mais longos que o necessario). Em outras palavras, identificando as malhas de pontosinternos com os grafos direcionados da matriz de discretizacao, de modo que existe um arco direcionado entredois pontos da malha se e somente se eles sao vizinhos, os esquemas de discretizacao considerados garantemque estes grafos sao fortemente conexos.

As matrizes obtidas atraves de diferencas finitas em geral sao irredutıveis, pois elas satisfazem a pro-priedade FC. E difıcil imaginar um esquema de diferencas finitas para uma malha sobre um domınio conexoem que nao houvesse um caminho direcionado entre pontos vizinhos (isto e, em que tivessemos aij = 0para dois pontos vizinhos Pi e Pj). Outra maneira de pensar sobre isso e observar que se uma matriz dediscretizacao fosse (apos permutacao de linhas e colunas) da forma

[B C0 D

],

isso implicaria que um conjunto de pontos da malha (os correspondentes ao bloco D) teriam diferencasfinitas independentes do conjunto dos pontos restantes da malha (os correspondentes ao bloco D); piorainda, estes ultimos poderiam ter diferencas finitas dependentes dos primeiros (ja que o bloco C poderiaser nao-nulo). Em ultima analise, seria possıvel reduzir o problema de resolver o sistema linear associado adiscretizacao a dois problemas mais simples. E difıcil imaginar um esquema de diferencas finitas com estapropriedade, embora talvez possa ocorrer em algum domınio com geometria altamente irregular em que amalha de pontos interiores se dividisse em essencialmente duas malhas independentes. Tal situacao deve serevitada com cuidado na hora de discretizar tais regioes.

2.6.2 Esquema de Coordenadas Polares

As mesmas observacoes anteriores valem para a matriz de discretizacao obtida atraves do esquema de coorde-nadas polares do capıtulo anterior, isto e, ela satisfaz a propriedade FC. Para verificar que ela e diagonalmentedominante, note que para todas as linhas, exceto a primeira que deve ser tratada separadamente, temos

|aii| = γi =1ri

ri+1/2 + ri−1/2

∆r2+

2r2i

1∆θ2

.

Alem disso, para todas as linhas, excetuando a primeira e as linhas correspondentes a pontos adjacentes afronteira do disco temos

n∑

j=1j 6=i

|aij | = αi + βi + 2δi =1

∆r2

ri−1/2

ri+

1∆r2

ri+1/2

ri+

2r2i

1∆θ2

= |aii| .


Nestas linhas existe dominancia diagonal, enquanto que nas linhas correspondentes a pontos adjacentes afronteira do disco temos

(n−1)×m+1∑

j=1j 6=i

|aij | = αi + 2δi < |aii| ,

isto e, temos dominancia diagonal estrita. Finalmente, para a primeira linha tambem temos dominanciadiagonal, pois

|a00| = 4∆r2

,

(n−1)×m+1∑

j=1j 6=0

|a0j | = m2π

∆θ

∆r2= 4

m

2π

∆θ

∆r2=

4∆r2

= |a00| .

2.6.3 Esquema de Shortley-Weller

Se a geometria e razoavelmente regular, o esquema de Shortley-Weller para o problema de Dirichlet devesatisfazer a propriedade FC : aij 6= 0 sempre que Pi e Pj sao pontos internos vizinhos, e se a geometria nao ealtamente irregular (por exemplo, se o domınio e “razoavelmente” convexo) existe um caminho direcionado deum ponto interno arbitrario a qualquer outro ponto interno da malha passando apenas por pontos internos dodomınio. Caso contrario, a matriz de discretizacao obtida pode deixar de ser irredutıvel, mas isso deve ocorrerapenas devido a quebra da malha de pontos internos em varias submalhas desconexas, e cada submalha porsi so deve ser fortemente conexa. Portanto, a matriz de discretizacao total deve ser uma matriz em blocos,cada bloco satisfazendo a propriedade FC, logo a matriz e invertıvel.

Capıtulo 3

Metodos Iterativos para a Resolucaode Sistemas Lineares

Neste capıtulo investigaremos metodos iterativos para a resolucao de sistemas lineares

Ax = b.

Embora a matriz A que temos em mente e em geral uma matriz grande e esparsa, do tipo que apareceem esquemas de diferencas finitas, os metodos considerados aqui requerem apenas que A seja uma matrizinvertıvel com todas as entradas diagonais aii nao-nulas.

Metodos iterativos requerem um chute inicial x0, um vetor inicial que aproxima a solucao exata x (senao ha nenhuma informacao disponıvel sobre a solucao exata, de modo que nao temos como construir ochute inicial de forma inteligente, x0 pode ser uma aproximacao muito ruim de x). Uma vez que x0 e dado,o metodo iterativo gera a partir de x0 uma nova aproximacao x1, que esperamos deve aproximar melhor asolucao exata. Em seguida, x1 e usada para gerar uma nova melhor aproximacao x2 e assim por diante.Desta forma, gera-se uma sequencia de vetores

(xk

)que espera-se convergir para x. Como na pratica nao

podemos iterar para sempre, algum criterio de parada deve ser estabelecido a priori. Uma vez que xk estejasuficientemente proximo da solucao exata quanto se precise, de acordo com uma margem de tolerancia aceita,para-se o processo de iteracao e aceita-se xk como a solucao aproximada adequada para o problema. Porexemplo, o criterio de parada pode ser estabelecido atraves de uma cota de tolerancia τ : quando

∥∥b−Axk∥∥ < τ

ou quando ∥∥xk+1 − xk∥∥ < τ

as iteracoes sao interrompidas e o ultimo valor aproximado obtido e aceito como a melhor aproximacao dasolucao dentro das circunstancias.

Os metodos discutidos neste capıtulo nao necessitam de um bom chute inicial (embora, e claro, quantomelhor o chute inicial, menor o numero de iteracoes necessarias para se chegar a solucao aproximada com aprecisao especificada).

3.1 Metodos Iterativos Lineares

Nesta secao apresentamos alguns exemplos classicos de metodos iterativos lineares. Na proxima secao dare-mos condicoes necessarias e suficientes para estabelecer a sua convergencia.

50


3.1.1 Metodo de Jacobi

O primeiro metodo iterativo (que ja foi descrito como o mais lento para convergir, embora isso realmentedepende da matriz A do sistema) e o algoritmo de Jacobi. Escrevendo o sistema Ax = b na forma

n∑j=1

a1jxj = b1

...n∑

j=1

anjxj = bn

,

se aii 6= 0 para todo i, cada xi pode ser isolado na i-esima equacao e escrito na forma

xi =1aii

bi −

n∑

j=1j 6=i

aijxj

.

Isso sugere definir um metodo iterativo da seguinte forma: suposto xk =(xk

1 , . . . , xkn

)obtido no passo

anterior, obtemos xk+1 =(xk+1

1 , . . . , xk+1n

)por

xk+1i =

1aii

bi −

n∑

j=1j 6=i

aijxkj

. (3.1)

No caso da formula de cinco pontos para o problema de Poisson com ∆x = ∆y, como a equacao paracada ponto (i, j) e dada por

−ui,j−1 − ui,j+1 + 4ui,j − ui−1,j − ui+1,j = ∆x2fi,j

o metodo de Jacobi euk+1

i,j =14

(uk

i,j−1 + uki,j+1 + uk

i−1,j + uki+1,j + ∆x2fi,j

). (3.2)

No caso especial da equacao de Laplace (f = 0) com condicao de fronteira de Dirichlet nao-nula, o metodode Jacobi e simplesmente a propriedade do valor medio discreta

uk+1i,j =

14

(uk


i−1,j + uki+1,j

). (3.3)

Em outras palavras, calculados os valores de u em todos os pontos da malha na iteracao anterior, o novovalor de u em um ponto interior da malha nesta iteracao e calculado atraves da media dos seus quatropontos vizinhos. Os valores iniciais de u nos pontos interiores da malha para a primeira iteracao (isto e, ochute inicial) podem ser atribuidos arbitrariamente ou atraves de algum argumento razoavel; por exemplo,podemos utilizar uma media ponderada dos valores de fronteira para o valor inicial em cada ponto interiorda malha, de acordo com a posicao do ponto em relacao aos pontos das quatro fronteiras discretizadas.

Em forma matricial, o algoritmo de Jacobi pode ser descrito da seguinte forma. Denotando por D = diag(a11, . . . , ann) a matriz diagonal cujas entradas sao as entradas diagonais de A, temos que

xk+1 = D−1[(D −A) xk + b

](3.4)

ouxk+1 = D−1

(Cxk + b

)(3.5)

onde C = D −A e a matriz consistindo dos elementos restantes de A fora da diagonal principal.


3.1.2 Metodo de Gauss-Seidel

Um metodo iterativo que converge cerca de duas vezes mais rapido que o metodo de Jacobi (pelo menos emvarias aplicacoes) e o metodo de Gauss-Seidel, onde os valores de x sao atualizados dentro de cada iteracao,sem esperar pela proxima. Em outras palavras, obtido o valor de xk+1

l este e usado no lugar de xkl no calculo

seguinte. No sistema Ax = b em que aii 6= 0 para todo i, como antes isolamos cada xi na i-esima equacaomas desta vez escrevemos

xi =1aii

bi −

i−1∑

j=1

aijxj +n∑

j=i+1

aijxj

.

Entao definimos

xk+1i =

1aii

bi −

i−1∑

j=1

aijxk+1j +

n∑

j=i+1

aijxkj

(3.6)

pois os valores xk+11 , . . . , xk+1

i−1 ja foram computados nesta iteracao, enquanto que os valores xki+1, . . . , x

kn sao

fornecidos pela iteracao anterior.Por exemplo, no caso da equacao de Laplace, poderıamos utilizar a formula

uk+1i,j =

14

(uk+1

i,j−1 + uki,j+1 + uk+1

i−1,j + uki+1,j

)(3.7)

assumindo que os pontos da malha sao percorridos na ordem lexicografica, de modo que quando vamoscalcular o valor de u no ponto i, j na iteracao k + 1, nesta mesma iteracao ja calculamos os valores de u emi − 1, j e em i, j − 1, e usamos estes valores para calcular uk+1

i,j ao inves dos valores uki,j−1 e uk

i−1,j obtidosna iteracao anterior.

Em forma matricial, o algoritmo de Jacobi pode ser descrito da seguinte forma. Dada uma matriz A,existe uma unica decomposicao

A = D − L− U (3.8)

onde D e uma matriz diagonal, L e uma matriz estritamente triangular inferior e U e uma matriz estritamentetriangular superior; de fato, D = diag (a11, . . . , ann) e a parte diagonal de A, −L e a parte estritamentetriangular inferior de A e −U e a parte estritamente triangular superior de A. Entao o algoritmo de Jacobipode ser definido por

xk+1 = D−1(Lxk+1 + Uxk + b

)(3.9)

ou(D − L)xk+1 = Uxk + b,

dondexk+1 = (D − L)−1 (

Uxk + b). (3.10)

E importante ressaltar que existem matrizes para as quais o metodo de Jacobi converge e o metodo deGauss-Seidel diverge, e vice-versa. Veja a proxima secao sobre a convergencia dos metodos.

3.1.3 Metodo SOR

O processo de corrigir uma equacao atraves da modificacao de uma variavel e as vezes chamado de relax-amento. Antes da correcao, a equacao nao e verdadeira; como um conjunto de partes que nao se ajustam,ela esta em estado de tensao. A correcao de uma variavel relaxa a tensao. O metodo de Gauss-Seidel efetuarelaxamento sucessivo, ou seja, passa de equacao para equacao, relaxando uma depois da outra. [Watkins]

Por este motivo, os metodos de Jacobi e de Gauss-Seidel sao tambem chamados metodos de relaxamento.Em muitos casos, a convergencia pode ser substancialmente acelerada atraves de sobrerelaxamento. Issosignifica que ao inves de fazer uma correcao para a qual a equacao e satisfeita exatamente, nos fazemos umacorrecao maior. No caso mais simples, escolhe-se um fator de relaxamento ω > 1 que sobrecorrige por aquele


fator em cada passo (se mover um passo na direcao de xk para xk+1 e bom, mover naquela direcao ω > 1passos e melhor). Este e o chamado metodo de sobrerelaxamento sucessivo (SOR, successive overrelaxation):usando o metodo de Gauss-Seidel obtemos

xk+1i =

1aii

bi −

i−1∑

j=1

aijxk+1j +

n∑

j=i+1

aijxkj

;

daı tomamosxk+1

i = xki + ω

(xk+1

i − xki

).

Isso pode ser resumido em

xk+1i = xk

i + ω

1

aii

bi −

i−1∑

j=1

aijxk+1j −

n∑

j=i+1

aijxkj

− xk

i

. (3.11)

Quando ω = 1, o metodo SOR e exatamente o metodo de Gauss-Seidel. Um fator ω < 1 (subrelaxamento)normalmente diminui a velocidade de convergencia.

Para a maioria dos problemas, o melhor valor para o fator de relaxamento e desconhecido. Para a matrizde discretizacao obtida a partir da formula de cinco pontos, e sabido que o valor otimo de ω e, como veremosna proxima secao,

ω =2

1 + sen (π∆x). (3.12)

Em forma matricial, o metodo SOR pode ser descrito da seguinte forma. Como antes, dada uma matrizA escrevemos

A = D − L− U (3.13)

onde D e uma matriz diagonal, L e uma matriz estritamente triangular inferior e U e uma matriz estritamentetriangular superior. Entao, escrevendo o algoritmo SOR na forma

aiixk+1i = aiix

ki + ω

bi −

i−1∑

j=1

aijxk+1j −

n∑

j=i

aijxkj

,

temosDxk+1 = Dxk + ω

[Lxk+1 + (U −D) xk + b

](3.14)

ou (1ω

D − L

)xk+1 =

(1− ω

ωD + U

)xk + b,

donde

xk+1 =(

1ω

D − L

)−1 [(1− ω

ωD + U

)xk + b

]. (3.15)

3.1.4 Comparacao da Velocidade de Convergencia dos Tres Metodos

A tabela a seguir foi extraıda de [Watkins], pags. 533 e 542. Os metodos introduzidos acima foram usadospara resolver o sistema linear Ax = b onde A e a matriz de discretizacao obtida a partir da formula doscinco pontos do laplaciano no quadrado unitario Ω = (0, 1)2 e b e estabelecido pela condicao de fronteira deDirichlet dada por

g (x, y) =

0 se x = 0,y se x = 1,

(x− 1) sen x se y = 0,x (2− x) se y = 1,


ou seja, para resolver o problema discretizado −∆dud = 0 em Ωd,

ud = gd sobre ∂Ωd.

As iteracoes foram interrompidas quando∣∣uk+1 − uk

∣∣2

|uk+1|2< 10−8.

O numero de iteracoes necessarias para convergir de acordo com esta margem de tolerancia, para tres refina-mentos possıveis da malha (correspondentes a matrizes de dimensoes n = 81, 361 e 1521, respectivamente),de acordo com cada metodo e para diferentes valores de ω no caso do metodo SOR e apresentado na tabelaabaixo.

∆x = 0.1 ∆x = 0.05 ∆x = 0.025Jacobi 299 1090 3908SOR (ω = 0.8) 235 845 3018Gauss-Seidel 160 581 2082SOR (ω = 1.4) 67 262 955SOR (ω = 1.6) 42 151 577SOR (ω = 1.7) 57 96 412SOR (ω = 1.8) 86 89 252SOR (ω = 1.9) 176 180 179SOR (ω = 2.0) ∞ ∞ ∞

Vemos que o metodo de Gauss-Seidel e cerca de duas vezes mais rapido para convergir que o metodo deJacobi e que dependendo da escolha de ω, o metodo SOR pode ser ate dez vezes mais rapido que o metodode Gauss-Seidel para a malha mais refinada. Subrelaxamento nao ajuda e para ω = 2 o metodo SOR edivergente.

3.1.5 Metodo de Jacobi Amortecido

O metodo de Gauss-Seidel pode ser sobrerelaxado atraves de um parametro ω > 1 para obter um metodoque converge mais rapido.Ja o metodo de Jacobi nao pode em geral ser sobrerelaxado, porque o metodoobtido nao converge. Ele pode no entanto ser subrelaxado atraves de um parametro ω < 1 para obter ummetodo convergente, se bem que mais vagaroso. A vantagem de se utilizar um tal metodo e que para certosvalores de ω ele e um otimo suavizador de erro (em um sentido que sera explicado no proximo capıtulo),enquanto que o metodo de Jacobi usual nao possui esta propriedade. Assim, o metodo de Jacobi amortecidopode ser usado em metodos multigrid (veja o proximo capıtulo).

Pelo metodo de Jacobi usual obtemos

xk+1i =

1aii

bi −

n∑

j=1j 6=i

aijxkj

,

e tomamosxk+1

i = xki + ω

(xk+1

i − xki

),

ou seja,

xk+1i = xk

i + ω

1aii

bi −

n∑

j=1j 6=i

aijxkj

− xk

i

. (3.16)


Este metodo e conhecido como metodo de Jacobi amortecido, metodo de Jacobi ponderado ou aindametodo de relaxamento simultaneo (diferente do metodo de relaxamento sucessivo, baseado no metodo deGauss-Seidel, em que cada variavel e substituıda sucessivamente dentro da mesma iteracao a medida queela e atualizada; no metodo de Jacobi, as variaveis sao todas substituıdas simultameamente na proximaiteracao).

Em forma matricial, o metodo de Jacobi amortecido pode ser descrito da seguinte forma. Denotando porD a parte diagonal de A, temos

aiixk+1i = aiix

ki + ω

bi −

n∑

j=1

aijxkj

,

temosDxk+1 = Dxk + ω

[b−Axk

](3.17)

ou (1ω

D

)xk+1 =

(1ω

D −A

)xk + ωb,

donde

xk+1 =(

1ω

D

)−1 [(1ω

D −A

)xk + b

]. (3.18)

Em contraste com o metodo SOR, que converge em geral para 0 < ω < 2, o metodo de Jacobi amortecidoconverge para 0 < ω 6 1 (veja a proxima secao).

3.2 Analise de Convergencia dos Metodos Iterativos Lineares

Os metodos descritos na secao anterior sao casos especiais de uma classe geral de metodos chamados metodositerativos lineares ou metodos de correcao residual. Um metodo iterativo linear para resolver o sistemalinear

Ax = b

envolve a decomposicao da matriz A na forma

A = B − C, (3.19)

onde B e necessariamente uma matriz invertıvel, e entao a resolucao iterativa do sistema de equacoes

Bxk+1 = Cxk + b (3.20)

ou, mais explicitamente,xk+1 = B−1

(Cxk + b

).

Se xk → x, entao Bx = Cx + b, donde Ax = b. Do ponto de vista pratico, e importante que a matriz Bseja “facil de resolver” (mesmo que a inversa de B nao seja efetivamente calculada), como nos exemplos dasecao anterior:

B C

Jacobi D D −A

Gauss-Seidel D − L U

SOR1ω

D − L1− ω

ωD + U

Para obter uma convergencia rapida, tambem gostarıamos que B ≈ A e C ≈ 0. Deste ponto de vista, o idealseria B = A e C = 0 (convergencia em uma iteracao), mas isso viola em geral o criterio que B seja “facilde resolver”. Um compromisso e necessario: B deve aproximar A o melhor possıvel sem se tornar muitocomplicada.


3.2.1 Convergencia dos Metodos Iterativos Lineares

Para metodos iterativos em geral, definimos o erro algebrico por

ek = x− xk, (3.21)

enquanto que o erro residual e dado por

rk = Ax−Axk = f −Axk. (3.22)

O erro algebrico tem interesse puramente teorico (para provar que determinado metodo iterativo converge,precisamos mostrar que o erro algebrico tende a zero), ja que ele so pode ser calculado uma vez que seconhece a solucao exata, e se este for o caso obviamente nao ha necessidade de resolver o sistema. Ja o erroresidual pode ser usado como criterio de parada para o metodo iterativo. Como

Bek+1 = Bx−Bxk+1 = Ax + Cx− Cxk − b = C(x− xk

)= Cek,

segue queek+1 = B−1Cek.

Observe queB−1C = B−1 (B −A) = I −B−1A.

A matrizR = I −B−1A = B−1C (3.23)

e chamada a matriz de iteracao ou matriz de propagacao do erro do algoritmo considerado, porque

xk+1 = Rxk + B−1b. (3.24)

e o erro e dado porek+1 = Rek. (3.25)

Em particular,ek = Rke0 (3.26)

de modo que o erro converge para 0, independentemente do chute inicial x0, se e somente se Rk → 0. Issoocorre se e somente se existe alguma norma matricial ‖·‖ tal que ‖R‖ < 1. Obter uma norma matricialque satisfaz esta propriedade, no entanto, e difıcil. Vamos obter uma condicao necessaria e suficiente paraRk → 0 em termos do raio espectral da matriz de iteracao (Corolario 3.5 a seguir), que e em geral um poucomais facil de calcular. Antes, para motivar o resultado, suponha que A seja uma matriz diagonalizavel comλ1, . . . , λn os seus autovalores e v1, . . . , vn uma correspondente base de autovetores. Escrevendo o erroinicial como uma combinacao linear dos autovetores, temos

e0 =n∑

i=1

aivi.

Logo,

ek = Rke0 =n∑

i=1

aiλki vi,

de modo que∣∣ek

∣∣ 6n∑

i=1

|ai| |λi|k |vi| .

Como |λi|k → 0 se e somente se |λi| < 1, concluımos que ek → 0 qualquer que seja o erro inicial (isto e,qualquer que seja o chute inicial), se e somente se ρ (R) = max16i6n |λi| < 1 .


3.1 Lema. Se A ∈ Mn (C) e ‖·‖ e qualquer norma matricial, entao

ρ (A) 6 ‖A‖ .

Prova. Seja λ um autovalor qualquer de A e x um autovetor nao-nulo correspondente a λ, de modo que

Ax = λx.

Considere a matriz X ∈ Mn (C) cujas colunas sao todas iguais ao vetor x. Temos tambem

AX = λX

de modo que|λ| ‖X‖ = ‖AX‖ 6 ‖A‖ ‖X‖ ,

donde|λ| 6 ‖A‖

para todo autovalor λ de A. Como existe um autovalor λ de A tal que ρ (A) = λ, isso prova o resultado. ¥

3.2 Lema. Seja A ∈ Mn (C) e ε > 0 dado. Entao existe uma norma matricial ‖·‖ tal que

ρ (A) 6 ‖A‖ 6 ρ (A) + ε. (3.27)

Prova. Toda matriz complexa e triangularizavel atraves de uma matriz unitaria (isto e, uma matriz U quesatisfaz U∗U = UU∗ = I; sua inversa e a sua adjunta ou transposta conjugada). Sejam entao

T =

λ1 a12 a22 . . . a1n

λ2 a23 . . . a2n

λ3 . . . a3n

. . ....

λn

uma matriz triangular e U uma matriz unitaria tais que

A = U∗TU.

Considere a matriz diagonal

Dt =

tt2

. . .tn

.

Temos

DtTD−1t =

λ1 a12t−1 a22t

−2 . . . . . . a1nt−n+1

λ2 a23t−1 . . . . . . a2nt−n+2

λ3 . . . . . . a3nt−n+3

. . ....

λn−1 an−1,nt−1

λn

.

Logo, para t > 0 suficientemente grande, a matriz DtTD−1t tem a propriedade que a soma dos valores

absolutos de elementos fora da diagonal principal e menor que ε. Em particular, se ‖·‖L denota a norma domaximo das somas das linhas, podemos garantir que

∥∥DtTD−1t

∥∥L

6 ρ (A) + ε


para t suficientemente grande. Portanto, fixado um tal t, se definirmos uma norma por

‖A‖ :=∥∥DtUAU∗D−1

t

∥∥L

=∥∥∥(U∗D−1

t

)−1AU∗D−1

t

∥∥∥L

,

teremos‖A‖ =

∥∥DtUAU∗D−1t

∥∥L

=∥∥DtTD−1

t

∥∥L

6 ρ (A) + ε.

Pelo lema anterior, ρ (A) 6 ‖A‖. ¥

3.3 Lema. Seja A ∈ Mn (C). Se existe alguma norma matricial ‖·‖ tal que ‖A‖ < 1, entao

Ak → 0.

Prova. Se ‖A‖ < 1, entao ∥∥Ak∥∥ 6 ‖A‖k → 0.

¥

3.4 Proposicao. Seja A ∈ Mn (C). EntaoAk → 0

se e somente seρ (A) < 1.

Prova. Se existe algum autovalor λ de A tal que |λ| > 1 e x e um autovetor nao-nulo correspondente, entao

Akx = λkx

nao converge para 0. Reciprocamente, se ρ (A) < 1, entao pelo Lema 3.2 existe uma norma matricial ‖·‖ talque ‖A‖ < 1, logo Ak → 0 pelo lema anterior. ¥

3.5 Corolario. Seja R a matriz de iteracao de um metodo iterativo linear. Entao

ek → 0

se e somente seρ (R) < 1.

Em outras palavras, um metodo iterativo linear e convergente independentemente da escolha do chuteinicial se e somente se todos os autovalores da matriz de iteracao tem valor absoluto menor que 1.

3.2.2 Velocidade de Convergencia dos Metodos Iterativos Lineares

O raio espectral tambem da informacao sobre a velocidade de convergencia. Se nos tivermos dois metodositerativos lineares diferentes, isto e, duas maneiras diferentes de decompor a matriz A:

A = B1 − C1 = B2 − C2,

entao o segundo metodo convergira mais rapido se e somente se

ρ (R2) < ρ (R1) .

Vamos analisar a velocidade de convergencia dos metodos iterativos com maior precisao. Novamente atıtulo de motivacao, suponha que R e uma matriz diagonalizavel com seu maior autovalor sendo um autovalorsimples. Ordene os autovalores de R na forma

|λ1| > |λ2| > . . . > |λn|


e seja v1, . . . , vn uma correspondente base de autovetores. Escrevendo de novo

e0 =n∑

i=1

aivi,

donde

ek = Rke0 =n∑

i=1

aiλki vi,

segue que

ek = λk1

[a1x1 +

n∑

i=2

ai

(λi

λ1

)k

vi

].

Como (λi

λ1

)k

→ 0,

a taxa de convergencia e determinada por |λ1|k. Para k grande, temos

ek ≈ λk1a1v1.

Portanto, ∣∣ek+1∣∣

|ek| = |λ1| = ρ (R) . (3.28)

Em outras palavras, a convergencia e linear com taxa de convergencia igual ao raio espectral. Se a1 =0 a convergencia sera mais rapida, pois dependera do modulo do segundo autovalor, mas e obviamenteextremamente raro que o chute inicial satisfaca esta condicao. Para o caso geral, precisamos do seguinteresultado:

3.6 Proposicao. Seja A ∈ Mn (C) e ‖·‖ uma norma matricial. Entao

ρ (A) = lim∥∥Ak

∥∥1/k.

Prova. Como os autovalores da matriz Ak sao as k-esimas potencias dos autovalores de A, temos que

ρ (A)k = ρ(Ak

)6

∥∥Ak∥∥ ,

dondeρ (A) 6

∥∥Ak∥∥1/k

.

Dado ε > 0, a matriz

B =1

ρ (A) + εA

tem raio espectral menor que 1, logo Bk → 0. Portanto, existe algum N = N (ε,A) tal que∥∥Bk

∥∥ < 1

ou seja, ∥∥Ak∥∥1/k

< ρ (A) + ε

para todo k > N . ¥Definimos a taxa media de convergencia de um metodo iterativo linear com matriz de iteracao R por

Rk (R) = − log10

∥∥Rk∥∥1/k

= −1k

log10

∥∥Rk∥∥ (3.29)

e a taxa assintotica de convergencia por

R∞ (R) = limk→∞

Rk (R) . (3.30)


3.7 Corolario. Seja R a matriz de iteracao de um metodo iterativo linear. Entao a taxa assintotica deconvergencia do metodo e dada por

R∞ (R) = − log10 ρ (R) . (3.31)

Prova. PoisR∞ (R) = − lim

k→∞log10

∥∥Rk∥∥1/k

= − log10 limk→∞

∥∥Rk∥∥1/k

= − log10 ρ (R) .

¥A taxa assintotica de convergencia mede o aumento no numero de casas decimais corretas na solucao poriteracao. De fato, usando a norma matricial do Lema 3.2 e medindo as normas dos vetores de acordo, temos

∣∣ek+1∣∣

|ek| =

∣∣Rk+1e0∣∣

|Rke0| 6 ‖R‖ = ρ (R) + ε,

donde

− log10

∣∣ek+1∣∣

|ek| = − log10 ρ (R) + O (ε)

oulog10

∣∣ek∣∣− log10

∣∣ek+1∣∣ = R∞ (R) + O (ε) . (3.32)

Assim, se∣∣ek

∣∣ = O(10−p

),∣∣ek+1

∣∣ = O(10−q

),

teremosq − p ≈ R∞ (R) ,

isto e, reduzimos R∞ (R) ≈ q − p casas decimais no erro. Visto de outra forma, como∣∣ek+m

∣∣|ek| =

∣∣Rk+me0∣∣

|Rke0| 6 ‖Rm‖ = ρ (R)m + O (ε) ,

donde

− log10

∣∣ek+m∣∣

|ek| ≈ −m log10 ρ (R) ,

ou

m =log10

(∣∣ek+m∣∣ /

∣∣ek∣∣)

log10 ρ (R)(3.33)

e o numero de iteracoes necessarias para diminuir o erro de um numero prescrito de casas decimais.

3.2.3 Convergencia para Matrizes Simetricas Positivas Definidas

Para matrizes reais simetricas positivas definidas e mais facil provar a convergencia dos metodos iterativoslineares. Temos o seguinte resultado basico a seguir. Antes precisamos da seguinte definicao:

Definicao. Introduzimos uma ordenacao parcial em Mn (C) definindo

A 6 B

se〈Ax, x〉 6 〈Bx, x〉

para todo x ∈ Cn.


Em particular, se A e uma matriz positiva definida, segue que A > εI para algum ε (o menor autovalor deA) e denotamos este fato por

A > 0.

3.8 Teorema. Seja A uma matriz simetrica positiva definida e seja A = B − C com B invertıvel. Entaoo metodo iterativo linear com matriz de iteracao R = B−1C converge se e somente se Bt + C e umamatriz simetrica positiva definida.

Prova. Medimos a norma do erro atraves da norma induzida por A

|x|A := 〈Ax, x〉1/2

e consideraremos a norma matricial ‖·‖A induzida por esta norma. Se provarmos que

‖R‖A < 1,

o metodo convergira. Temos

‖R‖2A =∥∥B−1C

∥∥2

A= sup

x6=0

∣∣B−1Cx∣∣2A

|x|2A= sup

x 6=0

⟨AB−1Cx,B−1Cx

⟩

〈Ax, x〉 = supx 6=0

⟨CtB−tAB−1Cx, x

⟩

〈Ax, x〉 . (3.34)

Suponha que Bt + C e uma matriz simetrica, positiva definida. Temos

CtB−tAB−1C =(Bt −A

)B−tAB−1 (B −A) =

(I −AB−t

)A

(I −B−1A

)

= A− (AB−tA + AB−1A−AB−tAB−1A

)

= A−AB−t(B + Bt −A

)B−1A

= A− (B−1A

)t (B + Bt −A

)B−1A

ouCtB−tAB−1C = A− (

B−1A)t (

Bt + C)B−1A, (3.35)

de modo que CtB−tAB−1C e uma matriz simetrica, positiva definida. Logo, por (4.8), mostrar que ‖R‖A < 1e equivalente a provar que

CtB−tAB−1C < A,

e por (4.16) CtB−tAB−1C < A se e somente se(B−1A

)t (Bt + C) B−1A > 0, o que e verdade porque Bt+Ce positiva definida. ¥

3.3 Convergencia dos Metodos Iterativos Lineares para as Ma-trizes de Discretizacao

3.3.1 Convergencia do Metodo de Jacobi

3.9 Teorema. Se A e uma matriz irredutıvel, diagonalmente dominante tal que |aii| >n∑

j=1j 6=i

|aij | para pelo

menos alguma linha i, entao o metodo de Jacobi converge.

Prova. Seja D a parte diagonal da matriz A e R = D−1 (D −A) = I − D−1A a matriz de iteracao dometodo de Jacobi para A. Suponha por absurdo que exista um autovalor λ de R tal que |λ| > 1. Comoλ det

(λ−1R− I

)= det (R− λI) = 0, temos

det(I − λ−1R

)= 0.


Por outro lado, observe que I − λ−1R tambem e irredutıvel, pois

Rij =(I −D−1A

)ij

=

0 se i = j,

−aij

aiise i 6= j,

(I − λ−1R

)ij

=

1 se i = j,

λ−1 aij

aiise i 6= j,

de modo que, onde A se anula, I−λ−1R tambem se anula. Alem disso, I−λ−1R e diagonalmente dominantee estritamente dominante nas linhas onde A e, pois |λ|−1 6 1,

(I − λ−1R

)ii

= 1 e

n∑

j=1j 6=i

∣∣∣(I − λ−1R

)ij

∣∣∣ =|λ|−1

|aii|n∑

j=1j 6=i

|aij | 6 1|aii|

n∑

j=1j 6=i

|aij | .

Mas, pela Proposicao 2.16, isso implica que I − λ−1R e invertıvel, uma contradicao. ¥O Teorema 3.8 mostra que o metodo de Jacobi converge para as matrizes de discretizacao obtidas atravesdos esquemas de diferencas finitas do Capıtulo 2.

Atraves do Teorema 3.9, fomos capazes de provar a convergencia do metodo de Jacobi para as matrizes dediscretizacao sem calcular explicitamente os seus raios espectrais. Para analizar a velocidade de convergenciado metodo de Jacobi, no entanto, e necessario obter os raios espectrais destas matrizes. Vamos fazer issopara as matrizes de discretizacao obtidas a partir da formula de tres pontos unidimensional e a partir daformula de cinco pontos bidimensional.

3.10 Teorema. Seja A a matriz de discretizacao obtida a partir da formula de tres pontos unidimensionalou a partir da formula de cinco pontos bidimensional com ∆x = ∆y. Seja R = D−1 (D −A) a matrizde iteracao do metodo de Jacobi. Entao

ρ (R) = cosπ

n. (3.36)

Prova. Para o metodo de Jacobi, a matriz de discretizacao xk+1 = Rxk +D−1b e obtida atraves da formula:

uk+1i,j =

14

(uk


i−1,j + uki+1,j

).

Ja vimos no Lema 1.2 que

−ukli−1,j − ukl

i+1,j + 4ukli,j − ukl

i,j−1 − ukli,j+1 =

(λkl∆x2

)ukl

i,j

com

λkl =2

∆x2

(2− cos

kπ

n− cos

lπ

n

).

Daı segue queukl

i,j−1 + ukli,j+1 + ukl

i−1,j + ukli+1,j =

(4− λkl∆x2

)ukl

i,j

Logo14

(ukl

i,j−1 + ukli,j+1 + ukl

i−1,j + ukli+1,j

)= µlkukl

i,j

para

µlk = 1− 14λkl∆x2 = 1− 1

2

(2− cos

kπ

n− cos

lπ

n

)=

12

(cos

kπ

n+ cos

lπ

n

).


Estes sao os autovalores da matriz de iteracao de Jacobi para a matriz de discretizacao obtida a partir daformula de cinco pontos (observe que elas possuem os mesmos autovetores; no entanto R possui autovaloresnulos). Segue que o maximo autovalor ocorre quando k = l = 1, logo

ρ (R) = cosπ

n.

O argumento para a formula de tres pontos e analogo. ¥Para o quadrado unitario temos

ρ (R) = cos (π∆x) . (3.37)

Vemos em particular que ρ (R) → 1 quando ∆x → 0, de modo que a velocidade de convergencia do metodode Jacobi vai ficando cada vez menor para malhas mais refinadas. Podemos dizer mais usando a expansaoda funcao cosseno em torno da origem

cosx = 1− 12x2 + O

(x4

);

se ∆x e pequeno podemos aproximar

cos (π∆x) ≈ 1− π2

2∆x2,

de modo que ρ (R) → 1 quadraticamente quando ∆x → 0. Em outras palavras, para uma malha duas vezesmais refinada (isto e, ∆x reduzido pela metade), o metodo de Jacobi e cerca de quatro vezes mais vagarosoem media (consulte novamente a tabela no final da secao anterior). A tabela abaixo mostra os valores doraio espectral para alguns valores de ∆x:

∆x 0.1 0.05 0.025

ρ (R) 0.9511 0.9877 0.9969

Para ∆x = 0.025 (correspondente a uma matriz de tamanho n = 39× 39 = 1521), temos

R∞ (R) = − log10 (0.9969) = 0.0013484,

de modo que para reduzir o erro pelo fator de uma casa decimal precisamos de

m =log10 0.1

log10 ρ (R)= − 1

log10 ρ (R)=

10.00135

≈ 742

iteracoes.

3.3.2 Convergencia do Metodo de Gauss-Seidel


j=1j 6=i

|aij | para pelo

menos alguma linha i, entao o metodo de Gauss-Seidel converge.

Prova. Sejam D a parte diagonal, −L a parte triangular inferior estrita e −U a parte triangular superiorestrita da matriz A, e seja R = (D − L)−1

U a matriz de iteracao do metodo de Gauss-Seidel para A.Escrevemos

R = (D − L)−1U =

[D

(I −D−1L

)]−1U

ouR =

(I −D−1L

)−1D−1U. (3.38)


Suponha por absurdo que exista um autovalor λ de R tal que |λ| > 1; como na demonstracao do Teorema3.9, temos

det(I − λ−1R

)= det

(I − λ−1

[(I −D−1L

)−1D−1U

])= 0.

Agora, observando quedet

(I −D−1L

)= 1

porque I −D−1L e uma matriz triangular inferior com apenas 1’s na diagonal principal, escrevemos

0 = det(I − λ−1

[(I −D−1L

)−1D−1U

])

= det(I −D−1L

)det

(I − λ−1

[(I −D−1L

)−1D−1U

])

= det(

I −D−1L) (

I − λ−1[(

I −D−1L)−1

D−1U])

= det(I −D−1L− λ−1D−1U

).

Por outro lado,D−1A = I −D−1L−D−1U

e irredutıvel, diagonalmente dominante e estritamente dominante nas linhas onde A e porque

(D−1A

)ij

=

1 se i = j,aij

aiise i 6= j.

Logo, a matriz I − D−1L − λ−1D−1U tambem satisfaz estas propriedades, pois I, −D−1L e −D−1U saorespectivamente a parte diagonal, a parte triangular inferior estrita e a parte triangular superior estrita damatriz D−1A, e multiplicar a parte triangular inferior estrita pelo numero λ−1 cujo modulo e menor que ouigual a 1 nao alterara a dominancia diagonal (na verdade so tende a melhora-la) nem acrescentara zeros amatriz. A Proposicao 2.16 implica entao que I −D−1L− λ−1D−1U e invertıvel, um absurdo. ¥Usando o Teorema 3.11, concluımos que o metodo de Gauss-Seidel converge para as matrizes de discretizacaoobtidas atraves dos esquemas de diferencas finitas do Capıtulo 1. Para analizar a velocidade de convergenciado metodo de Gauss-Seidel, vamos obter os raios espectrais para as matrizes de discretizacao obtidas a partirda formula de tres pontos unidimensional e a partir da formula de cinco pontos bidimensional.

3.12 Teorema. Seja A a matriz de discretizacao obtida a partir da formula de tres pontos unidimensionalou a partir da formula de cinco pontos bidimensional com ∆x = ∆y. Seja R = (D − L)−1

U a matrizde iteracao do metodo de Gauss-Seidel. Entao

ρ (R) = cos2π

n. (3.39)

Prova. Para obter o raio espectral da matriz de iteracao R, queremos encontrar os autovalores µ de R:

Ru = (D − L)−1Uu = µu,

ou seja,Uu = µ (D − L)u

(um problema de autovalor generalizado). No caso da matriz de discretizacao da formula de cinco pontos,isso significa encontrar µ tal que

ui,j+1 + ui+1,j = µ (4ui,j − ui,j−1 − ui−1,j) . (3.40)

Para os autovalores nao-nulos, podemos fazer a substituicao

ui,j = µi+j2 vi,j (3.41)


para transformar a equacao de autovalor naquela que aparece no metodo de Jacobi. Temos

µi+j+1

2 vi,j + µi+j+1

2 vi+1,j = µ(4µ

i+j2 vi,j − µ

i+j−12 vi,j−1 − µ

i+j−12 vi−1,j

)

= 4µi+j+2

2 vi,j − µi+j+1

2 vi,j−1 − µi+j+1

2 vi−1,j ,

de modo que, dividindo por µi+j+1

2 , obtemos

vi−1,j + vi+1,j + vi,j−1 + vi,j+1 = µ1/24vi,j .

Portanto os autovalores da matriz de iteracao de Gauss-Seidel para esta matriz sao exatamente os quadradosdos autovalores da matriz de iteracao de Jacobi (e os autovetores sao os mesmos):

µlk =14

(cos

kπ

n+ cos

lπ

n

)2

.

Portanto, o maximo autovalor ocorre quando k = l = 1 e

ρ (R) = cos2π

n.

O argumento para a formula de tres pontos e analogo. ¥Para o quadrado unitario temos

ρ (R) = cos2 (π∆x) ,

e usando

cos2 x =[1− 1

2x2 + O

(x4

)]2

= 1− x2 + O(x4

),


cos2 (π∆x) ≈ 1− π2∆x2.

No metodo de Gauss-Seidel ainda temos ρ (R) → 1 quadraticamente quando ∆x → 0, mas a sua velocidadede convergencia para a matriz de discretizacao de cinco pontos do quadrado unitario e duas vezes maior quea do metodo de Jacobi. Para ver isso, faca a expansao do logaritmo em torno do ponto x = 1:

log (1 + x) = x + O(∆x2

).

Segue que

R∞ (RJacobi) =π2

2∆x2 + O

(∆x4

), (3.42)

R∞ (RGauss-Seidel) = π2∆x2 + O(∆x4

). (3.43)

3.3.3 Convergencia do Metodo SOR

3.13 Teorema. Se o metodo SOR converge, entao

0 < ω < 2.

Prova. A matriz de iteracao do metodo SOR e

R =(

1ω

D − L

)−1 (1− ω

ωD + U

)=

[1ω

D(I − ωD−1L

)]−1 (1− ω

ωD + U

)

=(I − ωD−1L

)−1ωD−1

(1− ω

ωD + U

)


ouR =

(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

]. (3.44)

Se λ1, . . . , λn sao os autovalores de R, entao

detR = λ1 . . . λn.

Mas,

detR = det(

I − ωD−1L)−1 [

(1− ω) I + ωD−1U]

= det(I − ωD−1L

)−1det

[(1− ω) I + ωD−1U

]

= (1− ω)n,

ja que I −ωD−1L e uma matriz triangular inferior com apenas 1 na diagonal principal e (1− ω) I +ωD−1Ue uma matriz triangular superior com apenas 1− ω na diagonal principal. Logo

λ1 . . . λn = (1− ω)n.

Em particular, pelo menos um dos autovalores λj de R deve satisfazer

|λj | > |1− ω| .

Mas, se o metodo SOR converge, devemos ter tambem |λ| < 1 para todo autovalor λ de R. Logo

|1− ω| < 1,

donde0 < ω < 2.

¥

3.14 Corolario. Se R e a matriz de iteracao n× n para o metodo SOR, entao

detR = (1− ω)n.

Em particular, diferente das matrizes de iteracao dos metodos de Jacobi e de Gauss-Seidel (para a matriz dediscretizacao de cinco pontos), zero nao e um autovalor para a matriz de iteracao do metodo SOR se ω 6= 1(para nenhuma matriz).


j=1j 6=i

|aij | para pelo

menos alguma linha i, entao o metodo SOR converge se 0 < ω 6 1.

Prova. A demonstracao e analoga a do Teorema 3.11. A matriz de iteracao do metodo SOR e

R =(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

].

Suponha por absurdo que exista um autovalor λ de R tal que |λ| > 1; temos

det(I − λ−1R

)= det

(I − λ−1

(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

])= 0.

Agora, observando quedet

(I − ωD−1L

)= 1


porque I − ωD−1L e uma matriz triangular inferior com apenas 1’s na diagonal principal, escrevemos

0 = det(I − λ−1

(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

])

= det(I − ωD−1L

)det

(I − λ−1

(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

])

= det[(

I − ωD−1L) (

I − λ−1(

I − ωD−1L)−1 [

(1− ω) I + ωD−1U])]

= det(I − ωD−1L− λ−1

[(1− ω) I + ωD−1U

])

= det[

1− λ−1 (1− ω)]I − ωD−1L− λ−1ωD−1U

.

Por outro lado, como vimos na demonstracao do Teorema 3.11, a matriz

D−1A = I −D−1L−D−1U

e irredutıvel, diagonalmente dominante e estritamente dominante nas linhas onde A e, logo a matriz

S =[1− λ−1 (1− ω)

]I − ωD−1L− λ−1ωD−1U

tambem satisfaz estas propriedades. De fato, S tem zeros nas mesmas posicoes que I −D−1L−D−1U , logoa sua irredutibilidade nao e afetada. Alem disso, pela dominancia diagonal de D−1A, sabemos que se

bij =(D−1L

)ij

,

cij =(D−1U

)ij

.

entao

1 >i−1∑

j=1

|bij |+n∑

j=i+1

|cij | .

Para provar a dominancia diagonal de S, observamos que os valores que S possui na diagonal principal sao

1− λ−1 (1− ω) = 1− 1− ω

λ=

λ + ω − 1λ

,

de modo que precisamos provar que

∣∣∣∣λ + ω − 1

λ

∣∣∣∣ > ω

i−1∑

j=1

|bij |+ ω

|λ|n∑

j=i+1

|cij |

se 0 < ω 6 1 e |λ| > 1. Provaremos que∣∣∣∣λ + ω − 1

λ

∣∣∣∣ > ω,

∣∣∣∣λ + ω − 1

λ

∣∣∣∣ > ω

|λ| .

Para isso, observe que como |λ| > 1 basta provar a primeira desigualdade, a qual por sua vez e equivalente a

|λ + ω − 1| > |λ|ω.

E facil ver que esta desigualdade e valida quando λ ∈ R, pois

|λ + ω − 1| = λ + ω − 1 > λω porque λ− 1 > λω − ω = ω (λ− 1) .


Para o caso geral em que λ ∈ C, fazemos cair no caso real escrevendo

|λ + ω − 1|2 = |λ− (1− ω)|2 = |λ|2 − 2 (Re λ) (1− ω) + (1− ω)2

> |λ|2 − 2 |λ| (1− ω) + (1− ω)2 = [|λ| − (1− ω)]2

= [|λ|+ ω − 1]2 > |λ|2 ω2.

O resultado acima continua valendo com desigualdade estrita nas linhas onde a desigualdade e estrita. AProposicao 2.16 implica entao que S e invertıvel, contradizendo det S = 0. ¥

3.16 Teorema. Seja A uma matriz simetrica positiva definida. Entao o metodo SOR converge se 0 < ω < 2.

Prova. Usaremos o Teorema 3.8. Escrevendo A = D − L − U , temos Lt = U porque A e simetrica e asentradas diagonais de D positivas porque A e positiva definida. Para o metodo SOR temos

B =1ω

D − L e C =1− ω

ωD + U,

logo

Bt + C =1ω

D − Lt +1− ω

ωD + U =

2− ω

ωD

e uma matriz simetrica positiva definida se 0 < ω < 2. ¥Na verdade, se as entradas diagonais de uma matriz simetrica sao positivas, a condicao de ser definidapositiva e equivalente a convergencia do metodo SOR para 0 < ω < 2, como o proximo resultado mostra.

3.17 Teorema. Seja A uma matriz simetrica com entradas diagonais positivas. Entao o metodo SORconverge se e somente se A e positiva definida e 0 < ω < 2.

Prova. Assuma que A e positiva definida e que 0 < ω < 2. Seja

R =(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

]

a matriz de iteracao do metodo SOR. Se λ e um autovalor de R e x um autovetor associado, temos Rx = λx,donde [

(1− ω) I + ωD−1U]x = λ

(I − ωD−1L

)x.

Fazendo o produto interno canonico (hermitiano) de Cn de ambos os lados com o vetor x, segue que

(1− ω) 〈x, x〉+ ω⟨x,D−1Ux

⟩= λ

(〈x, x〉 − ω⟨x,D−1Lx

⟩)

Isolando λ,

λ =(1− ω) 〈x, x〉+ ω

⟨x,D−1Ux

⟩

〈x, x〉 − ω 〈x,D−1Lx〉 . (3.45)

Como A e simetrica, o produto de matrizes simetricas D−1A = I − D−1U − D−1L tambem e; comoD−1U,D−1L sao respectivamente a parte estritamente triangular superior e estritamente triangular infe-rior de uma matriz simetrica, temos (

D−1U)t

= D−1L.

Logo ⟨x, D−1Ux

⟩=

⟨(D−1U

)tx, x

⟩=

⟨(D−1L

)x, x

⟩= 〈x, (D−1L) x〉,

e definindo

z =

⟨x,

(D−1L

)x⟩

〈x, x〉 ,


podemos escrever

λ =(1− ω) + ωz

1− ωz. (3.46)

Os argumentos acima assumem que o denominador e nao-nulo. E, de fato, temos

Re z =12

(z + z) =12

(⟨x,

(D−1L

)x⟩

〈x, x〉 +

⟨x,

(D−1U

)x⟩

〈x, x〉

)=

12

⟨x,

(D−1L + D−1U

)x⟩

〈x, x〉

=12

⟨x,

(I −D−1A

)x⟩

〈x, x〉 =12

(1−

⟨x,

(D−1A

)x⟩

〈x, x〉

).

e como A e positiva definida, D−1A tambem e, o que implica⟨x,

(D−1A

)x⟩

〈x, x〉 > 0

dondeRe z <

12.

de modo que a parte real do denominador 1− ωz de λ e nao-nula para 0 < ω < 2. Segue que

|λ|2 = λλ =[(1− ω) + ωz] [(1− ω) + ωz]

(1− ωz) (1− ωz)=

(1− ω)2 + 2ω (1− ω)Re z + ω2 |z|21− 2ω Re z + ω2 |z|2

=ω2 − 2ω2 Re z − 2ω + 4ω Re z + 1− 2ω Re z + ω2 |z|2

1− 2ω Re z + ω2 |z|2

= 1− ω (2− ω) (1− 2 Re z)1− 2ω Re z + ω2 |z|2 .

Como 0 < ω < 2 e Re z <12, temos

ω (2− ω) (1− 2Re z) > 0,

e concluımos que|λ| < 1

para todo autovalor λ de R, logo o metodo SOR converge. A demonstracao da recıproca (assim como umademonstracao alternativa, variacional, deste teorema) pode ser vista em [Young]. ¥Usando o Teorema 3.15, concluımos que o metodo SOR converge para as matrizes de discretizacao obtidasatraves dos esquemas de diferencas finitas do Capıtulo 1 se 0 < ω 6 1. Isso permite apenas subrelaxamentodo metodo de Gauss-Seidel, o que em geral reduz a velocidade de convergencia. Por outro lado, usando oTeorema 3.16 ou o Teorema 3.17, concluımos que o metodo SOR converge para as matrizes de discretizacaoobtidas a partir da formula de tres pontos unidimensional e a partir da formula de cinco pontos bidimensionalse 0 < ω < 2, ja que estas sao matrizes simetricas, positivas definidas (ja as matrizes de discretizacao obtidasatraves de coordenadas polares ou pelo esquema de Shortley-Weller nao sao simetricas, em geral, comovimos).

Em seguida fazemos uma analise da velocidade de convergencia do metodo SOR para a matriz de dis-cretizacao da formula de cinco pontos, bem como obtemos o melhor valor do fator de relaxamento ω paraeste caso.

3.18 Lema. Seja A a matriz de discretizacao obtida a partir da formula de tres pontos unidimensional oua partir da formula de cinco pontos bidimensional com ∆x = ∆y. Se λ 6= 0 e um autovalor de RSOR,entao existe um autovalor λJ de RJ tal que

λJ =1− ω − λ

λ1/2ω2. (3.47)


Reciprocamente, se λJ e um autovalor de RJ e λ ∈ C satisfaz a equacao acima, entao λ e um autovalorde RSOR.

Prova. Argumentamos como na demonstracao do Teorema 3.12. Para obter o raio espectral da matriz deiteracao RSOR, queremos encontrar os autovalores λ de RSOR:

RSORu =(I − ωD−1L

)−1 [(1− ω) I + ωD−1U

]u = λu,

ou seja, [(1− ω) I + ωD−1U

]u = λ

(I − ωD−1L

)u

No caso da matriz de discretizacao da formula de cinco pontos, isso significa encontrar λ tal que

(1− ω) ui,j +ω

4ui,j+1 +

ω

4ui+1,j = λ

(ui,j − ω

4ui,j−1 − ω

4ui−1,j

)

ou1− ω − λ

ωui,j =

14

(ui,j+1 + ui+1,j + λui,j−1 + λui−1,j) . (3.48)

Fazendo a substituicaoui,j = λ

i+j2 vi,j

e dividindo por µi+j+1

2 , segue que

vi−1,j + vi+1,j + vi,j−1 + vi,j+1 =1− ω − λ

λ1/2ω4vi,j

e daı o resultado. ¥Resolvendo a equacao (3.47) como uma equacao quadratica em

√λ, vemos que as duas raızes λ± =

(√λ±

)2

podem ser escritas na forma

λ± =14

[−ωλJ ±

√ω2λ2

J − 4 (ω − 1)]2

. (3.49)

DenotaremosΛω,λJ = max (|λ+| , |λ−|) (3.50)

e por λJ = ρ (RJ ) o maior autovalor do metodo de Jacobi.

3.19 Proposicao. Seja A a matriz de discretizacao obtida a partir da formula de tres pontos unidimensionalou a partir da formula de cinco pontos bidimensional com ∆x = ∆y. Entao

ρ (RSOR,ω) = Λω,λJ(3.51)

Prova. Por definicao,ρ (RSOR,ω) = max

λJ

Λω,λJ.

De (3.49) segue que

Λω,λJ =14

∣∣∣∣ωλJ +√

ω2λ2

J − 4 (ω − 1)∣∣∣∣2

.

Se 0 < ω 6 1, ω2λ2

J − 4 (ω − 1) > 0 e Λω,λJe uma funcao crescente de λJ , logo o maximo e atingido em λJ .

Se ω > 1, defina

λc =

√4 (ω − 1)

ω2.


Se λJ > λc, ω2λ2

J − 4 (ω − 1) > 0 e segue a conclusao como no caso anterior. Se λJ 6 λc, entao ω2λ2

J −4 (ω − 1) 6 0 e √

ω2λ2

J − 4 (ω − 1) =√

4 (ω − 1)− ω2λ2

J i,

onde i =√−1, logo

Λω,λJ=

∣∣∣∣ωλJ +√

ω2λ2

J − 4 (ω − 1)∣∣∣∣2

=

∣∣∣∣∣

√ω2λ2

J +[4 (ω − 1)− ω2λ

2

J

]∣∣∣∣∣

2

= ω − 1,

e novamente Λω,λJe uma funcao crescente de λJ . ¥

Definaωotimo =

2

1 +√

1− λ2

J

. (3.52)

Note que 1 < ωotimo < 2. Mostraremos que ωotimo e de fato o melhor valor para o fator de relaxamento nometodo SOR. Antes precisamos do seguinte resultado:

3.20 Proposicao. Seja A a matriz de discretizacao obtida a partir da formula de tres pontos unidimensionalou a partir da formula de cinco pontos bidimensional com ∆x = ∆y. Entao

ρ (RSOR,ω) =

14

(ωλJ +

√ω2λ

2

J − 4 (ω − 1))2

se 0 < ω 6 ωotimo,

ω − 1 se ωotimo 6 ω < 2.

(3.53)

Prova. Temos ω2λ2

J − 4 (ω − 1) > 0 para 0 < ω < 2 se e somente se ω 6 ωotimo. De fato, as raızes def (ω) = ω2λ

2

J − 4ω + 4 sao

ω± =4± 4

√1− λ

2

J

2λ2

J

=2

λ2

J

(1±

√1− λ

2

J

)

de modo que a raiz positiva de f e maior que 2, logo para que f (ω) > 0 se 0 < ω < 2, devemos ter

ω 6 2

λ2

J

(1−

√1− λ

2

J

)=

2

λ2

J

1−(1− λ

2

J

)

1 +√

1− λ2

J

=2

1 +√

1− λ2

J

.

O resultado segue entao como na demonstracao da proposicao anterior. ¥

3.21 Teorema. Seja A a matriz de discretizacao obtida a partir da formula de tres pontos unidimensionalou a partir da formula de cinco pontos bidimensional com ∆x = ∆y. Entao o fator de relaxamentootimo para o metodo SOR e dado por

ωotimo =2

1 + senπ

n

(3.54)

e o fator de relaxamento otimo para o metodo SOR.

Prova. Se 0 < ω 6 ωotimo, entao ω2λ2

J − 4 (ω − 1) > 0 e

d

dω

(ωλJ +

√ω2λ

2

J − 4 (ω − 1))

=λJ

√ω2λ

2

J − 4 (ω − 1) + ωλ2

J − 2√ω2λ

2

J − 4 (ω − 1).


Temos ωλ2

J − 2 < 0, porque 0 < ω < 2 e λJ < 1, e∣∣∣ωλ

2

J − 2∣∣∣ > λJ

√ω2λ

2

J − 4 (ω − 1),

pois∣∣∣ωλ

2

J − 2∣∣∣2

= ω2λ4

J − 4λ2

Jω + 4 > ω2λ4

J − 4λ2

Jω + 4λ2

J > ω2λ4

J − 4λ2

J (ω − 1)

=[λJ

√ω2λ

2

J − 4 (ω − 1)]2

.

Isso implicad

dω

(ωλJ +

√ω2λ

2

J − 4 (ω − 1))

< 0,

logo ρ (RSOR,ω) e decrescente de 0 ate ωotimo. Para ωotimo 6 ω < 2, ρ (RSOR,ω) = ω − 1 e claramentecrescente. Portanto, ρ (RSOR,ω) atinge o seu mınimo em ωotimo.

Pelo Teorema 3.10, temosλJ = cos

π

n,

logo

ωotimo =2

1 +√

1− λ2

J

=2

1 +√

1− cos2π

n

=2

1 + senπ

n

.

¥Para o quadrado unitario temos

ωotimo =2

1 + sen (π∆x)e consequentemente

ρ (RSOR,ω) =2

1 + sen (π∆x)− 1 =

1− sen (π∆x)1 + sen (π∆x)

.

e usando1− x

1 + x= 1− 2x + O

(x2

),

sen x = x + O(x3

),


1− sen (π∆x)1 + sen (π∆x)

≈ 1− 2π∆x + O(∆x2

).

Portanto, usando o valor otimo de ω no metodo SOR, temos ρ (R) → 1 linearmente quando ∆x → 0, umresultado muito melhor que o obtido nos metodos de Jacobi e de Gauss-Seidel. Para uma comparacao maisprecisa, usando

log (1 + x) = x + O(∆x2

)

temos queR∞ (RSOR) = 2π∆x + O

(∆x2

). (3.55)

Segue queR∞ (RSOR)

R∞ (RGauss-Seidel)≈ 2π∆x

π2∆x2=

2π∆x

.

Em particular, se ∆x = 0.025, temos ωotimo = 1. 8545 e R∞ (RSOR) /R∞ (RGauss-Seidel) = 25.5, isto e, ometodo SOR e 25 vezes mais rapido que o metodo de Gauss-Seidel. Quanto mais refinada a malha, maior ea diferenca na velocidade de convergencia entre os dois metodos.


3.3.4 Convergencia do Metodo de Jacobi Amortecido

3.22 Teorema. Se o metodo de Jacobi converge, entao o metodo de Jacobi amortecido converge para

0 < ω 6 1.

Prova. Vamos escrever a matriz de iteracao RJ,ω do metodo de Jacobi amortecido em funcao da matriz deiteracao do metodo de Jacobi RJ . Temos

RJ = D−1 (D −A)

de modo que

RJ,ω =(

1ω

D

)−1 (1ω

D −A

)= ωD−1

(1ω

D −D + D −A

)= ωD−1

(1ω

D −D

)+ ωD−1 (D −A)

dondeRJ,ω = (1− ω) I + ωRJ . (3.56)

Em particular,RJv = λv

se e somente se[RJ,ω − (1− ω) I] v = ωλv.

Portanto, λJ e um autovalor de RJ se e somente se

λJ,ω = ωλJ + 1− ω (3.57)

e um autovalor de RJ,ω. Logo, se todo autovalor de RJ satisfaz |λJ | < 1 (isto e, ρ (RJ) < 1 equivalente aometodo de Jacobi convergir) e ω < 1, entao

|λJ,ω|2 = (ωλJ + 1− ω)(ωλJ + 1− ω

)

= ω2 |λJ |2 + 2 Re λJω (1− ω) + (1− ω)2

6 ω2 |λJ |2 + 2 |λJ |ω (1− ω) + (1− ω)2

= (ω |λJ |+ 1− ω)2

< 1.

¥Segue do Teorema 3.8 que o metodo de Jacobi amortecido converge para as matrizes de discretizacao doCapıtulo 1 se 0 < ω 6 1.

3.23 Corolario.ρ (RJ,ω) = ω [ρ (RJ)− 1] + 1. (3.58)

Para o quadrado unitario temosρ (RJ,ω) = ω [cos (π∆x)− 1] + 1. (3.59)

Usando

cos x = 1− 12x2 + O

(x4

),

log (1 + x) = x + O(∆x2

),



ρ (RJ,ω) ≈ 1− ωπ2

2∆x2 + O

(∆x4

),

R∞ (RJ,ω) ≈ ωπ2

2∆x2.

Vemos que a velocidade de convergencia do metodo de Jacobi amortecido e da mesma ordem que a do metodode Jacobi, um pouco pior para valores de ω proximos de 1 e muito pior para valores de ω proximos de 0.

3.3.5 Resumo

Metodo ρ (R) R∞ (R)

Jacobi cos (π∆x)π2

2∆x2 + O

(∆x4

)

Gauss-Seidel cos2 (π∆x) π2∆x2 + O(∆x4

)

SOR otimo 1− 2π∆x + O(∆x2

)2π∆x + O

(∆x2

)

Jacobi amortecido 1− ωπ2

2∆x2 + O

(∆x4

)ω

π2

2∆x2 + O

(∆x4

)

3.4 Metodo do Gradiente Conjugado

Nesta secao, A sera sempre uma matriz real simetrica, positiva definida. Neste caso, a resolucao do sistemaAx = b e equivalente a resolucao de um problema de minimizacao de um funcional quadratico:

3.24 Teorema. (Metodo Variacional para a Resolucao de Sistemas Lineares) Seja A ∈ Mn (R) uma matrizsimetrica positiva definida e b ∈ Rn. Entao a solucao do sistema

Ax = b

e o unico ponto x que minimiza o funcional quadratico

f (y) =12ytAy − ytb. (3.60)

Prova: Uma matriz simetrica positiva definida e invertıvel, logo existe uma unica solucao x para o sistemaAx = b. Para provar o teorema, comecamos observando que, como ytAx ∈ R e um escalar, temos

ytAx =(ytAx

)t = xtAty = xtAy.

Daı,

f (y)− f (x) =12ytAy − ytb− 1

2xtAx + xtb

=12ytAy − ytAx− 1

2xtAx + xtAx

=12ytAy − ytAx +

12xtAx

=12ytAy − 1

2ytAx− 1

2xtAy +

12xtAx

=12ytA (y − x)− 1

2xtA (y − x)


ouf (y)− f (x) =

12

(y − x)tA (y − x) . (3.61)

Como A e positiva definida, segue que

(y − x)tA (y − x) = 〈A (y − x) , (y − x)〉 > 0

e(y − x)t

A (y − x) = 0

se e somente se y = x. Portanto,f (y) > f (x)

para todo y 6= x e o mınimo de f ocorre em x. ¥Em muitos problemas, o funcional f tem significado fısico, correspondente a um funcional de energia quequando e minimizado corresponde a um estado de equilıbrio do sistema. Observe que definindo um produtointerno a partir da matriz simetrica positiva definida A da maneira usual por 〈v, w〉A = vtAw e considerandoa norma induzida ‖v‖A = 〈v, v〉1/2

A , o funcional f pode ser escrito na forma

f (y) =12〈y,Ay〉 − 〈y,Ax〉 (3.62)

ouf (y) =

12‖y‖2A − 〈y, x〉A . (3.63)

Outra maneira de enxergar o resultado do teorema anterior e observar que o gradiente do funcional f e

∇f (y) = Ay − b. (3.64)

Se x e um ponto de mınimo temos ∇f (x) = 0, ou seja,

Ax = b.

Este metodo variacional e a base dos metodos iterativos de descida em geral, e do metodo do gradienteconjugado em particular. A ideia e usar as ideias do calculo diferencial para encontrar o mınimo do funcionalquadratico f .

3.4.1 Metodos de Descida

A filosofia dos metodos de descida e comecar com um chute inicial x0 e gerar uma sequencia de iteradosx1, x2, . . . , xk, . . . que satisfazem

f(xk+1

)6 f

(xk

)

ou, melhor ainda,f

(xk+1

)< f

(xk

)

de tal modo que xk convirja para o minimizador de f . Em outras palavras, em um metodo de descidabuscamos encontrar uma sequencia minimizante

(xk

)que convirja para a solucao do sistema.

O passo de xk para xk+1 envolve dois ingredientes: (1) uma direcao de busca e (2) um avanco decomprimento especificado na direcao de busca. Uma direcao de busca significa a escolha de um vetor pk queindicara a direcao que avancaremos de xk para xk+1. O comprimento do avanco e equivalente a escolha deum escalar αk multiplicando o vetor pk. Assim,

xk+1 = xk + αkpk.


A escolha de αk e tambem chamada uma busca na reta, ja que queremos escolher um ponto na retaxk + αpk : α ∈ R

tal quef

(xk + αpk

)6 f

(xk

).

Idealmente, gostarıamos de escolher αk de tal modo que

f(xk+1

)= f

(xk + αkpk

)= min

α∈Rf

(xk + αpk

)

Esta e chamada uma busca na reta exata. Para funcionais quadraticos, a busca na reta exata e trivial eobtemos uma formula para o valor de αk, como veremos a seguir. Denotaremos o resıduo em cada iteracaopor

rk = b−Axk. (3.65)

3.25 Proposicao. Seja αk ∈ R tal que

f(xk + αkpk

)= min

α∈Rf

(xk + αpk

).

Entao

αk =

(pk

)trk

(pk)tApk

=

⟨pk, rk

⟩

〈pk, Apk〉 . (3.66)

Prova: Considere o funcionalg (α) = f

(xk + αpk

).

g e um polinomio quadratico em α, pois

g (α) =12

(xk + αpk

)tA

(xk + αpk

)− (xk + αpk

)tb

=12

(xk

)tAxk − (

xk)t

b +α

2(xk

)tApk +

α

2(pk

)tAxk +

α2

2(pk

)tApk − α

(pk

)tb

= f(xk

)+ α

[12

(pk

)tAxk +

12

(pk

)tAxk − (

pk)t

b

]+

α2

2(pk

)tApk

= f(xk

)− α(pk

)tArk +

α2

2(pk

)tApk,

portanto o mınimo de g e atingido no vertice −B/2A da parabola Y = AX2 + BX + C. ¥Observe que αk = 0 se e somente se

(pk

)trk = 0, isto e, a direcao de busca e ortogonal ao resıduo. Como

gostarıamos sempre que possıvel de ter xk+1 6= xk, devemos sempre escolher a direcao de busca de forma anao ser ortogonal a rk. Se esta escolha e feita, entao teremos sempre f

(xk+1

)< f

(xk

).

Exemplo 1. (Metodo de Gauss-Seidel) Considere o metodo de descida em que as primeiras n direcoes debusca p1, . . . , pn sao os vetores e1, . . . , en da base canonica de Rn, e isso e repetido a cada n iteracoes,de modo que pk+n = ek para todo k = 1, . . . , n, com uma busca na reta exata executada em cadaiteracao. Entao cada grupo de n iteracoes corresponde a uma iteracao do metodo de Gauss-Seidel.

Exemplo 2. (Metodo SOR) Usando as mesmas direcoes de busca do exemplo anterior, mas com xk+1 =xk + ωαkpk, ω 6= 1, obtemos um metodo de descida em que as buscas nas retas sao inexatas. Cadagrupo de n iteracoes corresponde a uma iteracao do metodo SOR.


3.4.2 Metodo da Descida Mais Acentuada

Do Calculo Diferencial, sabemos que a direcao em que a funcao cresce a uma taxa mais rapida a partir deum ponto e a direcao do gradiente neste ponto. Esta observacao e a base da escolha da direcao de busca nometodo da descida mais acentuada. Em outras palavras, escolhemos

pk = −∇f(xk

)= b−Axk

oupk = rk. (3.67)

Buscar na direcao da descida mais acentuada e uma ideia natural, mas que na pratica nao funciona semmodificacoes. De fato, em alguns casos o metodo e de velocidade comparavel a do metodo de Jacobi, comona matriz de discretizacao da formula de cinco pontos aplicada ao problema descrito na primeira secao destecapıtulo [Watkins]:

∆x = 0.1 ∆x = 0.05 ∆x = 0.025Jacobi 299 1090 3908Descida Mais Acentuada 304 1114 4010

De fato, como as iteracoes do metodo de descida mais acentuada sao bem mais custosas que as do metodode Jacobi, o primeiro e muito pior que este ultimo.

Para entender melhor o metodo da descida mais acentuada, porque ele pode ser lento e as modificacoes quevamos fazer para torna-lo mais rapido levando ao metodo do gradiente conjugado, vamos entender o processodo ponto de vista geometrico. Como vimos na demonstracao do Teorema 3.24, o funcional quadratico f eda forma

f (y) =12

(y − x)tA (y − x) + c (3.68)

onde c = f (x) = 12xtAx − xtb e uma constante. Ja que A e uma matriz simetrica, existe uma matriz

ortogonal P tal que P tAP e uma matriz diagonal D , cujos valores na diagonal principal sao exatamente osautovalores positivos de A. Nas coordenadas

z = P t (y − x) ,

o funcional f tem a forma

f (z) =12ztDz + c =

12

n∑

i=1

λiz2i + c. (3.69)

As curvas de nıvel do funcional f neste sistema de coordenadas sao elipses (em R2, elipsoides em R3 ehiperelipsoides em Rn) centradas na origem com eixos paralelos aos eixos coordenados e f (0) = c e nıvelmınimo de f ; elipses correspondentes a menores valores de f estao dentro de elipses correspondentes amaiores valores de f . Como P e uma aplicacao ortogonal, as curvas de nıvel de f no sistema de coordenadasoriginal tambem sao elipses, centradas em x, e uma reta de um ponto y ate o ponto x corta elipses de nıveiscada vez menores ate chegar ao mınimo da funcao f em x, centro de todas as elipses. O vetor gradiente eperpendicular as curvas de nıvel, logo e perpendicular as elipses. Seguir a direcao de descida mais acentuadaequivale a cortar a elipse que contem xk ortogonalmente na direcao do interior da elipse ate encontrar umponto xk+1 situado em uma elipse que a reta tangencie, pois a partir daı a reta ira na direcao de elipses comnıveis maiores, portanto este e o ponto da reta onde f atinge o seu mınimo. Em particular, vemos que aproxima direcao pk+1 e ortogonal a direcao anterior pk, tangente a esta elipse. Em geral, a direcao de descidamais acentuada nao e a direcao de x (quando bastaria uma iteracao para atingir a solucao exata) a nao serque A seja um multiplo escalar da identidade, de modo que todos os autovalores de A sao iguais e as elipsessao cırculos. Por outro lado, se os autovalores de A tem valores muito diferentes uns dos outros, com algunsmuito pequenos e alguns muito grandes, as elipses serao bastante excentricas e, dependendo do chute inicial,


a convergencia pode ser muito lenta (matrizes com estas propriedades sao chamadas mal-condicionadas; paraque o metodo de descida acentuada seja lento, a matriz A nao precisa ser muito mal-condicionada).

Como vimos na secao anterior, os algoritmos de Gauss-Seidel e SOR podem ser encarados como algoritmosde descida. A discussao no paragrafo anterior tambem pode ser usada para entender a relativa lentidao destesalgoritmos.

3.4.3 Metodo do Gradiente Conjugado

Todos os metodos iterativos que vimos neste capıtulo sao limitados pela sua falta de memoria, no sentido deque apenas informacao sobre xk e usada para obter xk+1. Toda a informacao sobre as iteracoes anteriores edeletada. O metodo do gradiente conjugado e uma variacao simples do metodo da descida mais acentuadaque funciona melhor porque a informacao obtida atraves das iteracoes anteriores e utilizada.

Para entender brevemente como isso funciona, observe que depois de j iteracoes xk+1 = xk + αkpk deum metodo de descida temos

xj = x0 + α0p0 + α1p

1 + . . . + αj−1pj−1,

de modo que xj esta no subespaco afim gerado pelo chute inicial x0 e pelos vetoresp0, p1, . . . , pj−1

.

Enquanto o metodo da descida mais acentuada minimiza o funcional de energia f apenas ao longo das jretas xk + αkpk, cuja uniao constitui apenas um pequeno subconjunto de x0 +

⟨p0, p1, . . . , pj−1

⟩, o metodo

do gradiente conjugado minimiza f sobre todo o subespaco afim x0 +⟨p0, p1, . . . , pj−1

⟩.

Para definir as direcoes de busca do metodo do gradiente conjugado (que e, antes de mais nada, ummetodo de descida), lembramos que o funcional f foi escrito na forma

f (y) =12‖y‖2A − 〈y, x〉A .

Defina o erroe = x− y. (3.70)

Pela regra do paralelogramo, temos

‖x + y‖2A + ‖x− y‖2A = 2 ‖x‖2A + 2 ‖y‖2A ,

donde

2 ‖y‖2A = ‖x− y‖2A + ‖x‖2A + 2 〈y, x〉A + ‖y‖2A − 2 ‖x‖2A= ‖x− y‖2A + 2 〈y, x〉A − ‖x‖2A + ‖y‖2A ,

ou‖y‖2A − 2 〈y, x〉A = ‖x− y‖2A − ‖x‖2A .

Logo, podemos escrever

f (y) =12‖e‖2A −

12‖x‖2A . (3.71)

Consequentemente, minimizar o funcional f e equivalente a minimizar a A-norma do erro.Agora, em um metodo de descida, depois de j iteracoes temos:

ej = x− xj = x− x0 − (α0p

0 + α1p1 + . . . + αj−1p

j−1)

= e0 − (α0p

0 + α1p1 + . . . + αj−1p

j−1).

Logo, minimizar∥∥ej

∥∥2

Ae equivalente a minimizar

∥∥e0 − (α0p

0 + α1p1 + . . . + αj−1p

j−1)∥∥

A,

o que por sua vez e equivalente a encontrar a melhor aproximacao do vetor e0 no subespaco Wj =⟨p0, p1, . . . , pj−1

⟩.

Esta e dada pelo lema da melhor aproximacao:


3.26 Proposicao. Sejam A ∈ Mn (R) uma matriz simetrica positiva definida, v ∈ Rn e W um subsespacode Rn. Entao existe um unico w ∈ W tal que

‖v − w‖A = minz∈W

‖v − z‖A .

O vetor w e caracterizado pela condicao v − w ⊥A W .

Segue deste resultado que∥∥ej

∥∥A

e minimizado quando escolhemos p = α0p0 + α1p

1 + . . . + αj−1pj−1 ∈ Wj

tal que ej = e0 − p satisfazej ⊥A pi para i = 1, . . . , j − 1. (3.72)

Definicao. Dois vetores y, z que sao ortogonais com respeito ao produto interno 〈·, ·〉A, isto e, tais que

〈y, z〉A = 0

sao chamados conjugados.

Nosso objetivo entao e desenvolver um metodo em que o erro a cada passo e conjugado com todas as direcoesde busca anteriores. O proximo resultado, que e basicamente uma reafirmacao da Proposicao 3.25, mostraque em qualquer metodo de descida em que a busca na reta e exata satisfaz automaticamente ej ⊥A pj−1,isto e, (3.72) e valido para a ultima iteracao (o erro da iteracao presente e A-ortogonal a direcao de buscada iteracao anterior).

3.27 Proposicao. Seja xk+1 = xk + αkpk obtido atraves de uma busca na reta exata. Entao

rk+1 ⊥ pk

eek+1 ⊥A pk.

Prova: Temosb−Axk+1 = b−Axk − αkApk,

de modo que a sequencia dos resıduos e dada pela formula

rk+1 = rk − αkApk. (3.73)

Logo,⟨rk+1, pk

⟩=

⟨rk+1, pk

⟩− αk

⟨Apk, pk

⟩=

⟨rk, pk

⟩−⟨pk, rk

⟩

〈pk, Apk〉⟨Apk, pk

⟩= 0.

Alem disso, comoAek+1 = rk+1,

segue que ⟨ek+1, pk

⟩A

=⟨Aek+1, pk

⟩=

⟨rk+1, pk

⟩= 0.

¥O significado geometrico deste resultado e que o mınimo do funcional f na reta xk + αkpk ocorre quando aderivada direcional de f na direcao de busca e zero, ou seja,

0 =∂f

∂pk

(xk+1

)=

⟨∇f(xk+1

), pk

⟩=

⟨rk+1, pk

⟩.

De acordo com a Proposicao 3.27, depois do primeiro passo temos e1 ⊥A p0. Para manter os errossubsequentes conjugados a p0, como

ek+1 = x− xk+1 = x− xk − αkpk


ouek+1 = ek − αkpk, (3.74)

basta escolher as direcoes de busca subsequentes conjugadas a p0. Se escolhemos p1 conjugado a p0, obtemosx2 para o qual o erro satisfaz e2 ⊥A p1; como p1 ⊥A p0, segue de (3.74) que e2 ⊥A p0 tambem. Para manteros erros subsequentes conjugados a p0 e p1, basta escolher as direcoes de busca subsequentes conjugadas ap0 e p1. Assim, vemos que para obter a condicao (3.72) basta escolher as direcoes de busca de tal forma que

pi ⊥A pj para todos i 6= j.

Um metodo com estas caracterısticas e chamado um metodo de direcoes conjugadas. Estes resultadossao resumidos na proposicao a seguir:

3.28 Teorema. Se um metodo emprega direcoes de busca conjugadas e performa buscas na reta exatas,entao

ej ⊥A pi para i = 1, . . . , j − 1,

para todo j. Consequentemente ∥∥ej∥∥

A= min

p∈Wj

∥∥e0 − p∥∥

A,

onde Wj =⟨p0, p1, . . . , pj−1

⟩.

Prova: A demonstracao e por inducao. Para j = 1, temos e1 ⊥A p0 pela Proposicao 3.27 porque a buscana reta e exata. Em seguida, assuma ej ⊥A pi para i = 1, . . . , j − 1; queremos mostrar que ej+1 ⊥A pi

para i = 1, . . . , j. Comoej+1 = ej − αjp

j ,

para i = 1, . . . , j − 1 temos⟨ej+1, pi

⟩A

=⟨ej − αjp

j , pi⟩

A=

⟨ej , pi

⟩A− αj

⟨pj , pi

⟩A

= 0− 0 = 0

porque as direcoes de busca sao conjugadas. ej+1 ⊥A pj segue novamente da Proposicao 3.27. ¥Quando a direcao inicial e dada pelo vetor gradiente de f , como na primeira iteracao do metodo da descidamais acentuada, obtemos o metodo do gradiente conjugado. As direcoes subsequentes sao escolhidasatraves de A-ortogonalizar o resıduo (ou vetor gradiente de f , que e a direcao de busca em cada iteracaodo metodo da descida mais acentuada) com todas as direcoes de busca anteriores, para isso utilizando oalgoritmo de Gram-Schmidt. Assim, dado um chute inicial p0, a primeira direcao e

p0 = −∇f(x0

)= b−Ax0 = r0

ou seja, a direcao inicial e o primeiro resıduo:

p0 = r0. (3.75)

Depois de k passos com direcoes de busca conjugadas p0, . . . , pk, escolhemos

pk+1 = rk+1 −k∑

i=0

ckipi (3.76)

onde os cki sao dados pelo algoritmo de Gram-Schmidt:

cki =

⟨rk+1, pi

⟩A

〈pi, pi〉A. (3.77)

de forma que pk+1 ⊥A pi para todos i = 1, . . . , k. Felizmente, como veremos a seguir depois de algum trabalhopreliminar (Corolario 3.32), cki = 0 para todo i exceto i = k, o que torna necessario que apenas a direcao


de busca mais recente pk seja armazenada na memoria do computador, o que garante que a implementacaodo gradiente conjugado e eficiente:

pk+1 = rk+1 −⟨rk+1, pk

⟩A

〈pk, pk〉Apk = rk+1 −

⟨rk+1, Apk

⟩

〈pk, Apk〉 pk (3.78)

ou, definindo

βk = −⟨rk+1, Apk

⟩

〈pk, Apk〉 , (3.79)

temos quepk+1 = rk+1 + βkpk. (3.80)

Esta e a modificacao do metodo do gradiente conjugado em relacao ao metodo da descida mais acentuada,no qual tomamos pk+1 = rk+1.

Podemos obter uma expressao mais simples para o escalar βk, em funcao apenas dos resıduos. Comefeito, temos ⟨

rk+1, rk+1⟩

=⟨rk+1, rk

⟩− αk

⟨rk+1, Apk

⟩= −αk

⟨rk+1, Apk

⟩

porque os resıduos obtidos atraves do metodo do gradiente conjugado sao mutualmente ortogonais (vejaCorolario 3.31), logo

β = −⟨rk+1, Apk

⟩

〈pk, Apk〉 =

⟨rk+1, rk+1

⟩

αk 〈pk, Apk〉 .

Temos

αk =

⟨pk, rk

⟩

〈pk, Apk〉 =

⟨rk + βpk−1, rk

⟩

〈pk, Apk〉 =

⟨rk, rk

⟩

〈pk, Apk〉 ,

porque⟨pk−1, rk

⟩= 0 pela Proposicao 3.27, logo

αk =

⟨rk, rk

⟩

〈pk, Apk〉 . (3.81)

Portanto

β =

⟨rk+1, rk+1

⟩

〈rk, rk〉 . (3.82)

Podemos obter um algoritmo ainda mais eficiente para o metodo do gradiente conjugado se observarmos quepara calcular o resıduo rk+1 = b−Axk+1 em cada iteracao nao e necessario calcular Axk+1 explicitamente;de fato, como vimos na demonstracao da Proposicao 3.27, temos rk+1 = rk − αkApk. Assim, um algoritmo


eficiente para o metodo do gradiente conjugado poderia ser escrito da seguinte forma:

initialize x;set b;r ← b−Ax;rScalarR ← 〈r, r〉 ;set M ; //maximumNumberOfIterationsnumberOfIterations = 0;do until numberOfIterations > M

Ap ← Ap;pScalarAp ← 〈p,Ap〉 ;α ← rScalarR/pScalarAp;x ← x + αp;r ← r − αAp;rNewScalarRNew ← 〈r, r〉 ;β ← rNewScalarRNew/rScalarR;p ← r + βp;rScalarR ← rNewScalarRNew;numberOfIterations + +;

3.5 Convergencia do Metodo do Gradiente Conjugado

Vamos agora provar uma serie de resultados com o objetivo principal de demonstrar o fato mencionadoacima que cki = 0 para todo i = 1, . . . , k − 1 e tambem que o metodo do gradiente conjugado converge emaritmetica exata em precisas n iteracoes para uma matriz de tamanho n.

Definicao. Dada uma matriz A ∈ Mn (C) e um vetor v ∈ Cn, o espaco de Krylov Kj (A, v) e o subespaco⟨v,Av, . . . , Aj−1v

⟩.

3.29 Teorema. Depois de j iteracoes do algoritmo do gradiente conjugado (com rk 6= 0 em cada iteracao),temos ⟨

p0, p1, . . . , pj−1⟩

=⟨r0, r1, . . . , rj−1

⟩= Kj

(A, r0

).

Prova: A demonstracao e por inducao. O resultado e trivial para j = 0, pois p0 = r0. Assuma o resultadovalido para j − 1. Em primeiro lugar, mostraremos que

⟨r0, r1, . . . , rj

⟩ ⊂ Kj+1

(A, r0

). (3.83)

Em vista da hipotese de inducao, basta mostrar que rj ∈ Kj+1

(A, r0

). Como rj = rj−1 − αj−1Apj−1 e

rj−1 ∈ Kj

(A, r0

) ⊂ Kj+1

(A, r0

)por hipotese de inducao, basta provar que Apj−1 ∈ Kj+1

(A, r0

). Mas,

tambem por hipotese de inducao, pj−1 ∈ Kj+1

(A, r0

), logo

Apj−1 ∈ Kj

(A,Ar0

)=

⟨Ar0, A2r0, . . . , Ajr0

⟩ ⊂ ⟨r0, Ar0, A2r0, . . . , Ajr0

⟩= Kj+1

(A, r0

).

Em seguida, mostraremos que⟨p0, p1, . . . , pj

⟩ ⊂ ⟨r0, r1, . . . , rj

⟩. (3.84)

Por hipotese de inducao, basta provar que pj ∈ ⟨r0, r1, . . . , rj

⟩. Isso segue de (3.76) e da hipotese de inducao.

Ate aqui provamos que⟨p0, p1, . . . , pj

⟩ ⊂ ⟨r0, r1, . . . , rj

⟩ ⊂ Kj+1

(A, r0

). (3.85)


Para provar que eles sao iguais, basta mostrar que eles tem a mesma dimensao. Isso decorre de

dim⟨r0, r1, . . . , rj

⟩6 j + 1,

dimKj+1

(A, r0

)6 j + 1

edim

⟨p0, p1, . . . , pj

⟩= j + 1,

o ultimo porque os vetores p0, p1, . . . , pj sao vetores nao-nulos A-ortogonais. ¥

3.30 Corolario. Depois de j iteracoes do algoritmo do gradiente conjugado, temos

ej ⊥A Kj

(A, r0

)

para todo j.

Prova: Segue imediatamente do teorema anterior e do Teorema 3.28. ¥

3.31 Corolario. Depois de j iteracoes do algoritmo do gradiente conjugado, temos

rj ⊥ Kj

(A, r0

)

para todo j.

Prova: Em vista do Teorema 3.29, basta provar que rj ⊥ p0, p1, . . . , pj−1 para todo j. Como Aej+1 = rj+1,⟨rj+1, pi

⟩=

⟨Aej+1, pi

⟩=

⟨ej+1, pi

⟩A

= 0

para todo i = 1, . . . , j − 1, como vimos na demonstracao do Teorema 3.28. ¥

3.32 Corolario. cki = 0 para todo i = 1, . . . , k − 1.

Prova: Temos que provar que ⟨rk+1, pi

⟩A

=⟨rk+1, Api

⟩= 0

para todos i = 1, . . . , k − 1. Pelo Teorema 3.29, pi ∈ ⟨p0, p1, . . . , pi

⟩=

⟨r0, Ar0, . . . , Air

⟩= Ki+1

(A, r0

),

logoApi ∈ ⟨

Ar0, A2r0, . . . , Ai+1r⟩ ⊂ Ki+2

(A, r0

) ⊂ Kk+1

(A, r0

)

e o resultado segue do corolario anterior. ¥

3.33 Teorema. Seja A uma matriz simetrica positiva definida n×n. Entao o metodo do gradiente conjugadoconverge em n iteracoes.

Prova: Se fizemos n − 1 iteracoes em obter x, pelo Corolario 3.32 os vetores r0, r1, . . . , rn−1 formam umabase ortogonal para Rn. Depois de mais uma iteracao, de acordo com este mesmo corolario o resıduo rn

satisfaz rn ⊥ ⟨r0, r1, . . . , rn−1

⟩= Rn, logo rn = 0. ¥

De fato, na maioria das aplicacoes o metodo do gradiente conjugado converge ainda mais rapido, se apenasuma boa aproximacao e requerida. Defina o numero de condicao de uma matriz simetrica positiva definidapor

κ (A) =max λ : λ e um autovalor de Amin λ : λ e um autovalor de A ; (3.86)

assim, quanto maior o numero de condicao de uma matriz, ela e mais mal-condicionada e a convergenciade metodos de descida e mais vagarosa. Pode-se provar a seguinte estimativa de erro para o metodo dogradiente conjugado (veja [Strikwerda]):

∥∥ek∥∥

A6 2

∥∥e0∥∥

A

(√κ (A)− 1√κ (A) + 1

)k

. (3.87)


Esta estimativa e uma estimativa grosseira, mas mostra que o metodo do gradiente conjugado convergemais rapidamente para matrizes bem-condicionadas (κ (A) ∼ 1). Uma comparacao entre a velocidade deconvergencia dos dois metodos para a matriz de discretizacao da formula de cinco pontos aplicada ao problemadescrito na primeira secao deste capıtulo, desta vez com o tamanho das matrizes indicado na linha superiorda tabela, e dada a seguir [Watkins].

n = 81 n = 361 n = 1521Descida Mais Acentuada 304 1114 4010Gradiente Conjugado 29 60 118

No caso desta matriz de discretizacao no quadrado unitario temos

κ (A) =sen2 (n− 1) π

2n

sen2π

2n

= cot2π

2n= cot2

π∆x

2≈ 4

π2∆x2

de modo que √κ (A)− 1√κ (A) + 1

≈ 1− π∆x/21 + π∆x/2

≈ 1− π∆x,

o que da uma velocidade de convergencia para o metodo do gradiente conjugado duas vezes maior que ado metodo SOR com o fator de relaxamento otimo. No entanto, deve-se ter em mente que enquanto que ataxa de covergencia que obtivemos para o metodo SOR e precisa, a estimativa de erro (3.87) para o metododo gradiente conjugado e apenas um limitante superior grosseiro (veja [Watkins] para algumas estimativasmelhoradas).

Capıtulo 4

Metodos Multigrid

Neste capıtulo consideraremos o metodo multigrid, que e o metodo mais rapido para resolver equacoeselıpticas em geral. Embora o metodo possa ser empregado em malhas de elementos finitos e volumes fini-tos tambem, neste capıtulo consideraremos o seu emprego apenas em malhas de diferencas finitas para aequacao de Poisson no quadrado. A tabela a seguir (adaptada de [TOS]) compara o custo de processamentoem uma maquina serial de alguns dos metodos mais populares para resolver sistemas lineares que surgem nadiscretizacao do problema de Poisson (a excecao do metodo de eliminacao gaussiana cujo custo de armazena-mento e O

(n2

), todos os demais metodos tem custo de armazenamento O (n)). Como estamos comparando

metodos diretos (eliminacao gaussiana e transformada de Fourier rapida (FFT) ) com metodos iterativos(todos os demais), assumimos um unico criterio de parada para os varios metodos iterativos; se o criterio deparada for escolhido da ordem do erro de discretizacao da malha, um fator O (log n) deve ser multiplicadopara todos os metodos iterativos, a excecao do multigrid completo.

Metodo numero de operacoes (2D; n = N2)Eliminacao Gaussiana O

(n3

)Jacobi O

(n2

)Gauss-Seidel O

(n2

)SOR O

(n3/2

)Gradiente Conjugado O

(n3/2

)FFT O (n log n)Multigrid iterativo O (n)Multigrid completo O (n)

A ideia do metodo multigrid e baseada em dois princıpios: suavizacao do erro e a sua correcao emum grid mais grosseiro (menos refinado). Estes princıpios serao explicados em detalhes nas proximassecoes.

Em linhas gerais, a ideia basica e eliminar os componentes de alta frequencia do erro em uma malharefinada. Para que isso ocorra, e necessario que estes componentes de alta frequencia correspondam aosmenores autovalores da matriz de iteracao porque, como vimos na Secao 3.2.2, estes sao eliminados rapi-damente pelos metodos iterativos lineares (a velocidade de convergencia de cada metodo e dada pelo raioespectral da matriz de iteracao, que corresponde ao valor absoluto do maior autovalor |λ1| < 1, enquantoque as componentes do erro correspondentes aos menores autovalores λj convergem para zero muito maisrapidamente (|λj/λ1| ¿ 1); isso significa que este metodo iterativo suaviza o erro, pois quanto maior ainfluencia das componentes de maior frequencia (maior oscilacao), menos suave e a funcao. Aqui e util fazeruma analogia com a serie de Fourier: e exatamente a presenca de componentes de oscilacao arbitrariamentemaior que permite que a serie convirja para uma funcao nao diferenciavel, ou mesmo descontınua; se a seriefor truncada a qualquer momento o resultado e sempre uma funcao suave, pois e a combinacao linear finitade autofuncoes suaves. Esta visualizacao tambem permanece verdade para funcoes discretizadas em mal-

85


has de diferencas finitas escritas como uma combinacao linear das autofuncoes da matriz de iteracao nestamalha: mesmo que o numero de componentes da funcao seja finito, porque a malha e discreta a presenca decomponentes de alta oscilacao dao origem a um grafico com um aspecto escarpado, nao suave.

Assim, como o nosso objetivo e eliminar apenas as componentes de alta frequencia do erro, e nao todo oerro, poucas iteracoes do metodo iterativo sao necessarias nesta malha refinada, onde o custo computacional ealto (malhas muito refinadas significa que elas possuem muitos pontos, o que por sua vez implica em matrizesde discretizacao muito grandes). Ocorre que algumas autofuncoes de frequencia baixa em uma malha maisrefinada correspondem a autofuncoes de frequencia alta em uma malha mais grosseira (como veremos). Umavez tendo eliminado as componentes de alta frequencia do erro na malha mais refinada, tendo deixado ascomponentes de baixa frequencia praticamente intocadas, transferimos o problema para uma malha maisgrosseira, cujos componentes de alta frequencia do erro correspondem a alguns dos componentes de baixafrequencia do erro na malha mais refinada anterior, que nao puderam ser eliminados com as poucas iteracoesdo metodo iterativo permitidas na malha mais refinada. Com poucas iteracoes do metodo iterativo nestamalha mais grosseira, estes erros tambem sao rapidamente eliminados, a um custo computacional mais baixodo que se tivessemos tentado elimina-los ainda na malha mais refinada. Este processo e a correcao do erroem uma malha mais grosseira. Ele e repetido em malhas cada vez mais grosseiras ate que todo o erro eeliminado, a um custo computacional muito mais baixo do que se tivessemos trabalhado sempre na malhamais refinada original.

4.1 A Malha de Multigrid

A discretizacao uniforme do problema de Poisson −∆u = f em Ω,

u = 0 sobre ∂Ω,

onde Ω = (0, 1)2 ⊂ R2 e o quadrado unitario, sera denotada por −∆huh = fh em Ωh,

uh = 0 sobre ∂Ωh,(4.1)

onde uh como usual denota a solucao do problema discretizado (aproximacao da solucao exata), fh a dis-cretizacao da funcao f em Ωh,

h =1n

, (4.2)

Ωh = (x, y) ∈ Ω : (x, y) = (ih, jh) , 1 6 i, j 6 n− 1 ,

∂Ωh = (x, y) ∈ ∂Ω : (x, y) = (ih, jh) , i = 0 ou i = n e 0 6 j 6 n; j = 0 ou j = n e 0 6 i 6 ne

−∆huh =1h2

−1−1 4 −1

−1

(4.3)

ou, em outras palavras,

−∆huh =−uh (xi−1, yj)− uh (xi+1, yj) + 4uh (xi, yj)− uh (xi, yj−1)− uh (xi, yj+1)

h2,

com (xi, yj) = (ih, jh), e o operador de discretizacao dado pela formula dos cinco pontos. Denotaremosusualmente a solucao aproximada uh na iteracao k (ou seja, uma aproximacao da solucao discretizada, deacordo com o metodo iterativo utilizado) por

umh (4.4)


de modo que o erro do metodo iterativo na iteracao m e dado por

emh (xi, yj) = uh (xi, yj)− um

h (xi, yj) . (4.5)

Em geral, tomaremos n par, ou mesmo n = 2p para algum p. Assim, uma malha Ωh e mais refinada queuma malha Ω2h (esta e mais grosseira que a primeira). Temos uma sequencia de malhas progressivamentemais grosseiras:

Ωh ⊂ Ω2h ⊂ Ω4h ⊂ . . . ⊂ Ω2ph = Ω1,

onde Ω1 possui apenas uma celula.

4.2 Frequencias Altas e Baixas

Para analizar as propriedades de suavizacao de um metodo iterativo de maneira rigorosa, precisamos dis-tinguir de maneira precisa entre as frequencias baixas e altas. Estas devem ser definidas de acordo com amalha usada.

As autofuncoes dos metodos iterativos lineares considerados no capıtulo anterior sao exatamente asautofuncoes do laplaciano discretizado −∆h na malha discretizada Ωh, dadas por

ϕklh (x, y) = sen kπx sen lπy, 1 6 k, l 6 n− 1 (4.6)

onde x, y denotam as variaveis discretizadas (isto e, x = ih e y = jh para 0 6 i, j 6 n). Assim, o erro nam-esima iteracao pode ser escrito na forma

emh (x, y) =

n−1∑

k,l=1

αmk,lϕ

klh (x, y) =

n−1∑

k,l=1

αmk,l sen kπx sen lπy. (4.7)

O erro ser suavizado significa que apos algumas poucas iteracoes temos∣∣∣αm

k,l

∣∣∣ muito pequeno para k, l grandes,isto e, para

ϕklh (x, y) = sen kπx sen lπy de alta frequencia,

enquanto que o valor de∣∣∣αm

k,l

∣∣∣ para k, l pequenos, isto e, para

ϕklh (x, y) = sen kπx sen lπy de baixa frequencia,

pode ter mudado muito pouco. Como o fato de k, l serem grandes ou pequenos e definido relativamente deacordo com o valor de n (valores de k, l proximos de n sao considerados grandes, enquanto que valores de k, ldistantes de n sao considerados pequenos), segue que autofuncoes de baixa frequencia em uma malha maisrefinada (n maior) podem ser autofuncoes de alta frequencia em uma malha mais grosseira (n relativamentepequeno). Para propositos de analise, vamos dar uma definicao precisa a este conceito:

Definicao. Para 1 6 k, l 6 n− 1, dizemos que ϕklh e uma autofuncao (ou componente) de

baixa frequencia se max (k, l) <n

2,

alta frequencia sen

26 max (k, l) < n.

Alem disso, se considerarmos especialmente a passagem da malha mais refinada Ωh para a malha maisgrosseira Ω2h com o dobro do espacamento de malha, apenas as autofuncoes de frequencias mais baixas emΩh sao visıveis em Ω2h, pois todas as autofuncoes de frequencia alta em Ωh coincidem com as autofuncoesde frequencia baixa em Ω2h ou desaparecem em Ω2h. De fato, como

ϕk,lh (x, y) = −ϕn−k,l

h (x, y) = −ϕk,n−lh (x, y) = ϕn−k,n−l

h (x, y) para (x, y) ∈ Ω2h, (4.8)


estas quatro autofuncoes nao podem ser distinguidas umas das outras em Ω2h. Alem disso, se k = n/2 oul = n/2, temos

ϕk,lh (x, y) = 0 para (x, y) ∈ Ω2h. (4.9)

Para provar estas afirmacoes, escrevemos, por exemplo,

ϕn−k,lh (i (2h) , j (2h)) = sen

((n− k)π

2i

n

)sen

(lπ

2j

n

)= sen

(−kπ

2i

n+ 2iπ

)sen

(lπ

2j

n

)

= − sen(

kπ2i

n

)sen

(lπ

2j

n

)= −ϕk,l

h (i (2h) , j (2h))

e

ϕn/2,lh (i (2h) , j (2h)) = ϕ

n/2,lh

(2i

n,2j

n

)= sen

(n

2π

2i

n

)sen

(lπ

2j

n

)= sen iπ sen

2jlπy

n= 0.

Assim, podemos decompor o erro em duas somas, uma representando os componentes de baixa frequenciae a outras os componentes de alta frequencia:

emh (x, y) =

n/2−1∑

k,l=1

αmk,lϕ

klh (x, y) +

n−1∑

max(k,l)> n2

αmk,lϕ

klh (x, y) (4.10)

=∑baixa

αmk,lϕ

klh (x, y) +

∑altaαm

k,lϕklh (x, y) . (4.11)

4.3 Suavizacao do Erro

Os dois metodos iterativos classicos, o metodo de Jacobi amortecido e o metodo de Gauss-Seidel (incluindo ometodo SOR) sao metodos iterativos lineares suavizadores de erro. Como ja vimos acima, isso significa apenasque o erro torna-se mais suave com poucas iteracoes, mesmo que nao fique necessariamente menor (em outraspalavras, aqui a velocidade de convergencia nao e o fator principal). Componentes de alta frequencia do errosao eliminadas rapidamente, em comparacao com as componentes de baixa frequencia. As propriedades desuavizacao de cada um dos metodos dependem da escolha correta dos parametros de suavizacao e, no caso dometodo de Gauss-Seidel, tambem da ordenacao dos pontos da malha. Apesar do metodo de Gauss-Seidel serum melhor suavizador que o metodo de Jacobi amortecido, analisaremos rigorosamente apenas este ultimo(sua analise e mais simples porque as autofuncoes da sua matriz de iteracao sao as mesmas do laplacianodiscretizado; veja [TOS] para uma analise completa do poder de suavizacao do metodo de Gauss-Seidel).Uma comparacao entre os poderes suavizadores dos metodos e dada na seguinte tabela (adaptada de [TOS]):

Metodo Fator suavizante SuavizacaoJacobi amortecido, ω = 1 1 NenhumaJacobi amortecido, ω = 0.5 0.75 Nao satisfatoriaJacobi amortecido, ω = 0.8 0.6 AceitavelGauss-Seidel (ordem lexicografica) 0.5 BoaGauss-Seidel (ordem vermelho-negra) 0.25 Muito boa

4.3.1 Metodo de Jacobi Amortecido

Embora no que se refere a velocidade de convergencia, a escolha de ω = 1 no metodo de Jacobi amortecidoe a melhor possıvel (ou seja, correspondendo ao metodo de Jacobi), isso nao e verdade com respeito aspropriedades de suavizacao do erro, como veremos a seguir. A formula de iteracao para o metodo de Jacobipara o problema de Poisson discretizado e dada por

uk+1h (xi, yj) =

ukh (xi−1, yj) + uk

h (xi+1, yj) + ukh (xi, yj−1) + uk

h (xi, yj+1) + h2fh (xi, yj)4

. (4.12)


Em notacao de operadores, esta formula pode ser escrita como

uk+1h = Rhuk

h +h2

4fh, (4.13)

onde o operador de iteracao Rh e dado por

Rh = Ih − h2

4Lh, (4.14)

Ih sendo o operador identidade e Lh = −∆h. No metodo de Jacobi amortecido, introduzimos o parametrode relaxamento 0 < ω 6 1:

uk+1h (xi, yj) = uk

h (xi, yj)

+ ω

(uk

h (xi−1, yj) + ukh (xi+1, yj) + uk

h (xi, yj−1) + ukh (xi, yj+1) + h2fh (xi, yj)

4− uk

h (xi, yj))

.

Logo

uk+1h = Ihuk

h + ω

(Shuk

h +h2

4fh − Ihuk

h

)

= Ihukh + ω

(Ihuk

h −h2

4Lhuk

h +h2

4fh − Ihuk

h

)

= Ihukh −

ωh2

4Lhuk

h +ωh2

4fh,

ou

uk+1h = Rh (ω)uk

h +ωh2

4fh, (4.15)

onde

Rh (ω) = Ih − ωh2

4Lh. (4.16)

Em notacao estencil, o operador iteracao para o metodo de Jacobi amortecido pode ser escrito na forma

Rh (ω) =

1

− ωh2

41h2

−1−1 4 −1

−1

=

14

14 1− ω 1

414

ou tambem

Rh (ω) =ω

4

1

1 4(

1ω− 1

)1

1

.

Em particular, de (4.16) segue que

Rh (ω) = Ih +ωh2

4∆h,

logo os autovalores de Rh e −∆h estao relacionados da seguinte forma: λ e um autovalor de −∆h se esomente se

(Rh − Ih) v = −ωh2

4λv,


isto e, se e somente se

λh (ω) = 1− ωh2

4λ (4.17)

e um autovalor de Rh e as autofuncoes de Rh sao as mesmas autofuncoes de −∆h. As autofuncoes de −∆h

sao, como ja vimos,ϕkl

h (x, y) = sen kπx sen lπy, 1 6 k, l 6 n− 1,

enquanto que os correspondentes autovalores de −∆h sao (veja o Teoremas 3.10)

λkl =2h2

(2− cos kπh− cos lπh) .

Logo, os correspondentes autovalores de Rh sao

λklh (ω) = 1− ω

2(2− cos kπh− cos lπh) . (4.18)

[O raio espectral de Rh, correspondente ao maior autovalor em modulo, e

ρ (Rh) =∣∣∣λ1,1

h

∣∣∣ = |1− ω (1− cosπh)| = 1−O(ωh2

)

para 0 < ω 6 1, de modo que ω = 1 (metodo de Jacobi) oferece a melhor velocidade de convergencia,enquanto que ρ (Rh) > 1 para ω > 1 se h e suficientemente pequeno e o metodo nao converge.]

Para analisar as propriedades suavizadoras do metodo de Jacobi amortecido quantitativamente, intro-duzimos o fator suavizante de Rh:

Definicao. O fator suavizante µh (ω) de Rh e definido por

µh (ω) = max∣∣λkl

h (ω)∣∣ :

n

26 max (k, l) 6 n− 1

.

Definimos tambemµ∗ (ω) = sup

h∈Hµh (ω) ,

onde H = h = 1/n : n ∈ N e n > 4 denota o conjunto dos tamanhos de malha admissıveis.

Observe que µh (ω) e o maior autovalor dentre as maiores frequencias e representa o pior fator pelo qual oscomponentes de alta frequencia do erro sao reduzidos por passo de iteracao. Para entender isso, fixe umtamanho de malha h e escreva os autovalores de Rh (como na Secao 3.2.2) na forma

λ1 > λ2 > . . . > λq,

onde q = (n− 1)2, com ϕ1, . . . , ϕq sendo a correspondente base de autofuncoes. Escrevendo o erro inicialna forma

e0h =

q∑

i=1

αiϕi,

temos

ekh = Rk

he0h =

q∑

i=1

αiλki ϕi.

Como|λi|k → 0,


se |λi|k < 1, a taxa de eliminacao para o componente ϕi do erro e determinada por |λi|k e em cada iteracaoeste componente e reduzido por um fator exatamente igual a |λi|. Como

µh (ω) = max∣∣∣1− ω

2(2− cos kπh− cos lπh)

∣∣∣ :n

26 max (k, l) 6 n− 1

,

µ∗ (ω) = max∣∣∣1− ω

2

∣∣∣ , |1− 2ω|

,

segue que para 0 < ω < 1 o fator suavizante e menor que 1 e permanece longe de 1 por um limitanteindependente de h. Para ω = 1, o fator suavizante e da ordem de 1−O

(h2

)apenas; os menores autovalores

do metodo de Jacobiλkl =

12

(cos kπh + cos lπh)

estao associados as autofuncoes de frequencias medias, logo as autofuncoes de frequencias altas nao saorapidamente eliminadas e nao ha suavizacao. Por exemplo,

µh (ω) =

cosπh se ω = 1,

2 + cos πh

4se ω = 0.5,

1 + 2 cos πh

5se ω = 0.8,

µ∗ (ω) =

1 se ω = 1,

34

se ω = 0.5,

35

se ω = 0.8,

A escolha de ω = 0.8 e otima no sentido de que

inf0<ω61

µ∗ (ω) = µ∗ (0.8) = 3/5, (4.19)

enquanto que

inf0<ω61

µh (ω) = µh

(4

4 + cos πh

)=

3 cos πh

4 + cos πh=

35−

∣∣O (h2

)∣∣ . (4.20)

Isso significa que um passo do metodo de Jacobi amortecido com ω = 0.8 reduz todos os componentes doerro de alta frequencia por um fator de pelo menos 3/5, independente do tamanho h da malha.

4.4 O Ciclo de Duas Malhas

O segundo princıpio basico do metodo multigrid e a de que um termo de erro suave pode ser bem aproximadoem uma malha grosseira. Uma malha grosseira, por conter menos pontos, necessita de menos operacoes paraexecutar esta aproximacao (ela e muito mais barata que uma malha refinada). Introduzimos o ciclo de duasmalhas, que e a base para qualquer algoritmo de multigrid.

Enquanto que o erro na iteracao m e dado por

emh = uh − um

h ,

o resıduo (ou defeito) e definido porrmh = fh − Lhum

h . (4.21)

A equacao discretizada original Lhuh = fh e equivalente a equacao do resıduo

Lhemh = rm

h . (4.22)

Para transferir funcoes definidas em uma malha mais refinada Ωh para funcoes definidas em uma malha maisgrosseira Ω2h e vice-versa, precisamos definir dois operadores lineares de transferencia: um operador derestricao

I2hh : G (Ωh) −→ G (Ω2h) (4.23)


e um operador de interpolacao (ou de prolongamento)

Ih2h : G (Ω2h) −→ G (Ωh) . (4.24)

O operador de restricao sera usado para restringir o resıduo rmh obtido na malha mais refinada Ωh para a

malha mais grosseira Ω2h onde ele sera corrigido:

rm2h = I2h

h rmh . (4.25)

O operador de interpolacao sera usado para estender a correcao em2h obtida na malha mais grosseira Ω2h ate

a malha mais refinada Ωh:emh = Ih

2hem2h. (4.26)

4.1 Exemplo. Um operador de restricao particularmente simples de implementar computacionalmente e ooperador de injecao, definido por

(I2hh vh

)(x, y) = vh (x, y) para todo (x, y) ∈ Ω2h. (4.27)

Outro operador frequentemente usado e o operador peso total, que em notacao estencil e dado por

116

1 2 12 4 21 2 1

,

ou seja,

(I2hh vh

)(x, y) =

116

[4vh (x, y) + 2vh (x, y − h) + 2vh (x− h, y) + 2vh (x + h, y) + 2vh (x, y + h)

+vh (x− h, y − h) + vh (x + h, y − h) + vh (x− h, y + h) + vh (x + h, y + h)] .

Um terceiro operador de restricao e o operador metade peso:

18

0 1 01 4 10 1 0

.

¤

4.2 Exemplo. Um dos operadores de interpolacao mais simples de implementar e o operador de interpolacaobilinear :

(Ih2hv2h

)(x, y) =

v2h (x, y) se (x, y) = (2kh, 2lh) ,

12

(v2h (x, y − h) + v2h (x, y + h)) se (x, y) = (2kh, (2l − 1) h) ,

12

(v2h (x− h, y) + v2h (x, y + h)) se (x, y) = ((2k − 1)h, 2lh) ,

14

[vh (x− h, y − h) + vh (x + h, y − h)

+ vh (x− h, y + h) + vh (x + h, y + h)]se (x, y) = ((2k − 1)h, (2l − 1)h) .

para 1 6 k, l 6 n. Em notacao estencil, ele e denotado por

14

1 2 12 4 21 2 1

¤


Cada passo de iteracao (ciclo) de um metodo de duas malhas pode ser resumido no algoritmo seguinte(adaptado de [TOS]):

Ciclo de 2 Malhas

1. Pre-suavizacao

a) Calcule umh atraves de n1 passos de um suavizador aplicado a um

h :

umh = SUAVIZEn1 (um

h , Lh, fh).

2. Correcao na malha grosseira

a) Calcule o resıduo rmh = fh − Lhum

h .

b) Restrinja o resıduo a malha mais grosseira: rm2h = I2h

h rmh .

c) Calcule o erro na malha mais grosseira: L2hem2h = rm

2h.

d) Interpole a correcao para a malha mais refinada: emh = Ih

2hem2h.

e) Calcule a aproximacao corrigida: umh = um

h + emh .

3. Pos-suavizacao

a) Calcule um+1h atraves de n2 passos de um suavizador aplicado a um

h :

um+1h = SUAVIZEn2 (um

h , Lh, fh).

A necessidade da pos-suavizacao deve-se ao fato que as frequencias mais baixas na malha mais grosseiracorrespondem nao somente as frequencias mais baixas na malha mais refinada, como tambem as frequenciasmais altas, como vimos em (4.8) (em outras palavras, frequencias baixas em Ω2h sao mapeadas para a mesmafrequencia baixa em Ωh e para tres outras frequencias altas em Ωh); para evitar que estas componentesdo erro reaparecam, fazemos uma segunda suavizacao. Observe que varios componentes individuais dociclo de duas malhas devem ser especificados, e sua escolha pode ter uma forte influencia na eficiencia doalgoritmo: o procedimento suavizador SUAVIZE (um

h , Lh, fh); os numeros n1 e n2 de passos de suavizacao,a malha grosseira (aqui escolhemos Ω2h, mas outras escolhas sao possıveis) e os operadores de restricao e deinterpolacao.

4.5 O Ciclo Multigrid: Ciclos V

O ciclo de duas malhas per si e obviamente de pouco significado pratico, ja que o custo computacional namalha Ω2h ainda e relativamente alto. A ideia de um ciclo multigrid e nao resolver a equacao de correcaodo resıduo L2hem

2h = rm2h exatamente, mas suaviza-la e transferir o problema para uma malha ainda mais

grosseira Ω4h, onde o custo computacional e ainda menor. Esta ideia e repetida ate que podemos em princıpiochegar na malha Ω1, onde a correcao do resıduo pode entao ser calculada exatamente. Daı, voltamos paraa malha mais refinada original, formando um ciclo no formato da letra “V”.

Capıtulo 5

Metodo dos Volumes Finitos

A discretizacao do domınio no metodos dos volumes finitos difere da do metodo de diferencas finitas em quenesta o domınio e substituıdo por um conjunto de pontos, enquanto que na primeira o domınio e subdivididoem volumes de controle ou celulas. Os pontos nodais ou simplesmente nos, sao os centros das celulas.No metodo dos volumes finitos, ao inves de aproximarmos diretamente a equacao diferencial como no metodode diferencas finitas, ela e antes integrada sobre cada volume de controle. As integrais obtidas sao entaoaproximadas. As equacoes integrais estao na forma de leis de conservacao, o que assegura a conservacaodas grandezas fısicas tratadas em cada volume de controle (conservacao no nıvel discreto) e portanto estemetodo e bastante adequado para tratar de fenomenos fısicos que envolvem leis de conservacao. Muitasvezes pode-se trabalhar diretamente com as equacoes integrais, sem passar pelas equacoes diferenciais, o quetorna o metodo particularmente util para tratar de fenomenos descontınuos melhor modelados por equacoesintegrais, tais como fenomenos que envolvem ondas de choque.

5.1 Leis de Conservacao

Muitas das equacoes diferenciais parciais fundamentais sao obtidas atraves de leis de conservacao.Leis de conservacao sao essencialmente leis de balanceamento, expressando o fato de que alguma substancia

e balanceada. Aqui, o termo substancia pode indicar uma substancia realmente material, ou ate mesmo umconceito abstrato, tal como energia ou uma populacao de animais. Por exemplo, a primeira lei da ter-modinamica e a lei de conservacao da energia: a variacao de energia interna de um sistema e igual ao calortotal adicionado ao sistema mais o trabalho realizado sobre o sistema. Como outro exemplo, considere umfluido escoando em alguma regiao do espaco, consistindo de substancias sofrendo reacoes quımicas: paracada substancia quımica individual, a taxa de variacao da quantidade total da substancia na regiao e iguala taxa com que a substancia flui para dentro da regiao, menos a taxa com que ela flui para fora da regiao,mais a taxa com que ela e criada, ou consumida, pelas reacoes quımicas. Como ultimo exemplo, a taxa devariacao de uma dada populacao de animais em uma regiao e igual a taxa de nascimentos, menos a taxa demortes, mais a taxa de migracao para dentro ou fora da regiao.

Matematicamente, leis de conservacao traduzem-se em equacoes integrais, de onde podem ser deduzidasequacoes diferenciais, na maior parte dos casos. Estas equacoes descrevem como o processo evolui com otempo. Por este motivo, elas sao tambem chamadas de equacoes de evolucao. Vamos examinar primeiroo caso unidimensional.

5.1.1 Lei de Conservacao Unidimensional

Seja u = u(x, t) a densidade ou concentracao de alguma substancia, por unidade de volume, que dependeapenas de uma variavel espacial x ∈ R e do tempo t > 0. Novamente enfatizamos que a substancia cujadensidade estamos medindo pode ser massa, momento, energia, populacao, ou qualquer outra coisa, material

94


ou abstrata. Por exemplo, no caso da equacao do calor, a temperatura u e uma medida da densidade deenergia termica. De fato, se e(x, t) denota a densidade de energia termica, isto e, a quantidade de energiatermica por unidade de volume, entao a densidade de energia termica e a temperatura estao relacionadasatraves da equacao

e(x, t) = c(x)ρ(x)u(x, t),

cujo significado e: a energia termica por unidade de volume e igual a energia termica por unidade de massapor unidade de temperatura (i.e., o calor especıfico), vezes a temperatura, vezes a densidade volumetrica demassa.

Imaginamos que a substancia esta distribuıda em um tubo uniforme com secao transversal de areaconstante A. Por hipotese, u e constante em cada secao transversal do tubo, variando apenas na direcao x.Considere um segmento arbitrario do tubo, entre as secoes transversais localizadas em x = a e em x = b.Chamamos este segmento de volume de controle. A quantidade total da substancia dentro do volume decontrole no instante de tempo t e

Quantidade total da substanciadentro do volume de controle =

∫ b

a

u(x, t)Adx.

Assuma agora que existe movimento da substancia atraves do tubo na direcao axial. Definimos o fluxoφ(x, t) da substancia no tempo t como sendo a quantidade da substancia fluindo atraves da secao transversalem x no tempo t por unidade de area, por unidade de tempo. Assim as dimensoes de φ sao [φ] = quantidadeda substancia / (area × tempo). Por convencao, φ sera positivo se a substancia estiver se movendo na direcaopositiva do eixo x, e negativo se ela estiver se movendo na direcao negativa do eixo x. Portanto, no tempo t,a quantidade lıquida de substancia permanecendo no volume de controle sera a diferenca entre a quantidadeda substancia entrando em x = a e a quantidade da substancia saindo em x = b:

Taxa de transferencia lıquida da substanciapara dentro do volume de controle = φ(a, t)A− φ(b, t)A.

A substancia pode ser criada ou destruıda dentro do volume de controle por uma fonte interna ou externa.A taxa de criacao ou destruicao da substancia, que chamaremos de termo fonte e denotaremos por f(x, t, u),tem dimensoes [f ] = quantidade da substancia / (volume × tempo), tendo sinal positivo se a substancia ecriada dentro do volume de controle e negativa se a substancia for destruıda dentro do volume de controle.Observe que ela pode depender da propria quantidade da substancia disponıvel, medida pela densidade u.A taxa de criacao ou destruicao da substancia dentro do volume de controle e entao dada por

Taxa de criacao da substanciadentro do volume de controle =

∫ b

a

f(x, t, u)Adx.

A lei de conservacao para a substancia pode ser formulada da seguinte forma:

Taxa de variacao

da quantidade de substancia

dentro do volume de controle

=Taxa de transferencia lıquida de substancia

para dentro do volume de controle

atraves de sua fronteira

+Taxa de criacao da substancia

dentro do volume de controle

ou, em termos matematicos, apos cancelar o termo comum A,

d

dt

∫ b

a

u(x, t) dx = φ(a, t)− φ(b, t) +∫ b

a

f(x, t, u) dx. (5.1)

Esta e a lei de conservacao na forma integral, valendo mesmo se u, φ ou f nao forem funcoes diferenciaveis(o que pode ocorrer em certos fenomenos fısicos, como por exemplo naqueles que envolvem ondas de choque


ou outros tipos de descontinuidade). Se estas funcoes forem continuamente diferenciaveis, podemos derivarsob o sinal de integracao na primeira integral

d

dt

∫ b

a

u(x, t) dx =∫ b

a

ut(x, t) dx,

e usar o Teorema Fundamental do Calculo para escrever

φ(a, t)− φ(b, t) = −∫ b

a

φx(x, t) dx,

obtendo a equacao diferencial parcialut + φx = f(x, t, u) (5.2)

que e a lei de conservacao na forma diferencial.

5.1.2 Lei de Conservacao em Varias Dimensoes

Vamos formular a lei de conservacao nas formas integral e diferencial para os espacos Rn, n = 2 ou n = 3(na verdade, tudo o que deduzirmos aqui, vale para qualquer n > 2). Considere um volume de controle V emRn, em que a densidade ou concentracao u = u(x, t) de alguma substancia por unidade de volume dependede n variaveis espaciais x = (x1, . . . , xn) e do tempo t > 0. Temos

Quantidade total da substanciadentro do volume de controle =

∫

V

u(x, t) dV

e, se f(x, t, u) denota o termo fonte,

Taxa de criacao da substanciadentro do volume de controle =

∫

V

f(x, t, u) dV.

Em n dimensoes, o fluxo pode ser em qualquer direcao, logo ele e uma grandeza vetorial que denotaremospor φ(x, t). Se η(x) denota o vetor unitario normal apontando para fora da regiao V , a taxa de transferencialıquida da substancia para fora do volume de controle atraves de sua fronteira ∂V e dada por

Taxa de transferencia lıquida da substanciapara fora do volume de controle =

∫

∂V

φ(x, t) · η(x) dS.

A lei de conservacao e, portanto,

d

dt

∫

V

u(x, t) dV = −∫

∂V

φ(x, t) · η(x) dS +∫

V

f(x, t, u) dV. (5.3)

Se u, φ e f forem todas de classe C1 (assim como a regiao V ), podemos derivar sob o sinal de integracao eusar o Teorema da Divergencia

∫

∂V

φ(x, t) · η(x) dS =∫

V

div φ(x, t) dV,

para obter a lei de conservacao em forma diferencial

ut + div φ = f(x, t, u). (5.4)


5.1.3 Relacoes Constitutivas

A lei de conservacao na forma diferencial e uma equacao diferencial parcial em duas incognitas, u e φ.Precisamos, portanto, de uma segunda equacao para obter um sistema bem determinado. A equacao adicionale frequentemente baseada nas propriedades fısicas do meio, as quais frequentemente decorrem de observacoesempıricas. Tais equacoes sao chamadas de relacoes constitutivas ou equacoes de estado.

4.1 Exemplo. (Equacao do Calor) No caso da equacao do calor, a relacao constitutiva e a lei de Fourier:

φ(x, t) = −k (x) ux(x, t)

onde a constante de condutividade termica k = k (x) depende do material e muitas vezes pode serconsiderada constante.

Em dimensoes mais altas, a lei de Fourier assume a forma

φ(x, t) = −k (x)∇u(x, t).

De fato, para materiais isotropicos (isto e, materiais em que nao existem direcoes preferenciais) verifica-se experimentalmente que o calor flui de pontos quentes para pontos frios na direcao em que a diferencade temperatura e a maior. O fluxo de calor e proporcional a taxa de variacao da temperatura nestadirecao, com a constante de proporcionalidade k sendo por definicao a condutividade termica, comono caso unidimensional. Como sabemos, a direcao onde uma funcao cresce mais rapido e exatamenteaquela dada pelo vetor gradiente da funcao, e o modulo do gradiente fornece a magnitude da taxade variacao da funcao nesta direcao. O sinal negativo ocorre, como no caso unidimensional, porque ovetor gradiente aponta na direcao de crescimento da temperatura, enquanto que o fluxo do calor se dana direcao oposta (da temperatura maior para a temperatura menor). O fluxo do calor em uma regiaobi ou tridimensional pode ser facilmente visualizado quando se lembra que o gradiente de uma funcao eperpendicular as superfıcies de nıvel da funcao. No caso em que a funcao e a temperatura, as superfıciesde nıvel sao chamadas superfıcies isotermicas ou, simplesmente, isotermas. Assim, o calor flui dasisotermas mais quentes para as isotermas mais frias, e em cada ponto da isoterma perpendicularmentea isoterma. Em outras palavras, as linhas de corrente do fluxo de calor correspondem as linhas de fluxodo campo gradiente da temperatura.

Substituindo a relacao constitutiva na lei de conservacao, obtemos a equacao do calor: na formadivergente,

ut = div (k∇u) + f(x, t, u),

ou, quando k e constante, na forma usual envolvendo o laplaciano,

ut = k∆u + f(x, t, u).

¤

4.2 Exemplo. (Equacao da Difusao) Em muitos outros processos fısicos observa-se que a substancia fluia uma taxa diretamente proporcional ao gradiente de densidade, de regioes de maior densidade pararegioes de menor densidade. Esta relacao geral e chamada de lei de Fick :

φ(x, t) = −D (x)∇u(x, t),

onde D = D (x) e a constante de difusao. Assumindo D constante, se o termo fonte independe de u,obtemos a equacao da difusao

ut = D∆u + f(x, t),

caso contrario a equacao diferencial parcial obtida e chamada equacao da difusao-reacao

ut = D∆u + f(x, t, u),


que aparece na teoria de combustao e em biologia. Se D nao e constante, obtemos as respectivasequacoes na forma divergente. O nome difusao vem do fato de que a substancia difunde-se para regioesadjacentes por causa de gradientes (i.e., diferencas) de concentracao, e nao porque e transportada pelacorrente (i.e., nao atraves de conveccao). Por este motivo, o termo D∆u e chamado de termo difusivo.

Alem do calor, exemplos de outras substancias que se comportam assim sao substancias quımicasdissolvidas em algum fluido (neste caso, u representa a concentracao quımica) e ate mesmo populacoesde insetos. Alem de ser confirmada atraves de observacoes empıricas, a lei de Fick que governa estese varios outros fenomenos fısicos e biologicos pode ser justificada teoricamente atraves de argumentosbaseados em modelos probabilısticos e caminhos aleatorios. ¤

Neste texto sobre equacoes elıpticas, obviamente estamos interessados na equacao de estado estacionarioresultante da equacao da difusao ou de difusao-reacao, isto e, no caso em que ut = 0:

−∆u = f(x, t, u), (5.5)

ou, na forma divergente,− div (A (x)∇u) = f(x, t, u), (5.6)

onde no caso mais geral A e uma matriz n× n. Em termos da lei de conservacao, isto se escreve na forma

−∫

∂V

A (x)∇u(x, t) · η(x) dS =∫

V

f(x, t, u) dV. (5.7)

5.2 O Caso Unidimensional

Consideramos a seguinte equacao elıptica na forma divergente com condicao de Dirichlet:− d

dx

(a (x)

du

dx

)= f (x, u) em [0, L] ,

u (0) = u0, u (L) = uL.(5.8)

O primeiro passo e gerar a malha de volumes finitos no intervalo [0, L], isto e, discretizar o domınioem volumes de controle. Para isso, inserimos um numero n de pontos nodais ou nos P1, . . . , Pn entreos pontos 0 e L da fronteira do domınio. Os n volumes de controle V1, . . . , Vn serao centrados nestes nos.As faces (fronteiras) dos volumes de controle serao posicionadas no ponto medio entre dois nos. Em geral,posiciona-se os volumes de controle de modo que as fronteiras do domınio coincidem com faces dos volumesde controle, isto e, o ponto 0 esta na face esquerda do primeiro volume de controle e o ponto L esta na facedireita do ultimo volume de controle. Para simplificar a apresentacao, assumiremos que os pontos nodaisforam posicionados de modo a estarem igualmente espacados, de modo que todos os volumes de controletem mesma largura igual a ∆x.

Estabelecemos a seguinte notacao (esta convencao e frequentemente utilizada em dinamica dos fluidoscomputacional, onde o metodo dos volumes finitos e bastante popular): um ponto nodal arbitrario seradesignado simplesmente por P e os seus pontos nodais vizinhos serao designados por W (oeste, isto e, oponto nodal vizinho a esquerda) e E (leste, correspondendo ao vizinho a direita). A face esquerda (a oeste)do volume de controle sera designada por w e a face direita (a leste) por e. Assim, a distancia entre dois nosvizinhos, assim como a distancia entre as duas faces de um volume de controle e igual a ∆x.

Uma vez discretizado o domınio com a geracao da malha de volumes de controle, integrando a equacaodiferencial parcial em cada volume de controle para coloca-la na forma integral (reobtendo a lei de con-servacao; e claro que podemos desde o inıcio trabalhar diretamente com esta, se estiver disponıvel):

−∫

Vp

[d

dx

(a (x)

du

dx

)]dx =

∫

Vp

f (x, u) dx.


Segue pelo teorema fundamental do calculo que

a (xw)du

dx(xw)− a (xe)

du

dx(xe) = fV ∆x (5.9)

onde fV denota o valor medio de f sobre o volume de controle, isto e,

fVP=

1∆x

∫

Vp

f (x, u) dx.

Observe que a equacao integral obtida e uma equacao exata, ainda nao discretizada. Na linguagem de leisde conservacao, ela diz simplesmente que o fluxo de u deixando a face direita do volume de controle menoso fluxo deixando a face esquerda do mesmo (respeitando a nossa convencao de sinal para fluxos) e igual aquantidade de u gerada pela fonte dentro do volume de controle:

φw − φe = fVP∆x.

Agora procedemos a discretizacao da equacao integral. Valores nas faces devem ser dados em funcoes devalores nos pontos nodais. Consideremos primeiro os volumes de controle interiores V2, . . . , Vn−1. Usandointerpolacao linear, podemos obter valores aproximados para a (xe) , a (xw), calculados nas faces dos volumesde controle, em termos dos valores de a nos pontos nodais dos volumes de controle:

aw := a (xw) =aW + aP

2, (5.10)

ae := a (xe) =aP + aE

2. (5.11)

As derivadas primeiras, ou seja, os fluxos, podem ser aproximadas atraves de diferencas finitas apropriadas,por exemplo diferencas finitas centradas:

du

dx

∣∣∣∣w

:=du

dx(xw) =

uP − uW

∆x, (5.12)

du

dx

∣∣∣∣e

:=du

dx(xe) =

uE − uP

∆x. (5.13)

O termo fonte, que pode expressar uma dependencia nao linear do valor de u, pode ser linearizado e assumidoconstante ao longo do volume de controle, produzindo

fVP=

1∆x

∫

Vp

(f0

P + f1P up

)dx =

f0P + f1

P up

∆x

∫

Vp

dx = f0P + f1

P up. (5.14)

(Como queremos obter um sistema linear no final, nao e possıvel aproximar o termo fonte por uma aprox-imacao de ordem maior que 1. A linearizacao do termo linear sera discutida em maiores detalhes na secao4 deste capıtulo) Daı,

awuP − uW

∆x− ae

uE − uP

∆x=

(f0

P + f1P up

)∆x,

ouapuP + aW uW + aEuE = bp, (5.15)

onde

ap =aw

∆x2+

ae

∆x2− f1

P , (5.16)

aW = − aw

∆x2, aE = − ae

∆x2, (5.17)

bp = f0P . (5.18)


O tratamento dos volumes de controle adjacentes a fronteira e ligeiramente diferente. Para o volume decontrole V1 adjacente a fronteira esquerda (oeste) do domınio, temos

aw = a0, (5.19)

edu

dx

∣∣∣∣w

=uP − u0

∆x/2, (5.20)

porque a distancia entre P e 0 e apenas ∆x/2; neste caso somos forcados a utilizar uma diferenca finitaprogressiva para aproximar a derivada primeira em w. Assim, a equacao discretizada correspondente a estevolume de controle e

2a0uP − u0

∆x− ae

uE − uP

∆x=

(f0

P + f1P up

)∆x,

ouapuP + aEuE = bp, (5.21)

onde

ap =2a0

∆x2+

ae

∆x2− f1

P , (5.22)

aE = − ae

∆x2, (5.23)

bp = f0P +

2a0

∆x2u0. (5.24)

Para o volume de controle Vn adjacente a fronteira direita temos

ae = aL,

du

dx

∣∣∣∣e

=uL − uP

∆x/2,

utilizando uma diferenca finita regressiva para aproximar a derivada primeira em e, de modo que a equacaodiscretizada correspondente a este volume de controle e

awuP − uW

∆x− 2ae

uL − uP

∆x=

(f0

P + f1P up

)∆x,

ouapuP + aEuE = bp, (5.25)

onde

ap =aw

∆x2+

2aL

∆x2− f1

P , (5.26)

aW = − aw

∆x2, (5.27)

bp = f0P +

2aL

∆x2uL. (5.28)

Ordenando os volumes de controle (geralmente da esquerda para a direita), obtemos um sistema linear cujasolucao sera uma solucao aproximada para a equacao com as condicoes de fronteira dadas.

4.3 Exemplo. (Equacao de Poisson) Vamos aplicar o metodo de volumes finitos a equacao de Poisson comcondicao de fronteira de Dirichlet

−u′′ = f (x) em [0, L] ,u (0) = u0, u (L) = uL.

(5.29)


Aqui a (x) ≡ 1 e f (x, u) = f (x), de modo que f1P = 0. Se decidimos aproximar o valor medio de f no

volume de controle pelo valor de f em P , segue que

ap =2

∆x2, aW = − 1

∆x2, aE = − 1

∆x2, bp = fP

para os volumes de controle interiores V2, . . . , Vn−1. Para os volumes de controle adjacentes a fronteira,para o primeiro volume de controle V1 temos

ap =3

∆x2, aE = − 1

∆x2, bp = fP +

2∆x2

u0,

enquanto que para o ultimo volume de controle Vn temos

ap =3

∆x2, aW = − 1

∆x2, bp = fP +

2∆x2

uL.

O sistema discretizado e, portanto:

1∆x2

3 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 3

u1

u2

...

...un−1

un

=

f1 +2

∆x2u0

f2

...

...fn−1

fn +2

∆x2uL

.

Compare com o correspondente sistema discretizado obtido pelo metodo de diferencas finitas; a unicadiferenca esta na primeira e ultima linhas dos sistemas. ¤

4.4 Exemplo. (Equacao Elıptica Linear) Consideremos agora o seguinte problema linear elıptico comcondicao de fronteira de Dirichlet −u′′ = Au + B em [0, L] ,

u (0) = u0, u (L) = uL.(5.30)

Novamente a (x) ≡ 1, mas desta vez f (x, u) = f (u) = Au + B, de modo que f0P = B e f1

P = A. Segueque

ap =2

∆x2−A, aW = − 1

∆x2, aE = − 1

∆x2, bp = B

para os volumes de controle interiores V2, . . . , Vn−1. Para os volumes de controle adjacentes a fronteira,para o primeiro volume de controle V1 temos

ap =3

∆x2−A, aE = − 1

∆x2, bp = B +

2∆x2

u0,

enquanto que para o ultimo volume de controle Vn temos

ap =3

∆x2−A, aW = − 1

∆x2, bp = B +

2∆x2

uL.

O sistema discretizado e, portanto:

1∆x2

3−A −1−1 2−A −1

−1. . . . . .. . . . . . −1

−1 2−A −1−1 3−A

u1

u2

...

...un−1

un

=

B +2

∆x2u0

B......B

B +2

∆x2uL

.


Como e sabido, podemos assegurar que o problema linear elıptico possui solucao unica se A 6 0,utilizando o princıpio do maximo. Isso se traduz do ponto de vista numerico, no fato de que a matrizdiscretizada permanece diagonalmente dominante. No caso em que A > 0 e preciso ter cuidado, poispode haver infinitas solucoes exatas e nao existir solucao numerica e vice-versa, pois os autovaloresdo problema exato nao sao iguais aos autovalores da matriz de discretizacao (na maioria dos casosestes ultimos nao sao nem boas aproximacoes para os primeiros: usualmente as aproximacoes saorazoavelmente boas apenas para os primeiros autovalores e em malhas bastante refinadas, com umnumero enorme de pontos ou celulas). Para evitar este tipo de problema, e possıvel modificar alinearizacao; veja a secao 4 deste capıtulo. ¤

5.3 O Caso Bidimensional

Considere agora a seguinte equacao elıptica na forma divergente com condicao de Dirichlet em um domınioretangular: − div [A (x, y)∇u] = f (x, y, u) em Ω = [0, 1]× [0, 1] ,

u (x, y) = g (x, y) . sobre ∂Ω,(5.31)

onde A (x, y) e uma matriz 2× 2. Vamos considerar o caso mais simples em que A (x, y) = a (x, y) I.No caso bidimensional, os quatro pontos nodais vizinhos de um ponto nodal arbitrario P serao designados

por W (oeste), E (leste), S (sul) e N (norte), e as faces correspondentes do volume de controle por w, e,s e n. A distancia horizontal entre dois nos vizinhos (que e a largura de um volume de controle) sera iguala ∆x, enquanto que a distancia vertical entre dois nos vizinhos (altura do volume de controle) sera igual a∆y.

Integrando a equacao diferencial parcial em cada volume de controle, como

div [a (x, y)∇u] =∂

∂x

(a (x, y)

∂u

∂x

)+

∂

∂y

(a (x, y)

∂u

∂y

),

obtemos agora

−∫

Vp

[∂

∂x

(a (x, y)

∂u

∂x

)]dxdy −

∫

Vp

[∂

∂y

(a (x, y)

∂u

∂y

)]dxdy =

∫

Vp

f (x, y, u) dxdy,

ou

−∫ n

s

(∫ e

w

[∂

∂x

(a (x, y)

∂u

∂x

)]dx

)dy −

∫ e

w

(∫ n

s

[∂

∂y

(a (x, y)

∂u

∂y

)]dy

)dx =

∫

Vp

f (x, y, u) dxdy.

Atraves do teorema fundamental do calculo obtemos a equacao exata

−∫ n

s

[a (xe, y)

∂u

∂x(xe, y)− a (xw, y)

∂u

∂x(xw, y)

]dy −

∫ e

w

[a (x, yn)

∂u

∂y(x, yn)− a (x, ys)

∂u

∂y(x, ys)

]dx

=∫

Vp

f (x, y, u) dxdy.

Para continuar o processo de integracao, precisamos aproximar as integrais. Escolhemos a aproximacao dointegrando pelo ponto medio do intervalo:

∫ n

s

[a (xe, y)

∂u

∂x(xe, y)− a (xw, y)

∂u

∂x(xw, y)

]dy ≈

[a (xe, yp)

∂u

∂x(xe, yp)− a (xw, yp)

∂u

∂x(xw, yp)

]∆y,

∫ e

w

[a (x, yn)

∂u

∂y(x, yn)− a (x, ys)

∂u

∂y(x, ys)

]dx ≈

[a (xp, yn)

∂u

∂y(xp, yn)− a (xp, ys)

∂u

∂y(xp, ys)

]∆x


Obtemos, portanto, a seguinte equacao parcialmente discretizada (diferente do caso unidimensional, estaequacao nao e exata):

a (xw, yp)∂u

∂x(xw, yp)∆y − a (xe, yp)

∂u

∂x(xe, yp)∆y + a (xp, ys)

∂u

∂y(xp, ys)∆x− a (xp, yn)

∂u

∂y(xp, yn)∆x

= fV ∆x∆y,

ondefVP

=1

∆x∆y

∫

Vp

f (x, u) dxdy.

Em termos de fluxos discretizados,

φw − φe + φs − φn = fVP∆x∆y.

Usando interpolacao linear como antes, obtemos valores aproximados para a (xw) , a (xe) , a (xs) , a (xn), cal-culados nas faces dos volumes de controle, em termos dos valores de a nos pontos nodais dos volumes decontrole:

aw := a (xw, yp) =aW + aP

2, (5.32)

ae := a (xe, yp) =aP + aE

2, (5.33)

as := a (xp, xs) =aS + aP

2, (5.34)

an := a (xp, xn) =aP + aN

2. (5.35)

Os fluxos sao aproximadas atraves de diferencas finitas centradas:

∂u

∂x

∣∣∣∣w

:=∂u

∂x(xw, yp) =

uP − uW

∆x, (5.36)

∂u

∂x

∣∣∣∣e

:=∂u

∂x(xe, yp) =

uE − uP

∆x, (5.37)

∂u

∂y

∣∣∣∣s

:=∂u

∂y(xp, ys) =

uP − uS

∆y, (5.38)

∂u

∂y

∣∣∣∣n

:=∂u

∂y(xp, yn) =

uN − uP

∆y. (5.39)

O termo fonte e linearizado

fVP=

1∆x∆y

∫

Vp

(f0

P + f1P up

)dxdy =

f0P + f1

P up

∆x∆y

∫

Vp

dxdy = f0P + f1

P up. (5.40)

Daı,

awuP − uW

∆x∆y − ae

uE − uP

∆x∆y + as

uP − uS

∆y∆x− an

uN − uP

∆y∆x =

(f0

P + f1P up

)∆x∆y,

ouapuP + aW uW + aEuE + aSuS + aNuN = bp. (5.41)

com

ap =aw

∆x2+

ae

∆x2+

as

∆y2+

an

∆y2− f1

P , (5.42)

aW = − aw

∆x2, aE = − ae

∆x2, aS = − as

∆y2, aN = − an

∆y2, (5.43)

bp = f0P . (5.44)


O tratamento dos volumes de controle adjacentes a fronteira e diferente. Por exemplo, para volumes decontrole adjacentes a fronteira esquerda (oeste), que nao sejam os dois volumes de controle dos cantos, temos

aw = a (0, yp) , (5.45)

e∂u

∂x

∣∣∣∣w

=uP − u (0, yp)

∆x/2, (5.46)

porque a distancia horizontal entre P e 0 e ∆x/2. Assim, a equacao discretizada correspondente a estevolume de controle e

2a (0, yp)uP − u (0, yp)

∆x∆y − ae

uE − uP

∆x∆y + as

uP − uS

∆y∆x− an

uN − uP

∆y∆x =

(f0

P + f1P up

)∆x∆y,

ouapuP + aEuE + aSuS + aNuN = bp, (5.47)

com

ap =2a (0, yp)

∆x2+

ae

∆x2+

as

∆y2+

an

∆y2− f1

P , (5.48)

aE = − ae

∆x2, aS = − as

∆y2, aN = − an

∆y2, (5.49)

bp = f0P +

2a (0, yp)∆x2

g (0, yp) . (5.50)

Formulas semelhantes sao obtidas para volumes de controle adjacentes as demais fronteiras que nao estejamem um dos quatro cantos do domınio retangular. Para os volumes de controle nos cantos do retangulo,precisamos fazer mais uma modificacao. Por exemplo, para o volume de controle no canto superior esquerdotemos

aw = a (0, yp) , (5.51)an = a (xp, 1) , (5.52)

e∂u

∂x

∣∣∣∣w

=uP − u (0, yp)

∆x/2, (5.53)

∂u

∂y

∣∣∣∣n

=u (xp, 1)− uP

∆y/2, (5.54)

e a equacao discretizada correspondente a este volume de controle e

2a (0, yp)uP − u (0, yp)

∆x∆y−ae

uE − uP

∆x∆y+as

uP − uS

∆y∆x−2a (xp, 1)

u (xp, 1)− uP

∆y∆x =

(f0

P + f1P up

)∆x∆y,

ouapuP + aEuE + aSuS = bp, (5.55)

com

ap =2a (0, yp)

∆x2+

ae

∆x2+

as

∆y2+

2a (xp, 1)∆y2

− f1P , (5.56)

aE = − ae

∆x2, aS = − as

∆y2, (5.57)

bp = f0P +

2a (0, yp)∆x2

g (0, yp) +2a (xp, 1)

∆x2g (xp, 1) . (5.58)

Ordenando os volumes de controle (por exemplo, usando a ordem lexicografica), obtemos um sistema linearcuja solucao sera uma solucao aproximada para a equacao com as condicoes de fronteira dadas.


4.5 Exemplo. (Equacao de Poisson) Vamos aplicar o metodo de volumes finitos a equacao de Poisson comcondicao de fronteira de Dirichlet

−∆u = f (x, y) em [0, 1]2 ,

u = g (x, y) sobre ∂ [0, 1]2 .(5.59)

Temos a (x) ≡ 1, f1P = 0, f0

P = fP , e optamos por discretizar a malha por volumes de controlequadrados, isto e, satisfazendo ∆x = ∆y. Segue que a linha do sistema discretizado corresponde a umvolume de controle interior tem a forma (multiplicamos todas as linhas do sistema por ∆x2)

elemento na diagonal: ap = 4,

elementos fora da diagonal: a∗ = −1 (4 elementos),

elemento constante: bP = fP ∆x2.

Para volumes de controle adjacentes a fronteira, nao localizados nos cantos, a linha correspondente nosistema discretizado e



elemento constante: bP = fP ∆x2 + 2g (∗) .

Finalmente, para volumes de controle localizados nos cantos, temos



elemento constante: bP = fP ∆x2 + 2g (∗) + 2g (∗∗) .

Compare com o correspondente sistema discretizado obtido pelo metodo de diferencas finitas; como nocaso unidimensional, as diferencas surgem apenas para as linhas correspondentes a celulas e pontos nafronteira do domınio. ¤

5.4 Linearizacao do Termo Fonte

Ao linearizar o termo fontef (u) = f0

P + f1P up

devemos ter cuidado para esolher a linearizacao de tal forma a obter

f1P 6 0. (5.60)

A necessidade matematica desta escolha ja foi discutida no Exemplo 4.4. Fisicamente, esta exigencia tambemfaz sentido: a maioria dos termos fontes em fenomenos transientes que tendem a um estado estacionario emgeral tem derivada primeira negativa, caso contrario o sistema nao tenderia a um regime permanente. Porexemplo, na difusao do calor, a existencia de um termo linear com derivada positiva implicaria na acumulacaode energia termica dentro do domınio, a nao ser que o calor pudesse ser rapidamente dissipado atraves dafronteira, o que geral nao ocorre, pois mesmo o calor perdido por um objeto quente atraves da sua imersaoem um recipiente cheio de lıquido frio e transferido para o lıquido a uma taxa linear. Isso tende a geraruma situacao instavel que eventualmente leva ao colapso termico do sistema (explosao ou derretimento doobjeto).


5.4.1 Termo Fonte do Tipo f (u) = Au + B com A < 0

A linearizacao obvia neste caso e tomarf0

P = B, f1P = A, (5.61)

o que aumentara a dominancia diagonal da matriz, como ja vimos no Exemplo 4.4. Outra possibilidade eusar um processo iterativo, definindo

f0P = Auk−1

P + B e f1P = 0, (5.62)

usando o valor Auk−1P +B obtido na iteracao anterior no lado direito do sistema a ser resolvido nesta iteracao.

Como o termo fonte neste caso e linear, a primeira sugestao e mais aconselhada neste caso.

5.4.2 Termo Fonte do Tipo f (u) = Au + B com A > 0

Neste caso, como ja mencionado varias vezes, nao e aconselhavel tomar a linearizacao obvia (especialmentese existirem outras nao-linearidades, e um processo nao-iterativo se fazer necessario para resolver o sistema,isso pode levar o processo iterativo a divergir). A segunda sugestao da subsecao anterior e a mais adequadaneste caso, isto e, tomar

f0P = Auk−1

P + B e f1P = 0, (5.63)

e usar um processo iterativo.

5.4.3 Termo Fonte do Tipo f (u) com f ′ (u) < 0

A maneira mais simples de lidar com um termo fonte nao-linear e usar um processo iterativo simples,definindo

f0P = f

(uk−1

P

)e f1

P = 0. (5.64)

A desvantagem deste metodo e que ele nao toma conhecimento da dependencia de f em u na iteracaocorrente. Uma linearizacao que leva isto em conta e a seguinte: escrevendo

f(uk

P

)= f

(uk−1

P

)+

df

du

(uk−1

P

) (uk

P − uk−1P

), (5.65)

tomamosf0

P = f(uk−1

P

)− df

du

(uk−1

P

)uk−1

P e f1P =

df

du

(uk−1

P

). (5.66)

Por exemplo, se f (u) = 4− 5u3, terıamos

f(uk

P

)= 4− 5

(uk−1

P

)3 − 15(uk−1

P

)2 (uk

P − uk−1P

)

= 4 + 10(uk−1

P

)3 − 15(uk−1

P

)2uk

P .

Referencias Bibliograficas

[Asmar] Nakhle ASMAR, Partial Differential Equations and Boundary Value Problems, Pren-tice Hall, 2000.

[Biezuner] Rodney Josue BIEZUNER, Notas de Aula: Equacoes Diferenciais Parciais, UFMG,2005.

[BHM] William L. BRIGGS, Van Emden HENSON e Steve F. McCORMICK, A MultigridTutorial, SIAM, 2000.

[Demmel] James W. DEMMEL, Applied Numerical Linear Algebra, SIAM, 1997.

[Hackbusch] W. HACKBUSCH, Elliptic Differential Equations: Theory and Numerical Treatment,Springer Series in Computational Mathematics 18, Springer, 1992.

[Heuveline] Vincent HEUVELINE, On the computation of a very large number of eigenvalues forselfadjoint elliptic operators by means of multigrid methods, Journal of ComputationalPhysics 184 (2003), 321–337.

[Horn-Johnson] Roger A. HORN e Charles R. JOHNSON, Matrix Analysis, Cambridge UniversityPress, 1985.

[Maliska] CLOVIS R. MALISKA, Transferencia de Calor e Mecanica dos Fluidos Computa-cional, 2a. Edicao, LTC, 2004.

[Patankar] S. V. PATANKAR, Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980.

[Rosser1] J. Barkley ROSSER, Nine point difference solutions for Poisson’s equation, Comp.Math. Appl. 1 (1975), 351–360.

[Rosser2] J. Barkley ROSSER, Finite-difference solution of Poisson’s equation in rectangles ofarbitrary proportions, Zeitschrift fur Angewandte Mathematik und Physik (ZAMP)28 (1977), no.2, 185–196.

[Strang] Gilbert STRANG, Linear Algebra and its Applications, 3rd Ed., Harcourt Brace Jo-vanovich, 1988.

[Strikwerda] John C. STRIKWERDA, Finite Difference Schemes and Partial Differential Equa-tions, 2nd Ed., SIAM, 2004.

[Thomas1] J. W. THOMAS, Numerical Partial Differential Equations: Finite Difference Meth-ods, Texts in Applied Mathematics 22, Springer, 1995.

[Thomas2] J. W. THOMAS, Numerical Partial Differential Equations: Conservation Laws andElliptic Equations, Texts in Applied Mathematics 33, Springer, 1999.

107


[TOS] Ulrich TROTTENBERG, Cornelis OOSTERLEE e Anton SCHULLER, Multigrid,Elsevier, 2001.

[Versteeg-Malalasekera] H. K. VERSTEEG e W. MALALASEKERA, An introduction to computational fluiddynamics: The finite volume method, Prentice Hall, 1995.

[Watkins] David S. WATKINS, Fundamentals of Matrix Computations, 2nd Ed., John Wiley &Sons, 2002.

[Young] David M. YOUNG, Iterative Solutions of Large Linear Systems, Academic Press,1971.

Documents

Notas de Aula M¶etodos Num¶ericos para Equa»c~oes ...rodney/notas_de_aula/metodos_numericos_EDPs... · encial parcial, procedemos µa discretiza»c~ao desta. Usando diferen»cas