Método das Diferenças Finitas

Capıtulo 2

Metodo de DiferencasFinitas

O metodo de diferencas finitas e uma das varias tecnicas para a diferenciacaode uma funcao discreta, i.e., um conjunto discreto de valores da variavel depen-dente em pontos conhecidos da variavel independente. O calculo do declive naoe mais do que o calculo discreto de dy/dx num dado ponto x0. O seu valordy/dx = (y1 − y0)/(x1 − x0) corresponde a uma formula de diferencas finitaspara o calculo da primeira derivada. Neste capıtulo serao apresentadas algumastecnicas para gerar formulas de diferencas finitas.

2.1 Expansao em serie de Taylor

Seja φ(x)uma funcao definida no intervalo a ≤ x ≤ b, e suponha que o inter-valo [a, b] contem o conjunto x0 = a, x1 . . . , xi . . . , xN+1 = b. A representacaodiscreta de φ(x) e [φi] = [φ(a), φ(x1), . . . , φ(xi), . . . , φ(b)] onde o valor deφ(xi) se representa por φi. O valor de xi+1 − xi, representa o espacamentoda malha e para uma malha uniforme xi = a + i∆x com i = 0, . . . , N + 1e ∆x = (b − a)/(N + 1). A figura 2.1 representa uma malha computacionalpara um problema uni-dimensional e a figura 2.2 uma malha estruturada carte-siana ortogonal para um problema bidimensional. Qualquer ponto (xi, yi) ficarepresentado na malha por (i, j) e os vizinhos relativamente a esse ponto vemrepresentados por (i± 1, j ± 1).

Figura 2.1: Malha Computacional Uni-dimensional

2-1

2-2 CAPITULO 2. METODO DE DIFERENCAS FINITAS

Figura 2.2: Malha computacional bidimensional

Formulas de diferencas finitas podem ser obtidas atraves do desenvolvimentoem serie de Taylor (ver por exemplo Collatz (1966), Ames (1977),Roache (1976),e Lagidus e Piner (1982) ).

Por exemplo, aproximar ∂φ/∂x no ponto (xi, yi) por uma diferenca discretapara um valor finito de ∆x, expande-se em serie de Taylor:

φ(xi + ∆x, yi) =∞∑m=0

(∆x)m

m!∂mφ(xi, yi)

∂xm(2.1)

que para as equacoes a direita e a esquerda do ponto (xi, yi):

φ(xi + ∆x, yi) = φ(xi, yi) +∂φ

∂x

∣∣∣∣i

∆x+∂2φ

∂x2

∣∣∣∣i

∆x2

2!+∂3φ

∂x3

∣∣∣∣i

∆x3

3!+ . . . (2.2)

φ(xi −∆x, yi) = φ(xi, yi)−∂φ

∂x

∣∣∣∣i

∆x+∂2φ

∂x2

∣∣∣∣i

∆x2

2!+∂3φ

∂x3

∣∣∣∣i

∆x3

3!+ . . . (2.3)

Se pretendemos obter no ponto (i, j) a primeira derivada atraves de diferencas a jusanteou, progressivas temos, por (2.2):

∂φ

∂x

∣∣∣∣i

=(φi+1,j − φi,j)

∆x+O

(∆x,

∂2φ

∂x2

)(2.4)

ou, atraves de diferencas a montante ou regressivas ou ainda apelidadas de”UPWIND”, por (2.3):

∂φ

∂x

∣∣∣∣i

=(φi,j − φi−1,j)

∆x+O

(∆x,

∂2φ

∂x2

)(2.5)

2.1. EXPANSAO EM SERIE DE TAYLOR 2-3

Para obter uma formula de diferencas centrais subtrai-se (2.3) de (2.2) e aserie e truncada:

∂φ

∂x

∣∣∣∣i

=(φi+1,j − φi−1,j)

2∆xi+O

(∆x2,

∂3φ

∂x3

)(2.6)

Logo ha um erro induzido pela truncatura da serie. Se mais termos da seriede Taylor fossem considerados seria necessario aproximar ∂3φ

∂x3 ,∂5φ∂x5 , . . ., etc... que

levava a considerar mais pontos e mais series de Taylor.Note-se que diferencas a jusante ou montante sao aproximacoes de 1ª ordem

dado que, por (2.4) e (2.5) vem que o(∆xn, ∂mφ∂xm ) para n = 1 significa que e

proporcional a ∆x e para m = 2( que a aproximacao e exacta para φ = a+ bxou seja a aproximacao e de 1ª ordem ou linear). Da analise do erro de truncaturapara diferencas centrais, n = 2 e m = 3 (e exacto para φ = a + bx + cx2 que euma parabola, ou seja, de 2ª ordem). Assim, pode-se designar a ordem ou nooperador ∆xn, por ordem n, ou na derivada ∂mφ/∂xm, por ordem (m− 1).

Igualmente obtem-se a aproximacao da 2ª derivada, ∂2φ/∂x2:

∂2φ

∂x2=

(φi+1,j + φi−1,j − 2φi,j)∆x2

+O(

∆x2,∂4φ

∂x4

)(2.7)

que e de 2ª ordem de precisao no operador ∆x e de 3ª ordem na derivada(entendido no grau maximo do polinomio em que a aproximacao e exacta).Para ilustrar o uso da expansao em serie de Taylor para aproximar expressoescom derivadas parciais cruzadas consideramos ∂2

∂x∂y . As expansoes em serie deTaylor para as duas variaveis x e y,

φ(x+ ∆x, y + ∆y) =φ(x, y) + ∆x∂φ

∂x+ ∆y

∂φ

∂y+

(∆x)2

2!∂2φ

∂x2+

(∆y)2

2!∂2φ

∂y2

+ 2∆x∆y

2!∂2φ

∂x∂y+O

[(∆x)3, (∆y)3

](2.8)

φi+1,j+1 =φi,j + ∆x∂φ

∂x+ ∆y

∂φ

∂y+ ∆x∆y

∂2φ

∂x∂y+

(∆x)2

2∂2φ

∂x2

(∆y)2

2

+∂2φ

∂y2+O

[(∆x)3, (∆y)3

] (2.9)

Expressoes analogas podem ser obtidas para φi−1,j−1;φi+1,j−1;φi−1,j+1 ori-ginando:

∂2φ

∂x∂y=φi+1,j+1 − φi+1,j−1 − φi−1,j+1 + φi−1,j−1

4(∆x)(∆y)+O

[(∆x)2, (∆y)2

](2.10)

de um modo semelhante pode-se recorrer a:


∂2φ

∂x∂y=

∂

∂x

(∂φ

∂x

)(2.11)

∂2φ

∂x∂y=

∂

∂x

[φi,j+1 − φi,j−1

2∆y+O(∆y)2

]=

12∆y

[∂φ

∂x

∣∣∣∣i,j+1

− ∂φ

∂x

∣∣∣∣i,j−1

]+O(∆y)2

(2.12)

Usando diferencas centrais para ∂φ∂x conduz a:

∂2φ

∂x∂y=

12∆y

[φi+1,j+1 − φi−1,j+1

2∆x− φi+1,j−1 − φi−1,j−1

2∆x

]+O

[(∆y)2, (∆y)2

] (2.13)

O mesmo procedimento e valido para derivadas temporais, ou seja:

∂φ

∂t= (φn+1 − φn)/∆t+O

(∆t,

∂2φ

∂t2

)(2.14)

em que os indıces n e n + 1 designam dois nıveis temporais. o nıvel (n + 1)representa o futuro e o nıvel n, o presente, φn e conhecido.

