Notas de Aula Equa C Oes Diferenciais Num Ericas: Rodney Josu e Biezuner
Notas de Aula Equa C Oes Diferenciais Num Ericas: Rodney Josu e Biezuner
Notas de Aula Equa C Oes Diferenciais Num Ericas: Rodney Josu e Biezuner
1
Rodney Josué Biezuner
Departamento de Matemática
Instituto de Ciências Exatas (ICEx)
Universidade Federal de Minas Gerais (UFMG)
9 de novembro de 2015
1
E-mail: rodney@mat.ufmg.br; homepage: http://www.mat.ufmg.br/∼rodney.
Sumário
1
Rodney Josué Biezuner 2
Para complementar este capı́tulo, consulta às referências [Iserles], [HNW] e [Butcher] é fortemente recomen-
dada.
1.1.1 Discretização
Denote a aproximação numérica para a solução exata y (t) para t > t0 por y. Embora em princı́pio uma
aproximação para a solução exata pudesse estar definida em todo o domı́nio em que a solução exata y estiver
4
Rodney Josué Biezuner 5
definida, nossa aproximação numérica y, computacionalmente calculável, só estará definida em um número
finito de pontos. Fixe um intervalo [t0 , T ]. Se quisermos analisar o comportamento assintótico da solução
quando t → ∞, este intervalo poderá ser mais tarde aumentado, escolhendo valores sucessivamente maiores
para T . Dividimos o intervalo [t0 , T ] em n subintervalos de mesmo comprimento igual a h = (T − t0 ) /n
através de n − 1 pontos interiores uniformemente espaçados
t0 = t0
t1 = t0 + h
t2 = t0 + 2h
.. (1.2)
.
tn−1 = t0 + (n − 1) h
tn = t0 + nh = T
de modo que [t0 , T ] = [t0 , t1 ] ∪ [t1 , t2 ] ∪ . . . ∪ [tn−1 , tn ]. Esta partição uniforme do intervalo [t0 , T ] é chamada
uma discretização uniforme do intervalo [t0 , T ] de norma h. Ela produz um conjunto discreto de pontos
{t0 , t1 , . . . , tn−1 , tn } que chamamos uma malha discretizada para o intervalo [t0 , T ] (também chamada
uma malha de nós). Introduzimos a notação:
yi = y(ti ) (1.3)
para denotar o valor da solução exata no instante de tempo ti ; o valor da solução aproximada neste ponto
será denotado correspondentemente por yi . Assim, em princı́pio, a solução aproximada y existirá apenas
na malha discretizada; através de interpolação, valores para y podem ser produzidos em outros pontos
do intervalo. Uma vez discretizado o domı́nio da equação diferencial, procedemos à discretização desta
propriamente dita. Isso pode ser feito de várias maneiras. Veremos algumas destas maneiras nas seções a
seguir. Para analizar o quão bem estas diversas soluções aproximam a solução exata, introduzimos o seguinte
conceito:
1.3 Definição. Definimos o erro absoluto da solução aproximada no ponto i da malha por
ei = yi − y i , (1.4)
o erro de discretização local da solução aproximada no ponto i da malha por
|ei | = |yi − yi | (1.5)
e o erro de discretização global da solução aproximada no intervalo [t0 , tn ] discretizado por uma partição
de norma h por
E ([t0 , T ] , h) = max |ei | = max |yi − yi | . (1.6)
i=1,...,n i=1,...,n
Substituindo a derivada primeira pela expressão dada na EDO, podemos definir uma solução aproximada y
por
yi+1 = yi + hf (ti , yi ) . (1.9)
Este é o chamado método de Euler. Geometricamente, estamos aproximando yi+1 através de seguir a reta
tangente à solução em yi durante um intervalo de tempo h. Por este motivo, o método de Euler também é
chamado método da reta tangente. Vamos agora examinar o quão bem a solução numérica obtida através do
método de Euler aproxima a solução exata. Para isso, introduzimos o seguinte conceito:
1.4 Definição. O erro de truncamento local do método de Euler no ponto t é definido por
y (t + h) − y (t)
L (t, h) = − f (t, y) . (1.10)
h
O erro de truncamento global do método de Euler é definido por
O erro de truncamento local do método de Euler mede simplesmente o quanto o quociente de diferença que
aproxima a derivada y 0 (t) difere do valor exato para y 0 (t) dado por f (t, y) na equação diferencial. Além
disso, se yi = yi , isto é, se começássemos da solução exata, a diferença entre a solução exata yi+1 e a solução
aproximada yi+1 seria dada por
de modo que o erro de truncamento local vezes h é precisamente o erro de discretização absoluto local que
seria produzido pelo método de Euler em um único passo se começássemos da solução exata. Evidentemente,
embora no primeiro passo comecemos do valor exato da solução, dado pela condição inicial, este erro vai se
propagando à medida que calculamos os valores aproximados y1 , . . . , yn−1 , yn , com cada valor nesta sequência
dependendo do valor aproximado calculado anteriormente. Precisamos saber se o erro pode ser limitado e,
crucialmente, se o erro pode ser tornado arbitrariamente pequeno tomando h pequeno. Em outras palavras,
precisamos saber se a solução aproximada converge para a solução exata quando h → 0. Além disso,
estamos interessados na velocidade de convergência, isto é, o quão pequeno h precisa ser para um certo valor
especificado do erro.
1.5 Teorema (Convergência do Método de Euler). Considere o problema (1.1) com solução exata no in-
∂f
tervalo [t0 , T ]. Se f é de classe C 1 e é limitada em [t0 , T ] × R, então a solução aproximada obtida pelo
∂y
método de Euler converge para a solução exata quando h → 0 e o erro é de primeira ordem, isto é,
E ([t0 , T ] , h) = O (h) .
Rodney Josué Biezuner 7
Prova. Considere uma discretização uniforme do intervalo [t0 , T ] com norma h, de modo que tn = T , para
n = (tn − t0 ) /h. Seja
∂f
λ = sup (t, y) (1.12)
t∈[t0 ,T ] ∂y
y∈R
a constante de Lipschitz para f . Observe que como f é de classe C 1 , segue da regra da cadeia que y é uma
solução de classe C 2 , pois
∂f
|f (t, x) − f (t, y)| 6 sup (t, y) |x − y| = λ |x − y| ,
t∈[t0 ,T ] ∂y
y∈R
Denotando
A = 1 + λh,
B = N h2 ,
temos uma desigualdade iterativa para o erro de discretização local com a forma
|ei+1 | 6 A |ei | + B.
|ei+1 | 6 A |ei | + B
6 A [A |ei−1 | + B] + B
= A2 |ei−1 | + AB + B
6 A2 [A |ei−2 | + B] + AB + B
= A3 |ei−2 | + A2 B + AB + B
6 ...
= Ai+1 |e0 | + Ai B + . . . + A2 B + AB + B.
|ei+1 | 6 Ai + . . . + A2 + A + 1 B.
Observe que estes números são todos positivos, pois o erro tende a crescer à medida que i cresce, isto é,
à medida que executamos mais passos no método de Euler. De qualquer modo, o crescimento do erro é
controlado, pois para todo i temos
|ei | 6 An−1 + . . . + A2 + A + 1 B
An − 1
= B
A−1
n
(1 + λh) − 1
= N h2
λh
n
(1 + λh)
6 N h.
λ
Portanto,
n
(1 + λh)
E ([t0 , T ] , h) = max |ei | 6 N h. (1.18)
i=1,...,n λ
Usando a desigualdade
1 + x 6 ex ,
válida para todo x > 0 (que por sua vez segue diretamente da expansão em série da função exponencial),
escrevemos
eλnh eλ(T −t0 )
E ([t0 , T ] , h) 6 Nh 6 N h. (1.19)
λ λ
Em particular, E ([t0 , T ] , h) = O (h).
Portanto, a convergência do método de Euler é apenas linear. Além disso, quanto maior é o intervalo
considerado, isto é, quanto maior é T , a tendência do erro é aumentar, embora isso fosse esperado, pela
acumulação dos erros de aproximação. Se f não é de Lipschitz, a convergência não é assegurada, mas isso
vale também para a solução exata: sua existència não é assegurada quando f não é de Lipschitz. Por outro
lado, a estimativa (1.19) obtida no teorema não é a melhor possı́vel; na prática, o resultado obtido pode ser
muito melhor, dependendo de f .
Rodney Josué Biezuner 9
Como f (t, y) = −100y, segue que f é contı́nua, uniformemente de Lipschitz em R2 com constante de
Lipschitz 100. A solução exata deste problema é
y (t) = e−100t .
Logo,
y 00 (t) = 104 e−100t
e
sup |y 00 | 6 104 .
[0,+∞)
Pela estimativa (1.19) obtida no teorema anterior, o erro de discretização global no intervalo [0, T ] é limitado
por
E ([0, T ] , h) 6 100e100T h. (1.20)
Por outro lado, a fórmula iterativa de Euler neste caso se reduz a
de modo que
y0 = 1,
y1 = 1 − 100h,
2
y2 = (1 − 100h) ,
..
.
n
yn = (1 − 100h) .
Assim
i
E ([0, T ] , h) = max |yi − yi | = max e−100ih − (1 − 100h)
i=1,...,n i=1,...,n
i
= max e−100ih − (1 − 100h) .
i=1,...,n
Este erro efetivo é muito menor que a estimativa de erro dada pelo teorema em (1.20). Por exemplo, se
h = 10−2 , segue que
E ([0, T ] , h) = max e−i = e−1 ,
i=1,...,n
E ([0, T ] , h) 6 e100T ,
é Z t
y (t) = y0 + f (s) ds.
t0
Existe uma teoria rica e métodos robustos e eficientes para calcular integrais numericamente. Esta é a base
dos métodos de Runge-Kutta para resolver equações diferenciais ordinárias numericamente. Integrando de
ti até ti+1 = ti + h temos
Z ti +h
yi+1 = yi + f (s, y (s)) ds.
ti
Por um método de quadratura (veja [Iserles], p. 33 para uma breve introdução) a integral pode ser numeri-
camente integrada, fornecendo
m
X
yi+1 = yi + h aj f (ti + hbj , y (ti + hbj )) (1.23)
j=1
para alguns números aj , bj ∈ [0, 1], j = 1, . . . , m. Antes de fazermos algumas escolhas especı́ficas de métodos
de Runge-Kutta com diferentes velocidades de convergência, vamos analizar a convergência de métodos de
passo único de uma maneira mais geral. Na próxima subseção, definiremos o que queremos dizer exatamente
por métodos de Runge-Kutta e, além disso, forneceremos uma formulação mais conveniente para trabalhar.
1.7 Definição. Um método de passo único (explı́cito) para a resolução numérica de equações diferenciais
ordinárias é um método da forma
yi+1 = yi + hφ (ti , yi ; h) , (1.24)
onde φ = φ (t, y; h) é chamada a função de iteração.
A função de iteração representa uma aproximação para f (t, y) no intervalo [t0 , T ] onde buscamos a solução
aproximada y para a solução exata y (t). O método de Euler é um método de Runge-Kutta de primeira
ordem onde
φ (ti , yi ; h) = f (ti , yi ) (1.25)
é a função de iteração.
Rodney Josué Biezuner 11
1.8 Exemplo (Método de Heun). O método de Heun é um método de Runge-Kutta de segunda ordem
(veja Corolário 1.14 a seguir) onde
1
φ (ti , yi ; h) = [f (ti , yi ) + f (ti+1 , yi + hf (ti , yi ))] (1.26)
2
é a função de iteração. Neste caso, aproximamos f (ti+1 , yi+1 ) pela média dos valores f (ti , yi ) e f (ti+1 , y
bi+1 ),
onde y bi+1 é a aproximação de yi+1 que seria obtida pelo método de Euler. Geometricamente, ao invés de
seguirmos a direção da reta tangente à curva y (t) em ti , como fazemos no método de Euler, no método de
Heun seguimos a média das direções das retas tangentes em ti e ti+1 .
1.9 Exemplo (Método de Euler modificado). O método de Euler modificado ou método do ponto
médio é um método de Runge-Kutta de segunda ordem (veja Corolário 1.14 a seguir) que tem como função
de iteração
h h
φ (ti , yi ; h) = f ti + , yi + f (ti , yi ) (1.27)
2 2
Geometricamente, ao invés de seguirmos a direção da reta tangente à curva y (t) em ti como no método
de Euler, no método de Euler modificado seguimos a direção da reta tangente no ponto médio do intervalo
[ti , ti+1 ].
1.10 Exemplo (Método de Runge-Kutta de Quarta Ordem Clássico). O método de Runge-Kutta de
quarta ordem clássico (veja Corolário 1.21 a seguir), imensamente popular, tem como função de iteração
1
φ (ti , yi ; h) = (F1 + 2F2 + 2F3 + F4 ) (1.28)
6
onde
F1 = f (ti , yi ) ,
h h
F2 = f t i + , y i + F1 , (1.29)
2 2
h h
F3 = f t i + , y i + F2 ,
2 2
F4 = f (ti + h, yi + hF3 ) .
1.11 Definição. O erro de iteração local de um método de passo único no ponto (t + h, y (t + h)) é
definido por
η (t, h) = y (t + h) − [y (t) + hφ (t, y (t) ; h)] . (1.30)
Dizemos que um método de passo único tem ordem de consistência p se
η (t, h) = O hp+1 .
1.12 Teorema. Considere um método de passo único
yi+1 = yi + hφ (ti , yi ; h)
para o problema (1.1) no intervalo [t0 , T ]. Se φ é uniformemente de Lipschitz em relação à segunda variável,
independentemente de 0 < h 6 T − t, e tem ordem de consistência p, então o método tem ordem de
convergência p.
Rodney Josué Biezuner 12
para todo t ∈ [t0 , T ], para todos x, y ∈ R e para todo 0 < h 6 T − t. Seja também
ei = yi − yi ,
ηi = η (ti , h) .
yi+1 = yi + hφ (ti , yi ; h) + ηi ,
yi+1 = yi + hφ (ti , yi ; h) ,
donde
ei+1 = ei + h [φ (ti , yi ; h) − φ (ti , yi ; h)] + ηi . (1.33)
Logo,
Prova. Temos
φ (t, y; h) = a1 f (t, y) + a2 f (t + b1 h, y + b2 hf (t, y)) .
Denote
∂f
M= sup .
[t0 ,T ]×R ∂y
Como
segue que φ satisfaz a condição de Lipschitz requerida pelo Teorema 1.12, com constante de Lipschitz
λ = M [a1 + a2 + a2 b2 M (T − t0 )] .
φ (t, y; 0) = a1 f (t, y) + a2 f (t + 0, y + 0)
= (a1 + a2 ) f (t, y)
= f (t, y) ,
donde
η (t, h) = y (t + h) − [y (t) + hφ (t, y (t) ; h)] = O h3 .
1.14 Corolário. O método de Heun e o método de Euler modificado são métodos convergentes de segunda
ordem.
Prova. Os coeficientes da função de iteração do método de Heun satisfazem a1 = a2 = 12 e b1 = b2 = 1.
Os coeficientes da função de iteração do método de Euler modificado satisfazem a1 = 0, a2 = 1 e
b1 = b2 = 21 .
A convergência de métodos de Runge-Kutta de ordem mais alta será vista na próxima subseção.
a21 ,
a31 , a32 ,
a41 , a42 , a43 , (1.39)
..
.
as1 , as2 , as3 , . . . , as,s−1 ,
b1 , . . . , b s ,
c2 , . . . , cs ,
Rodney Josué Biezuner 15
coeficientes reais. Um método de Runge-Kutta explı́cito com s estágios é um método de passo único com
função de iteração da forma
φ (ti , yi ; h) = b1 k1 + . . . + bs ks , (1.40)
onde
k1 = f (ti , yi ) ,
k2 = f (ti + c2 h, yi + a21 k1 h) ,
k3 = f (ti + c3 h, yi + (a31 k1 + a32 k2 ) h) (1.41)
..
.
ks = f (ti + cs h, yi + (as1 k1 + . . . + as,s−1 ks−1 ) h) ,
com a condição
i−1
X
ci = ai,j . (1.42)
j=1
Observamos que a condição (1.42) expressa que todos os pontos onde f é calculada são aproximações de
primeira ordem para a solução e simplificam consideravelmente a obtenção de condições para métodos de
ordem de convergência alta; para ordem de convergência baixa, no entanto, ela não é necessária (veja [HNW],
p. 142, Exercı́cio 6).
1.16 Notação. Os coeficientes de um método de Runge-Kutta especificado são geralmente organizados em
uma tabela RK com o seguinte formato:
0
c2 a21
c3 a31 a32
.. .. .. ..
. . . .
cs as1 as2 ... as,s−1
b1 b2 ... bs−1 bs
1.17 Exemplo. Métodos RK de 2 estágios tem tabelas RK com formato
0
c2 a21
b1 b2
Nesta notação, as condições do Teorema 1.13 para que um método RK de 2 estágios seja de segunda ordem
são
b1 + b2 = 1
1
b2 c2 = (1.43)
2
c2 = a21
A última condição é também à condição (1.42). Os métodos de Heun e de Euler modificado tem as seguintes
tabelas RK, respectivamente,
0 0
1 1 1 1
e 2 2 .
1 1
2 2 0 1
Rodney Josué Biezuner 16
1.18 Teorema (Convergência de Métodos de Runge-Kutta de Terceira Ordem). Considere o problema (1.1)
∂f
com solução exata no intervalo [t0 , T ], tal que f é de classe C 3 e é limitada em [t0 , T ] × R. Então um
∂y
método RK com tabela RK
0
c2 a21
c3 a31 a32
b1 b2 b3
tem ordem de convergência 3 se os seus coeficientes satisfazem
b1 + b2 + b3 = 1
1
b2 c2 + b3 c3 =
2
1
b2 c22 + b3 c23 =
3 (1.44)
1
b3 a32 c2 =
6
a = c
21
2
a31 = c3 − a32
Prova. A expansão de Taylor da solução exata (como f é de classe C 3 , a solução exata é de classe C 4 ) é
1 1
y (t + h) = y (t) + y 0 (t) h + y 00 (t) h2 + y 000 (t) h3 + O h4 .
2 6
Segue de (1.1) e da regra da cadeia que
y 0 (t) = f (t, y (t)) = f,
∂f ∂f
y 00 (t) = (t, y (t)) + f (t, y (t)) (t, y (t))
∂t ∂y
∂f ∂f
= +f ,
∂t ∂y
∂2f ∂2f
y 000 (t) = 2 (t, y (t)) + (t, y (t)) y 0 (t)
∂t ∂t∂y
∂f ∂f ∂f
+ (t, y (t)) + f (t, y (t)) (t, y (t)) (t, y (t))
∂t ∂y ∂y
2
∂2f
∂ f 0
+ f (t, y (t)) (t, y (t)) + 2 (t, y (t)) y (t)
∂y∂t ∂y
2
∂2f ∂2f ∂2f ∂2f
∂f ∂f ∂f
= 2 +f + +f +f + f2 2
∂t ∂t∂y ∂t ∂y ∂y ∂y∂t ∂y
2 2
2 2
∂ f ∂ f ∂f ∂f ∂f ∂ f
= 2 + 2f + +f + f2 2 .
∂t ∂t∂y ∂t ∂y ∂y ∂y
Daı́,
y (t + h) = y (t)
( " 2 #)
h2 ∂2f ∂2f 2
h ∂f ∂f ∂f ∂f ∂f 2∂f
+h f + +f + 2
+ 2f + +f +f
2 ∂t ∂y 6 ∂t ∂t∂y ∂t ∂y ∂y ∂y 2
4
+O h .
Rodney Josué Biezuner 17
Se mostrarmos que
" 2 #
h2 ∂2f ∂2f 2
h ∂f ∂f ∂f ∂f ∂f 2∂
f
+ O h3 , (1.45)
φ (t, y; h) = f + +f + 2
+ 2f + +f +f
2 ∂t ∂y 6 ∂t ∂t∂y ∂t ∂y ∂y ∂y 2
obteremos
η (t, h) = y (t + h) − [y (t) + hφ (t, y (t) ; h)] = O h4
Comparando (1.52) com (1.45), vemos que as condições para que estas duas expressões sejam idènticas são
exatamente
b1 + b2 + b3 = 1,
1
b2 c2 + b3 c3 = ,
2
1
b2 c22 + b3 c23 = ,
3
1
b3 a32 c2 = .
6
1.19 Exemplo. O método de Runge-Kutta de terceira ordem clássico tem tabela RK
0
1 1
2 2
.
1 −1 2
1 2 1
6 3 6
Ou seja,
h
yi+1 = yi + (F1 + 4F2 + F3 ) , (1.53)
6
com
F1 = f (ti , yi ) ,
h h
F2 = f ti + , yi + F1 ,
2 2
F3 = f (ti + h, yi − hF1 + 2hF2 ) .
Prova. Uma demonstração análoga à demonstração do Teorema 1.19 é possı́vel, embora muito, muito
mais longa. Sugerimos consultar [HNW], Capı́tulo II-2, pp. 143–155, e [Butcher], Capı́tulo 3, para uma
demonstração bem mais elegante baseada em teoria dos grafos. Além disso, esta técnica é essencial na
obtenção e análise da convergência de métodos RK de ordem superior a 4.
1.21 Corolário. O método de Runge-Kutta de quarta ordem clássico é de fato um método convergente de
quarta ordem.
Prova. A tabela RK do método de Runge-Kutta de quarta ordem clássico é
0
1 1
2 2
1 1
0
2 2
1 0 0 1
1 1 1 1
6 3 3 6
É fácil, embora tedioso, verificar que todas as condições do teorema são satisfeitas.
O maior esforço computacional requerido pelos métodos de Runge-Kutta (abreviação: métodos RK) é
o número de avaliações de f necessários. O método de Euler requer apenas uma avaliação da função f ,
enquanto que os métodos de Runge-Kutta de segunda ordem requerem duas avaliações de f , métodos RK
de terceira ordem requerem três avaliações e métodos de RK de quarta ordem requerem quatro avaliações.
O maior número de avaliações é contrabalançado pelo maior tamanho do passo h necessário para se atingir
uma precisão determinada e, consequentemente, um menor número de cálculos. Além disso, o tamanho
do passo tem fundamental importância em computação prática, pois quanto menor o tamanho do passo,
maiores são os erros de arredondamento, exatamente porque quanto menor h mais cálculos são necessários
Rodney Josué Biezuner 21
Aproximando f (t, y (t)) por um polinômio p (t) segue que a solução aproximada será dada por
Z ti+1
yi+1 = yi + p (t) dt. (1.55)
ti
Os métodos numéricos vistos anteriormente se aplicam em forma vetorial. Por exemplo, o método de Euler
se escreve da mesma forma
yi+1 = yi + f (ti , yi ) h (1.62)
mas desta vez a solução aproximada y é um vetor com n componentes:
−u00 = f (x)
em [0, L] ,
u (0) = a, u (L) = b.
2.1.2 Discretização
Dividimos o intervalo [0, L] em n subintervalos de comprimento ∆x = L/n através de n − 1 pontos interiores
uniformemente espaçados:
ui = u(xi ),
fi = f (xi ) .
Esta é uma discretização uniforme do intervalo [0, L]. Uma vez discretizado o domı́nio da equação dife-
rencial parcial, procedemos à discretização desta. Usando diferenças centradas para cada ponto interior xi ,
1 6 i 6 n − 1, temos
−ui−1 + 2ui − ui+1
= fi . (2.7)
∆x2
Para os pontos de fronteira, a condição de Dirichlet implica
u0 = a e un = b. (2.8)
Rodney Josué Biezuner 26
Portanto, para encontrar a solução discretizada temos que resolver o sistema linear com n − 1 equações a
n − 1 incógnitas:
∆x−2 (2u1 − u2 ) = f1 + a∆x−2
−2
∆x (−u1 + 2u2 − u3 ) = f2
.. ,
.
−2
∆x (−un−3 + 2un−2 − un−1 ) = fn−2
∆x−2 (−un−2 + 2un−1 ) = fn−1 + b∆x−2
ou seja,
f1 + a∆x−2
2 −1 u1
−1 2 −1 u2 f2
.. ..
.. ..
1 −1 . .
.
= .
.
∆x2 .. .. .. ..
. . −1
.
.
−1 2 −1 un−2 fn−2
−1 2 un−1 fn−1 + b∆x−2
Esta é uma matriz tridiagonal simétrica, esparsa. Além disso, como veremos na próxima subseção, ela é
positiva definida (isto é, seus autovalores são positivos) e portanto possui uma inversa, o que garante a
existência e unicidade da solução. Dada sua simplicidade, ela pode ser resolvida por eliminação gaussiana
ou sua inversa pode ser efetivamente calculada. Por exemplo, para n = 4, 5, 6 temos
1
1 12 13
−1 2 0 0 1 0 0
2 −1 0 3 2 1
−1
2
2
1
1
2 −1 = 0 1 3 0 3 0 2 1 0 = 4 2 4 2 ,
0 −1 2 1 2 3
1 2
0 0 1 0 0 43 3 3 1
1
1 21 13 14
2 0 0 0 1 0 0 0
−1
2 −1 0 0 0 1 2 2 0 2 0 0 1 1 0 0
4 3 2 1
−1 2 −1 0 3 4 3 2 1 3 6 4 2
= =
0 −1 2 −1 0 0 1 3 0 0 3 0 1 2 1 0 5 2 4 6 3
0 0 −1 2
4
4
3 3
1 2 3 4
4 1 2 3
0 0 0 1 0 0 0 5 4 4 4 1
1 1 1 1 1
1 2 3 4 5 2 0 0 0 0 1 0 0 0 0
−1 2 2 2
2
1
2 −1 0 0 0 0
1 3 4 5
0
3 0 0 0
2 1 0 0 0
−1 2 −1 0 0
3
3 3
1 2
0
−1 2 −1 0
0
= 0 1 4 5 0 0 4 0 0
3 3 1 0 0
0 0 −1 2 −1
4 4 1 1 3
0 0 0 −1 2 0
0 0 1 5 0 0 0 5 0
4 2 4 1 0
5 1 2 3 4
0 0 0 0 1 0 0 0 0 6 5 5 5 5 1
5 4 3 2 1
4 8 6 4 2
1
= 3 6 9 6 3
.
6
2 4 6 8 4
1 2 3 4 5
jπx
Uj (x) = sen ,
L
este fato sugere que os autovetores uj da matriz A são 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) ,
Prova. Temos
jπ 2jπ
jπ 2 sen − sen
sen n n
n
2 −1
jπ 2jπ 3jπ
−1 2 −1 2jπ
− sen − sen
sen + 2 sen
.. .. n n n n
−1 . . ..
= ..
.. ..
. .
. . −1 − sen (n − 3) jπ + 2 sen (n − 2) jπ − sen (n − 1) jπ
(n − 2) jπ
−1 2 −1 sen
n n n
n
−1 2
(n − 1) jπ
sen (n − 2) jπ (n − 1) jπ
n − sen + 2 sen
n n
Rodney Josué Biezuner 28
jπ
sen
n
2jπ
sen
n
jπ
..
= 2 1 − cos ,
n
.
(n − 2) jπ
sen
n
(n − 1) jπ
sen
n
pois
jπ 2jπ jπ jπ jπ jπ jπ
2 sen − sen = 2 sen − 2 sen cos = 2 1 − cos sen ,
n n n n n n n
(n − k − 1) jπ (n − k) jπ (n − k + 1) jπ
− sen + 2 sen − sen
n n n
(n − k) jπ jπ (n − k) jπ (n − k) jπ jπ
= − sen − + 2 sen − sen +
n n n n n
(n − k) jπ jπ (n − k) jπ jπ (n − k) jπ
= − sen cos + cos sen + 2 sen
n n n n n
(n − k) jπ jπ (n − k) jπ jπ
− sen cos − cos sen
n n n n
jπ (n − k) jπ
= 2 1 − cos sen ,
n n
e
(n − 2) jπ (n − 1) jπ
− sen + 2 sen
n n
(n − 1) jπ jπ (n − 1) jπ
= − sen − + 2 sen
n n n
(n − 1) jπ jπ (n − 1) jπ jπ (n − 1) jπ
= − sen cos + cos sen + 2 sen
n n n n n
(n − 1) jπ jπ (n − 1) jπ jπ (n − 1) jπ
= − sen cos − sen cos + 2 sen
n n n n n
jπ (n − 1) jπ
= 2 1 − cos sen ,
n n
(n − 1) jπ jπ (n − 1) jπ jπ
cos sen = − sen cos
n n n n
porque
(n − 1) jπ jπ (n − 1) jπ jπ (n − 1) jπ jπ
0 = sen jπ = sen + = sen cos + cos sen .
n n n n n n
Os autovalores de A são positivos, portanto A é uma matriz positiva definida. Observe que, fixado j, se n é
arbitrariamente grande então
jπ j 2 π2
cos ≈1− ,
n 2n2
Rodney Josué Biezuner 29
2n2 4n2
2 (n − 1) π (n − 1) π
λn−1 = 1 − cos = 1 − cos ≈ ,
∆x2 n L2 n L2
que não é uma boa aproximação para um autovalor do laplaciano. Vemos que se aumentarmos o número de
pontos de discretização (malha mais refinada) obteremos melhores aproximações e uma quantidade maior de
autovalores próximos aos autovalores do laplaciano. Para comparar, veja a tabela a seguir para os autovalores
do laplaciano no intervalo [0, π]; na primeira coluna temos os
autovalores
exatos do laplaciano, enquanto que
2n2 jπ
na demais colunas os autovalores da matriz A, λj = 2 1 − cos , com a linha superior indicando o
π n
número n de subintervalos na malha
n = 11 n = 21 n = 31 n = 51 n = 101 n = 1001
1 0.993 221 21 0.998 136 38 0.999 144 44 0.999 683 82 0.999 919 37 0.999 999 18
4 3.892 419 95 3.970 248 82 3.986 325 21 3.994 943 16 3.998 710 15 3.999 986 87
9 8.462 720 39 8.849 945 24 8.930 889 79 8.974 415 97 8.993 471 18 8.999 933 51
16 14.333 863 96 15.528 221 28 15.782 100 25 15.919 213 41 15.979 370 36 15.999 789 87
25 21.030 205 54 23.855 895 28 24.469 653 89 24.802 991 47 24.949 649 29 24.999 486 99
36 28.009 247 34 33.646 940 78 34.904 404 68 35.592 050 94 35.895 629 79 35.998 936 22
49 34.705 588 92 44.682 641 99 46.979 277 93 48.245 465 23 48.806 722 35 48.998 029 23
64 40.576 732 50 56.716 479 58 60.570 369 11 62.715 235 6 63.670 436 30 63.996 637 97
81 45.147 032 93 69.479 637 52 75.538 215 24 78.946 473 26 80.472 391 97 80.994 614 71
100 48.046 231 68 82.687 007 94 91.729 225 95 96.877 607 56 99.196 334 56 99.991 792 02
onde
a b
∆x =
, ∆y = ,
n m
substituı́mos o domı́nio Ω pela malha (ou gride) uniforme
de forma que
Ωd = (x, y) ∈ Ω : x = i∆x, y = j∆y, 0 6 i 6 n, 0 6 j 6 m .
A equação 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 diferença centrada, obtendo
−ui−1,j + 2ui,j − ui+1,j
−uxx ≈ ,
∆x2
−ui,j−1 + 2ui,j − ui,j+1
−uyy ≈ .
∆y 2
Portanto, a equação de Poisson discretizada toma a forma
−ui−1,j + 2ui,j − ui+1,j −ui,j−1 + 2ui,j − ui,j+1
+ = fi,j . (2.11)
∆x2 ∆y 2
Como a função u é calculada em cinco pontos, esta equação é chamada a fórmula dos cinco pontos.
Para cada ponto interior da malha obtemos uma equação, logo temos um sistema linear de (n − 1) (m − 1)
equações com o mesmo número de incógnitas. Diferente do caso unidimensional, no entanto, não existe uma
maneira natural de ordenar os pontos da malha, logo não podemos obter imediatamente uma representação
matricial para o problema discretizado. Precisamos antes escolher uma ordenação para os pontos da malha,
e como existem várias ordenações possı́veis, existem várias matrizes associadas.
Talvez a mais simples ordenação é a ordem lexicográfica induzida de Z2 . Nesta ordem, os pontos da
malha são percorridos linha por linha, da esquerda para a direita, de baixo para cima:
Neste caso, a matriz associada ao sistema linear é uma matriz (n − 1) (m − 1) × (n − 1) (m − 1) que pode
Rodney Josué Biezuner 31
Observe que
1 1
aii = 2 +
∆x2 ∆y 2
para todo 1 6 i 6 (n − 1) (m − 1), enquanto que
1
aij = −
∆y 2
se o ponto j é vizinho à esquerda ou à direita do ponto i e
1
aij = −
∆x2
se o ponto j é vizinho acima ou abaixo do ponto i. Por exemplo, no caso especial ∆x = ∆y, se n = 4 e m = 6
Rodney Josué Biezuner 32
é o problema
−∆d ud = fd em Ωd ,
(2.13)
ud = 0 sobre ∂Ωd .
Para estabelecer a existência e unicidade da solução discreta, provaremos que a matriz de discretização A,
que é uma matriz simétrica, é também uma matriz positiva definida, pois isso implica em particular que A
é invertı́vel.
Lembrando que as autofunções de Dirichlet do laplaciano no retângulo [0, a] × [0, b] são as funções
kπx lπy
Ukl (x, y) = sen sen ,
a b
este fato sugere que os autovetores ukl da matriz A na ordem lexicográfica são os vetores de coordenadas
= 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) ,
Em particular, este método não depende da maneira como os pontos da malha são ordenados (não depende
da matriz A usada para representar o laplaciano discreto). Como no método de separação de variáveis
contı́nuo, assumimos que as soluções da equação discreta acima são produtos da forma
onde F e G são funções de uma variável inteira. Substituindo esta expressão na equação de Helmholtz
discreta, obtemos
F (i − 1) G (j) − 2F (i) G (j) + F (i + 1) G (j) F (i) G (j − 1) − 2F (i) G (j) + F (i) G (j + 1)
+ = −λF (i) G (j) .
∆x2 ∆y 2
Rodney Josué Biezuner 34
F (i + 1) − (A + 2) F (i) + F (i − 1) = 0,
G (j − 1) − (B + 2) G (j) + G (j + 1) = 0.
2α = A + 2, 2β = B + 2.
Observe que
1−α 1−β
λ=2 + . (2.23)
∆x2 ∆y 2
Vamos resolver a equação para F , já que a equação para G é completamente análoga. Substituindo em
(2.21) uma solução da forma
F (i) = z i (2.24)
obtemos
z i−1 − 2αz i + z i+1 = 0,
donde, dividindo por z i−1 extraı́mos a equação quadrática (análoga à equação indicial)
z 2 − 2αz + 1 = 0. (2.25)
para algumas constantes c1 , c2 . Para determinarmos estas constantes e também α, aplicamos as condições
de fronteira, que implicam
F (0) = F (n) = 0.
Rodney Josué Biezuner 35
Como a equação para F é homogênea, a constante c é arbitrária. Aplicando a segunda, segue que
n n
z+ = z− ,
ou, como z+ z− = 1,
2n
z+ =1
Conseqüentemente, z+ é uma 2n-ésima raiz complexa de 1:
z+ = eijπ/n (2.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
(2.26) produz todas as soluções não-triviais F de (2.21).
Portanto,
z+ + z− eiπk/n + e−iπk/n kπ
α= = = cos , 0 6 k 6 n − 1,
2 2 n
e, escolhendo c = 1/2,
ikπ
Fk (i) = eiπki/n − e−iπki/n = sen .
n
Analogamente,
lπ
β = cos , 0 6 l 6 m − 1,
m
e
jlπ
Gl (j) = sen .
m
Segue que os autovalores são
1 kπ 1 lπ
λkl = 2 1 − cos + 1 − cos
∆x2 n ∆y 2 m
2.3 Teorema (Existência e Unicidade da Solução Discreta). Seja Ω = (0, a) × (0, b). Então o problema
discretizado
−∆d ud = fd em Ωd ,
ud = 0 sobre ∂Ωd ,
possui uma única solução.
Prova. Pelo lema anterior, os autovalores da matriz simétrica A são positivos, logo ela é uma matriz
invertı́vel.
Rodney Josué Biezuner 36
2.4 Lema (Propriedade do Valor Médio). Se ∆d ud = 0, então para pontos interiores vale
(ou seja, é um máximo em relação aos seus quatro vizinhos), somente se cada um dos seus quatro vizinhos
assume este mesmo valor máximo, e a desigualdade torna-se uma identidade. Aplicando este argumento a
todos os pontos da malha, concluı́mos que ou não existe um máximo interior, e portanto o máximo é atingido
na fronteira, ou existe um máximo interior e todos os pontos da malha assumem o mesmo valor, isto é, ud é
constante.
No caso geral ∆x 6= ∆y, se ∆d ud > 0 temos
1 1 1 ui,j−1 + ui,j+1 ui−1,j + ui+1,j
+ u i,j 6 + .
∆x2 ∆y 2 2 ∆y 2 ∆x2
donde
1 1 1 ui,j + ui,j ui,j + ui,j
+ ui,j = + ,
∆x2 ∆y 2 2 ∆y 2 ∆x2
logo nenhum dos seus quatro vizinhos pode assumir um valor menor que ui,j , isto é, cada um dos quatro
vizinhos assume o mesmo valor máximo e o argumento segue como no caso anterior. O caso ∆d ud 6 0 é
provado considerando-se −ud .
Rodney Josué Biezuner 37
Em primeiro lugar, obtemos uma estimativa a priori discreta (que também pode ser visto como um resultado
de regularidade discreto) para soluções da equação de Poisson discreta com condição de Dirichlet homogênea:
2
2.6 Lema (Estimativa a Priori). Seja Ω = (0, 1) . Seja ud uma solução de
−∆d ud = fd em Ωd ,
ud = 0 sobre ∂Ωd .
Então
1
kud k∞ 6 k∆d ud k∞ . (2.28)
8
Prova. Considere a função " 2 2 #
1 1 1
w (x, y) = x− + y−
4 2 2
e sua versão discretizada wd definida por
" 2 2 #
1 1 1
wi,j = xi − + yj − . (2.29)
4 2 2
Então
w>0 e ∆w = 1,
e também
wd > 0 e ∆d wd = 1, (2.30)
pois
wi−1,j − 2wi,j + wi+1,j wi,j−1 − 2wi,j + wi,j+1
∆d wd = +
∆x2 ∆y 2
2 2 2 2 2 2
"
1 xi−1 − 12 + yj − 12 − 2 xi − 12 − 2 yj − 12 + xi+1 − 21 + yj − 12
=
4 ∆x2
2 2 2 2 2 2 #
xi − 21 + yj−1 − 12 − 2 xi − 12 − 2 yj − 12 + xi − 12 + yj+1 − 12
+
∆y 2
2 2 2 2 2 2 #
"
1 xi−1 − 12 − 2 xi − 12 + xi+1 − 12 yj−1 − 12 − 2 yj − 12 + yj+1 − 21
= +
4 ∆x2 ∆y 2
" 2 2 2 2 2 2 #
1 xi − ∆x − 12 − 2 xi − 12 + xi + ∆x − 12 yj − ∆y − 12 − 2 yj − 12 + yj + ∆y − 12
= +
4 ∆x2 ∆y 2
"
1 x2i + ∆x2 + 14 − 2xi ∆x − xi + ∆x − 2 x2i − xi + 14 + x2i + ∆x2 + 14 + 2xi ∆x − xi − ∆x
=
4 ∆x2
#
yj2 + ∆y 2 + 14 − 2yj ∆y − yj + ∆y − 2 yj2 − yj + 14 + yj2 + ∆y 2 + 14 + 2yj ∆y − yj − ∆y
+
∆y 2
1 2∆x2 2∆y 2
= + = 1.
4 ∆x2 ∆y 2
Rodney Josué Biezuner 38
Segue do Princı́pio do Máximo Discreto que a função ud − k∆d ud k∞ wd assume o seu mı́nimo na fronteira.
Este último é igual a − k∆d ud k∞ max∂Ωd wd . Por sua vez, o máximo de wd na fronteira é menor ou igual ao
máximo de w em ∂Ω, dado por
2 2
1 1 1 1 1
max x− = max y− = .
06x61 4 2 06x61 4 2 8
1
ui,j 6 ui,j − k∆d ud k∞ wi,j 6 k∆d ud k∞ (2.33)
8
para todos i, j. Reunindo as duas desigualdades, segue que
1
|ui,j | 6 k∆d ud k∞
8
para todos i, j, o que conclui a demonstração.
2
2.7 Teorema. Seja Ω = (0, 1) . Sejam u ∈ C 4 Ω uma solução clássica para o problema de Dirichlet
−∆u = f em Ω,
u=0 sobre ∂Ω,
kud − vd k∞ 6 C D4 u ∆x2 + ∆y 2 .
L∞ (Ω)
(2.34)
Prova. A hipótese f ∈ C 2,α Ω garante que u ∈ C 4 Ω . Lembre-se que
∂4u
D4 u L∞ (Ω)
= sup (x, y) .
(x,y)∈Ω ∂xp ∂y q
p+q=4
Rodney Josué Biezuner 39
∂4u ∂4u
1 2
(xi , yj )∆y 2 + O ∆x4 , ∆y 4 .
− (∆d ud )i,j = (fd )i,j − (x i , yj )∆x + (2.36)
3! ∂x4 ∂y 4
Rodney Josué Biezuner 40
Para uma demonstração destes resultados, veja [Hackbusch], págs. 60-61. Se quisermos uma melhor ordem
de convergência para as soluções discretizadas, é necessário considerar outras forma de discretizar o laplaciano
através de diferenças finitas. Isto será feito na próxima seção.
Como procuramos um esquema de diferenças finitas com ordem de convergência maior que 2, queremos obter
uma solução não-nula para o sistema
c + c2 + c3 + c4 + c5 = 0
1
−2c 1 − c2 + c 4 + 2c 5 = 0
2c1 + 1 c2 + 1 c4 + 2c5 1
=
2 2 ∆x2 ;
4 1 1 4
− c1 − c2 + c4 + c5 =
0
3 6 6 3
2c + 1 c + 1 c + 2c =
1 2 4 5 0
3 24 24 3
isso implicaria em princı́pio em um esquema com ordem de convergência pelo menos igual a 3:
Como a matriz
1 1 1 1 1
−2 −1 0 1 2
1 1
2 0 2
2 2
4 1 1 4
−
3 −6 0
6 3
2 1 1 2
0
3 24 24 3
tem determinante igual a 1, ela é invertı́vel e o sistema possui a solução única
1 1
c1 = − ,
12 ∆x2
4 1
c2 = ,
3 ∆x2
5 1
c3 =−
2 ∆x2
4 1
c4 = ,
3 ∆x2
1 1
c5 =− .
12 ∆x2
Rodney Josué Biezuner 42
Observe a distribuição dos nove pontos. Além dos cinco usuais, foram acrescentados os quatro pontos que
ocupam as posições diagonais. Para os quatro pontos vizinhos horizontais ou verticais do ponto central, a
fórmula de Taylor produz
5
1 1 1 1 1 1 ∂ u
+ ∆x5 − c1 + c3 − c4 + c6 − c7 + c9 (xi , yj )
120 120 120 120 120 120 ∂x5
5
1 1 1 1 ∂ u
+ ∆x4 ∆y − c1 − c3 + c7 + c9 (xi , yj )
24 24 24 24 ∂x4 ∂y
∂5u
1 1 1 1
+ ∆x3 ∆y 2 − c1 + c3 + c7 + c9 (xi , yj )
12 12 12 12 ∂x3 ∂y 2
∂5u
1 1 1 1
+ ∆x2 ∆y 3 − c1 − c3 − c7 + c9 (xi , yj )
12 12 12 12 ∂x2 ∂y 3
5
4 1 1 1 1 ∂ u
+ ∆x∆y − c1 + c3 − c7 + c9 (xi , yj )
24 24 24 24 ∂x∂y 4
5
5 1 1 1 1 1 1 ∂ u
+ ∆y − c1 − c2 − c3 + c7 + c8 + c9 (xi , yj )
120 120 120 120 120 120 ∂y 5
Rodney Josué Biezuner 46
Para obter um esquema com ordem de convergência pelo menos igual a 3, precisarı́amos obter uma solução
não-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
1
c1 + c3 + c4 + c6 + c7 + c9 =
∆x2
c1 − c3 − c7 + c9 = 0
1
c1 + c2 + c3 + c7 + c8 + c9 =
∆y 2
−c1 + c3 − c4 + c6 − c7 + c9 = 0
−c1 − c3 + c7 + c9 = 0
−c1 + c3 − c7 + c9 = 0
−c1 − c2 − c3 + c7 + c8 + c9 = 0
c1 + c3 + c4 + c6 + c7 + c9 = 0
c1 − c3 − c7 + c9 = 0
c1 + c3 + c7 + c9 = 0
c1 − c3 − c7 + c9 = 0
c1 + c2 + c3 + c7 + c8 + c9 = 0
Infelizmente este sistema não tem solução pois ele é inconsistente: a sexta e a última equação são incom-
patı́veis, assim como a quarta e a décima primeira. Portanto, não existe uma fórmula de nove pontos
compacta tal que
−∆d ud = −∆u + O ∆x3 , ∆y 3 .
No entanto, em 1975 o matemático e lógico Rosser introduziu a seguinte fórmula de nove pontos compacta
no caso especial ∆x = ∆y (em [Rosser1]; veja também [Rosser2])
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
∆d ud = , (2.42)
6∆x2
que pode ser resumida na forma
−1 −4 −1
1
− ∆d ud = −4 20 −4 , (2.43)
6∆x2
−1 −4 −1
a qual produz um esquema convergente de quarta ordem se a solução u ∈ C 6 Ω (ou mesmo se u ∈ C 5,1 Ω
apenas) dependendo
de como a função f é discretizada. Para entender como isso ocorre, observe que se
u ∈ C 8 Ω a fórmula de Taylor produz
∆x2 2 ∆x4 ∂ 4 ∂4 ∂4
+ 4 2 2 + 4 ∆u + O ∆x6
−∆d ud = −∆u − ∆ u− (2.44)
12 360 ∂x4 ∂x ∂y ∂y
∆x2 ∆x4 ∂ 4 ∂4 ∂4
f + O ∆x6 .
= −∆u + ∆f + 4
+ 4 2 2
+ 4
(2.45)
12 360 ∂x ∂x ∂y ∂y
O ponto crucial aqui é que o erro é expresso em termos de −∆u e, conseqüentemente, por f . Ainda é
necessário escolher uma discretização especial para f :
fi,,j−1 + fi−1,j + 8fi,j + fi+1,j + fi,j+1
fd = (2.46)
12
ou
1
1
fd = 1 8 1 . (2.47)
12
1
Rodney Josué Biezuner 47
Usando a fórmula de Taylor para f , obtemos que esta discretização especial para f satisfaz
∆x2
∆f + O ∆x4 .
fd = f + (2.48)
12
Somando esta estimativa com (2.45), e usando −∆d ud = fd , −∆u = f , obtemos
Para este esquema, pode-se provar (veja [Hackbusch], pág. 64) que existe uma constante C > 0 tal que
kud − vd k∞ 6 C∆x4 kukC 6 (Ω) ou kud − vd k∞ 6 C∆x4 kukC 5,1 (Ω) (2.49)
O esquema de Rosser também satisfaz o princı́pio do máximo. Concluindo, vemos que uma maior regularidade
da solução permite obter métodos de diferenças finitas com maior ordem de convergência, embora esta não
seja uma tarefa simples.
Embora esta condição não seja uma condição de fronteira e aparece apenas por causa do sistema de coor-
denadas utilizado, ela acaba funcionando como uma condição de fronteira em muitos métodos numéricos (e
mesmo analı́ticos), pois não deixa de ser uma condição na fronteira do retângulo (0, R) × (0, 2π).
∆r
∆θ
onde
R 2π
∆r = , ∆θ = .
n m
Rodney Josué Biezuner 48
ui,j = u (ri , θj ) ,
fi,j = f (ri , θj ) ,
para todo 0 6 j 6 m, já que existe apenas um ponto associado com i = 0 (a origem, correspondente a r = 0).
Além disso, pela condição de continuidade, devemos ter também
para todo 0 6 i 6 n. Usando uma diferença centrada usual para derivadas segundas, o terceiro termo do
laplaciano em coordenadas polares pode ser aproximado para pontos interiores do disco por
1 1 ui,j−1 − 2ui,j − ui,j+1
2
uθθ (ri , θj ) ≈ 2 . (2.52)
r ri ∆θ2
para 1 6 i 6 n − 1. Como este esquema de diferenças finitas foi obtido através de diferenças centradas,
ele deve ser de segunda ordem. No entanto, devemos ter cuidado ao discretizar a equação de Poisson na
Rodney Josué Biezuner 49
origem para preservar esta ordem de convergência. Para isso, multiplicamos a equação de Poisson por r e
integramos o resultado sobre um pequeno disco Dε centrado na origem de raio ε:
Z 2π Z ε Z 2π Z ε
1 1
f r drdθ = r (rur )r + 2 uθθ drdθ
0 0 0 0 r r
Z 2π Z ε Z ε Z 2π
1
= (rur )r drdθ + uθθ drdθ
0 0 0 r 0
Z 2π Z ε
ε 1 2π
= [rur ]0 dθ + [uθ ]0 drdθ
0 0 r
Z 2π
=ε ur (ε, θ) dθ,
0
aproximando a derivada primeira ur (∆r/2, θ) = (ur )i+1/2,j por diferenças centradas e f por f (0) (pois ∆r
é suposto pequeno), de modo que
u1,j − u0,j
ur (∆r/2, θj ) ≈ ,
∆r
Z 2π Z ∆r/2 Z 2π Z ∆r/2 ∆r/2
r2 π
f r drdθ ≈ f (0) r drdθ = 2πf (0) = f (0) ∆r2 ,
0 0 0 0 2 0 4
e assim
m−1
∆r X u1,j − u0,j π
∆θ = f (0) ∆r2 ,
2 j=0 ∆r 4
donde, como u0 := u0,j independe de j, segue que o valor de u na origem será dado por
m−1
∆θ ∆θ X π
m u0 = u1,j − f (0) ∆r2 ,
2 2 j=0 4
Au = f ,
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 ) . (2.57)
Rodney Josué Biezuner 50
Observe que existem (n − 1) × m + 1 incógnitas. Nesta ordenação, segue que A tem a forma em blocos
α0 b
a
B1 −β1 I
..
−α2 I B2 −β2 I .
A= −α3 I B3 −β3 I ,
(2.58)
. .. . .. . ..
−αn−2 I Bn−2 −βn−2 I
−αn−1 I Bn−1
onde
4
α0 = ,
∆r2
−α1
a = ... ,
−α1 m×1
1 ri−1/2
αi = , i = 1, . . . , n − 1,
∆r2 ri
1 ri+1/2
βi = , i = 1, . . . , n − 2,
∆r2 ri
b = −β0 . . . −β0 1×m ,
2 ∆θ
β0 = ,
π ∆r2
I = Im ,
γi −δi 0 −δi
−δi γi −δi
−δi γi −δi
Bi = ,
.. .. ..
. . .
−δi γi −δi
−δi −δi γi m×m
onde
1 ri+1/2 + ri−1/2 2 1
γi = 2
+ 2 ,
ri ∆r ri ∆θ2
1 1
δi = 2 .
ri ∆θ2
Rodney Josué Biezuner 51
A matriz A em geral não é simétrica. 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 0
0 0 −α2 0 0 0 −δ2 γ2 −δ2 0 0 0 −β2 0 0 0
0 0 0 −α2 0 0 0 −δ2 γ2 −δ2 0 0 0 −β2 0 0
0 0 0 0 −α2 0 0 0 −δ2 γ2 −δ2 0 0 0 −β2 0
0 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 0
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 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 são diferentes porque os pontos (0, j), j = 0, . . . , m, são realmente um
único ponto e este ponto é vizinho a todos os pontos (1, j), j = 0, . . . , m.
A matriz de discretização A no caso do anel será um pouco mais simples, já que ela será igual à matriz
de discretização no caso do disco menos a primeira linha e a primeira coluna.
De fato,
1 2 1 000 3
u(x− ) = u(x) − u0 (x) (x − x− ) + u00 (x) (x − x− ) − u (ξ− ) (x − x− ) ,
2 3!
1 2 1 000 3
u(x+ ) = u(x) + u0 (x) (x+ − x) + u00 (x) (x+ − x) + u (ξ+ ) (x+ − x) ,
2 3!
para alguns ξ− ∈ [x− , x] , ξ+ ∈ [x, x+ ], de modo que
u (x) − u (x− ) 1 1 2
− = −u0 (x) + u00 (x) (x − x− ) − u000 (ξ− ) (x − x− ) ,
x − x− 2 6
u (x+ ) − u (x) 1 1 2
= u0 (x) + u00 (x) (x+ − x) + u000 (ξ+ ) (x+ − x) ,
x+ − x 2 6
(xi − tW ∆x, yj ) , (xi + tE ∆x, yj ) , (xi , yj − tS ∆y) , (xi , yj + tN ∆y) , com t∗ ∈ (0, 1]
Rodney Josué Biezuner 53
Se (xi , yj ) é um ponto interior distante da fronteira (isto é, não adjacente à fronteira), então t∗ = 1 e para este
ponto vale a fórmula dos cinco pontos usual. Observe que a matriz obtida pelo esquema de Shortley-Weller
não é simétrica, em geral.
Embora a ordem de aproximação do laplaciano para pontos próximos à fronteira é apenas 1, o esquema
de Shortley-Weller é convergente de segunda ordem. No próximo capı́tulo, provaremos que o problema
discretizado possui solução única.
2.6 Exercı́cios
1. Implemente os métodos discutidos neste capı́tulo computacionalmente, verifique a precisão comparando
com a solução exata e também a velocidade de convergência.
2. Discretize o problema de Poisson com valor de fronteira de Dirichlet a seguir, usando a fórmula de
cinco 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 com
as soluções exatas.
3. Prove que a fórmula dos nove pontos compacta satisfaz o princı́pio do máximo discreto.
4. Prove resultados equivalentes ao Lema 1.5 e ao Teorema 1.6 para a fórmula dos nove pontos compacta.
5. Investigue a ordem de convergência do esquema de diferenças finitas misto: fórmula dos nove pontos nos
pontos interiores distantes da fronteira e fórmula dos cinco pontos para pontos adjacentes à fronteira.
6. Encontre um esquema de diferenças finitas de segunda ordem para a equação de laplace tridimensional
em um paralelepı́pedo reto. Escolha uma ordenação apropriada dos pontos da malha e descreva a
matriz de discretização obtida. Implemente o método no computador.
7. Mostre que o esquema de diferenças finitas em coordenadas polares introduzido neste capı́tulo satisfaz
o princı́pio do máximo discreto desde que o valor de u0 seja dado pela fórmula (2.56).
Rodney Josué Biezuner 54
8. Mostre que se ∆d denota o esquema de diferenças finitas em coordenadas polares introduzido neste
capı́tulo e Ω é o disco unitário, então vale a estimativa a priori: se ud é uma solução de
−∆d ud = fd em Ωd ,
ud = 0 sobre ∂Ωd ,
então
1
k∆d ud k∞
kud k∞ 6 (2.68)
4
desde que o valor de u0 seja dado pela fórmula (2.56). Conclua que este esquema tem ordem de
convergência 2.
Implemente alguns exemplos deste problema computacionalmente e compare os resultados obtidos com
as soluções exatas.
11. Mostre que tomando o “quadrado” da fórmula de três pontos para o laplaciano unidimensional (es-
quema de diferenças centradas para a derivada segunda) obtemos a seguinte fórmula de cinco pontos
para o operador biharmônico unidimensional (esquema de diferenças centradas para a derivada quarta):
ui−2 − 4ui−1 + 6ui − 4ui+1 + ui+2
δ 4 ui = (2.69)
∆x4
Usando a fórmula de Taylor, obtenha o expoente p tal que
12. O esquema de diferenças finitas mais simples para o operador biharmônico ∆2 em duas dimensões é a
seguinte fórmula de 13 pontos (para o caso ∆x = ∆y):
1
2 −8 2
2 1
∆ u= 1 −8 20 −8 1 . (2.70)
∆x4
2 −8 2
1
Mostre que esta fórmula pode ser obtida a partir do “quadrado” da fórmula de cinco pontos para
o laplaciano. Como a equação biharmônica não satisfaz o princı́pio do máximo, a demonstração da
ordem de convergência deste esquema necessita de argumentos diferentes dos usados neste capı́tulo
para o laplaciano. Na realidade, dependendo de como as duas
condições de fronteira são discretizadas,
a ordem de convergência deste método pode ser O ∆x3/2 ou O ∆x2 . Veja [Hackbusch], pág. 103 e
págs. 105-109, para detalhes e referências.
Capı́tulo 3
Determinar a existência e unicidade de soluções discretas para as matrizes de discretização obtidas via
esquemas de diferenças finitas através do cálculo de seus autovalores como fizemos no capı́tulo anterior para
diferenças centradas em uma dimensão e para a fórmula de cinco pontos é inviável em geral (tente calcular
os autovalores da matriz de discretização para a fórmula dos nove pontos, para o esquema em coordenadas
polares e para o esquema de Shortley-Weller). Neste capı́tulo, desenvolveremos métodos mais gerais e mais
fáceis de aplicar.
1. Norma l1
n
X
kAk1 = |aij | . (3.2)
i,j=1
De fato,
n
X n
X n
X n
X n
X n
X
kABk1 = aik bkj 6 |aik bkj | 6 |aik blj | = |aik | |blj | = kAk1 kBk1 .
i,j=1 k=1 i,j,k=1 i,j,k,l=1 i,j=1 k,l=1
2. Norma l2
1/2
n
X 2
kAk2 = |aij | . (3.3)
i,j=1
Com efeito,
n n 2 n n
! n
! n n
2
X X X X 2
X 2
X 2
X 2 2 2
kABk2 = aik bkj 6 |aik | |blj | = |aik | |blj | = kAk2 kBk2 .
i,j=1 k=1 i,j=1 k=1 l=1 i,k=1 j,l=1
55
Rodney Josué Biezuner 56
A norma l2 também é chamada norma euclidiana e, mais raramente e somente para matrizes, norma
de Schur, norma de Frobenius ou norma de Hilbert-Schmidt.
3. Norma l∞ modificada
A norma l∞
kAk∞ = max |aij | .
16i,j6n
é uma norma vetorial no espaço das matrizes complexas, mas não é uma norma matricial, pois se
1 1
A= ,
1 1
então
2 2
A2 =
2 2
e portanto
A2 ∞
= 2 > 1 = kAk∞ kAk∞ .
Mas um múltiplo escalar desta norma vetorial é uma norma matricial:
Com efeito,
n
X n
X n
X
kABkn∞ = n max aik bkj 6 n max |aik bkj | 6 n max kAk∞ kBk∞
16i,j6n 16i,j6n 16i,j6n
k=1 k=1 k=1
= n kAk∞ n kBk∞ = kABkn∞ .
4. Norma induzida
Dada uma norma vetorial |·| em Cn , ela induz uma norma matricial através da definição
|Ax|
kAk = max |Ax| = max . (3.5)
|x|=1 x6=0 |x|
De fato,
|ABx| |ABx| |Bx| |ABx| |Bx| |Ay| |Bx|
kABk = max = max 6 max max 6 max max = kAk kBk .
x6=0 |x| x6=0 |Bx| |x| x6=0 |Bx| x6=0 |x| y6=0 |y| x6=0 |x|
Esta norma também é chamada norma do operador. Ela satisfaz a propriedade muitas vezes útil
n
X n
X n
X
|Ax|∞ = max aij xj 6 max |aij xj | 6 max |aij | |x|∞ = kAkL |x|∞ ,
16i6n 16i6n 16i6n
j=1 j=1 j=1
Rodney Josué Biezuner 57
de modo que
max |Ax|∞ 6 kAkL .
|x|=1
n
X n
X n
X
max |Ax|∞ > |Ay|∞ = max aij yj > akj yj = |akj | .
|x|∞ =1 16i6n
j=1 j=1 j=1
Esta norma é induzida pela norma vetorial l1 . De fato, escrevendo A em termos de suas colunas
A = [A1 . . . An ]
segue que
kAkC = max |Aj |1 .
16j6n
donde
max |Ax|1 6 kAkC .
|x|1 =1
|Ay|1 = |Aj |1
7. p-normas
Este é o nome geral para as normas induzidas pela norma vetorial lp . O caso especial da norma induzida
pela norma vetorial l2 (a norma vetorial euclidiana) é também chamada a norma espectral e satisfaz
p n√ o
k|A|k2 = λmax = max λ : λ é um autovalor de A∗ A .
Rodney Josué Biezuner 58
De fato, A∗ A é uma matriz hermitiana e possui autovalores não-negativos, pois se A∗ Ay = λy, então
2 2
λ |y|2 = hy, λyi2 = hy, A∗ Ayi2 = hAy, Ayi2 = |Ay|2
e, além disso, pela caracterização variacional dos autovalores de uma matriz hermitiana temos
2
hA∗ Ax, xi2 |Ax|2
λmax = max 2 = max 2 .
x6=0 |x|2 x6=0 |x|2
Observe que a 2-norma é diferente da norma matricial l2 . Note também que se A é uma matriz
hermitiana, então A∗ A = A2 e k|A|k2 é portanto o módulo do maior autovalor de A, isto é, a norma
espectral de A é o raio espectral de A, definido como sendo o maior valor absoluto dos autovalores
de A:
ρ (A) = max |λi | ,
i=1,...,n
kAkS = S −1 AS (3.9)
Lembramos que todas as normas em um espaço vetorial são equivalentes, e isso vale em particular para
normas matriciais.
3.2 Proposição. Se A é uma matriz estritamente diagonalmente dominante, então A é invertı́vel.
Prova. Uma matriz A é invertı́vel se existe alguma norma matricial k·k tal que kI − Ak < 1. De fato, se
esta condição é satisfeita, então a inversa é dada explicitamente pela série
∞
X k
A−1 = (I − A) . (3.10)
k=0
Rodney Josué Biezuner 59
P∞
A condição kI − Ak < 1 garante a convergência desta série, pois a série geométrica k=0 rk tem raio de
convergência 1; como para todo N temos
N
X N
X N
X N
X +1
k k k k N +1
A (I − A) = [I − (I − A)] (I − A) = (I − A) − (I − A) = I − (I − A) ,
k=0 k=0 k=0 k=1
DR (a) = {z ∈ C : |z − a| 6 R} .
Rodney Josué Biezuner 60
denota a soma dos valores absolutos dos elementos da linha i de A excetuando o elemento da diagonal
principal, então todos os autovalores de A estão contidos na união dos n discos de Gershgorin
n
[
G (A) = DRi (A) (aii ) . (3.12)
i=1
Além disso, se uma união de k destes discos forma uma região que é disjunta dos n − k discos restantes,
então existem exatamente k autovalores de A nesta região.
Prova. Seja λ um autovalor de A e x = (x1 , . . . , xn ) 6= 0 um autovetor associado. Seja k um ı́ndice tal que
isto é, xk é a coordenada de x de maior valor absoluto. Denotando por (Ax)k a k-ésima coordenada do vetor
Ax = λx, temos
Xn
λxk = (Ax)k = akj xj
j=1
que é equivalente a
n
X
xk (λ − akk ) = akj xj .
j=1
j6=k
Daı́,
n
X n
X n
X
|xk | |λ − akk | 6 |akj xj | = |akj | |xj | 6 |xk | |akj | = |xk | Rk (A) ,
j=1 j=1 j=1
j6=k j6=k j6=k
ou seja,
|λ − akk | 6 Rk (A) .
Isso prova o resultado principal do Teorema de Gershgorin (como não sabemos qual k é apropriado para
cada autovalor λ, e um mesmo k pode servir para vários autovalores λ, tudo o que podemos afirmar é que
os autovalores estão na união dos discos).
Para provar a segunda afirmação, escreva A = D + B, onde D = diag (a11 , . . . , ann ) e defina
At = D + tB
DRi (At ) (aii ) = {z ∈ C : |z − aii | 6 Ri (At )} = {z ∈ C : |z − aii | 6 tRi (A)} ⊂ DRi (A) (aii ) ,
Rodney Josué Biezuner 61
logo
Gk (At ) ⊂ Gk (A)
e
Gk (A) ∩ [G (At ) \Gk (At )] = ∅
para 0 6 t 6 1. Porque os autovalores são funções contı́nuas das entradas de uma matriz, o caminho
λi (t) = λi (At )
é 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, não pode haver mais que
k autovalores de A em Gk (A), pois os n − k autovalores restantes de A0 = D começam fora do conjunto
Gk (A) e seguem caminhos contı́nuos que permanecem fora de Gk (A).
A união G (A) dos discos de Gershgorin é conhecida como a região de Gershgorin. Observe que enquanto
não podemos em geral afirmar com certeza que cada disco de Gershgorin possui um autovalor, a segunda
afirmação do teorema permite-nos fazer tal conclusão desde que os discos de Gershgorin sejam dois a dois
disjuntos.
O Teorema dos Discos de Gershgorin permite entender o resultado da Proposição 3.2: se uma matriz A é
estritamente diagonalmente dominante, então os discos de Gershgorin DRi (A) (aii ) não interceptam a origem,
logo 0 não pode ser um autovalor para a matriz A, o que implica que A é invertı́vel. Além disso, se todos
os elementos da diagonal principal de A são reais e positivos, então os autovalores de A estão localizados no
semiplano direito de C, de modo que se A é também simétrica, concluı́mos que todos os autovalores de A
são positivos.
A aplicação mais óbvia do Teorema dos Discos de Gershgorin é na estimativa dos autovalores de uma
matriz, o que é importante se vamos usar os autovalores de matrizes de discretização para aproximar os
autovalores do laplaciano:
Aplicação 1. Pelo Teorema dos Discos de Gershgorin, os autovalores da matriz de discretização do lapla-
ciano no intervalo (0, π) discretizado com n + 1 pontos (esquema de diferenças finitas centradas para
a derivada segunda unidimensional)
2 −1
−1 2 −1
2
. . . .
n −1 . .
A= 2
π .. ..
. . −1
−1 2 −1
−1 2
estão todos localizados no intervalo (A é simétrica, logo seusautovalores são todos reais) centrado em
x = 2n2 /π 2 de raio 2n2 /π 2 , ou seja, no intervalo 0, 4n2 /π 2 . Em particular o maior autovalor de A
não pode exceder 4n2 /π 2 . Como os autovalores do laplaciano neste intervalo são da forma λj = j 2 ,
para termos esperança em aproximar o autovalor λj por autovalores da matriz A precisamos que
j 2 6 4n2 /π 2 , isto é, precisamos discretizar o intervalo (0, π) com
π
n> j
2
pontos. Isso dá uma estimativa bastante grosseira do quão refinada a nossa malha precisa ser para
aproximar os autovalores do laplaciano. Na prática, vimos que apenas os primeiros autovalores de
A aproximam bem os primeiros autovalores do laplaciano e portanto precisamos de uma malha com
um número muito maior de pontos. Observe que uma estimativa semelhante vale para a matriz de
Rodney Josué Biezuner 62
2
discretização M fornecida pela fórmula de cinco pontos no quadrado (0, π) quando tomamos ∆x =
2 2
∆y = π/n: como os autovalores de M estão localizados no intervalo de centro em x = 4n /π de raio
4n2 /π 2 , isto é, em 0, 8n2 /π 2 , precisamos de
π p2
n> √ i + j2
2 2
pontos no eixos horizontal e vertical para aproximar o autovalor i2 + j 2 . Por outro lado, no caso
bidimensional isso implica em uma matriz de discretização da ordem de i2 + j 2 .
Usos mais refinados do Teorema de Gershgorin permitem obter conhecimento mais preciso sobre onde
os 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 discos
de Gershgorin equivalente para as colunas de uma matriz. Em particular, todos os autovalores de A estão
localizados na interseção destas duas regiões: G (A) ∩ G (At ). Isso implica a seguinte estimativa simples para
o raio espectral de uma matriz complexa:
3.4 Corolário. Se A ∈ Mn (C), então
n
X n
X
ρ (A) 6 min max |aij | , max |aij | = min (kAkL , kAkC ) .
i=1,...,n j=1,...,n
j=1 i=1
Prova. O ponto no i-ésimo disco de Gershgorin que é mais distante da origem tem módulo
n
X
|aii | + Ri (A) = |aij |
j=1
Em particular,
n n
1 X X 1
ρ (A) 6 min max pj |aij | , max pj |aij | . (3.14)
p1 ,...,pn >0 i=1,...,n pi j=1,...,n p
j=1 i=1 i
Rodney Josué Biezuner 63
3.4 Propriedade FC
Na nossa busca por propriedades para matrizes diagonalmente dominantes que garantirão a sua invertibili-
dade, uma observação fundamental é a de que se A é uma matriz diagonalmente dominante, então 0 não
pode ser um ponto interior de nenhum disco de Gershgorin. De fato, se λ é um autovalor de A interior a
algum disco de Gershgorin então devemos ter desigualdade estrita
n
X
|λ − aii | < Ri (A) = |aij |
j=1
j6=i
Tais pontos λ na região de Gershgorin G (A) (não necessariamente autovalores de A) constituem precisa-
mente a fronteira ∂G (A) da região de Gershgorin. Chamaremos a fronteira de um disco de Gershgorin
{z ∈ C : |z − aii | = Ri (A)} um cı́rculo de Gershgorin.
3.6 Lema. Seja A ∈ Mn (C) e λ um autovalor de A que não é um ponto interior de nenhum disco de
Gershgorin. Seja x = (x1 , . . . , xn ) 6= 0 um autovetor associado a λ e k um ı́ndice tal que
aij 6= 0,
então
|xj | = |xk |
e o j-ésimo cı́rculo de Gershgorin também passa por λ.
Prova. Como na demonstração do Teorema de Gershgorin, temos
n
X n
X n
X
|xi | |λ − aii | 6 |aij xj | = |aij | |xj | 6 |xk | |aij | = |xk | Ri (A) (3.15)
j=1 j=1 j=1
j6=k j6=k j6=k
|λ − aii | 6 Ri (A) .
donde
n
X
|aij | (|xi | − |xj |) = 0.
j=1
j6=k
Esta é uma soma de termos não-negativos, pois |xi | > |xj |, logo se aij 6= 0 necessariamente devemos ter
|xj | = |xi | = |xk |.
Este lema técnico tem as seguintes conseqüências úteis:
3.7 Teorema. Seja A ∈ Mn (C) uma matriz cujas entradas são todas não-nulas e seja λ um autovalor de A
que não é um ponto interior de nenhum disco de Gershgorin. Então todo cı́rculo de Gershgorin de A passa
por λ (isto é, λ está na interseção de todos os cı́rculos de Gershgorin de A) e se x = (x1 , . . . , xn ) 6= 0 é um
autovetor associado a λ então
|xi | = |xj | para todos i, j = 1, . . . , n.
Prova. Pois, como A é diagonalmente dominante, se 0 é um autovalor de A então 0 não pode ser um ponto
interior de nenhum disco de Gershgorin. Por outro lado, pelo teorema anterior, segue que todo cı́rculo de
Gershgorin passa por 0. Entretanto, o i-ésimo cı́rculo de Gershgorin centrado em aii e com raio Ri < |aii |
não pode passar por 0. Concluı́mos que 0 não é um autovalor de A, logo A é invertı́vel.
Na verdade, usando com maior cuidado a informação dada pelo Lema 3.6 podemos obter resultados ainda
melhores:
3.9 Definição. Dizemos que uma matriz A = (aij ) ∈ Mn (C) satisfaz a propriedade FC se para todo par de
inteiros distintos i, j existe uma seqüência de inteiros distintos i1 = i, i2 , i3 , . . . , im−1 , im = j, com 1 6 m 6 n,
tais que todas as entradas matriciais
ai1 i2 , ai2 i3 , . . . , aim−1 im
são não-nulas.
Por exemplo, a matriz diagonalmente dominante não-invertı́vel
4 2 1
0 1 1 ,
0 1 1
já vista anteriormente, não satisfaz a propriedade FC porque o par 2, 1 não admite tal seqüência (a única
seqüência possı́vel é a23 , a31 ). Já qualquer par de inteiros distintos i, j tal que aij 6= 0 admite a seqüência
trivial não-nula aij , de modo que uma matriz cujas entradas não-diagonais são todas não-nulas satisfaz a
propriedade FC. O significado da abreviatura “FC”, ou “fortemente conexo”, ficará claro mais adiante.
Rodney Josué Biezuner 65
3.10 Teorema. Seja A ∈ Mn (C) uma matriz que satisfaz a propriedade FC e seja λ um autovalor de A
que não é um ponto interior de nenhum disco de Gershgorin. Então todo cı́rculo de Gershgorin de A passa
por λ (isto é, λ está na interseção de todos os cı́rculos de Gershgorin de A) e se x = (x1 , . . . , xn ) 6= 0 é um
autovetor associado a λ então
|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
Em particular, segue novamente do Lema 3.6 que o j-ésimo cı́rculo de Gershgorin passa por λ. Como j é
arbitrário, isso prova o teorema.
3.11 Corolário. Se A ∈ Mn (C) é uma matriz que satisfaz a propriedade FC e diagonalmente dominante
n
P
tal que |aii | > |aij | para pelo menos alguma linha i, então A é invertı́vel.
j=1
j6=i
Prova. Segue do teorema anterior da mesma forma que o Corolário 3.8 segue do Teorema 3.7.
Vamos tentar entender melhor o significado da propriedade FC. Note que ela se refere apenas à localização
dos elementos não-nulos de A fora da diagonal principal – os elementos da diagonal principal e os valores
especı́ficos dos elementos fora da diagonal principal são irrelevantes. Isso motiva as seguintes definições:
3.12 Definição. Dada uma matriz A = (aij ) ∈ Mn (C) definimos o módulo da matriz A como sendo a
matriz
|A| = (|aij |)
cujos elementos são os módulos dos elementos da matriz A e a matriz indicadora de A como sendo a
matriz
M (A) = (µij ) ,
onde
1 se aij 6= 0,
µij =
0 se aij = 0.
O conceito de uma seqüência de entradas não-nulas da matriz A que aparece na definição da propriedade
FC pode ser visualizado em termos de caminhos em um grafo associado a A:
3.13 Definição. Dada uma matriz A ∈ Mn (C), o grafo direcionado de A é o grafo direcionado Γ (A)
com n nodos 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 Γ é uma seqüência de arcos Pi1 Pi2 , Pi2 Pi3 , . . . em Γ. O
comprimento de um caminho direcionado é o número de arcos sucessivos no caminho direcionado. Um ciclo
é um caminho direcionado que começa e termina no mesmo nó.
Dizemos que um grafo direcionado é fortemente conexo se entre qualquer par de nodos distintos
Pi , Pj ∈ Γ existir um caminho direcionado de comprimento finito que começa em Pi e termina em Pj .
Rodney Josué Biezuner 66
Observe que quando Γ é um grafo direcionado com n nodos, se existe um caminho direcionado entre dois
nodos de Γ, então sempre existe um caminho direcionado entre estes dois nodos de comprimento menor que
ou igual a n − 1.
3.14 Teorema. A ∈ Mn (C) satisfaz a propriedade FC se e somente se Γ (A) é fortemente conexo.
Prova. Provaremos o teorema por indução. Para m = 1 a afirmativa é trivial. Para m = 2, temos
n
X n
X
2
|A| = (|A|)ik (|A|)kj = |aik | |akj | ,
ij
k=1 k=1
2
de modo que |A| 6= 0 se e somente se aik , akj são ambos não-nulos para algum ı́ndice k. Mas isso é
ij
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
n
X n
X
m+1 m m
|A| = (|A| )ik (|A|)kj = (|A| )ik |akj | =
6 0
ij
k=1 k=1
m
se e somente se (|A| )ik , akj são ambos não-nulos para algum ı́ndice k. Por hipótese de indução, isso é
equivalente a existir um caminho direcionado de comprimento m em Γ (A) de Pi para Pk e um caminho
direcionado de comprimento 1 em Γ (A) de Pk para Pj , isto é, um caminho direcionado de comprimento
m + 1 em Γ (A) de Pi para Pj . O mesmo argumento vale para M (A).
3.16 Definição. Seja A = (aij ) ∈ Mn (C). Dizemos que A > 0 se aij > 0 para todos 1 6 i, j 6 n e que
A > 0 se aij > 0 para todos 1 6 i, j 6 n.
3.17 Corolário. Seja A ∈ Mn (C). Existe um caminho direcionado de comprimento m em Γ (A) de cada
nodo Pi para cada nodo Pj se e somente se
m
|A| > 0
ou, equivalentemente, se e somente se
m
M (A) > 0.
3.18 Corolário. Seja A ∈ Mn (C). A satisfaz a propriedade FC se e somente se
n−1
(I + |A|) >0
2 n−1
se e somente se para cada par de ı́ndices i, j com i 6= j pelo menos um dos termos |A| , |A| , . . . , |A|
tem uma entrada positiva em (i, j). Pelo Teorema 3.15, isso ocorre se e somente se existe algum caminho
direcionado em Γ (A) de Pi para Pj com comprimento 6 n−1. Isto é equivalente a A satisfazer a propriedade
FC. O mesmo argumento vale para M (A).
Em geral, a maneira como uma matriz foi obtida (como as nossas matrizes de discretização; veja a última
seção do capı́tulo) torna clara se elas são matrizes que satisfazem a propriedade FC ou não. Se isso
não é possı́vel, e pretende-se verificar a propriedade FC através do Corolário 3.18, é preferı́vel calcular
n−1
[I + M (A)] , já que M (A) é uma matriz composta apenas de 0’s e 1’s.
Matrizes de transposição são simétricas. O efeito de multiplicar uma matriz A por uma matriz de transposição
à esquerda é trocar a posição de duas linhas da matriz A (no caso acima, as linhas k e l), enquanto que a
multiplicação de A por uma matriz de transposição à direita muda a posição de duas colunas de A (no caso
acima, as colunas k e l).
1 0 0 0 a11 a12 a13 a14 a11 a12 a13 a14
0 0 1 0 a21 a22 a23 a24 a31 a32 a33 a34
TA = 0 1 0 0 a31 a32 a33 a34 = a21 a22 a23 a24 ,
P t AP = Tm . . . T1 AT1 . . . Tm
é portanto obtida através da permutação de linhas e colunas de A, de modo que nenhum novo elemento é
criado ou algum elemento existente de A destruı́do.
3.19 Definição. Dizemos que uma matriz A ∈ Mn (C) é redutı́vel se existe alguma matriz de permutação
P e algum inteiro 1 6 m 6 n − 1 tal que
B C
P t AP =
0 D
Da definição vemos que se |A| > 0, então A é irredutı́vel, e para que A seja redutı́vel, ela precisa ter pelo
menos n − 1 zeros (caso m = 1). A motivação para este nome é a seguinte. Suponha que queiramos resolver
o sistema Ax = b e que A seja redutı́vel. Então, se escrevermos
t B C
A = P AP = ,
0 D
n−1
h + |A|)
(I i possui entradas diagonais nulas, logo podemos assumir que para algum par i 6= j temos
não
n−1 m
(I + |A|) = 0, o que implica [|A| ]ij = 0 para todo 1 6 m 6 n − 1. Pelo Teorema 3.15 (e observação
ij
imediatamente posterior à Definição 3.13), não existe um caminho direcionado em Γ (A) de comprimento
finito entre Pi e Pj . Defina os conjuntos de nodos
Por definição destes conjuntos, não pode existir nenhum caminho de algum nodo de S2 para algum nodo de
m
S1 , logo [|A| ]lk = 0 se Pl ∈ S2 e Pk ∈ S1 . E ambos os conjuntos são não-vazios, pois Pj ∈ S1 e Pi ∈ S2 .
Renomeando os nodos de modo que
n o
S1 = Pe1 , . . . , Pem ,
n o
S2 = Pem+1 , . . . , Pen ,
De fato, P é justamente a matriz de permutação que troca as colunas de tal forma que as variáveis anteriores
correspondentes aos nodos Pe1 , . . . , Pem no sistema Ax = b são as novas m primeiras variáveis do sistema linear
Ax = b; como não existe nenhum caminho direcionado entre nenhum dos nodos Pem+1 , . . . , Pen e qualquer um
dos nodos Pe1 , . . . , Pem , temos aij = 0 para m + 1 6 i 6 n e 1 6 j 6 m pelo Teorema 3.15.
3.21 Corolário. Uma matriz A ∈ Mn (C) é irredutı́vel se e somente se ela satisfaz a propriedade FC.
n
P
3.22 Proposição. Se A é uma matriz irredutı́vel, diagonalmente dominante tal que |aii | > |aij | para
j=1
j6=i
pelo menos alguma linha i, então A é invertı́vel.
Além disso, se A é hermitiana e todos os elementos da diagonal principal de A são positivos, então todos
os autovalores de A são positivos.
Prova. O resultado segue do Teorema 3.20, do Corolário 3.11 e do Teorema dos Discos de Gershgorin (veja
comentários após o Teorema 3.3).
dados dois pontos distintos Pi , Pj é fácil encontrar uma seqüência de ı́ndices i1 = i, i2 , i3 , . . . , im−1 , im = j,
com 1 6 m 6 n, tais que todas as entradas matriciais
são não-nulas: no caso unidimensional, basta percorrer a malha diretamente de Pi até Pj (andando a partir
de Pi sempre para a direita ou sempre para a esquerda, conforme o caso, até encontrar Pj ), e no caso
bidimensional basta usar qualquer caminho interior de Pi até Pj (pode-se usar a ordem lexicográfica para
percorrer a malha, ou a ordem lexicográfica inversa, dependendo das posições relativas de Pi e Pj ; no entanto,
estes caminhos são mais longos que o necessário). Em outras palavras, identificando as malhas de pontos
internos com os grafos direcionados da matriz de discretização, de modo que existe um arco direcionado entre
dois pontos da malha se e somente se eles são vizinhos, os esquemas de discretização considerados garantem
que estes grafos são fortemente conexos.
As matrizes obtidas através de diferenças finitas em geral são irredutı́veis, pois elas satisfazem a proprie-
dade FC. É difı́cil imaginar um esquema de diferenças finitas para uma malha sobre um domı́nio conexo em
que não houvesse um caminho direcionado entre pontos vizinhos (isto é, em que tivéssemos aij = 0 para dois
pontos vizinhos Pi e Pj ). Outra maneira de pensar sobre isso é observar que se uma matriz de discretização
fôsse (após permutação de linhas e colunas) da forma
B C
,
0 D
isso implicaria que um conjunto de pontos da malha (os correspondentes ao bloco D) teriam diferenças
finitas independentes do conjunto dos pontos restantes da malha (os correspondentes ao bloco D); pior
ainda, estes últimos poderiam ter diferenças finitas dependentes dos primeiros (já que o bloco C poderia
ser não-nulo). Em última análise, seria possı́vel reduzir o problema de resolver o sistema linear associado à
discretização a dois problemas mais simples. É difı́cil imaginar um esquema de diferenças finitas com esta
propriedade, embora talvez possa ocorrer em algum domı́nio com geometria altamente irregular em que a
malha de pontos interiores se dividisse em essencialmente duas malhas independentes. Tal situação deve ser
evitada com cuidado na hora de discretizar tais regiões.
Além disso, para todas as linhas, excetuando a primeira e as linhas correspondentes a pontos adjacentes à
fronteira do disco temos
n
X 1 ri−1/2 1 ri+1/2 2 1
|aij | = αi + βi + 2δi = 2
+ 2
+ 2 = |aii | .
j=1
∆r ri ∆r ri ri ∆θ2
j6=i
Nestas linhas existe dominância diagonal, enquanto que nas linhas correspondentes a pontos adjacentes à
fronteira do disco temos
(n−1)×m+1
X
|aij | = αi + 2δi < |aii | ,
j=1
j6=i
Rodney Josué Biezuner 71
isto é, temos dominância diagonal estrita. Finalmente, para a primeira linha também temos dominância
diagonal, pois
4
|a00 | = ,
∆r2
(n−1)×m+1
X 2 ∆θ m ∆θ 4
|a0j | = m =4 = = |a00 | .
j=1
π ∆r2 2π ∆r2 ∆r2
j6=0
ut + [F (u)]x = 0, (4.1)
cujas soluções apresentam aspectos ondulatórios e a convecção domina a difusão. Pela regra da cadeia,
[F (u)]x = F 0 (u) vx ,
ut + F 0 (u) ux = 0, (4.2)
quando F é diferenciável, o que assumiremos de agora em diante. Uma condição inicial em t = 0 é assumida.
A esta classe pertence a equação da onda de primeira ordem (equação do transporte) linear com coeficientes
constantes, para a qual o problema de valor inicial homogêneo
ut + cux = 0 se x ∈ R e t ∈ R,
u(x, 0) = f (x) se x ∈ R,
Para equações do tipo (4.1) e mesmo (4.2), em geral temos resultados de existência local e a unicidade não é
garantida sem hipóteses adicionais. Em muitas situações temos a presença de ondas de choque, o que torna
bastante complexa até mesmo a definição do que se quer dizer por soluções da equação.
72
Rodney Josué Biezuner 73
Outro exemplo de (4.2) é a equação cinemática da onda, que descreve fenômenos não lineares em dinâmica
dos fluidos em que os efeitos dissipativos tais como viscosidade e difusão são ignorados. Ela é dada por
ut + c(u)ux = 0,
onde u representa densidade e c é uma função de classe C 1 dada. O termo c(u) significa que a velocidade
do fluido em um ponto depende exclusivamente de sua densidade naquele ponto.
Um caso especial importante é a equação de Burgers (invı́scida) que aparece no estudo da dinâmica dos
gases, dada por
ut + uux = 0,
quando, mais uma vez, ignora-se os efeitos dissipativos (note que F (u) = u2 /2). Neste caso, a velocidade do
fluido em cada ponto é diretamente proporcional à densidade do fluido no ponto.
Quando se leva em conta efeitos dissipativos, obtemos a equação de Burgers viscosa
ut + uux = εuxx ,
onde assume-se ε > 0 pequeno. A existência do pequeno elemento difusivo (dissipativo) permite a dissipação
de eventuais choques.
Referências para os métodos numéricos que trataremos neste capı́tulo são [Strikwerda], [Thomas1] e
[Thomas2].
xk = k∆x,
tn = n∆t,
vkn = v (xk , tn )
Esquema FTFS. No esquema FTFS (forward time, forward space), usamos diferenças progressivas para
aproximar ambas as derivadas primeiras:
un+1 − unk un − unk
k
+ c k+1 = 0,
∆t ∆x
ou seja,
∆t n
ukn+1 = unk − c uk+1 − unk .
(4.3)
∆x
Rodney Josué Biezuner 74
Esquema FTBS. No esquema FTBS (forward time, backward space), usamos uma diferença progressiva
para aproximar a derivada primeira em relação ao tempo e uma diferença regressiva para aproximar a
derivada primeira em relação ao espaço:
Todos os esquemas acima, com exceção do esquema Pulo do Sapo, são esquemas de passo único, isto é,
envolvem v apenas em dois nı́veis, n e n + 1. Dado o valor inicial u0k , que vem diretamente da condição
inicial do problema, todos os valores unk podem ser obtidos para todos os valores de n. O esquema Pulo do
Sapo é um esquema multipasso, pois envolve três nı́veis diferentes, n − 1, n e n + 1. Neste caso, o valor de
u1k também deve ser especificado ou calculado através de um esquema de passo único e então os valores unk
para n > 2 podem ser obtidos através do esquema. Algumas das propriedades do esquema Pulo do Sapo são
independentes da forma escolhida para inicializá-lo, outras não.
Definir esquemas de diferenças finitas é fácil. Já a sua análise, para determinar se eles são úteis para
obter aproximações numéricas da solução exata da equação requer ferramentas matemáticas sofisticadas,
Rodney Josué Biezuner 75
como veremos. Veja os Exemplos 1.3.1 e 1.3.2 em [Strikwerda], pp. 18–20, onde vemos entre outras coisas
que a convergência de um esquema de diferenças finitas depende crucialmente do valor de
∆t
λ= .
∆x
Dependendo do valor de λ um dado método pode ser convergente ou não, quando se faz ∆x, ∆t → 0.
u = (uk ) = (. . . , u−1 , u0 , u1 , . . .) .
Uma solução un para o esquema de diferenças finitas no instante de tempo discretizado n∆t é então denotada
por
un = . . . , un−1 , un0 , un1 , . . . .
e, similarmente,
Fn = . . . , F−1
n
, F0n , F1n , . . . .
Denotaremos um subespaço qualquer não especificado do espaço de sequências reais por `. Nestas notas,
consideraremos o espaço das sequências limitadas
`∞ = u : Z −→ R : sup |uk | < ∞ ,
Z
com norma
kuk∞ = sup |uk | .
Z
e v
u ∞
u X
2
kuk2,∆x = t |uk | ∆x.
k=−∞
Rodney Josué Biezuner 76
A necessidade de se considerar esta última norma é devida ao fato de que a norma `2 diverge para ∞ quando
∆x → 0 no seguinte sentido. Seja v ∈ L2 (R) uma função contı́nua e dado ∆x > 0 defina vk = v (k∆x);
denote a sequência resultante por
v∆x = (vk ) ;
em outras palavras, v∆x é a discretização da função v associada à malha uniforme sobre R com comprimento
∆x. É possı́vel provar que v∆x ∈ `2 , qualquer que seja o valor de ∆x. Como o vetor v∆x/2 é duas vezes mais
longo que o vetor v∆x (em um certo sentido, porque tem o dobro de pontos em qualquer intervalo finito;
não se esqueça que ambos são infinitamente longos), temos que
√
v∆x/2 2 ≈ 2 kv∆x k2 .
Iterando esta igualdade, temos que
v∆x/2p 2
≈ 2p/2 kv∆x k2 ,
o que implica que
kv∆x k2 → ∞ quando ∆x → 0.
Por outro lado, é fácil ver que
kvk2,∆x → kvkL2 (R) quando ∆x → 0.
Considere o problema de valor inicial
Lv = F se x ∈ R e t > 0,
v (x, 0) = f (x) se x ∈ R,
e sua correspondente aproximação através de um esquema de diferenças finitas
n n
Lk uk = Fkn k ∈ Z, n ∈ N,
u0k = fk k ∈ Z.
Várias noções de convergência da solução aproximada para a solução exata à medida que a norma da partição
da malha tende a 0 são possı́veis. Podemos considerar uma noção de convergência puntual ou uma noção de
convergência global. Esta última dependerá da norma adotada no espaço de sequências.
é
v (x, t) = f (x − ct) .
2 00 ∞
Consequentemente, se f ∈ C (R) e f ∈ L (R), então
onde
1 n ∆t n
Lnk unk = un+1 uk+1 + unk−1 − c uk+1 − unk−1
k −
2 ∆x
sujeito à restrição
∆t
c 6 1.
∆x
4.4 Teorema. Suponha que
∆t
ε0 6 c 61
∆x
para algum ε0 > 0. Se f ∈ C 2 (R) e f 00 ∈ L∞ (R), então o esquema de Lax-Friedrichs converge uniforme-
mente com ordem (1, 1).
Prova. Denotando
∆t
r=c ,
∆x
o esquema de Lax-Friedrichs pode ser escrito na forma
1 n ∆t n 1 n r n
un+1 n n
uk+1 + unk−1 − uk+1 − unk−1
k = uk+1 + uk−1 − c u k+1 − u k−1 =
2 ∆x 2 2
ou
1−r n 1+r n
un+1
k = uk+1 + uk−1 . (4.10)
2 2
Denotemos
zkn = unk − vkn . (4.11)
Rodney Josué Biezuner 78
A solução exata v é de classe C 1 em R × (0, +∞), como vimos no lema anterior. Podemos então usar a
fórmula de Taylor com resto dado pelo teorema do valor médio para escrever
n 1
vkn+1 = vkn + (vt )k ∆t + vtt (k∆x, t1 ) ∆t2 , (4.12)
2
para algum n∆t 6 t1 6 (n + 1) ∆t, e
n n 1
vk+1 = vkn + (vx )k ∆x + vxx (x1 , n∆t) ∆x2 ,
2
n 1
vk−1 = vk − (vx )k ∆x + vxx (x2 , n∆t) ∆x2 ,
n n
2
para alguns k∆x 6 x1 6 (k + 1) ∆x, (k − 1) ∆x 6 x2 6 k∆x, donde
n n 1
vk+1 + vk−1 = 2vkn + (vxx (x1 , n∆t) + vxx (x2 , n∆t)) ∆x2 , (4.13)
2
n n n 1
vk+1 − vk−1 = 2 (vx )k ∆x + (vxx (x1 , n∆t) − vxx (x2 , n∆t)) ∆x2 , (4.14)
2
Segue de (4.12) e (4.13) que
1 n n 1 1
vkn+1 = n
+ (vt )k ∆t + vtt (k∆x, t1 ) ∆t2 − (vxx (x1 , n∆t) + vxx (x2 , n∆t)) ∆x2 .
vk+1 + vk−1
2 2 4
Daı́, o fato de que v é a solução exata para a equação da advecção e (4.14) implicam que
1 n n 1 1
vkn+1 = n
− c (vx )k ∆t + vtt (k∆x, t1 ) ∆t2 − (vxx (x1 , n∆t) + vxx (x2 , n∆t)) ∆x2
vk+1 + vk−1
2 2 4
1 n n c 1
vk+1 − vk−1 − (vxx (x1 , n∆t) − vxx (x2 , n∆t)) ∆x2 ∆t
n n
= vk+1 + vk−1 −
2 2∆x 2
1 1
+ vtt (k∆x, t1 ) ∆t2 − (vxx (x1 , n∆t) + vxx (x2 , n∆t)) ∆x2
2 4
1 n n
c∆t n n
c
= vk+1 + vk−1 − vk+1 − vk−1 − (vxx (x1 , n∆t) − vxx (x2 , n∆t)) ∆x∆t
2 2∆x 4
1 1
+ vtt (k∆x, t1 ) ∆t − (vxx (x1 , n∆t) + vxx (x2 , n∆t)) ∆x2 ,
2
2 4
ou seja,
1 n r n c
vkn+1 = v n
+ vk−1 − v n
− vk−1 − (vxx (x1 , n∆t) − vxx (x2 , n∆t)) ∆x∆t
2 k+1 2 k+1 4
1 1
+ vtt (k∆x, t1 ) ∆t2 − (vxx (x1 , n∆t) + vxx (x2 , n∆t)) ∆x2 .
2 4
Usando o lema anterior, podemos escrever de forma mais compacta
1−r n 1+r n
vkn+1 = vk−1 + O ∆t2 + O (∆t∆x) + O ∆x2 ,
vk+1 + (4.15)
2 2
onde as constantes nos termos de ordem O independem de k e n se f 00 ∈ L∞ (R). Segue que
1−r n 1+r n
zkn+1 = zk−1 + O ∆t2 + O (∆t∆x) + O ∆x2 .
zk+1 + (4.16)
2 2
Como |r| 6 1 por hipótese, segue que 1 ± r > 0 e podemos escrever
1−r n 1+r n
zkn+1 6 zk−1 + O ∆t2 + O (∆t∆x) + O ∆x2
zk+1 +
2 2
1−r n 1+r n
kz k∞ + O ∆t2 + O (∆t∆x) + O ∆x2
6 kz k∞ +
2 2
6 kzn k∞ + C ∆t2 + ∆t∆x + ∆x2 ,
Rodney Josué Biezuner 79
pois
zk0 = u0k − vk0 = f (k∆x) − v (k∆x, 0) = 0.
Portanto, para todos k, n, temos
∆x2
un+1
k − v (k∆x, (n + 1) ∆t) 6 (n + 1) ∆tC ∆t + ∆x +
∆t
c
= (n + 1) ∆tC ∆t + ∆x + ∆x .
r
→0
quando ∆t, ∆x → 0 e (n + 1) ∆t → t, porque
c 1
6 .
r ε0
Provar diretamente que um esquema de diferenças finitas é convergente é difı́cil em geral. Duas pro-
priedades que são mais fáceis de verificar, a consistência e a estabilidade de esquemas numéricos, podem
ser usadas para provar indiretamente que um determinado esquema é convergente. A seguir definiremos e
estudaremos estas duas propriedades.
4.5 Definição. Dizemos que o esquema de diferenças finitas Lnk unk = Fkn é pontualmente consistente
com a equação diferencial parcial Lv = F se para todo (x, t) ∈ R × (0, +∞) e para toda função ϕ ∈
C ∞ (R × (0, +∞)) temos
(P ϕ) (xk , tn ) − (Pkn ϕ) (xk , tn ) → 0
quando ∆x, ∆t → 0 e (k∆x, (n + 1) ∆t) → (x, t).
Ele é uniformemente consistente se
kP ϕ − Pkn ϕk∞ → 0
é puntualmente consistente.
Prova. O operador diferencial parcial é
P ϕ = ϕt + cϕx
enquanto que o operador de diferenças para o esquema FTFS é o operador
n 1 n
ϕn+1 = ϕnk + (ϕt )k ∆t + (ϕtt )k ∆t2 + O ∆t3 ,
k
2
n 1 n
ϕnk+1 n
= ϕk + (ϕx )k ∆x + (ϕxx )k ∆x2 + O ∆x3 .
2
Logo,
n 1 n n c n
(Pkn ϕ) (xk , tn ) = (ϕt )k + (ϕtt )k ∆t + O ∆t2 + c (ϕx )k + (ϕxx )k ∆x + O ∆x2
2 2
n n 1 n c n
= (ϕt )k + c (ϕx )k + (ϕtt )k ∆t + (ϕxx )k ∆x + O ∆t2 + O ∆x2 .
2 2
Daı́,
1 n c n
(P ϕ) (xk , tn ) − (Pkn ϕ) (xk , tn ) = (ϕtt )k ∆t + (ϕxx )k ∆x + O ∆t2 + O ∆x2 .
2 2
Portanto,
|(P ϕ) (xk , tn ) − (Pkn ϕ) (xk , tn )| = O (∆t) + O (∆x) .
Rodney Josué Biezuner 81
4.5 Estabilidade
Estabilidade em geral quer dizer que pequenas diferenças nas condições iniciais causam pequenos erros na
solução aproximada, um conceito análogo ao de problemas bem postos em EDPs.
4.8 Definição. Um esquema de diferenças finitas é estável se para qualquer T > 0 existe uma constante
C = C (T ) > 0 tais que
kun k 6 C u0
para todos 0 6 n∆t 6 T .
un+1
k = αunk + βunk+1
e
un+1
k = αunk+1 + βunk−1
onde α, β ∈ R são constantes tais que
|α| + |β| 6 1
é estável.
Prova. Temos
∞ ∞
2 X 2 X 2
un+1 = un+1
k = αunk + βunk+1
k=−∞ k=−∞
∞
X 2 2 2 2
6 |α| |unk | + 2 |α| |β| |unk | unk+1 + |β| unk+1
k=−∞
∞ h i
X 2 2 2 2 2 2
6 |α| |unk | + |α| |β| |unk | + unk+1 + |β| unk+1 ,
k=−∞
Portanto,
un+1 6 (|α| + |β|) kun k .
Iterando esta desigualdade, obtemos
n
un+1 6 (|α| + |β|) u0 .
Rodney Josué Biezuner 83
c<0
e
∆t
c > −1.
∆x
Prova. Denotando
∆t
λ=c
∆x
o esquema FTFS é dado por
un+1
k = (1 + λ) unk − λunk+1 .
A condição
|1 + λ| + |λ| 6 1
é equivalente a
−1 6 λ 6 0.
∆t
c 6 1.
∆x
Prova. Denotando
∆t
λ=c ,
∆x
o esquema de Lax-Friedrichs é dado por
1 n λ n
un+1 uk+1 + unk−1 − uk+1 − unk−1
k =
2
2
1−λ n 1+λ
= uk+1 + unk−1 .
2 2
A condição
1−λ 1+λ
+ 61
2 2
é equivalente a
|λ| 6 1.
Rodney Josué Biezuner 84
un+1
k = αunk + βunk+1
tem como domı́nio de dependência um triângulo retângulo de pontos com um cateto vertical, um cateto no
eixo x e a hipotenusa inclinada negativamente. O domı́nio de dependência de um método da forma
un+1
k = αunk−1 + βunk + γunk+1
é a união destes dois triângulos. Ultimamente, é a base destes triângulos que dá o domı́nio de dependência
global, já que todos os segmentos horizontais dos triângulos dependem no final das contas apenas deste último
onde está definida a condição inicial. No primeiro caso, o domı́nio de dependência do ponto (xk , tn+1 ) será
o intervalo
[xk−n−1 , xk ] ,
no segundo caso, o intervalo
[xk , xk+n+1 ]
e no último caso a união destes intervalos, isto é, o intervalo
[xk−n−1 , xk+n+1 ] .
Mais geralmente, dizemos que uma equação diferencial parcial e um esquema de diferenças finitas asso-
ciado satisfaz a condição de Courant-Friedrichs-Lewy quando o domı́nio de dependência analı́tico (isto é, da
Rodney Josué Biezuner 85
solução exata) está contido no domı́nio de dependência numérico. Claramente, uma condição necessária para
que em esquema convirja é que a condição CFL seja satisfeita; caso contrário, como a solução exata em um
ponto P depende do valor inicial em um ponto P0 , se este ponto P0 estiver fora do domı́nio de dependência
numérico, o esquema númérico vai ignorá-lo e não pode produzir uma solução próxima à solução exata em
P , principalmente levando-se em conta que quando ∆x → 0 o domı́nio de dependência numérica ficará cada
vez mais distante do ponto P0 .
No nosso caso, considerando o ponto (xk , tn+1 ), esta condição CFL geral se traduz em que o ponto
x0 = xk − ctn+1
= k∆x − c (n + 1) ∆t
= [k − λ (n + 1)] ∆x,
onde definimos
∆t
λ=c ,
∆x
esteja dentro do intervalo
ou seja,
(k − n − 1) ∆x 6 [k − λ (n + 1)] ∆x 6 (k + n + 1) ∆x.
Daı́, segue que
k − n − 1 6 k − λ (n + 1) 6 k + n + 1,
donde
−n − 1 6 −λ (n + 1) 6 n + 1,
o que implica
−1 6 −λ 6 1,
isto é,
|λ| 6 1.
Uma consequência particular do resultado acima é que qualquer esquema de diferenças finitas explı́cito
para a equação da onda, e para equações hiperbólicas em geral, terá que satisfazer uma condição da forma
∆t
C1 6 6 C2
∆x
e, portanto, não temos a liberdade de escolher ∆t e ∆x independentemente.
Utilizando uma malha uniforme, dividimos o intervalo [0, L] em n subintervalos de mesmo comprimento
∆x = L/n escolhendo os pontos x0 = 0, x1 = ∆x, x2 = 2∆x, ..., xn−1 = (n − 1)∆x, xn = L; similarmente,
dividimos o intervalo [0, T ] em m subintervalos de mesmo comprimento ∆t = T /m escolhendo os pontos
t0 = 0, t1 = ∆t, t2 = 2∆t, ..., tm−1 = (m − 1)∆t, tm = T . Assim,
Em seguida, substituı́mos a equação diferencial parcial parabólica por uma equação parcial de diferenças
finitas. Por exemplo, escolhendo o esquema FTCS, isto é, uma diferença finita progressiva para ut e uma
diferença centrada para uxx , uma aproximação uji para a solução exata u do problema satisfaz a equação de
diferenças parciais
uj+1 − uji uj − 2uji + uji−1
i
= c2 i+1 (4.19)
∆t ∆x2
com condição inicial
u0i = f (i∆x) = f (xi ) = fi para i = 0, . . . , n,
e condição de fronteira
uj0 = ujn = 0 para j = 0, . . . , m.
A equação (4.19) pode ser resolvida em termos de uj+1
i :
∆t j
uj+1
i = uji + c2 ui+1 − 2uj
i + uj
i−1 . (4.20)
∆x2
Denotando
∆t
s = c2 , (4.21)
∆x2
ela também pode ser escrita na forma
uj+1
i = (1 − 2s)uj
i + s uj
i+1 + uj
i−1 . (4.22)
Este esquema é um esquema explı́cito, pois podemos obter o valor de u no instante de tempo j + 1 explici-
tamente em função do valor de u no instante de tempo j.
Na verdade, não é preciso considerar condições de fronteira homogêneas. Se as extremidades têm tem-
peraturas controladas por funções dadas, o método de diferenças finitas funciona do mesmo jeito. De fato,
considere o problema geral 2
ut = c uxx
(x, t) ∈ [0, L] × [0, T ],
u(0, t) = a(t) t ∈ [0, T ],
(4.23)
u(L, t) = b(t) t ∈ [0, T ],
u(x, 0) = f (x) x ∈ [0, L].
Então,
uj0 = a(j∆t) = a(tj ) = aj para j = 0, . . . , m,
(4.24)
ujn = b(j∆t) = b(tj ) = bj para j = 0, . . . , m.
Rodney Josué Biezuner 87
4.8.1 Consistência
Temos
j j j
uj+1 − uji u − 2ui + ui−1
ut (i∆x, j∆t) − c2 uxx (i∆x, j∆t) = i
− c2 i+1 + O(∆t) + O(∆x2 ),
∆t ∆x2
o que significa que a equação em diferenças finitas (4.19) aproxima a equação diferencial do calor na primeira
ordem em t e na segunda ordem em x. Isso não diz nada sobre quão bem a solução da equação em diferenças
finitas aproxima a solução da equação diferencial parcial.
É possı́vel obter um esquema de diferenças finitas para a equação do calor que é de segunda ordem
também em t. Basta usar uma diferença centrada também para ut :
j j j
ui,j+1 − ui,j−1 u − 2ui + ui−1
= c2 i+1
2k ∆x2
de modo que
∆t j
uj+1
i = uj−1
i + 2c2 u i+1 − 2uj
i + u j
i−1 .
∆x2
Este esquema é chamado o esquema leapfrog. No entanto, ele é um esquema de três nı́veis no tempo. Para
iniciá-lo, pode-se usar no primeiro passo o esquema de primeira ordem obtido anteriormente, e os valores de
primeira ordem do primeiro passo não contaminarão a solução, de modo que a aproximação ainda será de
segunda ordem, como veremos mais tarde.
Usaremos a seguinte notação:
segue que
j j j
!
j+1 j
j u − u u − 2u + u
ut − c2 uxx i − i i
− c2 i+1 i i−1
∆t ∆x2
∆t2 2
2∆x4
∆t 4 j j 2 2∆x j j
= − c (uxxxx )i − (uttt )i − ... + c (uxxxx )i + (uxxxxxx )i + ...
2 3! 4! 6!
2
∆t2 4
2 2∆x ∆t 4 j j 2 2∆x j
= c − c (uxxxx )i − (uttt )i − ... + c (uxxxxxx )i + ...
4! 2 3! 6!
∆x2
∆t 4 j
= c2 − c (uxxxx )i + O(∆t2 ) + O(∆x4 ),
4! 2
de modo que se escolhermos
∆x2
∆t = ,
6c2
obteremos !
j j j
j uj+1 − uji u − 2ui + ui−1
2
ut − c uxx − i
− c2 i+1 = O(∆t2 ) + O(∆x4 ).
i ∆t ∆x2
j uj1 − uj−1
0= (ux )0 = ,
2∆x
Rodney Josué Biezuner 89
donde
uj−1 = uj1 , (4.29)
e daı́, usando o nosso esquema de diferenças finitas,
∆t j
uj+1
0 = uj0 + c2 u1 − 2uj
0 + uj
−1
∆x2
segue que
∆t j
uj+1
0 = uj0 + c2 u 1 − u j
0 . (4.30)
∆x2
com u00 obviamente dado pela condição inicial.
Neste capı́tulo, estudaremos esquemas de diferenças finitas para a equação da difusão não-homogênea
como modelo para equações parabólicas unidimensionais. Consideraremos o problema de valor inicial
vt − Kvxx = F (x, t) se x ∈ R e t > 0,
(4.32)
v (x, 0) = f (x) se x ∈ R,
e o problema de valor inicial e de valor de fronteira
vt − Kvxx = 0 se x ∈ [0, L] e t > 0,
v (x, 0) = f (x) se x ∈ [0, L] ,
(4.33)
a1 (t) v (0, t) + a2 (t) vx (0, t) = 0 se t > 0,
b1 (t) v (L, t) + b2 (t) vx (L, t) = 0 se t > 0.
L = ∂t − K∂xx (4.34)
Lv = F (4.35)
onde
∆t
Lnk unk = un+1 n n n n
− uk + K u − 2u + u
k
∆x2 k−1 k k+1
sujeito à restrição
∆t 1
K 6 .
∆x2 2
Precisaremos antes do seguinte resultado teórico:
possui uma solução v ∈ C 0 (R × [0, +∞)) ∩ C ∞ (R × (0, +∞)), e existe uma constante C > 0 tal que para
todo x ∈ R temos
C
|vtt (x, t)| 6 kf kL∞ (R) ,
t3/2
C
|vxxxx (x, t)| 6 kf kL∞ (R) .
t3/2
Se, além disso, f ∈ C 4 (R) e f (4) ∈ L∞ (R), então existe uma constante C > 0 tal que para todo x ∈ R e
para todo t > 0 temos
Prova. Através da teoria analı́tica para a equação do calor na reta, sabe-se que a solução para este problema
é dada pela convolução com o núcleo do calor
1 x2
Φ (x, t) = √ e− 4Kt ,
4πKt
isto é, Z
v (x, t) = Φ (x − s, t) f (s) ds.
R
para todo x ∈ R. Por simplicidade, e sem perda de generalidade, vamos tomar K = 1. Por indução, as
derivadas parciais de Φ são dadas por
P (x, t) − x2
Φxp tq (x, t) = e 4t
Q (t)
Rodney Josué Biezuner 91
√
aij xi tj e Q (t) = ak πtk/2 . Por exemplo,
P
para alguns polinômios algébricos P (x, t) =
−x − x2
Φx = √ e 4t ,
2t 4πt
x2 − 2t − x2
Φt = √ e 4t ,
8 πt5/2
x2 − 2t − x2
Φxx = √ e 4t ,
8 πt5/2
x4 − 12x2 t + 12t2 − x2
Φtt = √ e 4t ,
32 πt9/2
−x3 + 6xt − x2
Φxt = √ e 4t .
16 πt7/2
e assim provar que v é de classe C ∞ em R × (0, +∞); usando o fato que f é contı́nua, também pode-se provar
que u é contı́nua até a fronteira (para detalhes, veja algum livro-texto de EDPs analı́ticas).
Agora, observe que dado β > 0 existe uma constante C = C (β) tal que
z
z β e−z 6 Ce− 2 para todo z > 0;
z z z
basta escrever z β e−z = z β e− 2 e− 2 e usar o fato que a exponencial negativa e− 2 tende a zero mais
z
rapidamente que o polinômio z β quando z → +∞, logo a função z β e− 2 assume um máximo em [0, +∞).
Usaremos esta propriedade da seguinte forma. Para cada t > 0, temos
x2
Pela propriedade mencionada, tomando z = − e β = 2, 1, 0 em cada um dos três termos do lado direito
4t
desta equação, obtemos
C x2 C x2
|Φtt (x, t)| 6 √ 5/2 e− 8t = √ e− 8t ,
2 πt t 3/2 8πt
de modo que
C
|Φtt (x, t)| 6 Φ (x, 2t) (4.38)
t3/2
para todo x ∈ R. Daı́,
Z Z Z
C
|vtt (x, t)| = Φtt (x − s, t) f (s) ds 6 kf kL∞ (R) |Φtt (x − s, t)| ds = 3/2 kf kL∞ (R) Φ (x, 2t) ds
R R t R
C
6 3/2 kf kL∞ (R) .
t
Rodney Josué Biezuner 92
Similarmente, obtemos
x4 − 12x2 t + 12t2 − x2
Φxxxx (x, t) = √ e 4t = Φtt (x, t) , (4.39)
32 πt9/2
logo o mesmo tipo de argumento produz a estimativa para vxxxx .
Suponha agora que f ∈ C 4 (R) e f (4) ∈ L∞ (R). Usando a mudança de variáveis r = x − s, podemos
escrever Z Z
v (x, t) = Φ (x − s, t) f (s) ds = Φ (r, t) f (x − r) dr.
R R
Derivando sob o sinal de integração, segue que
Z
vxxxx (x, t) = Φ (r, t) f (4) (x − r) dr,
R
logo Z
|vxxxx (x, t)| 6 f (4) Φ (r, t) dr = f (4) .
L∞ (R) R L∞ (R)
4.15 Teorema. Suponha que
∆t 1
2
K
6 .
∆x 2
Se f ∈ L∞ (R), então o esquema F T CS converge espacialmente uniformemente com precisão de ordem
(2, 1). Se f ∈ C 4 (R) ∩ L∞ (R) e f (4) ∈ L∞ (R), então o esquema F T CS converge uniformemente com
precisão de ordem (2, 1).
Prova. Denotando
∆t
r=K ,
∆x2
o esquema F T CS pode ser escrito na forma
Denotemos
zkn = unk − vkn . (4.41)
∞
A solução exata v é de classe C em R × (0, +∞), como vimos no lema anterior. Podemos então usar a
fórmula de Taylor com resto dado pelo teorema do valor médio e o fato que v é a solução exata para a
equação da difusão para escrever
n 1
vkn+1 = vkn + (vt )k ∆t + vtt (k∆x, t1 ) ∆t2
2
n 1
n
= vk + K (vxx )k ∆t + vtt (k∆x, t1 ) ∆t2 ,
2
para algum n∆t 6 t1 6 (n + 1) ∆t. Usando novamente a fórmula de Taylor com resto dado pelo teorema do
valor médio, segue também que
n n 1 n 1 n 1
vk+1 = vkn + (vx )k ∆x + (vxx )k ∆x2 + (vxxx )k ∆x3 + vxxxx (x1 , n∆t) ∆x4 ,
2 3! 4!
n n 1 n 1 n 1
vk−1 = vkn − (vx )k ∆x + (vxx )k ∆x2 − (vxxx )k ∆x3 + vxxxx (x2 , n∆t) ∆x4 ,
2 3! 4!
Rodney Josué Biezuner 93
n n n 1
vk−1 + vk+1 = 2vkn + (vxx )k ∆x2 + (vxxxx (x1 , n∆t) + vxxxx (x2 , n∆t)) ∆x4 ,
4!
e daı́,
n 1 1
(vxx )k = 2
−2vkn + vk+1
n n
+ vk−1 − (vxxxx (x1 , n∆t) + vxxxx (x2 , n∆t)) ∆x2 .
∆x 4!
Portanto,
1 1
vkn+1 = (1 − 2r) vkn + r vk−1
n n
+ vk+1 + vtt (k∆x, t1 ) ∆t2 − (vxxxx (x1 , n∆t) + vxxxx (x2 , n∆t)) ∆t∆x2 .
2 4!
Usando o lema anterior, podemos escrever de forma mais compacta
onde as constantes nos termos de ordem O independem de k e n se f ∈ C 4 (R) ∩ L∞ (R) e f (4) ∈ L∞ (R).
Segue que
zkn+1 = (1 − 2r) zkn + r zk−1
n n
+ O ∆t2 + O ∆t∆x2 .
+ zk+1 (4.43)
1
Como 0 < r 6 2 por hipótese, segue que 1 − 2r > 0 e podemos escrever
pois
zk0 = u0k − vk0 = f (k∆x) − v (k∆x, 0) = 0.
Portanto, para todos k, n, temos
quando ∆t, ∆x → 0 e (n + 1) ∆t → t.
Capı́tulo 5
Ax = b.
Embora a matriz A que temos em mente é em geral uma matriz grande e esparsa, do tipo que aparece
em esquemas de diferenças finitas, os métodos considerados aqui requerem apenas que A seja uma matriz
invertı́vel com todas as entradas diagonais aii não-nulas.
Métodos iterativos requerem um chute inicial x0 , um vetor inicial que aproxima a solução exata x (se
não há nenhuma informação disponı́vel sobre a solução exata, de modo que não temos como construir o
chute inicial de forma inteligente, x0 pode ser uma aproximação muito ruim de x). Uma vez que x0 é dado,
o método iterativo gera a partir de x0 uma nova aproximação x1 , que esperamos deve aproximar melhor a
solução exata. Em seguida, x1 é usada para geraruma nova melhor aproximação x2 e assim por diante.
Desta forma, gera-se uma seqüência de vetores xk que espera-se convergir para x. Como na prática não
podemos iterar para sempre, algum critério de parada deve ser estabelecido a priori. Uma vez que xk esteja
suficientemente próximo da solução exata quanto se precise, de acordo com uma margem de tolerância aceita,
pára-se o processo de iteração e aceita-se xk como a solução aproximada adequada para o problema. Por
exemplo, o critério de parada pode ser estabelecido através de uma cota de tolerância τ : quando
b − Axk < τ
ou quando
xk+1 − xk < τ
as iterações são interrompidas e o último valor aproximado obtido é aceito como a melhor aproximação da
solução dentro das circunstâncias.
Os métodos discutidos neste capı́tulo não necessitam de um bom chute inicial (embora, é claro, quanto
melhor o chute inicial, menor o número de iterações necessárias para se chegar à solução aproximada com a
precisão especificada).
94
Rodney Josué Biezuner 95
se aii 6= 0 para todo i, cada xi pode ser isolado na i-ésima equação e escrito na forma
n
1 bi −
X
xi = aij xj
.
aii
j=1
j6=i
Isso sugere definir um método iterativo da seguinte forma: suposto xk = xk1 , . . . , xkn obtido no passo
anterior, obtemos xk+1 = xk+1
1 , . . . , xk+1
n por
n
1 X
xk+1 aij xkj
= bi −
i
aii . (5.1)
j=1
j6=i
No caso da fórmula de cinco pontos para o problema de Poisson com ∆x = ∆y, como a equação para
cada ponto (i, j) é dada por
o método de Jacobi é
1 k
uk+1 ui,j−1 + uki,j+1 + uki−1,j + uki+1,j + ∆x2 fi,j .
i,j = (5.2)
4
No caso especial da equação de Laplace (f = 0) com condição de fronteira de Dirichlet não-nula, o método
de Jacobi é simplesmente a propriedade do valor médio discreta
1 k
uk+1 ui,j−1 + uki,j+1 + uki−1,j + uki+1,j .
i,j = (5.3)
4
Em outras palavras, calculados os valores de u em todos os pontos da malha na iteração anterior, o novo
valor de u em um ponto interior da malha nesta iteração é calculado através da média dos seus quatro
pontos vizinhos. Os valores iniciais de u nos pontos interiores da malha para a primeira iteração (isto é, o
chute inicial) podem ser atribuidos arbitrariamente ou através de algum argumento razoável; por exemplo,
podemos utilizar uma média ponderada dos valores de fronteira para o valor inicial em cada ponto interior
da malha, de acordo com a posição do ponto em relação 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 são as entradas diagonais de A, temos que
xk+1 = D−1 (D − A) xk + b
(5.4)
ou
xk+1 = D−1 Cxk + b
(5.5)
onde C = D − A é a matriz consistindo dos elementos restantes de A fora da diagonal principal.
Rodney Josué Biezuner 96
Então definimos
i−1 n
1 X X
xk+1
i = bi − aij xk+1
j + aij xkj (5.6)
aii j=1 j=i+1
fator em cada passo (se mover um passo na direção de xk para xk+1 é bom, mover naquela direção ω > 1
passos é melhor). Este é o chamado método de sobrerelaxamento sucessivo (SOR, successive overrelaxation):
usando o método de Gauss-Seidel obtemos
i−1 n
1 X X
bk+1
x i = bi − aij xk+1
j + aij xkj ;
aii j=1 j=i+1
daı́ tomamos
xk+1 bk+1
= xki + ω x − xki .
i i
temos
Dxk+1 = Dxk + ω Lxk+1 + (U − D) xk + b
(5.14)
ou
1 k+1 1−ω
D−L x = D + U xk + b,
ω ω
donde −1
k+1 1 1−ω k
x = D−L D+U x +b . (5.15)
ω ω
e tomamos
xk+1 bk+1
= xki + ω x − xki ,
i i
ou seja,
n
1 X
xk+1 k
aij xkj k
bi −
= x + ω − xi . (5.16)
i
i aii
j=1
j6=i
Rodney Josué Biezuner 99
Este método é conhecido como método de Jacobi amortecido, método de Jacobi ponderado ou ainda
método de relaxamento simultâneo (diferente do método de relaxamento sucessivo, baseado no método de
Gauss-Seidel, em que cada variável é substituı́da sucessivamente dentro da mesma iteração à medida que
ela é atualizada; no método de Jacobi, as variáveis são todas substituı́das simultameamente na próxima
iteração).
Em forma matricial, o método de Jacobi amortecido pode ser descrito da seguinte forma. Denotando por
D a parte diagonal de A, temos
Xn
aii xk+1
i = aii xki + ω bi − aij xkj ,
j=1
temos
Dxk+1 = Dxk + ω b − Axk
(5.17)
ou
1 1
D xk+1 = D − A xk + ωb,
ω ω
donde −1
k+1 1 1 k
x = D D−A x +b . (5.18)
ω ω
Em contraste com o método SOR, que converge em geral para 0 < ω < 2, o método de Jacobi amortecido
converge para 0 < ω 6 1 (veja a próxima seção).
ek = x − xk , (5.21)
O erro algébrico tem interesse puramente teórico (para provar que determinado método iterativo converge,
precisamos mostrar que o erro algébrico tende a zero), já que ele só pode ser calculado uma vez que se
conhece a solução exata, e se este for o caso obviamente não há necessidade de resolver o sistema. Já o erro
residual pode ser usado como critério de parada para o método iterativo. Como
segue que
ek+1 = B −1 Cek .
Observe que
B −1 C = B −1 (B − A) = I − B −1 A.
A matriz
R = I − B −1 A = B −1 C (5.23)
é chamada a matriz de iteração ou matriz de propagação do erro do algoritmo considerado, porque
Logo,
n
X
ek = R k e0 = ai λki vi ,
i=1
de modo que
n
X k
ek 6 |ai | |λi | |vi | .
i=1
k
Como |λi | → 0 se e somente se |λi | < 1, concluı́mos que ek → 0 qualquer que seja o erro inicial (isto é,
qualquer que seja o chute inicial), se e somente se ρ (R) = max16i6n |λi | < 1 .
Rodney Josué Biezuner 101
ρ (A) 6 kAk .
Ax = λx.
Considere a matriz X ∈ Mn (C) cujas colunas são todas iguais ao vetor x. Temos também
AX = λX
de modo que
|λ| kXk = kAXk 6 kAk kXk ,
donde
|λ| 6 kAk
para todo autovalor λ de A. Como existe um autovalor λ de A tal que ρ (A) = λ, isso prova o resultado.
5.2 Lema. Seja A ∈ Mn (C) e ε > 0 dado. Então existe uma norma matricial k·k tal que
Prova. Toda matriz complexa é triangularizável através de uma matriz unitária (isto é, uma matriz U que
satisfaz U ∗ U = U U ∗ = I; sua inversa é a sua adjunta ou transposta conjugada). Sejam então
λ1 a12 a22 . . . a1n
λ2 a23 . . . a2n
T =
λ3 . . . a3n
. . ..
. .
λn
A = U ∗ T U.
Dt T Dt−1 L
6 ρ (A) + ε
Rodney Josué Biezuner 102
para t suficientemente grande. Portanto, fixado um tal t, se definirmos uma norma por
−1
kAk := Dt U AU ∗ Dt−1 L
= U ∗ Dt−1 AU ∗ Dt−1 ,
L
teremos
kAk = Dt U AU ∗ Dt−1 L
= Dt T Dt−1 L
6 ρ (A) + ε.
Pelo lema anterior, ρ (A) 6 kAk.
5.3 Lema. Seja A ∈ Mn (C). Se existe alguma norma matricial k·k tal que kAk < 1, então
Ak → 0.
Ak x = λ k x
não converge para 0. Reciprocamente, se ρ (A) < 1, então pelo Lema 5.2 existe uma norma matricial k·k tal
que kAk < 1, logo Ak → 0 pelo lema anterior.
5.5 Corolário. Seja R a matriz de iteração de um método iterativo linear. Então
ek → 0
se e somente se
ρ (R) < 1.
Em outras palavras, um método iterativo linear é convergente independentemente da escolha do chute inicial
se e somente se todos os autovalores da matriz de iteração têm valor absoluto menor que 1.
A = B1 − C1 = B2 − C2 ,
Vamos analisar a velocidade de convergência dos métodos iterativos com maior precisão. Novamente à
tı́tulo de motivação, suponha que R é uma matriz diagonalizável com seu maior autovalor sendo um autovalor
simples. Ordene os autovalores de R na forma
donde
n
X
ek = R k e0 = ai λki vi ,
i=1
segue que " #
n k
k
X λi
e = λk1 a 1 x1 + ai vi .
i=2
λ1
Como k
λi
→ 0,
λ1
k
a taxa de convergência é determinada por |λ1 | . Para k grande, temos
ek ≈ λk1 a1 v1 .
Portanto,
ek+1
= |λ1 | = ρ (R) . (5.28)
|ek |
Em outras palavras, a convergência é linear com taxa de convergência igual ao raio espectral. Se a1 =
0 a convergência será mais rápida, pois dependerá do módulo do segundo autovalor, mas é obviamente
extremamente raro que o chute inicial satisfaça esta condição. Para o caso geral, precisamos do seguinte
resultado:
5.6 Proposição. Seja A ∈ Mn (C) e k·k uma norma matricial. Então
1/k
ρ (A) = lim Ak .
Prova. Como os autovalores da matriz Ak são as k-ésimas potências dos autovalores de A, temos que
k
ρ (A) = ρ Ak 6 Ak ,
donde
1/k
ρ (A) 6 Ak .
Dado ε > 0, a matriz
1
B= A
ρ (A) + ε
tem raio espectral menor que 1, logo B k → 0. Portanto, existe algum N = N (ε, A) tal que
Bk < 1
ou seja,
1/k
Ak < ρ (A) + ε
para todo k > N .
Definimos a taxa média de convergência de um método iterativo linear com matriz de iteração R por
1/k 1
Rk (R) = − log10 Rk = − log10 Rk (5.29)
k
e a taxa assintótica de convergência por
R∞ (R) = lim Rk (R) . (5.30)
k→∞
Rodney Josué Biezuner 104
5.7 Corolário. Seja R a matriz de iteração de um método iterativo linear. Então a taxa assintótica de
convergência do método é dada por
R∞ (R) = − log10 ρ (R) . (5.31)
Prova. Pois
1/k 1/k
R∞ (R) = − lim log10 Rk = − log10 lim Rk = − log10 ρ (R) .
k→∞ k→∞
A taxa assintótica de convergência mede o aumento no número de casas decimais corretas na solução por
iteração. De fato, usando a norma matricial do Lema 5.2 e medindo as normas dos vetores de acordo, temos
ek+1 Rk+1 e0
= 6 kRk = ρ (R) + ε,
|ek | |Rk e0 |
donde
ek+1
− log10 = − log10 ρ (R) + O (ε)
|ek |
ou
log10 ek − log10 ek+1 = R∞ (R) + O (ε) . (5.32)
Assim, se
ek = O 10−p ,
ek+1 = O 10−q ,
teremos
q − p ≈ R∞ (R) ,
isto é, reduzimos R∞ (R) ≈ q − p casas decimais no erro. Visto de outra forma, como
ek+m Rk+m e0 m
k
= 6 kRm k = ρ (R) + O (ε) ,
|e | |Rk e0 |
donde
ek+m
− log10 ≈ −m log10 ρ (R) ,
|ek |
ou
log10 ek+m / ek
m= (5.33)
log10 ρ (R)
é o número de iterações necessárias para diminuir o erro de um número prescrito de casas decimais.
A6B
se
hAx, xi 6 hBx, xi
para todo x ∈ Cn .
Rodney Josué Biezuner 105
Em particular, se A é uma matriz positiva definida, segue que A > εI para algum ε (o menor autovalor de
A) e denotamos este fato por
A > 0.
5.9 Teorema. Seja A uma matriz simétrica positiva definida e seja A = B − C com B invertı́vel. Então o
método iterativo linear com matriz de iteração R = B −1 C converge se e somente se B t + C é uma matriz
simétrica positiva definida.
Prova. Medimos a norma do erro através da norma induzida por A
1/2
|x|A := hAx, xi
e consideraremos a norma matricial k·kA induzida por esta norma. Se provarmos que
kRkA < 1,
C t B −t AB −1 C = B t − A B −t AB −1 (B − A) = I − AB −t A I − B −1 A
= A − AB −t A + AB −1 A − AB −t AB −1 A
= A − AB −t B + B t − A B −1 A
t
= A − B −1 A B + B t − A B −1 A
ou t
C t B −t AB −1 C = A − B −1 A B t + C B −1 A,
(5.35)
de modo que C t B −t AB −1 C é uma matriz simétrica, positiva definida. Logo, por (5.34), mostrar que
kRkA < 1 é equivalente a provar que
C t B −t AB −1 C < A,
t
e por (5.35) C t B −t AB −1 C < A se e somente se B −1 A (B t + C) B −1 A > 0, o que é verdade porque B t +C
é positiva definida.
det I − λ−1 R = 0.
Rodney Josué Biezuner 106
π2
cos (π∆x) ≈ 1 − ∆x2 ,
2
de modo que ρ (R) → 1 quadraticamente quando ∆x → 0. Em outras palavras, para uma malha duas vezes
mais refinada (isto é, ∆x reduzido pela metade), o método de Jacobi é cerca de quatro vezes mais vagaroso
em média (consulte novamente a tabela no final da seção anterior). A tabela abaixo mostra os valores do
raio espectral para alguns valores de ∆x:
de modo que para reduzir o erro pelo fator de uma casa decimal precisamos de
log10 0.1 1 1
m= =− = ≈ 742
log10 ρ (R) log10 ρ (R) 0.00135
iterações.
porque I − D−1 L é uma matriz triangular inferior com apenas 1’s na diagonal principal, escrevemos
h −1 −1 i
0 = det I − λ−1 I − D−1 L D U
h −1 −1 i
= det I − D−1 L det I − λ−1 I − D−1 L
D U
n h −1
io
= det I − D−1 L I − λ−1 I − D−1 L D−1 U
ou seja,
U u = µ (D − L) u
(um problema de autovalor generalizado). No caso da matriz de discretização da fórmula de cinco pontos,
isso significa encontrar µ tal que
para transformar a equação de autovalor naquela que aparece no método de Jacobi. Temos
i+j+1 i+j+1
i+j i+j−1 i+j−1
µ 2 vi,j + µ 2 vi+1,j = µ 4µ 2 vi,j − µ 2 vi,j−1 − µ 2 vi−1,j
i+j+2 i+j+1 i+j+1
= 4µ 2 vi,j − µ 2 vi,j−1 − µ 2 vi−1,j ,
Rodney Josué Biezuner 109
i+j+1
de modo que, dividindo por µ 2 , obtemos
Portanto os autovalores da matriz de iteração de Gauss-Seidel para esta matriz são exatamente os quadrados
dos autovalores da matriz de iteração de Jacobi (e os autovetores são os mesmos):
2
1 kπ lπ
µlk = cos + cos .
4 n n
Portanto, o máximo autovalor ocorre quando k = l = 1 e
π
ρ (R) = cos2 .
n
O argumento para a fórmula de três pontos é análogo.
Para o quadrado unitário temos
ρ (R) = cos2 (π∆x) ,
e usando
2
1
cos2 x = 1 − x2 + O x4 = 1 − x2 + O x4 ,
2
se ∆x é pequeno podemos aproximar
No método de Gauss-Seidel ainda temos ρ (R) → 1 quadraticamente quando ∆x → 0, mas a sua velocidade
de convergência para a matriz de discretização de cinco pontos do quadrado unitário é duas vezes maior que
a do método de Jacobi. Para ver isso, faça a expansão do logaritmo em torno do ponto x = 1:
log (1 + x) = x + O ∆x2 .
Segue que
π2
∆x2 + O ∆x4 ,
R∞ (RJacobi ) = (5.42)
2
R∞ (RGauss-Seidel ) = π 2 ∆x2 + O ∆x4 .
(5.43)
0 < ω < 2.
det R = λ1 . . . λn .
Rodney Josué Biezuner 110
Mas,
n −1 o
det R = det I − ωD−1 L (1 − ω) I + ωD−1 U
−1
= det I − ωD−1 L det (1 − ω) I + ωD−1 U
n
= (1 − ω) ,
já que I − ωD−1 L é uma matriz triangular inferior com apenas 1 na diagonal principal e (1 − ω) I + ωD−1 U
é uma matriz triangular superior com apenas 1 − ω na diagonal principal. Logo
n
λ1 . . . λn = (1 − ω) .
Em particular, pelo menos um dos autovalores λj de R deve satisfazer
|λj | > |1 − ω| .
Mas, se o método SOR converge, devemos ter também |λ| < 1 para todo autovalor λ de R. Logo
|1 − ω| < 1,
donde
0 < ω < 2.
5.15 Corolário. Se R é a matriz de iteração n × n para o método SOR, então
n
det R = (1 − ω) .
Em particular, diferente das matrizes de iteração dos métodos de Jacobi e de Gauss-Seidel (para a matriz de
discretização de cinco pontos), zero não é um autovalor para a matriz de iteração do método SOR se ω 6= 1
(para nenhuma matriz).
n
P
5.16 Teorema. Se A é uma matriz irredutı́vel, diagonalmente dominante tal que |aii | > |aij | para pelo
j=1
j6=i
menos alguma linha i, então o método SOR converge se 0 < ω 6 1.
Prova. A demonstração é análoga à do Teorema 5.12. A matriz de iteração do método SOR é
−1
R = I − ωD−1 L (1 − ω) I + ωD−1 U .
Suponha por absurdo que exista um autovalor λ de R tal que |λ| > 1; temos
n −1 o
det I − λ−1 R = det I − λ−1 I − ωD−1 L (1 − ω) I + ωD−1 U
= 0.
porque I − ωD−1 L é uma matriz triangular inferior com apenas 1’s na diagonal principal, escrevemos
n −1 o
0 = det I − λ−1 I − ωD−1 L (1 − ω) I + ωD−1 U
n −1 o
= det I − ωD−1 L det I − λ−1 I − ωD−1 L (1 − ω) I + ωD−1 U
h n −1 oi
= det I − ωD−1 L I − λ−1 I − ωD−1 L (1 − ω) I + ωD−1 U
= det I − ωD−1 L − λ−1 (1 − ω) I + ωD−1 U
é irredutı́vel, diagonalmente dominante e estritamente dominante nas linhas onde A é, logo a matriz
também satisfaz estas propriedades. De fato, S tem zeros nas mesmas posições que I − D−1 L − D−1 U , logo
a sua irredutibilidade não é afetada. Além disso, pela dominância diagonal de D−1 A, sabemos que se
bij = D−1 L ij ,
cij = D−1 U ij .
então
i−1
X n
X
1> |bij | + |cij | .
j=1 j=i+1
Para provar a dominância diagonal de S, observamos que os valores que S possui na diagonal principal são
1−ω λ+ω−1
1 − λ−1 (1 − ω) = 1 − = ,
λ λ
de modo que precisamos provar que
i−1 n
λ+ω−1 X ω X
>ω |bij | + |cij |
λ j=1
|λ| j=i+1
λ+ω−1
> ω,
λ
λ+ω−1 ω
> .
λ |λ|
Para isso, observe que como |λ| > 1 basta provar a primeira desigualdade, a qual por sua vez é equivalente a
|λ + ω − 1| > |λ| ω.
O resultado acima continua valendo com desigualdade estrita nas linhas onde a desigualdade é estrita. Isso
implica então que S é invertı́vel, contradizendo det S = 0.
5.17 Teorema. Seja A uma matriz simétrica positiva definida. Então o método SOR converge se 0 < ω < 2.
Rodney Josué Biezuner 112
Fazendo o produto interno canônico (hermitiano) de Cn de ambos os lados com o vetor x, segue que
Isolando λ,
(1 − ω) hx, xi + ω x, D−1 U x
λ= . (5.45)
hx, xi − ω hx, D−1 Lxi
Como A é simétrica, o produto de matrizes simétricas D−1 A = I − D−1 U − D−1 L também é; como
D−1 U, D−1 L são respectivamente a parte estritamente triangular superior e estritamente triangular infe-
rior de uma matriz simétrica, temos
t
D−1 U = D−1 L.
Logo D E
t
x, D−1 U x = D−1 U D−1 L x, x = hx, (D−1 L) xi,
x, x =
e definindo
x, D−1 L x
z= ,
hx, xi
podemos escrever
(1 − ω) + ωz
λ= . (5.46)
1 − ωz
Os argumentos acima assumem que o denominador é não-nulo. E, de fato, temos
!
x, D−1 L x x, D−1 U x 1 x, D−1 L + D−1 U x
1 1
Re z = (z + z) = + =
2 2 hx, xi hx, xi 2 hx, xi
!
1 x, I − D−1 A x x, D−1 A x
1
= = 1− .
2 hx, xi 2 hx, xi
Rodney Josué Biezuner 113
x, D−1 A x
>0
hx, xi
donde
1
.
Re z <
2
de modo que a parte real do denominador 1 − ωz de λ é não-nula para 0 < ω < 2. Segue que
2 2
2 [(1 − ω) + ωz] [(1 − ω) + ωz] (1 − ω) + 2ω (1 − ω) Re z + ω 2 |z|
|λ| = λλ = = 2
(1 − ωz) (1 − ωz) 1 − 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 Re z)
=1− 2 .
1 − 2ω Re z + ω 2 |z|
1
Como 0 < ω < 2 e Re z < , temos
2
ω (2 − ω) (1 − 2 Re z) > 0,
e concluı́mos que
|λ| < 1
para todo autovalor λ de R, logo o método SOR converge. A demonstração da recı́proca (assim como uma
demonstração alternativa, variacional, deste teorema) pode ser vista em [Young].
Usando o Teorema 5.16, concluı́mos que o método SOR converge para as matrizes de discretização obtidas
através dos esquemas de diferenças finitas do Capı́tulo 2 se 0 < ω 6 1. Isso permite apenas subrelaxamento
do método de Gauss-Seidel, o que em geral reduz a velocidade de convergência. Por outro lado, usando o
Teorema 5.17 ou o Teorema 5.18, concluı́mos que o método SOR converge para as matrizes de discretização
obtidas a partir da fórmula de três pontos unidimensional e a partir da fórmula de cinco pontos bidimensional
se 0 < ω < 2, já que estas são matrizes simétricas, positivas definidas (já as matrizes de discretização obtidas
através de coordenadas polares ou pelo esquema de Shortley-Weller não são simétricas, em geral, como
vimos).
Em seguida fazemos uma análise da velocidade de convergência do método SOR para a matriz de discre-
tização da fórmula de cinco pontos, bem como obtemos o melhor valor do fator de relaxamento ω para este
caso.
5.19 Lema. Seja A a matriz de discretização obtida a partir da fórmula de três pontos unidimensional ou
a partir da fórmula de cinco pontos bidimensional com ∆x = ∆y. Se λ 6= 0 é um autovalor de RSOR , então
existe um autovalor λJ de RJ tal que
1−ω−λ
λJ = . (5.47)
λ1/2 ω 2
Reciprocamente, se λJ é um autovalor de RJ e λ ∈ C satisfaz a equação acima, então λ é um autovalor de
RSOR .
Prova. Argumentamos como na demonstração do Teorema 5.13. Para obter o raio espectral da matriz de
iteração RSOR , queremos encontrar os autovalores λ de RSOR :
−1
RSOR u = I − ωD−1 L (1 − ω) I + ωD−1 U u = λu,
ou seja,
(1 − ω) I + ωD−1 U u = λ I − ωD−1 L u
Rodney Josué Biezuner 114
No caso da matriz de discretização da fórmula de cinco pontos, isso significa encontrar λ tal que
ω ω ω ω
(1 − ω) ui,j + ui,j+1 + ui+1,j = λ ui,j − ui,j−1 − ui−1,j
4 4 4 4
ou
1−ω−λ 1
ui,j = (ui,j+1 + ui+1,j + λui,j−1 + λui−1,j ) . (5.48)
ω 4
Fazendo a substituição
i+j
ui,j = λ 2 vi,j
i+j+1
e dividindo por µ 2 , segue que
1−ω−λ
vi−1,j + vi+1,j + vi,j−1 + vi,j+1 = 4vi,j
λ1/2 ω
e daı́ o resultado. √ p 2
Resolvendo a equação (5.47) como uma equação quadrática em λ, vemos que as duas raı́zes λ± = λ±
podem ser escritas na forma
2
1
q
λ± = −ωλJ ± ω 2 λ2J − 4 (ω − 1) . (5.49)
4
Denotaremos
Λω,λJ = max (|λ+ | , |λ− |) (5.50)
e por λJ = ρ (RJ ) o maior autovalor do método de Jacobi.
5.20 Proposição. Seja A a matriz de discretização obtida a partir da fórmula de três pontos unidimensional
ou a partir da fórmula de cinco pontos bidimensional com ∆x = ∆y. Então
√
onde i = −1, logo
q 2 r 2
2
h 2
i
Λω,λJ = ωλJ + ω 2 λJ − 4 (ω − 1) = ω 2 λ2J + 4 (ω − 1) − ω 2 λJ
= ω − 1,
Defina
2
ωótimo = q . (5.52)
2
1 + 1 − λJ
Note que 1 < ωótimo < 2. Mostraremos que ωótimo é de fato o melhor valor para o fator de relaxamento no
método SOR. Antes precisamos do seguinte resultado:
5.21 Proposição. Seja A a matriz de discretização obtida a partir da fórmula de três pontos unidimensional
ou a partir da fórmula de cinco pontos bidimensional com ∆x = ∆y. Então
2
1
q
2
ωλ + ω 2 λ − 4 (ω − 1) se 0 < ω 6 ωótimo ,
ρ (RSOR,ω ) = J J (5.53)
4
ω−1 se ωótimo 6 ω < 2.
2
Prova. Temos ω 2 λJ − 4 (ω − 1) > 0 para 0 < ω < 2 se e somente se ω 6 ωótimo . De fato, as raı́zes de
2
f (ω) = ω 2 λJ − 4ω + 4 são q
2
4 ± 4 1 − λJ
q
2 2
ω± = 2 = 2 1 ± 1 − λJ
2λJ λJ
de modo que a raiz positiva de f é maior que 2, logo para que f (ω) > 0 se 0 < ω < 2, devemos ter
2
2
q
2
2 1 − 1 − λ J 2
ω 6 2 1 − 1 − λJ = 2 q = q .
λJ λJ 1 + 1 − λ 2 2
J 1 + 1 − λJ
5.22 Teorema. Seja A a matriz de discretização obtida a partir da fórmula de três pontos unidimensional
ou a partir da fórmula de cinco pontos bidimensional com ∆x = ∆y. Então o fator de relaxamento ótimo
para o método SOR é dado por
2
ωótimo = π (5.54)
1 + sen
n
é o fator de relaxamento ótimo para o método SOR.
2
Prova. Se 0 < ω 6 ωótimo , então ω 2 λJ − 4 (ω − 1) > 0 e
q
2 2
λJ ω 2 λJ − 4 (ω − 1) + ωλJ − 2
q
d 2
ωλJ + ω 2 λJ − 4 (ω − 1) = q .
dω 2
2
ω λJ − 4 (ω − 1)
2
Temos ωλJ − 2 < 0, porque 0 < ω < 2 e λJ < 1, e
q
2 2
ωλJ − 2 > λJ ω 2 λJ − 4 (ω − 1),
pois
2 2 4 2 4 2 2 4 2
ωλJ − 2 = ω 2 λJ − 4λJ ω + 4 > ω 2 λJ − 4λJ ω + 4λJ > ω 2 λJ − 4λJ (ω − 1)
q 2
2
= λJ ω 2 λJ − 4 (ω − 1) .
Rodney Josué Biezuner 116
Isso implica q
d 2
ωλJ + ω 2 λJ − 4 (ω − 1) < 0,
dω
logo ρ (RSOR,ω ) é decrescente de 0 até ωótimo . Para ωótimo 6 ω < 2, ρ (RSOR,ω ) = ω − 1 é claramente
crescente. Portanto, ρ (RSOR,ω ) atinge o seu mı́nimo em ωótimo .
Pelo Teorema 5.11, temos
π
λJ = cos ,
n
logo
2 2 2
ωótimo = q = r = π.
2 π 1 + sen
1 + 1 − λJ 1 + 1 − cos 2
n
n
Para o quadrado unitário temos
2
ωótimo =
1 + sen (π∆x)
e conseqüentemente
2 1 − sen (π∆x)
ρ (RSOR,ω ) = −1= .
1 + sen (π∆x) 1 + sen (π∆x)
e usando
1−x
= 1 − 2x + O x2 ,
1+x
sen x = x + O x3 ,
1 − sen (π∆x)
≈ 1 − 2π∆x + O ∆x2 .
1 + sen (π∆x)
Portanto, usando o valor ótimo de ω no método SOR, temos ρ (R) → 1 linearmente quando ∆x → 0, um
resultado muito melhor que o obtido nos métodos de Jacobi e de Gauss-Seidel. Para uma comparação mais
precisa, usando
log (1 + x) = x + O ∆x2
temos que
R∞ (RSOR ) = 2π∆x + O ∆x2 .
(5.55)
Segue que
R∞ (RSOR ) 2π∆x 2
≈ 2 2 = .
R∞ (RGauss-Seidel ) π ∆x π∆x
Em particular, se ∆x = 0.025, temos ωótimo = 1. 8545 e R∞ (RSOR ) /R∞ (RGauss-Seidel ) = 25.5, isto é, o
método SOR é 25 vezes mais rápido que o método de Gauss-Seidel. Quanto mais refinada a malha, maior é
a diferença na velocidade de convergência entre os dois métodos.
0 < ω 6 1.
Rodney Josué Biezuner 117
Prova. Vamos escrever a matriz de iteração RJ,ω do método de Jacobi amortecido em função da matriz de
iteração do método de Jacobi RJ . Temos
RJ = D−1 (D − A)
de modo que
−1
1 1 1 1
RJ,ω = D D − A = ωD−1 D − D + D − A = ωD−1 D − D + ωD−1 (D − A)
ω ω ω ω
donde
RJ,ω = (1 − ω) I + ωRJ . (5.56)
Em particular,
RJ v = λv
se e somente se
[RJ,ω − (1 − ω) I] v = ωλv.
Portanto, λJ é um autovalor de RJ se e somente se
é um autovalor de RJ,ω . Logo, se todo autovalor de RJ satisfaz |λJ | < 1 (isto é, ρ (RJ ) < 1 equivalente ao
método de Jacobi convergir) e ω < 1, então
2
|λJ,ω | = (ωλJ + 1 − ω) ωλJ + 1 − ω
2 2
= ω 2 |λJ | + 2 Re λJ ω (1 − ω) + (1 − ω)
2 2
6 ω 2 |λJ | + 2 |λJ | ω (1 − ω) + (1 − ω)
2
= (ω |λJ | + 1 − ω)
< 1.
Segue do Teorema 5.9 que o método de Jacobi amortecido converge para as matrizes de discretização do
Capı́tulo 1 se 0 < ω 6 1.
5.24 Corolário.
ρ (RJ,ω ) = ω [ρ (RJ ) − 1] + 1. (5.58)
Para o quadrado unitário temos
ρ (RJ,ω ) = ω [cos (π∆x) − 1] + 1. (5.59)
Usando
1
cos x = 1 − x2 + O x4 ,
2
log (1 + x) = x + O ∆x2 ,
5.3.5 Resumo
Método ρ (R) R∞ (R)
π2
Jacobi cos (π∆x) ∆x2 + O ∆x4
2
Gauss-Seidel cos2 (π∆x) π 2 ∆x2 + O ∆x4
SOR ótimo 1 − 2π∆x + O ∆x2 2π∆x + O ∆x2
π2 π2
Jacobi amortecido 1−ω ∆x2 + O ∆x4 ω ∆x2 + O ∆x4
2 2
Ax = b
Daı́,
1 t 1
f (y) − f (x) = y Ay − y t b − xt Ax + xt b
2 2
1 t 1
= y Ay − y t Ax − xt Ax + xt Ax
2 2
1 t 1
= y Ay − y t Ax + xt Ax
2 2
1 t 1 1 1
= y Ay − y t Ax − xt Ay + xt Ax
2 2 2 2
1 t 1 t
= y A (y − x) − x A (y − x)
2 2
ou
1 t
f (y) − f (x) = (y − x) A (y − x) . (5.61)
2
Como A é positiva definida, segue que
t
(y − x) A (y − x) = hA (y − x) , (y − x)i > 0
e
t
(y − x) A (y − x) = 0
Rodney Josué Biezuner 119
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 que
quando é minimizado corresponde a um estado de equilı́brio do sistema. Observe que definindo um produto
interno a partir da matriz simétrica positiva definida A da maneira usual por hv, wiA = v t Aw e considerando
1/2
a norma induzida kvkA = hv, viA , o funcional f pode ser escrito na forma
1
f (y) = hy, Ayi − hy, Axi (5.62)
2
ou
1 2
kykA − hy, xiA .
f (y) = (5.63)
2
Outra maneira de enxergar o resultado do teorema anterior é observar que o gradiente do funcional f é
∇f (y) = Ay − b. (5.64)
Se x é um ponto de mı́nimo temos ∇f (x) = 0, ou seja,
Ax = b.
Este método variacional é a base dos métodos iterativos de descida em geral, e do método do gradiente
conjugado em particular. A idéia é usar as idéias do cálculo diferencial para encontrar o mı́nimo do funcional
quadrático f .
de tal modo que xk convirja para o minimizador de f . Em outras palavras, em um método de descida
buscamos encontrar uma seqüência minimizante xk que convirja para a solução do sistema.
O passo de xk para xk+1 envolve dois ingredientes: (1) uma direção de busca e (2) um avanço de
comprimento especificado na direção de busca. Uma direção de busca significa a escolha de um vetor pk que
indicará a direção que avançaremos de xk para xk+1 . O comprimento do avanço é equivalente à escolha de
um escalar αk multiplicando o vetor pk . Assim,
xk+1 = xk + αk pk .
A escolha de αk é também chamada uma busca na reta, já que queremos escolher um ponto na reta
k
x + αpk : α ∈ R
tal que
f xk + αpk 6 f xk .
Esta é chamada uma busca na reta exata. Para funcionais quadráticos, a busca na reta exata é trivial e
obtemos uma fórmula para o valor de αk , como veremos a seguir. Denotaremos o resı́duo em cada iteração
por
rk = b − Axk . (5.65)
Rodney Josué Biezuner 120
f xk + αk pk = min f xk + αpk .
α∈R
Então t
pk rk pk , rk
αk = t = . (5.66)
(pk ) Apk hpk , Apk i
Prova: Considere o funcional
g (α) = f xk + αpk .
pk = −∇f xk = b − Axk
ou
pk = rk . (5.67)
Buscar na direção da descida mais acentuada é uma idéia natural, mas que na prática não funciona sem
modificações. De fato, em alguns casos o método é de velocidade comparável à do método de Jacobi, como
na matriz de discretização da fórmula de cinco pontos aplicada ao problema descrito na primeira seção deste
capı́tulo [Watkins]:
De fato, como as iterações do método de descida mais acentuada são bem mais custosas que as do método
de Jacobi, o primeiro é muito pior que este último.
Para entender melhor o método da descida mais acentuada, porque ele pode ser lento e as modificações que
vamos fazer para torná-lo mais rápido levando ao método do gradiente conjugado, vamos entender o processo
do ponto de vista geométrico. Como vimos na demonstração do Teorema 5.25, o funcional quadrático f é
da forma
1 t
f (y) = (y − x) A (y − x) + c (5.68)
2
onde c = f (x) = 21 xt Ax − xt b é uma constante. Já que A é uma matriz simétrica, existe uma matriz
ortogonal P tal que P t AP é uma matriz diagonal D , cujos valores na diagonal principal são exatamente os
autovalores positivos de A. Nas coordenadas
z = P t (y − x) ,
As curvas de nı́vel do funcional f neste sistema de coordenadas são elipses (em R2 , elipsóides em R3 e
hiperelipsóides em Rn ) centradas na origem com eixos paralelos aos eixos coordenados e f (0) = c é nı́vel
mı́nimo de f ; elipses correspondentes a menores valores de f estão dentro de elipses correspondentes a
maiores valores de f . Como P é uma aplicação ortogonal, as curvas de nı́vel de f no sistema de coordenadas
original também são elipses, centradas em x, e uma reta de um ponto y até o ponto x corta elipses de nı́veis
cada vez menores até chegar ao mı́nimo da função f em x, centro de todas as elipses. O vetor gradiente é
perpendicular às curvas de nı́vel, logo é perpendicular às elipses. Seguir a direção de descida mais acentuada
equivale a cortar a elipse que contém xk ortogonalmente na direção do interior da elipse até encontrar um
ponto xk+1 situado em uma elipse que a reta tangencie, pois a partir daı́ a reta irá na direção de elipses com
nı́veis maiores, portanto este é o ponto da reta onde f atinge o seu mı́nimo. Em particular, vemos que a
próxima direção pk+1 é ortogonal à direção anterior pk , tangente a esta elipse. Em geral, a direção de descida
mais acentuada não é a direção de x (quando bastaria uma iteração para atingir a solução exata) a não ser
que A seja um múltiplo escalar da identidade, de modo que todos os autovalores de A são iguais e as elipses
são cı́rculos. Por outro lado, se os autovalores de A têm valores muito diferentes uns dos outros, com alguns
muito pequenos e alguns muito grandes, as elipses serão bastante excêntricas e, dependendo do chute inicial,
a convergência pode ser muito lenta (matrizes com estas propriedades são chamadas mal-condicionadas; para
que o método de descida acentuada seja lento, a matriz A não precisa ser muito mal-condicionada).
Como vimos na seção anterior, os algoritmos de Gauss-Seidel e SOR podem ser encarados como algoritmos
de descida. A discussão no parágrafo anterior também pode ser usada para entender a relativa lentidão destes
algoritmos.
xj = x0 + α0 p0 + α1 p1 + . . . + αj−1 pj−1 ,
de modo que xj está no subespaço afim gerado pelo chute inicial x0 e pelos vetores p0 , p1 , . . . , pj−1 .
Enquanto o método da descida mais acentuada minimiza o funcional de energia f apenas ao longo das j
Rodney Josué Biezuner 122
donde
2 2 2 2 2
2 kykA = kx − ykA + kxkA + 2 hy, xiA + kykA − 2 kxkA
2 2 2
= kx − ykA + 2 hy, xiA − kxkA + kykA ,
ou
2 2 2
kykA − 2 hy, xiA = kx − ykA − kxkA .
Logo, podemos escrever
1 2 1 2
kekA − kxkA .
f (y) = (5.71)
2 2
Conseqüentemente, minimizar o funcional f é equivalente a minimizar a A-norma do erro.
Agora, em um método de descida, depois de j iterações temos:
ej = x − xj = x − x0 − α0 p0 + α1 p1 + . . . + αj−1 pj−1
= e0 − α0 p0 + α1 p1 + . . . + αj−1 pj−1 .
2
Logo, minimizar ej A
é equivalente a minimizar
e0 − α0 p0 + α1 p1 + . . . + αj−1 pj−1
A
,
o que por sua vez é equivalente a encontrar a melhor aproximação do vetor e0 no subespaço Wj = p0 , p1 , . . . , pj−1 .
Esta é dada pelo lema da melhor aproximação:
5.29 Proposição. Sejam A ∈ Mn (R) uma matriz simétrica positiva definida, v ∈ Rn e W um subsespaço
de Rn . Então existe um único w ∈ W tal que
hy, ziA = 0
Nosso objetivo então é desenvolver um método em que o erro a cada passo é conjugado com todas as direções
de busca anteriores. O próximo resultado, que é basicamente uma reafirmação da Proposição 5.26, mostra
que em qualquer método de descida em que a busca na reta é exata satisfaz automaticamente ej ⊥A pj−1 ,
isto é, (5.72) é válido para a última iteração (o erro da iteração presente é A-ortogonal à direção de busca
da iteração anterior).
5.31 Proposição. Seja xk+1 = xk + αk pk obtido através de uma busca na reta exata. Então
rk+1 ⊥ pk
e
ek+1 ⊥A pk .
Prova: Temos
b − Axk+1 = b − Axk − αk Apk ,
de modo que a seqüência dos resı́duos é dada pela fórmula
Logo,
pk , rk
rk+1 , pk = rk+1 , pk − αk Apk , pk = rk , pk − Apk , pk = 0.
hpk , Apk i
Além disso, como
Aek+1 = rk+1 ,
segue que
ek+1 , pk A
= Aek+1 , pk = rk+1 , pk = 0.
O significado geométrico deste resultado é que o mı́nimo do funcional f na reta xk + αk pk ocorre quando a
derivada direcional de f na direção de busca é zero, ou seja,
∂f
xk+1 = ∇f xk+1 , pk = rk+1 , pk .
0=
∂pk
De acordo com a Proposição 5.31, depois do primeiro passo temos e1 ⊥A p0 . Para manter os erros
subseqüentes conjugados a p0 , como
ek+1 = x − xk+1 = x − xk − αk pk
ou
ek+1 = ek − αk pk , (5.74)
0 1 0
basta escolher as direções de busca subseqüentes conjugadas a p . Se escolhemos p conjugado a p , obtemos
x2 para o qual o erro satisfaz e2 ⊥A p1 ; como p1 ⊥A p0 , segue de (5.74) que e2 ⊥A p0 também. Para manter
os erros subseqüentes conjugados a p0 e p1 , basta escolher as direções de busca subseqüentes conjugadas a
p0 e p1 . Assim, vemos que para obter a condição (5.72) basta escolher as direções de busca de tal forma que
pi ⊥A pj para todos i 6= j.
Um método com estas caracterı́sticas é chamado um método de direções conjugadas. Estes resultados
são resumidos na proposição a seguir:
Rodney Josué Biezuner 124
5.32 Teorema. Se um método emprega direções de busca conjugadas e performa buscas na reta exatas,
então
ej ⊥A pi para i = 1, . . . , j − 1,
para todo j. Conseqüentemente
ej A
= min e0 − p A
,
p∈Wj
onde Wj = p0 , p1 , . . . , pj−1 .
Prova: A demonstração é por indução. Para j = 1, temos e1 ⊥A p0 pela Proposição 5.31 porque a busca
na reta é exata. Em seguida, assuma ej ⊥A pi para i = 1, . . . , j − 1; queremos mostrar que ej+1 ⊥A pi
para i = 1, . . . , j. Como
ej+1 = ej − αj pj ,
para i = 1, . . . , j − 1 temos
ej+1 , pi A
= ej − αj pj , pi A
= ej , pi A
− αj pj , pi A
=0−0=0
porque as direções de busca são conjugadas. ej+1 ⊥A pj segue novamente da Proposição 5.31.
Quando a direção inicial é dada pelo vetor gradiente de f , como na primeira iteração do método da descida
mais acentuada, obtemos o método do gradiente conjugado. As direções subseqüentes são escolhidas
através de A-ortogonalizar o resı́duo (ou vetor gradiente de f , que é a direção de busca em cada iteração
do método da descida mais acentuada) com todas as direções de busca anteriores, para isso utilizando o
algoritmo de Gram-Schmidt. Assim, dado um chute inicial p0 , a primeira direção é
p0 = −∇f x0 = b − Ax0 = r0
p0 = r 0 . (5.75)
rk+1 , pi A
cki = . (5.77)
hpi , pi iA
de forma que pk+1 ⊥A pi para todos i = 1, . . . , k. Felizmente, como veremos a seguir depois de algum trabalho
preliminar (Corolário 5.37), cki = 0 para todo i exceto i = k, o que torna necessário que apenas a direção
de busca mais recente pk seja armazenada na memória do computador, o que garante que a implementação
do gradiente conjugado é eficiente:
ou, definindo
rk+1 , Apk
βk = − , (5.79)
hpk , Apk i
temos que
pk+1 = rk+1 + βk pk . (5.80)
Rodney Josué Biezuner 125
Esta é a modificação do método do gradiente conjugado em relação ao método da descida mais acentuada,
no qual tomamos pk+1 = rk+1 .
Podemos obter uma expressão mais simples para o escalar βk , em função apenas dos resı́duos. Com
efeito, temos
rk+1 , rk+1 = rk+1 , rk − αk rk+1 , Apk = −αk rk+1 , Apk
porque os resı́duos obtidos através do método do gradiente conjugado são mutualmente ortogonais (veja
Corolário 5.36), logo
rk+1 , Apk rk+1 , rk+1
β=− k k
= .
hp , Ap i αk hpk , Apk i
Temos
pk , r k rk + βpk−1 , rk rk , rk
αk = = = ,
hpk , Apk i hpk , Apk i hpk , Apk i
porque pk−1 , rk = 0 pela Proposição 5.31, logo
rk , rk
αk = . (5.81)
hpk , Apk i
Portanto
rk+1 , rk+1
β= . (5.82)
hrk , rk i
Podemos obter um algoritmo ainda mais eficiente para o método do gradiente conjugado se observarmos que
para calcular o resı́duo rk+1 = b − Axk+1 em cada iteração não é necessário calcular Axk+1 explicitamente;
de fato, como vimos na demonstração da Proposição 5.31, temos rk+1 = rk − αk Apk . Assim, um algoritmo
eficiente para o método do gradiente conjugado poderia ser escrito da seguinte forma:
initialize x;
set b;
r ← b − Ax;
rScalarR ← hr, ri ;
set M ; //maximumNumberOfIterations
numberOf Iterations = 0;
do
until numberOf Iterations > M
Ap ← Ap;
pScalarAp ← hp, Api ;
α ← rScalarR/pScalarAp;
x ← x + αp;
r ← r − αAp;
rN ewScalarRN ew ← hr, ri ;
β ← rN ewScalarRN ew/rScalarR;
p ← r + βp;
rScalarR ← rN ewScalarRN ew;
numberOf Iterations + +;
5.34 Teorema. Depois de j iterações do algoritmo do gradiente conjugado (com rk 6= 0 em cada iteração),
temos
p0 , p1 , . . . , pj−1 = r0 , r1 , . . . , rj−1 = Kj A, r0 .
Prova: A demonstração é por indução. O resultado é trivial para j = 0, pois p0 = r0 . Assuma o resultado
válido para j − 1. Em primeiro lugar, mostraremos que
r0 , r1 , . . . , rj ⊂ Kj+1 A, r0 .
(5.83)
p0 , p1 , . . . , pj ⊂ r0 , r1 , . . . , rj . (5.84)
Por hipótese de indução, basta provar que pj ∈ r0 , r1 , . . . , rj . Isso segue de (5.76) e da hipótese de indução.
Até aqui provamos que
p0 , p1 , . . . , pj ⊂ r0 , r1 , . . . , rj ⊂ Kj+1 A, r0 .
(5.85)
Para provar que eles são iguais, basta mostrar que eles têm a mesma dimensão. Isso decorre de
dim r0 , r1 , . . . , rj 6 j + 1,
dim Kj+1 A, r0 6 j + 1
e
dim p0 , p1 , . . . , pj = j + 1,
o último porque os vetores p0 , p1 , . . . , pj são vetores não-nulos A-ortogonais.
ej ⊥A Kj A, r0
para todo j.
rj ⊥ Kj A, r0
para todo j.
Prova: Em vista do Teorema 5.34, basta provar que rj ⊥ p0 , p1 , . . . , pj−1 para todo j. Como Aej+1 = rj+1 ,
logo
Api ∈ Ar0 , A2 r0 , . . . , Ai+1 r ⊂ Ki+2 A, r0 ⊂ Kk+1 A, r0
5.38 Teorema. Seja A uma matriz simétrica positiva definida n×n. Então o método do gradiente conjugado
converge em n iterações.
Prova: Se fizemos n − 1 iterações em obter x, pelo Corolário 5.37 os vetores r0 , r1 , . . . , rn−1 formam uma
base ortogonal para Rn . Depois de mais uma iteração, de acordo com este mesmo corolário o resı́duo rn
satisfaz rn ⊥ r0 , r1 , . . . , rn−1 = Rn , logo rn = 0.
De fato, na maioria das aplicações o método do gradiente conjugado converge ainda mais rápido, se apenas
uma boa aproximação é requerida. Defina o número de condição de uma matriz simétrica positiva definida
por
max {λ : λ é um autovalor de A}
κ (A) = ; (5.86)
min {λ : λ é um autovalor de A}
assim, quanto maior o número de condição de uma matriz, ela é mais mal-condicionada e a convergência
de métodos de descida é mais vagarosa. Pode-se provar a seguinte estimativa de erro para o método do
gradiente conjugado (veja [Strikwerda]):
p !k
k 0 κ (A) − 1
e A
62 e A
p . (5.87)
κ (A) + 1
Esta estimativa é uma estimativa grosseira, mas mostra que o método do gradiente conjugado converge
mais rapidamente para matrizes bem-condicionadas (κ (A) ∼ 1). Uma comparação entre a velocidade de
convergência dos dois métodos para a matriz de discretização da fórmula de cinco pontos aplicada ao problema
descrito na primeira seção deste capı́tulo, desta vez com o tamanho das matrizes indicado na linha superior
da tabela, é dada a seguir [Watkins].
n = 81 n = 361 n = 1521
Descida Mais Acentuada 304 1114 4010
Gradiente Conjugado 29 60 118
(n − 1) π
sen2 π π∆x 4
κ (A) = 2n = cot2 = cot2 ≈ 2 2
2
π 2n 2 π ∆x
sen
2n
de modo que p
κ (A) − 1 1 − π∆x/2
p ≈ ≈ 1 − π∆x,
κ (A) + 1 1 + π∆x/2
o que dá uma velocidade de convergência para o método do gradiente conjugado duas vezes maior que a
do método SOR com o fator de relaxamento ótimo. No entanto, deve-se ter em mente que enquanto que a
taxa de covergência que obtivemos para o método SOR é precisa, a estimativa de erro (5.87) para o método
do gradiente conjugado é apenas um limitante superior grosseiro (veja [Watkins] para algumas estimativas
melhoradas).
Capı́tulo 6
Métodos Multigrid
Neste capı́tulo consideraremos o método multigrid, que é o método mais rápido para resolver equações
elı́pticas em geral. Embora o método possa ser empregado em malhas de elementos finitos e volumes fini-
tos também, neste capı́tulo consideraremos o seu emprego apenas em malhas de diferenças finitas para a
equação de Poisson no quadrado. A tabela a seguir (adaptada de [TOS]) compara o custo de processamento
em uma máquina serial de alguns dos métodos mais populares para resolver sistemas lineares que surgem na
discretização do
problema de Poisson (à exceção do método de eliminação gaussiana cujo custo de armazena-
mento é O n2 , todos os demais métodos tem custo de armazenamento O (n)). Como estamos comparando
métodos diretos (eliminação gaussiana e transformada de Fourier rápida (FFT) ) com métodos iterativos
(todos os demais), assumimos um único critério de parada para os vários métodos iterativos; se o critério de
parada for escolhido da ordem do erro de discretização da malha, um fator O (log n) deve ser multiplicado
para todos os métodos iterativos, à exceção do multigrid completo.
A idéia do método multigrid é baseada em dois princı́pios: suavização do erro e a sua correção em
um grid mais grosseiro (menos refinado). Estes princı́pios serão explicados em detalhes nas próximas
seções.
Em linhas gerais, a idéia básica é eliminar os componentes de alta freqüência do erro em uma malha
refinada. Para que isso ocorra, é necessário que estes componentes de alta freqüência correspondam aos
menores autovalores da matriz de iteração porque, como vimos no capı́tulo anterior, estes são eliminados
rapidamente pelos métodos iterativos lineares (a velocidade de convergência de cada método é dada pelo raio
espectral da matriz de iteração, que corresponde ao valor absoluto do maior autovalor |λ1 | < 1, enquanto
que as componentes do erro correspondentes aos menores autovalores λj convergem para zero muito mais
rapidamente (|λj /λ1 | 1); isso significa que este método iterativo suaviza o erro, pois quanto maior a
influência das componentes de maior freqüência (maior oscilação), menos suave é a função. Aqui é útil fazer
uma analogia com a série de Fourier: é exatamente a presença de componentes de oscilação arbitrariamente
maior que permite que a série convirja para uma função não diferenciável, ou mesmo descontı́nua; se a
série for truncada a qualquer momento o resultado é sempre uma função suave, pois é a combinação linear
finita de autofunções suaves. Esta visualização também permanece verdade para funções discretizadas em
128
Rodney Josué Biezuner 129
malhas de diferenças finitas escritas como uma combinação linear das autofunções da matriz de iteração nesta
malha: mesmo que o número de componentes da função seja finito, porque a malha é discreta a presença de
componentes de alta oscilação dão origem a um gráfico com um aspecto escarpado, não suave.
Assim, como o nosso objetivo é eliminar apenas as componentes de alta freqüência do erro, e não todo o
erro, poucas iterações do método iterativo são necessárias nesta malha refinada, onde o custo computacional é
alto (malhas muito refinadas significa que elas possuem muitos pontos, o que por sua vez implica em matrizes
de discretização muito grandes). Ocorre que algumas autofunções de freqüência baixa em uma malha mais
refinada correspondem a autofunções de freqüência alta em uma malha mais grosseira (como veremos). Uma
vez tendo eliminado as componentes de alta freqüência do erro na malha mais refinada, tendo deixado as
componentes de baixa freqüência praticamente intocadas, transferimos o problema para uma malha mais
grosseira, cujos componentes de alta freqüência do erro correspondem a alguns dos componentes de baixa
freqüência do erro na malha mais refinada anterior, que não puderam ser eliminados com as poucas iterações
do método iterativo permitidas na malha mais refinada. Com poucas iterações do método iterativo nesta
malha mais grosseira, estes erros também são rapidamente eliminados, a um custo computacional mais baixo
do que se tivéssemos tentado eliminá-los ainda na malha mais refinada. Este processo é a correção do erro
em uma malha mais grosseira. Ele é repetido em malhas cada vez mais grosseiras até que todo o erro é
eliminado, a um custo computacional muito mais baixo do que se tivéssemos trabalhado sempre na malha
mais refinada original.
onde uh como usual denota a solução do problema discretizado (aproximação da solução exata), fh a discre-
tização da função f em Ωh ,
1
h= , (6.2)
n
Ω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 n}
e
−1
1
− ∆h uh = 2 −1 4 −1 (6.3)
h
−1
ou, em outras palavras,
em m
h (xi , yj ) = uh (xi , yj ) − uh (xi , yj ) . (6.5)
Em geral, tomaremos n par, ou mesmo n = 2p para algum p. Assim, uma malha Ωh é mais refinada que
uma malha Ω2h (esta é mais grosseira que a primeira). Temos uma seqüência de malhas progressivamente
mais grosseiras:
Ωh ⊂ Ω2h ⊂ Ω4h ⊂ . . . ⊂ Ω2p h = Ω1 ,
onde Ω1 possui apenas uma célula.
ϕkl
h (x, y) = sen kπx sen lπy, 1 6 k, l 6 n − 1 (6.6)
onde x, y denotam as variáveis discretizadas (isto é, x = ih e y = jh para 0 6 i, j 6 n). Assim, o erro na
m-ésima iteração pode ser escrito na forma
n−1
X n−1
X
em
h (x, y) =
m kl
αk,l ϕh (x, y) = m
αk,l sen kπx sen lπy. (6.7)
k,l=1 k,l=1
m
O erro ser suavizado significa que após algumas poucas iterações temos αk,l muito pequeno para k, l grandes,
isto é, para
ϕkl
h (x, y) = sen kπx sen lπy de alta freqüência,
m
enquanto que o valor de αk,l para k, l pequenos, isto é, para
ϕkl
h (x, y) = sen kπx sen lπy de baixa freqüência,
pode ter mudado muito pouco. Como o fato de k, l serem grandes ou pequenos é definido relativamente de
acordo com o valor de n (valores de k, l próximos de n são considerados grandes, enquanto que valores de k, l
distantes de n são considerados pequenos), segue que autofunções de baixa freqüência em uma malha mais
refinada (n maior) podem ser autofunções de alta freqüência em uma malha mais grosseira (n relativamente
pequeno). Para propósitos de análise, vamos dar uma definição precisa a este conceito:
n
baixa freqüência se max (k, l) < ,
2
n
alta freqüência se 6 max (k, l) < n.
2
Além disso, se considerarmos especialmente a passagem da malha mais refinada Ωh para a malha mais
grosseira Ω2h com o dobro do espaçamento de malha, apenas as autofunções de freqüências mais baixas em
Rodney Josué Biezuner 131
Ωh são visı́veis em Ω2h , pois todas as autofunções de freqüência alta em Ωh coincidem com as autofunções
de freqüência baixa em Ω2h ou desaparecem em Ω2h . De fato, como
ϕk,l n−k,l
h (x, y) = −ϕh (x, y) = −ϕk,n−l
h (x, y) = ϕn−k,n−l
h (x, y) para (x, y) ∈ Ω2h , (6.8)
estas quatro autofunções não podem ser distingüidas umas das outras em Ω2h . Além disso, se k = n/2 ou
l = n/2, temos
ϕk,l
h (x, y) = 0 para (x, y) ∈ Ω2h . (6.9)
Para provar estas afirmações, escrevemos, por exemplo,
2i 2j 2i 2j
ϕhn−k,l (i (2h) , j (2h)) = sen (n − k) π sen lπ = sen −kπ + 2iπ sen lπ
n n n n
2i 2j
= − sen kπ sen lπ = −ϕk,l
h (i (2h) , j (2h))
n n
e
n/2,l n/2,l 2i 2j n 2i 2j 2jlπy
ϕh (i (2h) , j (2h)) = ϕh , = sen π sen lπ = sen iπ sen = 0.
n n 2 n n n
Assim, podemos decompor o erro em duas somas, uma representando os componentes de baixa freqüência
e a outras os componentes de alta freqüência:
n/2−1 n−1
X X
em
h (x, y) = m kl
αk,l ϕh (x, y) + m kl
αk,l ϕh (x, y) (6.10)
k,l=1 max(k,l)> n
2
Xbaixa Xalta
m kl m kl
= αk,l ϕh (x, y) + αk,l ϕh (x, y) . (6.11)
propriedades de suavização do erro, como veremos a seguir. A fórmula de iteração para o método de Jacobi
para o problema de Poisson discretizado é dada por
ukh (xi−1 , yj ) + ukh (xi+1 , yj ) + ukh (xi , yj−1 ) + ukh (xi , yj+1 ) + h2 fh (xi , yj )
uk+1
h (xi , yj ) = . (6.12)
4
Em notação de operadores, esta fórmula pode ser escrita como
h2
uk+1
h = Rh ukh + fh , (6.13)
4
onde o operador de iteração Rh é dado por
h2
Rh = Ih − Lh , (6.14)
4
Ih sendo o operador identidade e Lh = −∆h . No método de Jacobi amortecido, introduzimos o parâmetro
de relaxamento 0 < ω 6 1:
uk+1
h (xi , yj ) = ukh (xi , yj )
k
uh (xi−1 , yj ) + ukh (xi+1 , yj ) + ukh (xi , yj−1 ) + ukh (xi , yj+1 ) + h2 fh (xi , yj )
+ω − ukh (xi , yj ) .
4
Logo
h2
uk+1
h = Ih ukh+ω Sh ukh
+ fh − Ih uh k
4
2
h2
k k h k k
= Ih uh + ω Ih uh − Lh uh + fh − Ih uh
4 4
2 2
ωh ωh
= Ih ukh − Lh ukh + fh ,
4 4
ou
ωh2
uk+1
h = Rh (ω) ukh + fh , (6.15)
4
onde
ωh2
Rh (ω) = Ih −
Lh . (6.16)
4
Em notação estêncil, o operador iteração para o método de Jacobi amortecido pode ser escrito na forma
2 −1
ωh 1 −1
Rh (ω) = 1 − 4 −1
4 h2
−1
1
4
1 1
= 4 1−ω 4
1
4
ou também
1
ω 1
Rh (ω) =
1 4 −1 1
.
4
ω
1
Rodney Josué Biezuner 133
ωh2
Rh (ω) = Ih + ∆h ,
4
logo os autovalores de Rh e −∆h estão relacionados da seguinte forma: λ é um autovalor de −∆h se e
somente se
ωh2
(Rh − Ih ) v = − λv,
4
isto é, se e somente se
ωh2
λh (ω) = 1 − λ (6.17)
4
é um autovalor de Rh e as autofunções de Rh são as mesmas autofunções de −∆h . As autofunções de −∆h
são, como já vimos,
ϕkl
h (x, y) = sen kπx sen lπy, 1 6 k, l 6 n − 1,
para 0 < ω 6 1, de modo que ω = 1 (método de Jacobi) oferece a melhor velocidade de convergência,
enquanto que ρ (Rh ) > 1 para ω > 1 se h é suficientemente pequeno e o método não converge.]
Para analisar as propriedades suavizadoras do método de Jacobi amortecido quantitativamente, introdu-
zimos o fator suavizante de Rh :
Observe que µh (ω) é o maior autovalor dentre as maiores freqüências e representa o pior fator pelo qual os
componentes de alta freqüência do erro são reduzidos por passo de iteração. Para entender isso, fixe um
tamanho de malha h e escreva os autovalores de Rh (como no capı́tulo anterior) na forma
temos
q
X
ekh = Rhk e0h = αi λki ϕi .
i=1
Como
k
|λi | → 0,
k k
se |λi | < 1, a taxa de eliminação para o componente ϕi do erro é determinada por |λi | e em cada iteração
este componente é reduzido por um fator exatamente igual a |λi |. Como
n ω n o
µh (ω) = max 1 − (2 − cos kπh − cos lπh) : 6 max (k, l) 6 n − 1 ,
2 2
∗
n ω o
µ (ω) = max 1 − , |1 − 2ω| ,
2
segue que para 0 < ω < 1 o fator suavizante é menor que 1 e permanece longe de 1 por um limitante
independente de h. Para ω = 1, o fator suavizante é da ordem de 1 − O h2 apenas; os menores autovalores
do método de Jacobi
1
λkl = (cos kπh + cos lπh)
2
estão associados às autofunções de freqüências médias, logo as autofunções de freqüências altas não são
rapidamente eliminadas e não há suavização. Por exemplo,
cos πh se ω = 1,
1
se ω = 1,
2 + cos πh 3
se ω = 0.5, se ω = 0.5,
µh (ω) = 4 µ∗ (ω) = 4
1 + 2 cos πh 3
se ω = 0.8, se ω = 0.8,
5 5
A escolha de ω = 0.8 é ótima no sentido de que
enquanto que
4 3 cos πh 3
= − O h2 .
inf µh (ω) = µh = (6.20)
0<ω61 4 + cos πh 4 + cos πh 5
Isso significa que um passo do método de Jacobi amortecido com ω = 0.8 reduz todos os componentes do
erro de alta freqüência por um fator de pelo menos 3/5, independente do tamanho h da malha.
em m
h = uh − uh ,
Lh em m
h = rh . (6.22)
Rodney Josué Biezuner 135
Para transferir funções definidas em uma malha mais refinada Ωh para funções definidas em uma malha mais
grosseira Ω2h e vice-versa, precisamos definir dois operadores lineares de transferência: um operador de
restrição
Ih2h : G (Ωh ) −→ G (Ω2h ) (6.23)
e um operador de interpolação (ou de prolongamento)
h
I2h : G (Ω2h ) −→ G (Ωh ) . (6.24)
O operador de restrição será usado para restringir o resı́duo rhm obtido na malha mais refinada Ωh para a
malha mais grosseira Ω2h onde ele será corrigido:
m
r2h = Ih2h rhm . (6.25)
Outro operador freqüentemente usado é o operador peso total, que em notação estêncil é dado por
1 2 1
1
2 4 2 ,
16
1 2 1
ou seja,
1
Ih2h vh (x, y) =
[4vh (x, y) + 2vh (x, y − h) + 2vh (x − h, y) + 2vh (x + h, y) + 2vh (x, y + h)
16
+vh (x − h, y − h) + vh (x + h, y − h) + vh (x − h, y + h) + vh (x + h, y + h)] .
6.4 Exemplo. Um dos operadores de interpolação mais simples de implementar é o operador de interpolação
bilinear :
v2h (x, y)
se (x, y) = (2kh, 2lh) ,
1
(v2h (x, y − h) + v2h (x, y + h)) se (x, y) = (2kh, (2l − 1) h) ,
2
h 1
I2h v2h (x, y) = (v2h (x − h, y) + v2h (x, y + h)) se (x, y) = ((2k − 1) h, 2lh) ,
2
1
[vh (x − h, y − h) + vh (x + h, y − h)
se (x, y) = ((2k − 1) h, (2l − 1) h) .
4
+ vh (x − h, y + h) + vh (x + h, y + h)]
Rodney Josué Biezuner 136
Ciclo de 2 Malhas
1. Pré-suavização
a) Calcule um m
h através de n1 passos de um suavizador aplicado a uh :
um
h = SUAVIZE
n1
(um
h , Lh , fh ).
a) Calcule o resı́duo rm m
h = fh − Lh uh .
m m
em
e) Calcule a aproximação corrigida: uh = uh + eh .
3. Pós-suavização
a) Calcule um+1
h em
através de n2 passos de um suavizador aplicado a uh :
um+1
h = SUAVIZEn2 (e
um
h , Lh , fh ).
A necessidade da pós-suavização deve-se ao fato que as freqüências mais baixas na malha mais grosseira
correspondem não somente às freqüências mais baixas na malha mais refinada, como também às freqüências
mais altas, como vimos em (4.1) (em outras palavras, freqüências baixas em Ω2h são mapeadas para a mesma
freqüência baixa em Ωh e para três outras freqüências altas em Ωh ); para evitar que estas componentes
do erro reapareçam, fazemos uma segunda suavização. Observe que vários componentes individuais do
ciclo de duas malhas devem ser especificados, e sua escolha pode ter uma forte influência na eficiência do
algoritmo: o procedimento suavizador SUAVIZE (um h , Lh , fh ); os números n1 e n2 de passos de suavização,
a malha grosseira (aqui escolhemos Ω2h , mas outras escolhas são possı́veis) e os operadores de restrição e de
interpolação.
do resı́duo L2h em m
2h = r 2h exatamente, mas suavizá-la e transferir o problema para uma malha ainda mais
grosseira Ω4h , onde o custo computacional é ainda menor. Esta idéia é repetida até que podemos em princı́pio
chegar na malha Ω1 , onde a correção do resı́duo pode então ser calculada exatamente. Daı́, voltamos para
a malha mais refinada original, formando um ciclo no formato da letra “V”.
Capı́tulo 7
A discretização do domı́nio no métodos dos volumes finitos difere da do método de diferenças finitas em que
nesta o domı́nio é substituı́do por um conjunto de pontos, enquanto que na primeira o domı́nio é subdividido
em volumes de controle ou células. Os pontos nodais ou simplesmente nós, são os centros das células.
No método dos volumes finitos, ao invés de aproximarmos diretamente a equação diferencial como no método
de diferenças finitas, ela é antes integrada sobre cada volume de controle. As integrais obtidas são então
aproximadas. As equações integrais estão na forma de leis de conservação, o que assegura a conservação
das grandezas fı́sicas tratadas em cada volume de controle (conservação no nı́vel discreto) e portanto este
método é bastante adequado para tratar de fenômenos fı́sicos que envolvem leis de conservação. Muitas
vezes pode-se trabalhar diretamente com as equações integrais, sem passar pelas equações diferenciais, o que
torna o método particularmente útil para tratar de fenômenos descontı́nuos melhor modelados por equações
integrais, tais como fenômenos que envolvem ondas de choque.
138
Rodney Josué Biezuner 139
ou abstrata. Por exemplo, no caso da equação do calor, a temperatura u é uma medida da densidade de
energia térmica. De fato, se e(x, t) denota a densidade de energia térmica, isto é, a quantidade de energia
térmica por unidade de volume, então a densidade de energia térmica e a temperatura estão relacionadas
através da equação
e(x, t) = c(x)ρ(x)u(x, t),
cujo significado é: a energia térmica por unidade de volume é igual à energia térmica por unidade de massa
por unidade de temperatura (i.e., o calor especı́fico), vezes a temperatura, vezes a densidade volumétrica de
massa.
Imaginamos que a substância está distribuı́da em um tubo uniforme com seção transversal de área
constante A. Por hipótese, u é constante em cada seção transversal do tubo, variando apenas na direção x.
Considere um segmento arbitrário do tubo, entre as seções transversais localizadas em x = a e em x = b.
Chamamos este segmento de volume de controle. A quantidade total da substância dentro do volume de
controle no instante de tempo t é
Z b
Quantidade total da substância
= u(x, t)A dx.
dentro do volume de controle a
Assuma agora que existe movimento da substância através do tubo na direção axial. Definimos o fluxo
φ(x, t) da substância no tempo t como sendo a quantidade da substância fluindo através da seção transversal
em x no tempo t por unidade de área, por unidade de tempo. Assim as dimensões de φ são [φ] = quantidade
da substância / (área × tempo). Por convenção, φ será positivo se a substância estiver se movendo na direção
positiva do eixo x, e negativo se ela estiver se movendo na direção negativa do eixo x. Portanto, no tempo t,
a quantidade lı́quida de substância permanecendo no volume de controle será a diferença entre a quantidade
da substância entrando em x = a e a quantidade da substância saindo em x = b:
A substância pode ser criada ou destruı́da dentro do volume de controle por uma fonte interna ou externa.
A taxa de criação ou destruição da substância, que chamaremos de termo fonte e denotaremos por f (x, t, u),
tem dimensões [f ] = quantidade da substância / (volume × tempo), tendo sinal positivo se a substância é
criada dentro do volume de controle e negativa se a substância for destruı́da dentro do volume de controle.
Observe que ela pode depender da própria quantidade da substância disponı́vel, medida pela densidade u.
A taxa de criação ou destruição da substância dentro do volume de controle é então dada por
Z b
Taxa de criação da substância
= f (x, t, u)A dx.
dentro do volume de controle a
Esta é a lei de conservação na forma integral, valendo mesmo se u, φ ou f não forem funções diferenciáveis
(o que pode ocorrer em certos fenômenos fı́sicos, como por exemplo naqueles que envolvem ondas de choque
Rodney Josué Biezuner 140
ou outros tipos de descontinuidade). Se estas funções forem continuamente diferenciáveis, podemos derivar
sob o sinal de integração na primeira integral
d b
Z Z b
u(x, t) dx = ut (x, t) dx,
dt a a
Em n dimensões, o fluxo pode ser em qualquer direção, logo ele é uma grandeza vetorial que denotaremos
por φ(x, t). Se η(x) denota o vetor unitário normal apontando para fora da região V , a taxa de transferência
lı́quida da substância para fora do volume de controle através de sua fronteira ∂V é dada por
Z
Taxa de transferência lı́quida da substância
= φ(x, t) · η(x) dS.
para fora do volume de controle ∂V
Se u, φ e f forem todas de classe C 1 (assim como a região V ), podemos derivar sob o sinal de integração e
usar o Teorema da Divergência
Z Z
φ(x, t) · η(x) dS = div φ(x, t) dV,
∂V V
onde a constante de condutividade térmica k = k (x) depende do material e muitas vezes pode ser considerada
constante.
Em dimensões mais altas, a lei de Fourier assume a forma
De fato, para materiais isotrópicos (isto é, materiais em que não existem direções preferenciais) verifica-se
experimentalmente que o calor flui de pontos quentes para pontos frios na direção em que a diferença de
temperatura é a maior. O fluxo de calor é proporcional à taxa de variação da temperatura nesta direção, com a
constante de proporcionalidade k sendo por definição a condutividade térmica, como no caso unidimensional.
Como sabemos, a direção onde uma função cresce mais rápido é exatamente aquela dada pelo vetor gradiente
da função, e o módulo do gradiente fornece a magnitude da taxa de variação da função nesta direção. O sinal
negativo ocorre, como no caso unidimensional, porque o vetor gradiente aponta na direção de crescimento
da temperatura, enquanto que o fluxo do calor se dá na direção oposta (da temperatura maior para a
temperatura menor). O fluxo do calor em uma região bi ou tridimensional pode ser facilmente visualizado
quando se lembra que o gradiente de uma função é perpendicular às superfı́cies de nı́vel da função. No
caso em que a função é a temperatura, as superfı́cies de nı́vel são chamadas superfı́cies isotérmicas ou,
simplesmente, isotermas. Assim, o calor flui das isotermas mais quentes para as isotermas mais frias, e em
cada ponto da isoterma perpendicularmente à isoterma. Em outras palavras, as linhas de corrente do fluxo
de calor correspondem às linhas de fluxo do campo gradiente da temperatura.
Substituindo a relação constitutiva na lei de conservação, obtemos a equação do calor: na forma diver-
gente,
ut = div (k∇u) + f (x, t, u),
ou, quando k é constante, na forma usual envolvendo o laplaciano,
7.2 Exemplo (Equação da Difusão). Em muitos outros processos fı́sicos observa-se que a substância flui a
uma taxa diretamente proporcional ao gradiente de densidade, de regiões de maior densidade para regiões
de menor densidade. Esta relação geral é chamada de lei de Fick :
onde D = D (x) é a constante de difusão. Assumindo D constante, se o termo fonte independe de u, obtemos
a equação da difusão
ut = D∆u + f (x, t),
caso contrário a equação diferencial parcial obtida é chamada equação da difusão-reação
que aparece na teoria de combustão e em biologia. Se D não é constante, obtemos as respectivas equações
na forma divergente. O nome difusão vem do fato de que a substância difunde-se para regiões adjacentes por
causa de gradientes (i.e., diferenças) de concentração, e não porque é transportada pela corrente (i.e., não
através de convecção). Por este motivo, o termo D∆u é chamado de termo difusivo.
Além do calor, exemplos de outras substâncias que se comportam assim são substâncias quı́micas dissol-
vidas em algum fluido (neste caso, u representa a concentração quı́mica) e até mesmo populações de insetos.
Além de ser confirmada através de observações empı́ricas, a lei de Fick que governa estes e vários outros
fenômenos fı́sicos e biológicos pode ser justificada teoricamente através de argumentos baseados em modelos
probabilı́sticos e caminhos aleatórios.
Neste texto sobre equações elı́pticas, obviamente estamos interessados na equação de estado estacionário
resultante da equação da difusão ou de difusão-reação, isto é, no caso em que ut = 0:
O primeiro passo é gerar a malha de volumes finitos no intervalo [0, L], isto é, discretizar o domı́nio
em volumes de controle. Para isso, inserimos um número n de pontos nodais ou nós P1 , . . . , Pn entre
os pontos 0 e L da fronteira do domı́nio. Os n volumes de controle V1 , . . . , Vn serão centrados nestes nós.
As faces (fronteiras) dos volumes de controle serão posicionadas no ponto médio entre dois nós. Em geral,
posiciona-se os volumes de controle de modo que as fronteiras do domı́nio coincidem com faces dos volumes
de controle, isto é, o ponto 0 está na face esquerda do primeiro volume de controle e o ponto L está na face
direita do último volume de controle. Para simplificar a apresentação, assumiremos que os pontos nodais
foram posicionados de modo a estarem igualmente espaçados, de modo que todos os volumes de controle
têm mesma largura igual a ∆x.
Estabelecemos a seguinte notação (esta convenção é freqüentemente utilizada em dinâmica dos fluidos
computacional, onde o método dos volumes finitos é bastante popular): um ponto nodal arbitrário será
designado simplesmente por P e os seus pontos nodais vizinhos serão designados por W (oeste, isto é, o
ponto nodal vizinho à esquerda) e E (leste, correspondendo ao vizinho à direita). A face esquerda (à oeste)
do volume de controle será designada por w e a face direita (à leste) por e. Assim, a distância entre dois nós
vizinhos, assim como a distância entre as duas faces de um volume de controle é igual a ∆x.
Uma vez discretizado o domı́nio com a geração da malha de volumes de controle, integrando a equação di-
ferencial parcial em cada volume de controle para colocá-la na forma integral (reobtendo a lei de conservação;
é claro que podemos desde o inı́cio trabalhar diretamente com esta, se estiver disponı́vel):
Z Z
d du
− a (x) dx = f (x, u) dx.
Vp dx dx Vp
Rodney Josué Biezuner 143
Observe que a equação integral obtida é uma equação exata, ainda não discretizada. Na linguagem de leis
de conservação, ela diz simplesmente que o fluxo de u deixando a face direita do volume de controle menos
o fluxo deixando a face esquerda do mesmo (respeitando a nossa convenção de sinal para fluxos) é igual à
quantidade de u gerada pela fonte dentro do volume de controle:
φw − φe = f VP ∆x.
Agora procedemos à discretização da equação integral. Valores nas faces devem ser dados em funções de
valores nos pontos nodais. Consideremos primeiro os volumes de controle interiores V2 , . . . , Vn−1 . Usando
interpolação linear, podemos obter valores aproximados para a (xe ) , a (xw ), calculados nas faces dos volumes
de controle, em termos dos valores de a nos pontos nodais dos volumes de controle:
aW + aP
aw := a (xw ) = , (7.10)
2
aP + aE
ae := a (xe ) = . (7.11)
2
As derivadas primeiras, ou seja, os fluxos, podem ser aproximadas através de diferenças finitas apropriadas,
por exemplo diferenças finitas centradas:
du du uP − uW
:= (xw ) = , (7.12)
dx w dx ∆x
du du uE − uP
:= (xe ) = . (7.13)
dx e dx ∆x
O termo fonte, que pode expressar uma dependência não linear do valor de u, pode ser linearizado e assumido
constante ao longo do volume de controle, produzindo
fP0 + fP1 up
Z Z
1 0 1
dx = fP0 + fP1 up .
f VP = f + fP up dx = (7.14)
∆x Vp P ∆x Vp
(Como queremos obter um sistema linear no final, não é possı́vel aproximar o termo fonte por uma apro-
ximação de ordem maior que 1. A linearização do termo linear será discutida em maiores detalhes na seção
4 deste capı́tulo) Daı́,
uP − uW uE − uP
= fP0 + fP1 up ∆x,
aw − ae
∆x ∆x
ou
ap uP + aW uW + aE uE = bp , (7.15)
onde
aw ae
ap = + − fP1 , (7.16)
∆x2 ∆x2
aw ae
aW = − , aE = − , (7.17)
∆x2 ∆x2
bp = fP0 . (7.18)
Rodney Josué Biezuner 144
O tratamento dos volumes de controle adjacentes à fronteira é ligeiramente diferente. Para o volume de
controle V1 adjacente à fronteira esquerda (oeste) do domı́nio, temos
aw = a0 , (7.19)
e
du uP − u0
= , (7.20)
dx w ∆x/2
porque a distância entre P e 0 é apenas ∆x/2; neste caso somos forçados a utilizar uma diferença finita
progressiva para aproximar a derivada primeira em w. Assim, a equação discretizada correspondente a este
volume de controle é
uP − u0 uE − uP
= fP0 + fP1 up ∆x,
2a0 − ae
∆x ∆x
ou
ap uP + aE uE = bp , (7.21)
onde
2a0 ae
ap = + − fP1 , (7.22)
∆x2 ∆x2
ae
aE = − , (7.23)
∆x2
2a0
bp = fP0 + u0 . (7.24)
∆x2
Para o volume de controle Vn adjacente à fronteira direita temos
ae = aL ,
du uL − uP
= ,
dx e ∆x/2
utilizando uma diferença finita regressiva para aproximar a derivada primeira em e, de modo que a equação
discretizada correspondente a este volume de controle é
uP − uW uL − uP
= fP0 + fP1 up ∆x,
aw − 2ae
∆x ∆x
ou
ap uP + aE uE = bp , (7.25)
onde
aw 2aL
ap = + − fP1 , (7.26)
∆x2 ∆x2
aw
aW = − , (7.27)
∆x2
2aL
bp = fP0 + uL . (7.28)
∆x2
Ordenando os volumes de controle (geralmente da esquerda para a direita), obtemos um sistema linear cuja
solução será uma solução aproximada para a equação com as condições de fronteira dadas.
7.3 Exemplo (Equação de Poisson). Vamos aplicar o método de volumes finitos à equação de Poisson com
condição de fronteira de Dirichlet
−u00 = f (x)
em [0, L] ,
(7.29)
u (0) = u0 , u (L) = uL .
Rodney Josué Biezuner 145
Aqui a (x) ≡ 1 e f (x, u) = f (x), de modo que fP1 = 0. Se decidimos aproximar o valor médio de f no
volume de controle pelo valor de f em P , segue que
2 1 1
ap = 2
, aW = − 2
, aE = − , bp = fP
∆x ∆x ∆x2
para os volumes de controle interiores V2 , . . . , Vn−1 . Para os volumes de controle adjacentes à fronteira, para
o primeiro volume de controle V1 temos
3 1 2
ap = , aE = − , bp = fP + u0 ,
∆x2 ∆x2 ∆x2
enquanto que para o último volume de controle Vn temos
3 1 2
ap = , aW = − , bp = fP + uL .
∆x2 ∆x2 ∆x2
O sistema discretizado é, portanto:
2
3 −1
u1
f +
1 ∆x2 0 u
−1 2 −1 u2 f2
. .
. .. . .. ..
..
1 −1 =
.
∆x2 . . . .
. . . . −1 . . .
.
−1 2 −1 un−1 fn−1
−1 3 un
2
fn + 2
uL
∆x
Compare com o correspondente sistema discretizado obtido pelo método de diferenças finitas; a única dife-
rença está na primeira e última linhas dos sistemas.
7.4 Exemplo (Equação Elı́ptica Linear). Consideremos agora o seguinte problema linear elı́ptico com
condição de fronteira de Dirichlet
−u00 = Au + B
em [0, L] ,
(7.30)
u (0) = u0 , u (L) = uL .
Novamente a (x) ≡ 1, mas desta vez f (x, u) = f (u) = Au + B, de modo que fP0 = B e fP1 = A. Segue que
2 1 1
ap = 2
− A, aW = − 2
, aE = − , bp = B
∆x ∆x ∆x2
para os volumes de controle interiores V2 , . . . , Vn−1 . Para os volumes de controle adjacentes à fronteira, para
o primeiro volume de controle V1 temos
3 1 2
ap = − A, aE = − , bp = B + u0 ,
∆x2 ∆x2 ∆x2
enquanto que para o último volume de controle Vn temos
3 1 2
ap = − A, aW = − , bp = B + uL .
∆x2 ∆x2 ∆x2
O sistema discretizado é, portanto:
2
3−A −1
u1
B + u0
∆x2
−1 2 − A −1 u2 B
. .
. .. . .. ..
..
1 −1 =
.
∆x2 . . . .
. . . . −1 .
. .
.
−1 2 − A −1 un−1 B
−1 3 − A un
2
B+ 2
uL
∆x
Rodney Josué Biezuner 146
Como é sabido, podemos assegurar que o problema linear elı́ptico possui solução única se A 6 0, utilizando
o princı́pio do máximo. Isso se traduz do ponto de vista numérico, no fato de que a matriz discretizada
permanece diagonalmente dominante. No caso em que A > 0 é preciso ter cuidado, pois pode haver infinitas
soluções exatas e não existir solução numérica e vice-versa, pois os autovalores do problema exato não
são iguais aos autovalores da matriz de discretização (na maioria dos casos estes últimos não são nem
boas aproximações para os primeiros: usualmente as aproximações são razoavelmente boas apenas para os
primeiros autovalores e em malhas bastante refinadas, com um número enorme de pontos ou células). Para
evitar este tipo de problema, é possı́vel modificar a linearização; veja a seção 4 deste capı́tulo.
obtemos agora
Z Z Z
∂ ∂u ∂ ∂u
− a (x, y) dxdy − a (x, y) dxdy = f (x, y, u) dxdy,
Vp ∂x ∂x Vp ∂y ∂y Vp
ou
Z n Z e Z e Z n Z
∂ ∂u ∂ ∂u
− a (x, y) dx dy − a (x, y) dy dx = f (x, y, u) dxdy.
s w ∂x ∂x w s ∂y ∂y Vp
= f (x, y, u) dxdy.
Vp
Obtemos, portanto, a seguinte equação parcialmente discretizada (diferente do caso unidimensional, esta
equação não é exata):
∂u ∂u ∂u ∂u
a (xw , yp ) (xw , yp ) ∆y − a (xe , yp ) (xe , yp ) ∆y + a (xp , ys ) (xp , ys ) ∆x − a (xp , yn ) (xp , yn ) ∆x
∂x ∂x ∂y ∂y
= f V ∆x∆y,
onde Z
1
f VP = f (x, u) dxdy.
∆x∆y Vp
Em termos de fluxos discretizados,
φw − φe + φs − φn = f VP ∆x∆y.
Usando interpolação 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 de
controle:
aW + aP
aw := a (xw , yp ) = , (7.32)
2
aP + aE
ae := a (xe , yp ) = , (7.33)
2
aS + aP
as := a (xp , xs ) = , (7.34)
2
aP + aN
an := a (xp , xn ) = . (7.35)
2
Os fluxos são aproximadas através de diferenças finitas centradas:
∂u ∂u uP − uW
:= (xw , yp ) = , (7.36)
∂x w ∂x ∆x
∂u ∂u uE − uP
:= (xe , yp ) = , (7.37)
∂x e ∂x ∆x
∂u ∂u uP − uS
:= (xp , ys ) = , (7.38)
∂y s ∂y ∆y
∂u ∂u uN − uP
:= (xp , yn ) = . (7.39)
∂y n ∂y ∆y
O termo fonte é linearizado
f 0 + fP1 up
Z Z
1
fP0 + fP1 up dxdy = P dxdy = fP0 + fP1 up .
f VP = (7.40)
∆x∆y Vp ∆x∆y Vp
Daı́,
uP − uW uE − uP uP − uS uN − uP
∆x = fP0 + fP1 up ∆x∆y,
aw ∆y − ae ∆y + as ∆x − an
∆x ∆x ∆y ∆y
ou
ap uP + aW uW + aE uE + aS uS + aN uN = bp . (7.41)
com
aw ae as an
ap = 2
+ 2
+ 2
+ − fP1 , (7.42)
∆x ∆x ∆y ∆y 2
aw ae as an
aW = − , aE = − , aS = − 2 , aN = − 2 , (7.43)
∆x2 ∆x2 ∆y ∆y
bp = fP0 . (7.44)
Rodney Josué Biezuner 148
O tratamento dos volumes de controle adjacentes à fronteira é diferente. Por exemplo, para volumes de
controle adjacentes à fronteira esquerda (oeste), que não sejam os dois volumes de controle dos cantos, temos
aw = a (0, yp ) , (7.45)
e
∂u uP − u (0, yp )
= , (7.46)
∂x w ∆x/2
porque a distância horizontal entre P e 0 é ∆x/2. Assim, a equação discretizada correspondente a este
volume de controle é
uP − u (0, yp ) uE − uP uP − uS uN − uP
∆x = fP0 + fP1 up ∆x∆y,
2a (0, yp ) ∆y − ae ∆y + as ∆x − an
∆x ∆x ∆y ∆y
ou
ap uP + aE uE + aS uS + aN uN = bp , (7.47)
com
2a (0, yp ) ae as an
ap = 2
+ 2
+ 2
+ − fP1 , (7.48)
∆x ∆x ∆y ∆y 2
ae as an
aE = − 2
, aS = − 2 , aN = − 2 , (7.49)
∆x ∆y ∆y
2a (0, yp )
bp = fP0 + g (0, yp ) . (7.50)
∆x2
Fórmulas semelhantes são obtidas para volumes de controle adjacentes às demais fronteiras que não estejam
em um dos quatro cantos do domı́nio retangular. Para os volumes de controle nos cantos do retângulo,
precisamos fazer mais uma modificação. Por exemplo, para o volume de controle no canto superior esquerdo
temos
aw = a (0, yp ) , (7.51)
an = a (xp , 1) , (7.52)
e
∂u uP − u (0, yp )
= , (7.53)
∂x w ∆x/2
∂u u (xp , 1) − uP
= , (7.54)
∂y n ∆y/2
e a equação discretizada correspondente a este volume de controle é
uP − u (0, yp ) uE − uP uP − uS u (xp , 1) − uP
∆x = fP0 + fP1 up ∆x∆y,
2a (0, yp ) ∆y−ae ∆y+as ∆x−2a (xp , 1)
∆x ∆x ∆y ∆y
ou
ap uP + aE uE + aS uS = bp , (7.55)
com
2a (0, yp ) ae as 2a (xp , 1)
ap = 2
+ 2
+ 2
+ − fP1 , (7.56)
∆x ∆x ∆y ∆y 2
ae as
aE = − 2
, aS = − 2 , (7.57)
∆x ∆y
2a (0, yp ) 2a (xp , 1)
bp = fP0 + g (0, yp ) + g (xp , 1) . (7.58)
∆x2 ∆x2
Ordenando os volumes de controle (por exemplo, usando a ordem lexicográfica), obtemos um sistema linear
cuja solução será uma solução aproximada para a equação com as condições de fronteira dadas.
Rodney Josué Biezuner 149
7.5 Exemplo (Equação de Poisson). Vamos aplicar o método de volumes finitos à equação de Poisson com
condição de fronteira de Dirichlet
2
−∆u = f (x, y) em [0, 1] ,
2 (7.59)
u = g (x, y) sobre ∂ [0, 1] .
Temos a (x) ≡ 1, fP1 = 0, fP0 = fP , e optamos por discretizar a malha por volumes de controle quadrados,
isto é, satisfazendo ∆x = ∆y. Segue que a linha do sistema discretizado corresponde a um volume 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 à fronteira, não localizados nos cantos, a linha correspondente no sistema
discretizado é
elemento na diagonal: ap = 5,
elementos fora da diagonal: a∗ = −1 (3 elementos),
elemento constante: bP = fP ∆x2 + 2g (∗) .
elemento na diagonal: ap = 6,
elementos fora da diagonal: a∗ = −1 (2 elementos),
elemento constante: bP = fP ∆x2 + 2g (∗) + 2g (∗∗) .
Compare com o correspondente sistema discretizado obtido pelo método de diferenças finitas; como no caso
unidimensional, as diferenças surgem apenas para as linhas correspondentes a células e pontos na fronteira
do domı́nio.
fP1 6 0. (7.60)
A necessidade matemática desta escolha já foi discutida no Exemplo 4.4. Fisicamente, esta exigência também
faz sentido: a maioria dos termos fontes em fenômenos transientes que tendem a um estado estacionário em
geral têm derivada primeira negativa, caso contrário o sistema não tenderia a um regime permanente. Por
exemplo, na difusão do calor, a existência de um termo linear com derivada positiva implicaria na acumulação
de energia térmica dentro do domı́nio, a não ser que o calor pudesse ser rapidamente dissipado através da
fronteira, o que geral não ocorre, pois mesmo o calor perdido por um objeto quente através da sua imersão
em um recipiente cheio de lı́quido frio é transferido para o lı́quido a uma taxa linear. Isso tende a gerar
uma situação instável que eventualmente leva ao colapso térmico do sistema (explosão ou derretimento do
objeto).
Rodney Josué Biezuner 150
fP0 = Auk−1
P +B e fP1 = 0, (7.62)
[Asmar] Nakhlé ASMAR, Partial Differential Equations and Boundary Value Problems, Pren-
tice Hall, 2000.
[Biezuner] Rodney Josué BIEZUNER, Notas de Aula: Equações Diferenciais Parciais, UFMG,
2005.
[BHM] William L. BRIGGS, Van Emden HENSON e Steve F. McCORMICK, A Multigrid
Tutorial, SIAM, 2000.
[Butcher] J. C. BUTCHER, Numerical Methods for Ordinary Differential Equations, 2nd. Ed.
Wiley, 2008.
[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.
[HNW] E. HAIRER, S.P. NORSETT e G. WANNER, Solving Ordinary Differential Equations
I: Nonstiff problems, 2nd. Ed., Springer Series in Computational Mathematics 8,
Springer-Verlag, 1993.
[Heuveline] Vincent HEUVELINE, On the computation of a very large number of eigenvalues for
selfadjoint elliptic operators by means of multigrid methods, Journal of Computational
Physics 184 (2003), 321–337.
[Horn-Johnson] Roger A. HORN e Charles R. JOHNSON, Matrix Analysis, Cambridge University
Press, 1985.
[Iserles] A. ISERLES, A First Course in the Numerical Analysis of Differential Equations,
2nd. Ed., Cambridge Texts in Applied Mathematics, Cambridge University Press,
2008.
[Maliska] CLOVIS R. MALISKA, Transferência de Calor e Mecânica dos Fluidos Computaci-
onal, 2a. Edição, LTC, 2004.
[Patankar] S. V. PATANKAR, Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980.
[Plato] R. PLATO, Concise Numerical Mathematics, Graduate Studies in Mathematics 57,
American Mathematical Society, 2003.
[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 of
arbitrary proportions, Zeitschrift für Angewandte Mathematik und Physik (ZAMP)
28 (1977), no.2, 185–196.
151
Rodney Josué Biezuner 152
[Young] David M. YOUNG, Iterative Solutions of Large Linear Systems, Academic Press,
1971.