Para encontrar formulas de diferencas finitas que envolvam o valor da funcaoem pontos desejados e ao mesmo tempo maximizem a ordem de aproximacao,e.g.Moin(2000) (diminuicao do erro de truncatura), escreve-se a formula desejadana forma:

φi +n∑k=0

akfi+k = O(mınimo possıvel com pontos envolventes) (2.15)

O procedimento para determinar os coeficientes ak, para um dado numero depontos (n+1) da malha e o seguinte: fazem-se combinacoes lineares da serie deTaylor de modo que a funcao nos pontos considerados, seja expressa em funcaodo valores e suas derivadas, no ponto em que se pretende calcular a derivadaatraves da formula de diferencas. Por exemplo se pretender calcular φ

′

i usandoos pontos i, i+ 1, i+ 2 obtendo-se a tabela 2.1

Para o exemplo de considerar os tres pontos a equacao (2.15) e escrita coma ajuda da tabela 2.1 da seguinte forma:

2.1. EXPANSAO EM SERIE DE TAYLOR 2-5

φi φ′

i φ′′

i φ′′′

i

φ′

i 0 1 0 0a0φ

′

i a0 0 0 0a1φ

′

i+1 a1 a1∆x a1∆x2 a1

∆x3

6

a2φ′

i+2 a2 2a2∆x a2(2∆x)2

2 a2(2∆x)3

6

Tabela 2.1: Combinacao linear da serie de Taylor

φ′

i +2∑k=0

akφi+k =(a0 + a1 + a2)φ1

+(1 + a1∆x+ 2a2∆x)φ′

1

+(a1

∆x2

2+ a2

(2∆x)2

2

)φ′′

i

+(a1

∆x3

6+ a2

(2∆x)3

6

)φ′′′

i + . . .

(2.16)

Para obter a maior ordem de precisao e obviamente necessario diminuir oerro de truncatura o que se consegue eliminando os termos da equacao (2.16)que contem φi, φ

′

i e φ′′

i assim:

a0 + a1 + a2 = 0a1∆x+ 2a2∆x = −1a1∆x2/2 + 2a2∆x2 = 0

e obtem-se deste sistema de equacoes os coefiecientes:

a0 =3

2∆x; a1 = − 2

∆x; a2 =

12∆x

Substituindo estes valores em (2.16):

φ′

i =−3φi + 4φi+1 − φi+2

2∆x+O(∆x2) (2.17)

e o erro de truncatura e obtido substituindo a1 e a2 no ultimo termo daequacao (2.16). Com os pontos especificados nao e possıvel encontrar umaformula de diferencas que seja superior a segunda ordem de precisao.

Um segundo procedimento para aproximar uma derivada, consiste em re-presentar a funcao por um polinomio de ordem desejada. Os coeficientes dopolinomio sao calculados no referencial local centrado no ponto em que se de-seja obter a derivada discreta. Por exemplo, considere-se o polinomio de segundaordem

φ(x) = Ax2 +Bx+ C (2.18)


Seleccionando a origem em xi. Assim xi = 0, xi+1 = ∆x e xi+2 = 2∆x e osvalores da funcao nesses pontos sao

φ(xi) = φi;φ(xi+1) = φi+1 e φ(xi+2) = φi+2

assim

φi = Ax2i +Bxi + C = C

φi+1 = Ax2i+1 +Bxi+1 + C = A(∆x2) +B(∆x) + C

φi+2 = A2xi+2 +Bxi+2 + C = A(2∆x)2 +B(2∆x) + C

e os coeficientes sao dados por:

C = φiB = −φi+2+4φi+1−3φi

2∆x

A = φi+2−2φi+1+φi

2∆x

2

A primeira derivada de φ vem

∂φ

∂x= 2Ax+B

para xi = 0, ∂φ∂x∣∣∣i= B entao

∂φ

∂x=−φi+2 + 4φi+1 − 3φi

2∆x(2.19)

que e igual a aproximacao de diferencas finitas obtida por aproximacao progres-siva de serie de Taylor. A segunda derivada vem ∂2φ

∂x2 = 2A,

∂2φ

∂x2=φi+2 − 2φi+1 − φi

∆x

2

(2.20)

Se o espacamento entre os pontos i, i+1 e i+2, nao e igual, um procedimentoidentico e seguido mas agora e preciso ter em conta o espacamento nao uniformeda malha. Como sera analizado em 2.5.

2.2. CONSTRUCAO GERAL DE FORMULAS DE D. F. 2-7

2.2 Construcao Geral de Formulas de DiferencasFinitas

A ordem da aproximacao aumenta com o numero de pontos. O procedimentogeral para gerar formulas de diferencas finitas de qualquer ordem de precisaopode ser encontrada em Hirsch (1986). Definindo operadores que, por exemplo,para diferencas finitas a jusante:

∂φ

∂x

∣∣∣∣i,j

=φi+1,j−1 − φi,j

∆x+O(∆x) =

δ+φi∆x

+O(∆x) (2.21)

Assim o operador para diferencas a jusante

δ+φi = φi+1 − φi (2.22)

A tabela 2.2 mostra alguns dos operadores mais usados

Operador Sımbolo Expressao Relacoes

deslocamento E Eφi = φi+1

montante δ− δ−φi = φi−1 δ+ = E − 1central δ δφi = φi+1/2 − φi−1/2 δ− = 1− E−1

medio µ µφi = 12 (φi+1/2 − φi−1/2) µ = 1

2 (E1/2 + E−1/2)diferencial ∆ ∆φ = ∂φ

∂x

Tabela 2.2: Operadores e suas relacoes

As relacoes entre os operadores tambem estao listadas na tabela 2.2

En = φi+n (2.23)δ = E1/2 − E1/2 (2.24)

µ =12(E1/2 + E1/2) (2.25)

δ+δ+δ+ . . . δ+ = (E − 1)n (2.26)

Considerando a serie de Taylor:

φ(x+ ∆x) = φ(x) + ∆x∂φ(x)∂x

+∆x2

2!∂2φ(x)∂x2

+∆x3

3!∂3φ(x)∂x3

escrita com os operadores

Eφ(x) =(

1 + ∆xD +(∆xD)2

2!+

(∆xD)3

3!+ . . .

)φ(x) (2.27)

ou ainda atendendo ao desenvolvimento da exponencial

Eφ(x) = e∆xDφ(x) com E = e∆xD (2.28)


e tomando logaritmos∆xD = lnE (2.29)

Escolhendo diferencas progressivas e necessario substituir E = 1 + δ+,

δxD = lnE = ln(1 + δ+) = δ+ − (δ+)2

2+

(δ+)3

3− (δ+)4

4+ . . . (2.30)

Considerando diferencas regressivas para substituir E = 1 − δ na equacao(2.20) obtem-se:

∆xD = δ− +(δ−)2

2+

(δ−)3

3+

(δ−)4

4+ . . . (2.31)

e para diferencas centrais vira:

δφi = φi+ 1/2− φi−1/2 =(E1/2 − E−1/2

)φi

δ = e∆x/2 − e−∆x/2 = 2 sinh ∆x2

∆xD = δ − δ3

24+

3δ5

640− 7δ7

7168+ . . . (2.32)

Desprezando os termos na equacao (2.21) os termos da serie de potencias de(δ+)2 e superiores obtem-se ∆xD = δ+, ou seja:

∂φ

∂x

∣∣∣∣i,j

= Dφi =φi+1,j−1 − φi,j

∆x+O

(∆x,

∂2φ

∂x2

)(2.33)

se ∆xD = δ+ − (δ+)2

2 entao:

∂φ

∂x

∣∣∣∣i

=−3φi + 4φi+1 − φi+2

2∆x+O

(∆x2,

∂3φ

∂x3

)(2.34)

Para obter formulas de diferencas regressivas de realiza-se um procedimentosemelhante a partir da equacao (2.29).

Para obter formulas de diferencas centrais para a primeira derivada, a partirda equacao (2.32).

∆xD = δ =12(φi+1 − φi−1) (2.35)

∂φ

∂x=

12∆x

(φi+1 − φi−1) (2.36)

se forem introduzidos mais termos da equacao:

∆xD = δ − δ δ3

24(2.37)

∂φ

∂x=−φi+2 + 8φi+1 − 8φi−1 + φi−2

12∆x+O

(∆x4,

∂5

∂x5

)(2.38)


Derivadas de ordem elevada usando diferencas progressivas sao obtidas apartir de:

∂nφ∂xn = Dnφi = 1

∆xn [ln(1 + δ+)]nφi= 1

∆xn

[δ+n − n

2 δ+(n+1) + n(3n+5)

2n δ+(n+2)

−n(n+2)(n+3)48 δ+(n+3) + . . .

]φi

(2.39)

Diferencas regressivas

∂nφ∂xn = − 1

∆xn [ln(1− δ−)]nφi = 1∆xn

(δ− + δ−2

2 + δ−3

3 + . . .)n

φi

= 1∆xn

[δ−n + n

2 δ−(n+1) + n(3n+5)

24 δ−(n+2)

+ n(n+2)(n+3)48 + δ−(n+3)+...

]φi

(2.40)Ou centrais

Dnφi =(

2∆xsinh

−1 δ2

)nφi = 1

∆xn

[δ − δ3

24 + 3δ5

640 −5δ7

7168 + . . .]nφi

= δn

∆xn

[1− n

24δ2 + n

64

(22+5n

90

)δ4 − n

45

(57 + n−1

5 + (n−1)(n−2)35

)δ6 + . . .

]φi

(2.41)A tabela 2.3 lista os coeficientes da formula de diferencas finitas ate a 4ª

derivada e 6ª ou 8ª ordem de precisao para uma funcao discreta no ponto i.


i− 4 i− 3 i− 2 i− 1 i i+ 1 i+ 2 i+ 3 i+ 4

1derivada2 − 1

2 0 12

4 112 − 2

3 0 23 − 1

126 − 1

60320 − 3

4 0 34 − 3

20160

8 1290 − 4

10515 − 4

5 0 45 − 1

54

105 − 1280

2derivada2 1 -2 14 − 1

12 − 43 − 5

243 − 1

126 1

90 − 320

32 − 49

1832 − 3

20190

8 − 1560

8315 − 1

585 − 205

7285 − 1

58

315 − 1560

3derivada2 − 1

2 1 0 -1 12

4 18 −1 13

8 0 − 138 1 − 1

86 − 7

240310 − 169

1206130 0 − 61

30169120 − 3

107

2404derivada

2 1 −4 6 4 14 − 1

6 2 − 132

283 − 13

2 2 − 16

6 7240 − 2

516960 − 122

15918 − 122

1516960 − 2

57

240

Tabela 2.3: Coeficientes de formulas de diferencas finitas centrais em malhauniforme


Exemplo 2-1Obter formulas de fiferencas finitas regressivas, progressivas ou centrais para

a segunda derivada mantendo so o primeiro termo da expansao ou mantendo osdois primeiros termos.

(φxx)i =1

∆x2

(δ+2 − δ+3 +

1112δ+4 − 5

6δ+5 + . . .

)φi

(φxx)i =1

∆x2

(δ−2 − δ−3 +

1112δ−4 − 5

6δ+5 + . . .

)φi

(φxx)i =1

∆x2

(δ2 − δ4

12+δ6

90− δ8

560+ . . .

)φi

(φxx)i =1µ

∆x2

(δ2 − 5δ4

2x+

2595760

δ6 + o(∆x8))φi

Mantendo o primeiro termo da serie, vem:

(φxx)i =1

∆x2(φi+2 − 2φi+1 + φi)−∆xφxxx

(φxx)i =1

∆x2(φi − 2φi−1 + φi−2) + ∆xφxxx

(φxx)i =1

∆x2(φi+1 − 2φi + φi−1)−

∆x2

12

(∂4φ

∂x4

)(φxx)i =

12∆x2

(φi+3/2 − φi+1/2 − φi−1/2 + φi+3/2)

− 52x

∆x2

(∂4x

∂x4

)Mantendo dois termos da serie, vem para diferencas progressivas, regressivas

e centrais respectivamente:

(φxx)i =1

∆x2(2φi − 5φi+1 + 4φi+2 − φi+3) +

1112

∆x2

(∂4φ

∂x4

)(φxx)i =

1∆x2

(2φi − 5φi−1 + 4φi−2 − φi−3) +1112

∆x2

(∂4φ

∂x4

)(φxx)i =

112∆x2

(−φi+2 + 16φi+1 − 30φi + 16φi−1 − φi−2)

+∆x4

90

(∂6φ

∂x6

)Estas formulas sao importantes para a incorporacao de condicoes fronteira.


2.3 Ordem de Precisao das Formulas de Dife-rencas Finitas

O erro de truncatura da serie de Taylor e expresso pela ordem de grandeza doprimeiro termo desprezado da serie. A ordem da aproximacao e usualmente esta-belecida pela potencia, n, do intervalo discreto, ∆xn. Por exemplo Θ(∆x2, ∂

4φ∂x4 )

significa que a formula de diferencas finitas e de segunda ordem, ∆x2. contudoe importante ter em conta a ordem da derivada porque no erro de truncatura a4ª derivada significa que se a solucao exacta fosse dada por um polinomio do 3ºgrau o erro de truncatura seria nulo. Contudo a solucao analıtica e desconhecidae assim nao e quantificavel o erro de ∂4φ/∂x4. Da analise do erro de truncaturapode-se concluir o seguinte:

1. Se ha varias variaveis independentes, cada uma tera um erro de truncaturapor exemplo Θ(∆x2;∆t), ou seja, de segunda ordem em x e de primeiraordem em t.

2. A ordem do esquema numerico (formula de diferencas finitas), dependede propriedades locais da funcao. Por exemplo proximo de gradientes ele-vados de φ pode originar uma ordem de precisao menor do que a ordemteorica calculda pela expansao em serie de Taylor. Diferentes erros contri-buem para a solucao final podendo a ordem da precisao da solucao medidanuma dada norma, ser diferente da ordem de erro local num dado ponto.

3. E desejavel usar formulas de diferencas finitas de ordem elevada. Contudo,como veremos, o erro da formula de diferencas nao significa que se conhecao erro da solucao. As equacoes sao nao lineares e nao se conhece, em geral,como e que um erro local vai evoluir durante o processo iteractivo ate aobtencao da solucao final. A ordem de precisao estabelece como e que oerro vai variar com o refinamento da malha. Assim por exemplo para umaprecisao de segunda ordem, significa que o erro R obtido com a malha∆x sera de R/4 com a malha Θ(∆x/2)2. E sera de R/16 com a malha(∆x/4)2. De modo identico uma ordem de precisao de primeira ordemconduziria ao erro R com uma malha de espacamento ∆x e conduzira aoerro de R/2 com uma malha de espacamento ∆x/2. Para uma malha de∆x/n implicaria um erro R/n.

4. O erro da truncatura e cumulatativo, ou seja, para problemas por exemplonao estacionarios leva a que durante a evolucao temporal os erros vaoaumentando devido a erros em cada nıvel temporal.

5. O erro de truncatura em aplicacoes praticas de engenharia e quase sempremuito maior ao erro de truncatura da maquina, devido a um numero debits finito.

2.3. ORDEM DE PRECISAO DAS FORMULAS DE D. F. 2-13

Exemplo 2-2Para ilustrar os efeitos dos erros vamos considerar a funcao

f(x) =sinx

x3

num domınio discreto com incremento espacial h. Pretendemos obter o erro daformula de diferencas finitas para o calculo da 1ª derivada ∂t/∂x com:

1. Diferencas regressivas de 1ª ordem

2. Diferencas centrais de 2ª ordem

3. Diferencas centrais de 4ª ordem

Considere incrementos espaciais h = 10, 1, 0.1, 10−2, 10−3, 10−4, 2× 10−5. Afigura 2.3 mostra os resultados obtidos.

Figura 2.3: Ordem de precisao do exemplo 2.1

O decaimento do erro com a diminuicao de h esta de acordo com a ordemdo erro de truncatura da formula de diferencas finitas. A ordem de precisao dauma indicacao da precisao das formulas de diferencas finitas e diz como e que orefinamento da malha melhora a precisao.


2.4 Malhas nao uniformes

Quando a variavel dependente e uma funcao nao linear no domınio computa-cional podera ser um desperdıcio usar malhas muito finas para discretizar todoo domınio. Sera mais conveniente usar uma malha nao uniforme para refinarnas regioes onde se espera que a solucao contenha gradientes mais elevados. Seos pontos da malha nao estao uniformemente espacados vem por exemplo parao calculo de ∂φ/∂x, em pontos xi+1 − xi = Ax e xi − xi−1 = (1 + α)∆x

φ(x+ ∆x) = φ(x) + ∆x∂φ

∂x+

(∆x)2

2!∂2φ

∂x2+

(∆x)3

3!∂3φ

∂x3+O

((∆x)4

)(2.42)

φ[x+ (1 + α)∆x] = φ(x) + (1 + α)∆x+∂φ

∂x+

(1 + α)2(∆x)2

2!∂2φ

∂x2+

(1 + α)3(∆x)3

3!+∂3φ

∂x3+O

((∆x)4

)(2.43)

Multiplicando a equacao (2.42) por −(1+α) e adicionando a equacao (2.43)vem

φ′′

i = 2[

φi−1

∆x(1 + α) (∆x(1 + α) + ∆x)− φi

∆x(1 + α) + ∆x+

φi−1

∆x (∆x(1 + α) + ∆x)

](2.44)

Uma outra tecnica de ter em conta a nao uniformidade da malha consisteem transformar o domınio discreto.

Para ∂2φ∂x2 a formula de diferencas finitas e dada transformando 0 ≤ x ≤ 1 em

0 ≤ ξ ≤ π2 . Outra transformacao utilizada e, por exemplo, ξ = cos−1x. Para

um espacamento uniforme em ξ

ξj =π

2Nj j = 0, 1, 2, . . . , N (2.45)

corresponde a uma malha nao uniforme em x em que proximo de x = 1 amalha sera muito fina e proximo de x = 0, a malha sera muito espacada, verfigura 2.4

De um modo geral a transformacao

ξ = g(x) (2.46)

leva a que as derivadas no novo sistema de coordenadas sejam dadas por:

df

dx=dζ

dx

df

dζ= g′

df

dζ(2.47)

d2f

dx2=

d

dx

[g′df

dζ

]= g′′

df

dζ+ (g′)2

d2f

dζ2(2.48)

2.4. MALHAS NAO UNIFORMES 2-15

Figura 2.4: Transformacao de coordenadas

assim diferencas finitas para malhas uniformes sao usadas para ∂f∂ξ e ∂2f

∂ξ2 .


2.5 Condicoes Fronteiras de Dirichlet ou de vonNeumann

As condicoes fronteiras podem classificar-se como pertencendo a 3 classes:Dirichlet, von Neumann e Mistas. Diz-se que uma condicao fronteira e do tipode Dirichlet quando a variavel dependente e dada explıcitamente, mesmo queesteja a variar por exemplo com o tempo. ϕ(x, t) para ϕ(0, t) = C0(t) em que C0

e uma funcao conhecida de t. Quando a derivada normal a fronteira e conhecida,a condicao fronteira e do tipo de von Neumann. Por exemplo para a conducaode calor estacionaria com h, coeficiente de conveccao, e K, condutividade.

A equacao do calor unidimensional:

dφ1

dt= k

∂2φ

∂x2(2.49)

na equacao de diferencas finitas:

φn+1i = φni + r

(φni−1 − 2φni + φni+1

), r =

∆tK∆x2

(2.50)

para uma condicao fronteira de conveccao ou de Von Newmann:

−∂φ∂x

= − h

K(φ− φref ) (2.51)

com diferencas a jusante,

φ1 − φ0

∆x=

h

K(φ0 − φref ) (2.52)

Se fosse com diferencas centrais ha a necessidade de introduzir um pontovirtual i = −1.

φ1 − φ1

∆x=

h

K(φ0 − φref ) (2.53)

φ1 pode ser eliminado atraves da aplicacao da equacao de diferencas finitas parao ponto i = 0, ou seja:

φ−1 = −h∆xK

(φ0 − φref ) + φ1 (2.54)

φn+10 = φn0 + r

(φn1 −

h∆xK

(φn0 − φref )− 2φn0 + φn1

)(2.55)

O mesmo procedimento pode ser aplicado a fronteira para i = N. Por exem-plo:

∂φ

∂n

∣∣∣∣N

= CB(t) (2.56)

onde n e medido na direccao perpendicular a fronteira. Temos como exemplo

2.5. C. F. DE DIRICHLET OU DE VON NEUMANN 2-17

as condicoes fronteira de conveccao num problema de transmissao de calor. Acondicao fronteira do tipo misto pode ser expressa por:

aF (t)φF + bF (t)∂φ

∂n

∣∣∣∣F

= cF (t) (2.57)

onde as funcoes aF , bF , cF podem depender da posicao F na fronteira e aF ebF tem o mesmo sinal. Por exemplo, se a condicao fronteira e de:

∂φ

∂n

∣∣∣∣0,t

= C0(t) (2.58)

vem que a0(t) = 0 e b0(t) = −1 para um problema mais concreto como seja deconveccao nas fronteiras:

k∂φ

∂x= h(φ− φref ) (2.59)

∂φ

∂x=h

k(φ− φref ) (2.60)


2.6 Estrutura da Matriz de Coeficientes

A equacao de Laplace e de Poisson sao exemplos tıpicos de equacoes elıpticasque em coordenadas cartesianas na forma bidimensional sao expressas respecti-vamente por:

∂2φ

∂x2+∂2φ

∂y2= 0 (2.61)

∂2φ

∂x2+∂2φ

∂y2= f(x, y) (2.62)

Aproximando por diferencas centrais as segundas derivadas, vem que:

φi+1,j − 2φi,j + φi−1,j

(∆x)2+φi,j+1 − 2φi,j + φi,j−1

(∆y)2= 0 (2.63)

A molecula computacional de cinco pontos esta representada na figura 2.5.

Figura 2.5: Molecula computacional envolvendo cinco pontos

A equacao de Laplace tambem poderia ser discretizada por diferencas finitasde ordem mais elevada como por exemplo por diferencas centrais de 4ª. ordem

−φi−2,j + 16φi−1,j − 30φi,j + 16φi+1,j − φi+2,j

12(∆x)2+

+−φi,j−2 + 16φi,j−1 − 30φi,j + 16φi,j+1 − φi,j+2

12(∆y)2= 0

(2.64)

em que a molecula computacional esta representada na figura 2.6A discretizacao de segundas derivadas so raramente e que e feita com formulas

de diferencas de 4ª. ordem. Rearranjando os termos, vem que, para segundaordem de precisao:

2.6. ESTRUTURA DA MATRIZ DE COEFICIENTES 2-19

Figura 2.6: Molecula computacional envolvendo nove pontos

φi+1,j − 2φi,j + φi−1,j +(

∆x∆y

)2

(φi,j+1 − 2φi,j + φi,j−1) = 0 (2.65)

Definindo a razao do espacamento da malha por β, tal que, β = ∆x∆y a equacao

2.66 aparece como:

φi+1,j + φi−1,j + β2φi,j+1 + β2φi,j−1 − 2(1 + β2)φi,j = 0 (2.66)

Equacoes elıpticas requerem sempre a solucao de sistemas de equacoes. Eassim fundamental perceber a estrutura da matriz de coeficientes. Considera-mos uma malha de 6x6 como ilustrado na figura 2.7. As condicoes fronteira sao:

x = 0φ = u2, y = 0φ = u1

x = Lφ = u4, y = Hφ = u3

(2.67)

Figura 2.7: Malha computacional

Aplicando a equacao 2.66 a todos os pontos interiores do domınio computacional,


vem que:φ3,2 + φ1,2 + β2φ2,3 + β2φ2,1 − 2(1 + β2)φ2,2 = 0φ4,2 + φ2,2 + β2φ3,3 + β2φ3,1 − 2(1 + β2)φ3,2 = 0φ5,2 + φ3,2 + β2φ4,3 + β2φ4,1 − 2(1 + β2)φ4,2 = 0φ6,2 + φ4,2 + β2φ5,3 + β2φ5,1 − 2(1 + β2)φ5,2 = 0φ5,3 + φ1,3 + β2φ2,4 + β2φ2,2 − 2(1 + β2)φ2,3 = 0φ4,3 + φ2,3 + β2φ3,4 + β2φ3,2 − 2(1 + β2)φ3,3 = 0φ5,3 + φ3,3 + β2φ4,4 + β2φ4,2 − 2(1 + β2)φ4,3 = 0φ6,3 + φ4,3 + β2φ5,4 + β2φ5,2 − 2(1 + β2)φ5,3 = 0φ3,4 + φ1,4 + β2φ2,5 + β2φ2,3 − 2(1 + β2)φ2,4 = 0φ4,4 + φ2,4 + β2φ3,5 + β2φ3,3 − 2(1 + β2)φ3,4 = 0φ5,4 + φ3,4 + β2φ4,5 + β2φ4,3 − 2(1 + β2)φ4,4 = 0φ6,4 + φ4,4 + β2φ5,5 + β2φ5,3 − 2(1 + β2)φ5,4 = 0φ3,5 + φ1,5 + β2φ2,6 + β2φ2,4 − 2(1 + β2)φ2,5 = 0φ4,5 + φ2,5 + β2φ3,6 + β2φ3,4 − 2(1 + β2)φ3,5 = 0φ5,5 + φ3,5 + β2φ4,6 + β2φ4,4 − 2(1 + β2)φ4,5 = 0φ6,5 + φ4,5 + β2φ5,6 + β2φ5,4 − 2(1 + β2)φ5,5 = 0

ondeα = −2(1 + β2) (2.68)

A estrutura da matriz de coeficientes e pentadiagonal e e uma matriz defi-nida positiva. Nao existe na literatura uma notacao unica para o metodo dediferencas finitas. Ou se usa uma notacao indicial para localizar as variaveis nodomınio computacional, como a apresentada ate aqui, ou como varios autoresusam uma notacao de pontos cardeais, com Sul, Norte, Este e Oeste em tornodo ponto P . Em varios livros de texto a notacao para designar os pontos da ma-lha em volta de (i, j) e substituıda por pontos cardeais,N,S,E,W em torno doponto P como mostra a figura 2.7. Vamos seguidamente discretizar a equacaode conveccao difusao para um problema unidimensional.

Exemplo 2-3Obter a equacao de diferencas finitas e a matriz de coeficientes da equacao

de conveccao difusao:


∂

∂x(ρUφ) =

∂

∂x

(Γφ∂φ

∂x

)(2.69)

Para o problema unidimensional com e discretizacao por diferencas centrais[∂

∂x(ρUφ)

]i

=(ρUφ)i+1 − (ρUφ)i−1

xi+1 − xi−1

O termo difusivo tambem e discretizado por diferencas centrais

[∂

∂x(Γφ

∂φ

∂x)]i

=

(Γφ ∂φ∂x

)i+ 1

2

−(Γφ ∂φ∂x

)i− 1

212 (xi+1 − xi−1)

=2Γφ,i+ 1

2

φi+1−φi

xi+1−xi − Γφ,i− 12

φi−φi−1xi−xi−1

xi+1 − xi−1

O coeficiente difusivo so e conhecido nos pontos e e necessario o seu valorentre os nos:

Γφ,i+ 12

=12(Γφ,i + Γφ,i+1)

Γφ,i− 12

=12(Γφ,i + Γφ,i−1)

[∂

∂x(Γφ

∂φ

∂x)]i

=(Γφ,i + Γφ,i+1)φi+1 − φi

(xi+1 − xi)(xi+1 − xi− 1)−

−(Γφ,i + Γφ,i−1)φi − φi−1

(xi − xi− 1)(xi+1 − xi− 1)

Re-arranjando os termos, ficamos com a equacao de diferencas finitas paracada ponto

aipφi = aiEφi+1 + aiWφi−1

onde:

aiE = − (ρU)i+1

xi+1 − xi−1+

Γφ,i + Γφ,i+1

(xi+1 − xi)(xi+1 − xi−1)

aiW = − (ρU)i−1

xi+1 − xi−1+

Γφ,i + Γφ,i−1

(xi − xi−1)(xi+1 − xi−1)

aiP = − Γφ,i + Γφ,i+1

(xi+1 − xi)(xi+1 − xi−1)+

Γφ,i + Γφ,i−1

(xi − xi−1)(xi+1 − xi−1)

Pela equacao de continuidade temos


(ρU)i+1 = (ρU)i−1

que reservando o coeficiente da diagonal principal:

aiP = aiE + aiW

Na forma matricial o sistema de equacoes resultante da aplicacao da equacaode diferencas finitas a todos os pontos

12.

i, j.

N − 1N

1 2 · · · i · · · N × M26666666664

a1P −a1

E−a2

W a2P −a2

E 0. . .

−aiW ai

P −aiE

. . .0 . . .

−aNW aN

P

37777777775

266666664

φ1.

φi−1φi

φi+1.

φN

377777775=

2666666664

a1W0.0.0

aNE φN+1

3777777775

onde o valor de φ0 e φN+1 e conhecido atraves das condicoes fronteira Paraobter a solucao e necessario inverter a matriz triadiagonal.

Para um problema bidimensional de conveccao/difusao:

∂

∂x

(ρUφ− Γφ

∂φ

∂x

)+

∂

∂y

(ρV φ− Γφ

∂φ

∂y

)= qφ

Figura 2.8: Malha do Problema Bidimensional

Usando diferencas centrais para discretizar as primeiras e segundas derivadasobtem-se a equacao de diferencas finitas que se pode escrever na forma:

ai,jP φi,j = ai,jE φi+1,j + ai,jW φi−1,j + ai,jN φi,j+1 + ai,jS φi,j−1 +Qi,j

em que os coeficientes vem dados por:

aiE = − (ρU)i+1,j

xi+1 − xi−1+

Γφ,i,j + Γφ,i+1,j

(xi+1 − xi)(xi+1 − xi−1)

aiV = − (ρV )i,j+1

yi+1 − yi−1+

Γφ,i,j + Γφ,i,j+1

(yi+1 − yi)(yi+1 − yi−1)

aiW = +(ρU)i−1,j

xi+1 − xi−1+

Γφ,i,j + Γφ,i−1,j

(xi − xi−1)(xi+1 − xi−1)

aiS = +(ρV )i,j−1

yi+1 − yi−1+

Γφ,i,j + Γφ,i,j−1

(yi+1 − yi−1)(yi+1 − yi−1)


Discretizando a equacao da continuidade

∂(ρU)∂x

+∂(ρV )∂y

= 0

(ρU)i+1,j − (ρU)i−1,j

xi+1 − xi−1+

(ρV )i,j+1 − (ρV )i,j−1

yi+1 − yi−1= 0

De um modo analogo ao problema unidimensional e possıvel escrever o coe-ficiente aP como a soma dos coeficientes:

ai,jP =∑nb

ai,jnb, nb = E,W,N, S

Aφ = S e o sistema de equacoes na forma matricial e representado por:

12..

i, j..

N × M

1 2 · · · i · · · N × M26666666666664

a1,1P

a1,1N

a1,1E

−a1,2S

a1,2P

−a1,2N

a1,2E

. . .. . .

−ai,jW

−ai,jS

ai,jP

−ai,jN

ai,jE

. . . .. . . . .

−aN,MW

−aN,MS

aN,MP

37777777777775

266666666664

φ1,1φ1,2

.

.φi,j

.

.φN,M

377777777775=

266666666664

S1,1S1,2

.

.Si,j

.

.SN,M

377777777775

com:

S1,1 = Q1,1 + a1,1W φ0,1 + a1,1

S φ1,0

SN,M = QN,M + aN,ME φN+1,M + aN,MN φN,M+1


2.7 Aproximacao de Pade

O numero de pontos e sua localizacao ,que cada formula de diferencas finitasenvolve para o calculo de uma derivada, e conhecido como a molecula computa-cional. Pode-se perguntar qual sera a formula de calculo da derivada no pontoj, com um menor erro de truncatura envolvendo um dado numero, prescrito, depontos e tambem das derivadas f ′ em outros pontos. As aproximacoes Pade se-guem este raciocınio, fazendo envolver, por exemplo, f ′j , f

′j+1 e f ′j−1 e os valores

de funcao f,fj+1 e fj−1, ou seja dado

f ′j + a3f′j+1 + a4f

′j−1 = a0fj + a1fj−1 + a2fj−1 (2.70)

em que se pretende encontrar os coeficientes ai, i = 0, 1, . . . 4 de modo a maxi-mizar a ordem de aproximacao. A tabela 2.7 mostra os coeficinetes da expansaoem serie de Taylor, ou seja, cada linha apresenta os coeficientes das derivadasenvolvidas, por exemplo:

f ′ f ′′ f ′′ f ′′′ f iv fv

f ′j 0 1 0 0 0 0aof

′j a0 0 0 0 0 0

a1f′j+1 a1 a1h a1

h2

2 a1h3

6 a1h4

24 a1h5

120

a2f′j−1 a2 −a2h a2

h2

2 −a2h3

6 a2h4

24 −a2h5

120

a3f′j+1 0 a3 a3h a3

h2

2 a3h3

6 a3h4

24

a4f′j−1 0 a4 −a4h a4

h2

2 −a4h3

6 a4h4

24

Tabela 2.4: Coeficientes da serie de Taylor

a1fj+1 = +a1fj + a1hf′j + a1

h2

2f4j + a1

h3

6f ′′′ + . . . (2.71)

de modo a determinar a ordem de precisao mais elevada dever-se anular o maiornumero de termos de ordem pequena

a0 + a1 + a2 =0a1h− a2h+ a3 + a4 =− 1

a1h2

2+ a2

h2

2+ a3h− a4h =0

a1h

3− a2

h

3+ a3 + a4 =0

a1h

4+ a2

h

4+ a3 − a4 =0

(2.72)

a solucao do sistema (2.72) conduz a

a0 = 0; a1 = − 34h

; a2 =34h

; a3 = a4 =14

(2.73)

2.7. APROXIMACAO DE PADE 2-25

Substituındo na coluna e rearranjando conduz a formula de Pade

f ′j+1 + f ′j−1 + 4f ′j =3h

(fj+1 − fj−1) +h4

30f0j (2.74)

onde j = 1, 2, 3, . . . , n− 1. A equacao (2.74) conduz a um sistema de equacoeslineares para o calculo de f ′j . Ha (n − 1) equacoes para (n + 1) incognitas.Alteracao especial deve ser dada proximo das fronteiras. Usualmenete baixa-sea ordem e usam-se diferencas progressivas ou regressivas para f ′o e f ′n. Porexemplo:

f ′0 + 2f ′1 =1h

(−5

2f0 + 2f1 +

12f2

)f ′n + 2f ′n−1 =

1h

(52fn − 2fn−1 −

12fn−2

) (2.75)

Na forma matricial (2.74) e (2.75) aparecem como:

1 2 0 0 0 . . . 01 4 1 0 0 . . . 00 1 4 1 0 . . . 0...

.... . . . . . . . .

......

......

.... . . . . . . . .

...0 0 0 . . . 1 4 10 0 0 0 . . . 2 1

f ′0f ′1f ′2......

f ′n−1

f ′n

=

1h

− 52f0 + 2f1 + 1

23(f2 − f0)3(f3 − f1)

...

...3(fn − fn−2)

52fn − 2fn−1 − 1

2fn−2

(2.76)

A inversao da matriz triadiagonal e simples e rapida e permite o calculo daderivada em cada ponto. Este procedimento leva a uma ordem de precisaode 4ª ordem en que a aproximacao aumenta o domınio de numero de onda(frequencias) que se podem resolver, ver figura 2.9. Aproximacoes de Pade saoglobais porque fazem envolver todos os pontos na molecula computacional aolongo da direccao que se pretende a derivada. E possıvel estender a aproximacaoPade a qualquer outra derivada e por exemplo para a 2º derivada vem:

112f ′′i−1 +

1012f ′′i +

112f ′′i+1 =

fi+1 − 2fi + fi−1

h2(2.77)

Lele (1992) apresentou aproximacoes de Pade para a 1ª e 2ª derivadas ate a10ª ordem de aproximacao. Considerando a expressao (2.78):

β(φx)i−2,j + α(φx)i−1,j + (φx)i,j + α(φx)i+1,j + β(φx)i+2,j =

cφi+3,j − φi−3,j

6h+ b

φi+2,j − φi−2,j

4h+ a

φi+1 − φi−1

2h(2.78)

As tabelas 2.5 e 2.6 listam os coeficientes para as aproximacoes a 1ª e 2ªderivadas.


Nº

Dia

gona

isE

rro

daM

atri

zα

βa

bc

Ord

.de

Tru

ncat

ura

31 4

03 2(α

+2)

00

4ª4 5!(

3α−

1)h

4φ

5

31 3

03 2(α

+2)

1 4(4α−

1)0

6ª4 7!h

6φ

7

54 9

1 36

2 9(8−

3α)

1 18(5

7α−

17)

08ª

−16 9!h

8φ

9

51 2

1 36

1 6(1

2−

7α)

1150(5

86α−

183)

1 50(9α−

4)10

ª144

11!h

10φ

11

Tab

ela

2.5:

Coe

fient

esda

apro

xim

acao

dePad

epa

raa

1ªde

riva

da.

Nº

deD

iago

nais

Err

oM

atri

zα

βa

bc

Ord

.de

Tru

ncat

ura

31 10

04 3(1−α)

1 3(1

0k−

1)0

4ªh

4φ

6

32 11

012

11

3 11

06ª

h6φ

8

5344

1179

38α−

9214

696−

1191α

428

2454α−

294

535

08ª

h8φ

10

5344

899

43

1798

1065

1798

1038

1798

79

1798

10ª

h10φ

12

Tab

ela

2.6:

Coe

fient

esda

apro

xim

acao

dePad

epa

raa

2ªde

riva

da.

2.7. APROXIMACAO DE PADE 2-27

E de esperar que se a solucao for dada for uma funcao seno com umafrequencia elevada, obrigue a uma malha muito fina, ou seja, deve haver muitomais de que dois pontos por comprimento de onda caso contrario, a aproximacaode derivadas induzira um erro que nao e admissıvel.

O metodo de numero de onda aproximado permite verificar quais as frequenciasque a formula de diferencas finitas consegue resolver. E de esperar que sendo”φ”uma funcao sinosoidal, para baixas frequencias, ou seja grandes cumprimen-tos de onda, as formulas de diferencas finitas de baixa ordem tem suficienteresolucao para o calculo das derivadas. De modo oposto, para frequencias altas,so as formulas de alta ordem devem ser aplicadas e mesmo assim com forteslimitacoes de precisao. Para ilustrar o procedimento, considere uma funcao deuma harmonia de perıodo L.

f(x) = eikx (2.79)

onde K e o numero de onda (ou frequencia) que toma os seguintes valores

K =2πLn, n = 0, 1, 2, . . . , N/2 (2.80)

em que N e o numero de pontos em x. xi = LN j, j = 0, 1, 2, . . . , N − 1. A

derivada exacta de (2.79) e

df

dx= ikf (2.81)

Uma formula de diferencas finitas de 2ªordem central para a primeira deri-vada:

δf

δx

∣∣∣∣j

=fj+1 − fj−1

2h(2.82)

em que h = L/N e o incremento espacial do domınio de comprimento, L, como numero de pontos total N .

Substituindo (2.79) em (2.81) vem:

δf

δx

∣∣∣∣j

=ei2πn(j+1)/N − ei2πn(j−1)/N

2h=ei2πn/N − ei2πn/N

2hfj (2.83)

ou sejaδf

δx

∣∣∣∣j

= isin(2πn/N)

hfj = ik′fj (2.84)

k′ =sin(2πn/N)

h(2.85)

A equacao (2.84) mostra que a aproximacao a derivada conserva a mesmaforma excepto que k vem modificado por k′ (Numero de onda modificado). Deum modo analogo e possıvel obter valores de k′ para o presente caso da 1ªderivada. A figura 2.9 mostra a evolucao do numero de onda exacto, (linha


Figura 2.9: Numero de onda modificada

recta) com os valores de uma formula de diferencas finitas central de 2ª ordemou 4ª ordem de precisao e tambem com uma aproximacao Pade.

Para valores elevados de k a funcao varia rapidamente no domınio e asformulas de diferencas finitas fornecem uma precisao bastante fraca para a 1ªderivada. Da equacao 2.80 e de acordo com h = L/N vem que

k =2πh

n

Nn = 0, 1, 2, . . . , N/2

deste modo aumentando o numero de pontos N para, por exemplo, 2N permiteaumentar a gama de frequencia resoluveis pelo esquema numerico.

2.8. RESUMO DO CAPITULO 2-29

2.8 Resumo do Capıtulo

O aluno deve saber:

Gerar formulas de diferencas finitas usando Serie de Taylor ou aproximacaopolinomial ou formulas genericas.

Gerar formulas de diferencas finitas para a derivada de ordem m, pordiferencas progressivas ou regressivas ou centrais de ordem n

Gerar formulas de diferencas finitas de derivadas cruzadas

Gerar formulas de diferencas finitas em malhas nao uniformes

Saber o que significa a ordem de precisao e suas consequencias

Gerar os erros de truncatura de quaisquer formulas de diferencas

Saber obter aproximacoes a derivadas de 1ªordem ou 2ªordem por opera-dores compactos centrais

Saber obter uma equacao de diferencas finitas e incorporar as condicoesfronteira

Saber obter a matriz de coeficientes para a equacao de diferencas finitascom a incorporacao de condicoes fronteira


2.9 Referencias

[1] Ames, W. F. (1977). Numerical Methods for Partial Differential Equation,2nd ed., New York, Academic Press

[2] Collatz, L. (1966). The Numerical Treatment of Differential Equations, 3rded., Berlin, Springer Verlag

[3] Lapidus, L and Pinder G. (1982). Numerical Solution of Partial DifferentialEquation in Science and Engineering, Wiley

[4] Moin P.(2001). Fundamentals of Engineering Numerical Analysis, Cam-bridge Unviversity Press

[5] Roache, P. J. (1976). Computational Fluid Dynamics, Hermosa Pub., Al-buquerque, New Mexico

[6] Lele, S. K., J. Comput. Physics, 103, vol 16, 1992

[7] Hirsch, C. “Numerical Computation of Internal and External Flows Vol.1”, Wiley, 1986

2-31

2-32 2.9 REFERENCIAS

2.10 Problemas

P1. Aplique a serie de Taylor a expressao geral da forma:

(φx)i = aφi+2 + bφi+1 + cφi + dφi−1 + eφi−2

(a) Obtenha a aproximacao de precisao de 4ª. ordem para a primeiraderivada no ponto i.

(b) Repita o mesmo procedimento para a 2ª. derivada (uxx)i com omesmo numero de pontos. Calcule o erro de truncatura em ambosos casos.

P2. Obtenha os coeficientes a, b e c das formulas de diferencas centrais para:(∂φ

∂x

)i

=(−φi+2 + aφi+1 − cφi−1 + φi−2)

c+

(∆x4

30∂5φ

∂x5

)(∂2φ

∂x2

)i

=(−φi+2 + aφi+1 − bφi + cφi−1 − φi−2)

d∆x2+

190

∆x2

(∂6φ

∂x6

)(∂3φ

∂x3

)i

=(−φi+2 − aφi+1 + bφi−1 − φi−2)

c∆x2+

14∆x2

(∂5φ

∂x5

)(∂4φ

∂x4

)i

=(φi+2 − aφi+1 + bφi − cφi−1 + φi−2)

∆x4− ∆x2

6

(∂6φ

∂x6

)

P3. Obtenha os coeficientes a, b e c das formulas de diferencas finitas:

∂4φ

∂x4

∣∣∣∣i

=φi+2,j − aφi+1,j + bφi,j − cφi−1,j + φi−2,j

∆x4+O

(∆x2

)∂2φ

∂x2

∣∣∣∣i

=−φi+3,j + aφi+2,j − bφi+1,j + cφi,j

∆x2+O

(∆x2

)P4. Prove que as seguintes igualdades sao verdadeiras atraves da diferenciacao

numerica por diferencas finitas:

∂3φ

∂x3

∣∣∣∣i

=φi+2,j − 2φi+1,j + 2φi−1,j − φi−2,j

2x3+O

(∆x2

)∂2φ

∂x2

∣∣∣∣i

=2φi,j − 5φi−1,j + 4φi−2,j − φi−3,j

∆x2+O

(∆x2

)∂φ

∂x

∣∣∣∣i

=−φi+2,j + 8φi+1,j − 8φi−1,j + φi−2,j

12∆x+O

(∆x4

)∂2φ

∂x2

∣∣∣∣i

=−φi+2,j + 16φi+1,j − 30φi,j + 16φi−1,j − φi−2,j

12∆x2+O

(∆x4

)

2.10. PROBLEMAS 2-33

P5. Prove que as seguintes expressoes correspondem a aproximacoes de dife-rencas para derivadas parciais mistas:

∂2φ

∂x∂y

∣∣∣∣i,j

=1

∆x

(φi+1,j − φi+1,j−1

∆y− φi,j − φi,j−1

∆y

)+O (∆x∆y)

∂2φ

∂x∂y

∣∣∣∣i,j

=1

∆x

(φi,j+1 − φi,j

∆y− φi−1,j+1 − φi−1,j

∆y

)+O (∆x∆y)

∂2φ

∂x∂y

∣∣∣∣i,j

=1

∆x

(φi+1,j+1 − φi+1,j−1

2∆y− φi,j+1 − φi,j−1

2∆y

)+O

[∆x, (∆y)2

]∂2φ

∂x∂y

∣∣∣∣i,j

=1

∆x

(φi,j+1 − φi,j−1

2∆y− φi−1,j+1 − φi−1,j−1

2∆y

)+O

[∆x, (∆y)2

]P6. Obtenha um grafico semelhante ao da figura 2.9 mas correspondente a

segunda derivada de f = exp (ikx) em que a segunda derivada e −k2f .A aplicacao de uma formula de diferencas finitas para a segunda derivadaconfuz a −k′2f , em que k′2 se denomina pelo numero de onda modificadopaara a segunda derivada. Plote os valores de k′

2h2 em funcao de kh

no intervalo 0 ≤ kh ≤ π e os valores de k2h2 tambem em funcao de kh.Considere os operadores de finitas:

(a) Diferencas centrais: f ′′i = fi+1−2fi+2fi−1h2 .

(b) Formula de Pade para a segunda derivada.(c) Diferencas centrais de 4ª ordem de precisao.

P7. Considere as equacoes div (grad(φ)) = 0 e ∇2φ = 0 na forma unidi-mensional. O domınio discreto e formado por uma malha nao uniforme∆xi+1 = r∆xi. Pretende-se obter formulas de diferencas finitas centrais ede 2ª ordem (quando em malha uniforme), para as derivadas. Apresenteas formulas para malha nao uniforme e analise os erros de truncatura ob-tidos para ambas as formas da mesma equacao matematica,(Exame 2005,1ª chamada)

P8. Obtenha a equacao de diferencas finitas compacta de 4ª ou de 6ª ordempara a solucao de a∂φ∂x = bφ+ c no domınio Ω ∈ [0, 1] em que a, b e c saoconstantes. Supondo condicoes fronteira φ = φ1, φ = φ2 respectivamentepara os pontos x = 0 e x = 1, com o domınio com 7 pontos, apresente asmodificacoes para implementar as condicoes fronteira e a sua justificacao.(Exame 2005, 1ª chamada)

P9. Considere ∂∂xµ

∂U∂x e o seu desenvolvimento ∂

∂xµ∂U∂x = ∂2U

∂x2 + ∂µ∂x

∂U∂x em que

µ = µ(x).

(a) Justifique se ambas as expressoes, aproximadas por formulas de dife-rencas finitas centrais de 2ª ordem para as derivadas, tem o mesmoerro de truncatura numa malha uniforme ou nao uniforme, ∆xi+1 =r∆xi. (Exame 2005, 2ª chamada)


(b) Se aproximar

∂2U

∂x2= − 1

2∆x

[(∂U

∂x

)i+1

−(∂U

∂x

)i−1

]+

2∆x2

(Ui+1 − 2Ui + Ui−1) ,

qual seria a ordem de precisao, numa malha uniforme, da apro-ximacao se usar diferencas finitas centrais de 2ª ordem para a pri-meira derivada? Verifique para U = sinx.

(c) Considere a equacao ∂φ∂t + a∂φ∂x + b∂φ∂y = 0 em x ∈ [0, 1] e y ∈ [0, 1]

com condicoes periodicas em x = 0 e x = 1, e y = 0 e y = 1 emque a e b sao constantes. Considere um esquema temporal explıcitoe obtenha as equacoes de diferencas finitas compactas de 4ª ordemnum domınio com 4x3 pontos.

(d) Apresente a deducao da formula compacta de 4ª ordem central paraa segunda derivada numa malha uniforme.

P10. Considerando a funcao f(x) = sin xx3 num domınio discreto com incremento

espacial h, pretende-se:

(a) Derivar as formulas das diferencas finitas de 2ª e 4ª ordem.

(b) Estimar a ordem de grandeza do erro, de cada derivada, em relacaoa solucao analıtica.

(c) Obter o valor da segunda derivada em diferencas finitas compactas.

P11. Pretende-se calcular numericamente a rotacao do escoamento atmosfericodesde o vento geostroficos a 3000m ate a superfıcie da Terra. Este fenomenoe governado pelas equacoes de Ekman para a camada limite planetaria(PBL - Planetary Boundary Layer):

K∂2u

∂z2+ fv = 0

K∂2u

∂z2+ f (Ug − u) = 0

Para tal:

(a) Discretizar as equacoes com um esquema de diferencas finitas centraisde segunda ordem de precisao.

(b) Discretizar a equacao acoplada (4ª derivada) com um esquema dediferencas finitas centrais de segunda ordem de precisao.

(c) Comparar os resultados obtidos com a solucao exacta das equacoes.

2.10. PROBLEMAS 2-35

Para cada alınea, apresentar os resultados de perfil de u, de v e do anguloα (angulo entre a direccao do vento na camada limite e a direccao do ventogeostrofico) para os seguintes valores de viscosidades:

K = 0.001; 1.00; 10 [m2/s]

P12. Consideremos um escoamento inviscıdo bidimensional de um fluido invis-cido que entra numa camara, como representado na figura P2.12. Pretende-se calcular a funcao linha de corrente no interior da camara. Para umescoamento incompressıvel a equacao da continuidade vem expressa por:

∂u

∂x+∂v

∂y= 0

Afuncao linha de corrente pode ser definida de tal modo que:u = ∂ψ

∂y

v = ∂ψ∂x

Figura P2.12: Desenho da configuracao geometrica e condicoes fronteira doproblema 12

Substituindo pela funcao linha de corrente:

∂

∂x

(−∂ψ∂x

)− ∂

∂y

(∂ψ

∂y

)= 0⇔ ∂2ψ

∂x2+∂2ψ

∂y2= 0

Assim a solucao desta equacao, que e uma equacao elıptica, permite cal-cular a funcao linha de corrente no interior da camara, com as condicoesfronteira dadas na figura P2.12. Para a solucao do problema, divida odomınio num numero de pontos discretos atraves de uma malha computa-cional, em que, o espacamento ∆x e ∆y e igual a 0.2m. A condicao inicial


e de ψ = 0 no interior do domınio. Considere que a solucao numericaconvergiu quando o maximo erro no interior do domınio e menor do que∆ψ = 0.01, entre duas iteracoes consecutivas. Represente perfis de ψ aolongo de y, para x = 0.0, 1.0, 2.0, 3.0, 4.0 e 5.0. Represente valores daisofuncao de corrente.

P13. Considere o escoamento potencial na figura P2.13. O potencial de veloci-dades a entrada e dado por:

φ = Ax

onde A e uma constante igual a 1.

Usando o metodo de diferencas finitas, calcule os valores da funcao linhade corrente.

Figura P2.13: Figura do problema 13

Documents

Método das Diferenças Finitas