Baixe solução EDP eq onda e script e outras Provas em PDF para Engenharia Mecânica, somente na Docsity!
Solução numérica de equações
diferenciais parciais
13.1 – Advecção pura: a onda cinemática
Considere a equação
∂ u ∂ t
∂ u ∂ x
= 0, u ( x , 0) = g ( x ). (13.1)
A sua solução pode ser obtida pelo método das características, e é
u ( x , t ) = g ( x − c t ). (13.2)
Seja então o problema
∂ u ∂ t
∂ u ∂ x
u ( x , 0) = 2 x ( 1 − x ). (13.4)
A condição inicial, juntamente com u ( x , 1), u ( x , 2) e u ( x , 3) estão mostrados na figura 13.1. Observe que a solução da equação é uma simples onda cinemática.
0 2 4 6 8 10 0
1
2
3
4
u ( x ,^ t ) tempo
x
Figura 13.1: Condição inicial da equação 13.3.
Matemática Aplicada 234
Vamos adotar a notação
uni ≡ u ( xi , tn ), (13.5) xi = i ∆ x , (13.6) tn = n ∆ t , (13.7)
com
∆ x = L/Nx , (13.8) ∆ t = T /Nt (13.9)
onde L , T são os tamanhos de grade no espaço e no tempo, respectivamente, e Nx , Nt são os números de divisões no espaço e no tempo. Uma maneira simples de transformar as derivadas parciais em diferenças finitas na equação (13.3) é fazer
∂ u ∂ t (^) i , n
un i +^1 − uni ∆ t
∂ u ∂ x (^) i , n
uni + 1 − uni − 1 ∆ x
Substituindo na equação (13.3), obtemos o esquema de diferenças finitas explícito :
un i +^1 − uni ∆ t
= − c
Ç
uni + 1 − uni − 1 2 ∆ x
å ,
un i +^1 = uni −
c ∆ t 2 ∆ x
( uni + 1 − uni − 1 ), (13.12)
(com c = 2 no nosso caso). Esse é um esquema incondicionalmente instável , e vai fracassar. Vamos fazer uma primeira tentativa, já conformados com o fracasso an- tecipado. Ela vai servir para desenferrujar nossas habilidades de programação de métodos de diferenças finitas. O programa que implementa o esquema instável é o onda1d-ins.py, mostrado na listagem 13.1. Por motivos que ficarão mais claros na sequência, nós escolhemos ∆ x = 0,01, e ∆ t = 0,0005. O programa gera um arquivo de saída binário, que por sua vez é lindo pelo pró- ximo programa na sequência, surf1d-ins.py, mostrado na listagem 13.2. O único trabalho deste programa é selecionar algumas “linhas” da saída de onda1d-ins.py; no caso, nós o rodamos com o comando surf1d-ins.py 3 250
o que significa selecionar 3 saídas (além da condição inicial), de 250 em 250 interva- los de tempo ∆ t. Observe que para isto nós utilizamos uma lista (v), cujos elementos são arrays. O resultado dos primeiros 750 intervalos de tempo de simulação é mostrado na figura 13.2. Repare como a solução se torna rapidamente instável. Repare também como a solução numérica, em t = 750 ∆ t = 0,375, ainda está bastante distante dos tempos mostrados na solução analítica da figura 13.1 (que vão até t = 4). Clara- mente, o esquema explícito que nós programamos jamais nos levará a uma solução numérica satisfatória para tempos da ordem de t = 1!
Matemática Aplicada 236
Listagem 13.2: surf1d-ins.py — Seleciona alguns intervalos de tempo da solução numérica para plotagem 1 #!/usr/bin/python 2 # -- coding: iso-8859-1 -- 3 # ---------------------------------------------------------- 4 # surf1d-ins.py: imprime em +1 saídas de 5 # onda1d-ins a cada intervalos de tempo 6 # 7 # uso: ./surf1d-ins.py 8 # ---------------------------------------------------------- 9 from future import print_function 10 from future import division 11 from sys import argv 12 dx = 0. 13 dt = 0. 14 print('# dx = %9.4f' % dx) 15 print('# dy = %9.4f' % dt) 16 nx = int(10.0/dx) # número de pontos em x 17 print('# nx = %9d' % nx) 18 m = int(argv[1]) # m saídas 19 n = int(argv[2]) # a cada n intervalos de tempo 20 print('# m = %9d' % m) 21 print('# n = %9d' % n) 22 fin = open('onda1d-ins.dat', 23 'rb') # abre o arquivo com os dados 24 from numpy import fromfile 25 u = fromfile(fin,float,nx+1) # lê a condição inicial 26 v = [u] # inicializa a lista da "transposta" 27 for it in range(m): # para instantes: 28 for ir in range(n): # lê vezes, só guarda a última 29 u = fromfile(fin,float,nx+1) 30 v.append(u) # guarda a última 31 founam = 'surf1d-ins.dat' 32 print(founam) 33 fou = open(founam,'wt') # abre o arquivo de saída 34 for i in range(nx+1): 35 fou.write('%10.6f' % (i*dx)) # escreve o "x" 36 fou.write('%10.6f' % v[0][i]) # escreve a cond inicial 37 for k in range(1,m+1): 38 fou.write('%10.6f' % v[k][i])# escreve o k-ésimo 39 fou.write('\n') 40 fou.close()
237 13.1 – Advecção pura: a onda cinemática
0 2 4 6 8 10
u (
x ,
t )
tempo de simulação
x
Figura 13.2: Solução numérica produzida por onda1d-ins.py, para t = 250 ∆ t , 500 ∆ t e 750∆ t.
Por que o esquema utilizado em (13.12) fracassa? Uma forma de obter a resposta é fazer uma análise de estabilidade de von Neumann. A análise de estabilidade de von Neumann consiste primeiramente em observar que, em um computador real, (13.12) jamais será calculada com precisão infinita. O que o computador realmente calcula é um valor truncado u ˜ ni. Por enquanto, nós só vamos fazer essa distinção de notação, entre ˜ u e u , aqui, onde ela importa. O erro de truncamento é εni ≡ ˜ uni − uni. (13.13)
Note que (13.12) se aplica tanto para u quanto para ˜ u ; subtraindo as equações resul- tantes para ˜ un i +^1 e un i +^1 , obtém-se a mesma equação para a evolução de εni :
εn i +^1 = εni −
Co 2
( εni + 1 − εni − 1 ), (13.14)
onde
Co ≡
c ∆ t ∆ x
é o número de Courant. Isso só foi possível porque (13.12) é uma equação linear em u. Mesmo para equações não-lineares, entretanto, sempre será possível fazer pelo menos uma análise local de estabilidade. O próximo passo da análise de estabilidade de von Neumman é escrever uma série de Fourier para εni , na forma
tn = n ∆ t , xi = i ∆ x ,
εni =
∑^ N^ /^2
l = 1
ξl e atn^ ei kl^ xi^ , (13.16)
239 13.1 – Advecção pura: a onda cinemática
que é o critério de estabilidade de Courant-Friedrichs-Lewy.
A “mágica” de (13.21) é que ela introduz um pouco de difusão numérica ; de fato, podemos reescrevê-la na forma
un i +^1 − uni ∆ t
= − c
uni + 1 − uni − 1 2 ∆ x
uni + 1 − 2 uni + uni − 1 2 ∆ t
= − c
uni + 1 − uni − 1 2 ∆ x
Ç
∆ x^2 2 ∆ t
å uni + 1 − 2 uni + uni − 1 ∆ x^2
Não custa repetir: (13.24) é idêntica a (13.21). Porém, comparando-a com (13.12) (nosso esquema instável inicialmente empregado), nós vemos que ela também é equi- valente a esta última, com o termo adicional
∆ x^2 / 2 ∆ t
( uni + 1 − 2 uni + uni − 1 ) / ∆ x^2. O que este termo adicional significa? A resposta é uma derivada numérica de ordem 2. De fato, considere as expansões em série de Taylor
ui + 1 = ui +
du d x (^) i
∆ x +
d^2 u d x^2 i
ui − 1 = ui −
du d x (^) i
∆ x +
d^2 u d x^2 i
e some:
ui + 1 + ui − 1 = 2 ui +
d^2 u d x^2 i
∆ x^2 + O (∆ x^2 ),
d^2 u d x^2 i
ui + 1 − 2 ui + ui − 1 ∆ x^2
Portanto, a equação (13.24) — ou seja: o esquema de Lax (13.21) — pode ser inter- pretada também como uma solução aproximada da equação de advecção-difusão
∂ u ∂ t
∂ u ∂ c
= D
∂^2 u ∂ x^2
com
D =
Ç
∆ x^2 2 ∆ t
å .
Note que D tem dimensões de difusividade : J D K = L^2 T−^1. No entanto: não estamos
então resolvendo a equação errada? De certa forma, sim: estamos introduzindo um pouco de difusão na equação para amortecer as oscilações que aparecerão em decor- rência da amplificação dos erros de truncamento.
O quanto isto nos prejudica? Não muito, desde que o efeito da difusão seja muito menor que o da advecção que estamos tentando simular. Como a velocidade de advec- ção (“física”; “real”) que estamos simulando é c , precisamos comparar isto com (por exemplo) a magnitude das velocidades introduzidas pela difusão numérica; devemos
Matemática Aplicada 240
portanto verificar se
D ∂^
(^2) u ∂ x^2 c (^) ∂∂^ ux
D (^) ∆ ux 2 c (^) ∆ ux
D
∆ x
c ,
∆ x^2 2 ∆ t ∆ x
c , c ∆ t ∆ x
= Co
Em outras palavras, nós descobrimos que o critério para que o esquema seja acurado do ponto de vista físico é conflitante com o critério de estabilidade: enquanto que estabilidade demandava Co < 1, o critério de que a solução seja também fisicamente acurada demanda que Co 1 / 2. Na prática, isto significa que, para c = 2, ou o esquema é estável com muita difusão numérica, ou ele é instável. Isto praticamente elimina a possibilidade de qualquer uso sério de (13.21). Mesmo assim, vamos programá-lo! O programa onda1d-lax.py está mostrado na listagem 13.3. Ele usa os mesmos valores ∆ t = 0,0005 e ∆ x = 0,01, ou seja, Co = 0,10. O programa gera um arquivo de saída binário, que por sua vez é lido pelo pró- ximo programa na sequência, surf1d-lax.py, mostrado na listagem 13.4. O único trabalho deste programa é selecionar algumas “linhas” da saída de onda1d-lax.py; no caso, nós o rodamos com o comando surf1d-lax.py 3 500
o que significa selecionar 3 saídas (além da condição inicial), de 500 em 500 interva- los de tempo ∆ t. Com isto, nós conseguimos chegar até o instante 0,75 da simulação. O resultado dos primeiros 1500 intervalos de tempo de simulação é mostrado na figura 13.3. Observe que agora não há oscilações espúrias: o esquema é estável no tempo. No entanto, a solução está agora “amortecida” pela difusão numérica!
Upwind Um esquema que é conhecido na literatura como indicado por representar melhor o termo advectivo em (13.1) é o esquema de diferenças regressivas; neste esquema, chamado de esquema upwind — literalmente, “corrente acima” — na lite- ratura de língua inglesa, a discretização utilizada é
un i +^1 − uni ∆ t
= − c
uni − uni − 1 ∆ x
un i +^1 = uni − Co
î uni − uni − 1
ó
. (13.26)
Claramente, estamos utilizando um esquema de O (∆ x ) para a derivada espacial. Ele é um esquema menos acurado que os usados anteriormente, mas se ele ao mesmo tempo for condicionalmente estável e não introduzir difusão numérica, o resultado pode ser melhor para tratar a advecção. Antes de “colocarmos as mãos na massa”, sabemos que devemos analisar analiti- camente a estabilidade do esquema. Vamos a isto:
Matemática Aplicada 242
Listagem 13.4: surf1d-lax.py — Seleciona alguns intervalos de tempo da solução numérica para plotagem 1 #!/usr/bin/python 2 # -- coding: iso-8859-1 -- 3 # ---------------------------------------------------------- 4 # surf1d-lax.py: imprime em +1 saídas de 5 # onda1d-lax a cada intervalos de tempo 6 # 7 # uso: ./surf1d-lax.py 8 # ---------------------------------------------------------- 9 from future import print_function 10 from future import division 11 from sys import argv 12 dx = 0. 13 dt = 0. 14 print('# dx = %9.4f' % dx) 15 print('# dy = %9.4f' % dt) 16 nx = int(10.0/dx) # número de pontos em x 17 print('# nx = %9d' % nx) 18 m = int(argv[1]) # m saídas 19 n = int(argv[2]) # a cada n intervalos de tempo 20 print('# m = %9d' % m) 21 print('# n = %9d' % n) 22 fin = open('onda1d-lax.dat', 23 'rb') # abre o arquivo com os dados 24 from numpy import fromfile 25 u = fromfile(fin,float,nx+1) # lê a condição inicial 26 v = [u] # inicializa a lista da "transposta" 27 for it in range(m): # para instantes: 28 for ir in range(n): # lê vezes, só guarda a última 29 u = fromfile(fin,float,nx+1) 30 v.append(u) # guarda a última 31 founam = 'surf1d-lax.dat' 32 print(founam) 33 fou = open(founam,'wt') # abre o arquivo de saída 34 for i in range(nx+1): 35 fou.write('%10.6f' % (i*dx)) # escreve o "x" 36 fou.write('%10.6f' % v[0][i]) # escreve a cond inicial 37 for k in range(1,m+1): 38 fou.write('%10.6f' % v[k][i])# escreve o k-ésimo 39 fou.write('\n') 40 fou.close()
243 13.1 – Advecção pura: a onda cinemática
0 2 4 6 8 10
u (
x ,
t )
tempo de simulação
x
Figura 13.3: Solução numérica produzida por onda1d-lax.py, para t = 500 ∆ t , 1000 ∆ t e 1500∆ t.
ξl e a ( tn +∆ t )ei kl^ i ∆ x^ = ξl e atn^ ei kl^ i ∆ x^ − Co
î ξl e atn^ ei kl^ i ∆ x^ − ξl e atn^ ei kl^ ( i −^1 )∆ x^
ó
e a ∆ t^ ei kl^ i ∆ x^ = ei kl^ i ∆ x^ − Co
î ei kl^ i ∆ x^ − ei kl^ ( i −^1 )∆ x^
ó
e a ∆ t^ = 1 − Co
î 1 − e−i kl^ ∆ x^
ó
e a ∆ t^ = 1 − Co + Co cos( kl ∆ x ) − iCo sen( kl ∆ x ). (13.27) Desejamos que o módulo do fator de amplificação e a ∆ t^ seja menor que 1. O mó- dulo (ao quadrado) é
|e a ∆ t^ |^2 =
1 − Co + Co cos( kl ∆ x )
Co sen( kl ∆ x )
Para aliviar a notação, façamos
Ck ≡ cos( kl ∆ x ), Sk ≡ sen( kl ∆ x ).
Então,
|e a ∆ t^ |^2 = (Co Sk )^2 + (Co Ck − Co + 1 )^2 = Co^2 S k^2 + (Co^2 C k^2 + Co^2 + 1 ) + 2 (−Co^2 Ck + Co Ck − Co) = Co^2 ( S^2 k + C^2 k + 1 − 2 Ck ) + 2Co( Ck − 1 ) + 1 = 2Co^2 ( 1 − Ck ) + 2Co( Ck − 1 ) + 1.
A condição para que o esquema de diferenças finitas seja estável é, então,
2Co^2 ( 1 − Ck ) + 2Co( Ck − 1 ) + 1 ≤ 1, 2Co
Co( 1 − Ck ) + ( Ck − 1 )
1 − cos( kl ∆ x )
[Co − 1 ] ≤ 0, Co ≤ 1
245 13.2 – Difusão pura
Esta solução será vista no capítulo 17:
u ( x , t ) =
∑^ ∞
n = 1
An e−^
n^2 π^2 α^2 L^2 t^ sen nπx L
An =
L
∫ L
0
f ( x ) sen
nπx L
d x. (13.32)
Em particular, se
D = 2, L = 1, f ( x ) = 2 x ( 1 − x ),
An = 2
0
2 x ( 1 − x ) sen( nπx ) d x =
π^3 n^3
[ 1 − (− 1 ) n ].
Todos os An ’s pares se anulam. Fique então apenas com os ímpares:
A 2 n + 1 =
π^3 ( 2 n + 1 )^3
u ( x , t ) =
∑^ ∞
n = 0
π^3 ( 2 n + 1 )^3
e−(^2 (^2 n +^1 )
(^2) π (^2) ) t sen (( 2 n + 1 ) πx ) (13.33)
O programa difusao1d-ana.py, mostrado na listagem 13.5, implementa a solu- ção analítica para ∆ t = 0,0005 e ∆ x = 0,001. Da mesma maneira que os programas surf1d*.py, o programa divisao1d-ana.py, mostrado na listagem 13.6, seleciona alguns instantes de tempo da solução analítica para visualização: divisao1d-ana.py 3 100
A figura 13.5 mostra o resultado da solução numérica para t = 0, t = 0,05, t = 0,10 e t = 0,15. Este é praticamente o “fim” do processo difusivo, com a solução analítica tendendo rapidamente para zero.
Esquema explícito Talvez o esquema explícito mais óbvio para discretizar (13.28) seja un i +^1 − uni ∆ t
= D
uni + 1 − 2 uni + uni − 1 ∆ x^2
A derivada parcial em relação ao tempo é de O (∆ t ), enquanto que a derivada segunda parcial em relação ao espaço é, como vimos em (13.25), de O (∆ x^2 ). Mas não nos preocupemos muito, ainda, com a acurácia do esquema numérico. Nossa primeira preocupação, como você já sabe, é outra: o esquema (13.34) é estável? Explicitamos un i +^1 em (13.34):
un i +^1 = uni + Fo
î uni + 1 − 2 uni + uni − 1
ó , (13.35)
onde
Fo =
D ∆ t ∆ x^2
Matemática Aplicada 246
Listagem 13.5: difusao1d-ana.py — Solução analítica da equação da difusão 1 #!/usr/bin/python 2 # -- coding: iso-8859-1 -- 3 # ---------------------------------------------------------- 4 # difusao1d-ana: solução analítica de 5 # 6 # du/dt = D du^2/dx^ 7 # 8 # u(x,0) = 2x(1-x) 9 # u(0,t) = 0 10 # u(1,t) = 0 11 # 12 # uso: ./difusao1d-ana.py 13 # ---------------------------------------------------------- 14 from future import print_function 15 from future import division 16 fou = open('difusao1d-ana.dat','wb') 17 dx = 0. 18 dt = 0. 19 print('# dx = %9.4f' % dx) 20 print('# dy = %9.4f' % dt) 21 nx = int(1.0/dx) # número de pontos em x 22 nt = int(1.0/dt) # número de pontos em t 23 print('# nx = %9d' % nx) 24 print('# nt = %9d' % nt) 25 from math import pi, sin, exp 26 epsilon = 1.0e-6 # precisão da solução analítica 27 dpiq = 2pipi # 2pi^ 28 dzpic = 16/(pipipi) # 16/pi^ 29 def ana(x,t): 30 s = 0. 31 ds = epsilon 32 n = 0 33 while abs(ds) >= epsilon: 34 dnm1 = 2n + 1 # (2n+1) 35 dnm1q = dnm1dnm1 # (2n+1)^ 36 dnm1c = dnm1qdnm1 # (2n+1)^ 37 ds = exp(-dnm1qdpiqt) 38 ds = sin(dnm1pix) 39 ds /= dnm1c 40 s += ds 41 n += 1 42 return sdzpic 43 from numpy import zeros 44 u = zeros(nx+1,float) # um array para conter a solução 45 for n in range(nt+1): # loop no tempo 46 t = ndt 47 print(t) 48 for i in range(nx+1): # loop no espaço 49 xi = i*dx 50 u[i] = ana(xi,t) 51 u.tofile(fou) # imprime uma linha com os novos dados 52 fou.close()
Matemática Aplicada 248
0 0.2 0.4 0.6 0.8 1
u (
x ,
t )
x
t = 0,
t = 0,
t = 0, t = 0,
Figura 13.5: Solução analítica da equação de difusão para t = 0, t = 0,05, t = 0,10 e t = 0,15.
é o número de Fourier de grade (El-Kadi e Ling, 1993). A análise de estabiliade de von Neumann agora produz
ξl e a ( tn +∆ t )ei kl^ i ∆ x^ = ξl e atn^ ei kl^ i ∆ x^ + Fo
î ξl e atn^ ei kl^ ( i +^1 )∆ x^ − 2 ξl e atn^ ei kl^ i ∆ x^ + ξl e atn^ ei kl^ ( i −^1 )∆ x^
ó , e a ∆ t^ = 1 + Fo
î e+i kl^ ∆ x^ − 2 + e−i kl^ ∆ x^
ó
= 1 + 2Fo
cos( kl ∆ x ) − 1
= 1 − 4Fo sen^2
kl ∆ x 2
A análise de estabilidade requer que | ea ∆ t^ | < 1:
| ea ∆ t^ |^2 = 1 − 8Fo sen^2
kl ∆ x 2
kl ∆ x 2
ou
−8Fo sen^2
kl ∆ x 2
kl ∆ x 2
8Fo sen^2
kl ∆ x 2
− 1 + 2Fo sen^2
kl ∆ x 2
Fo <
Podemos agora calcular o número de Fourier que utilizamos para plotar a solução analítica (verifique nas listagens 13.5 e 13.6):
Fo =
2 × 0,
(0,001)^2
249 13.2 – Difusão pura
Listagem 13.7: difusao1d-exp.py — Solução numérica da equação da difusão: mé- todo explícito. 1 #!/usr/bin/python 2 # -- coding: iso-8859-1 -- 3 # ---------------------------------------------------------- 4 # difusao1d-exp resolve uma equação de difusão com um método 5 # explícito 6 # 7 # uso: ./difusao1d-exp.py 8 # ---------------------------------------------------------- 9 from future import print_function 10 from future import division 11 fou = open('difusao1d-exp.dat','wb') 12 dx = 0. 13 dt = 0. 14 print('# dx = %9.4f' % dx) 15 print('# dy = %9.4f' % dt) 16 from numpy import zeros 17 nx = int(round(1.0/dx,0)) # número de pontos em x 18 nt = int(round(1.0/dt,0)) # número de pontos em t 19 print('# nx = %9d' % nx) 20 print('# nt = %9d' % nt) 21 u = zeros((2,nx+1),float) # apenas 2 posições no tempo 22 # são necessárias! 23 def CI(x): # define a condição inicial 24 if 0 <= x <= 1.0: 25 return 2.0x(1.0-x) 26 else: 27 return 0. 28 for i in range(nx+1): # monta a condição inicial 29 xi = idx 30 u[0,i] = CI(xi) 31 u[0].tofile(fou) # imprime a condição inicial 32 old = False 33 new = True 34 D = 2.0 # celeridade da onda 35 Fon = Ddt/((dx)2) # número de Fourier 36 print("Fo = %10.6f" % Fon) 37 for n in range(nt): # loop no tempo 38 print(n) 39 for i in range(1,nx): # loop no espaço 40 u[new,i] = u[old,i] + Fon(u[old,i+1] - 2u[old,i] + u[old,i-1]) 41 u[new,0] = 0.0 # condição de contorno, x = 0 42 u[new,nx] = 0.0 # condição de contorno, x = 1 43 u[new].tofile(fou) # imprime uma linha com os novos dados 44 (old,new) = (new,old) # troca os índices 45 fou.close()
Utilizar os valores ∆ x = 0,0005 e ∆ x = 0,001 levaria a um esquema instável. Preci- samos diminuir ∆ t e/ou aumentar ∆ x. Com ∆ t = 0,00001 e ∆ x = 0,01,
Fo =
2 × 0,
(0,01)^2
= 0,2 < 0,5 (OK).
Repare que Fo < 1 / 2 é um critério de estabilidade muito mais exigente do que Co < 1 / 2 (para D = 2). Nós esperamos que nosso esquema explícito agora rode muito lentamente. Mas vamos implementá-lo. O programa que implementa o esquema é o difusao1d-exp.py, mostrado na listagem 13.7. O programa divisao1d-exp.py, mostrado na listagem 13.8, seleciona alguns instantes de tempo da solução analítica para visualização: divisao1d-exp 3 5000
251 13.2 – Difusão pura
0 0.2 0.4 0.6 0.8 1
u (
x ,
t )
x
t = 0,
t = 0,
t = 0, t = 0,
Figura 13.6: Solução numérica com o método explícito (13.35) (círculos) versus a solução analítica (linha cheia) da equação de difusão para t = 0, t = 0,05, t = 0,10 e t = 0,15. Apenas 1 a cada 5 pontos da grade numérica são mostrados, para facilitar a comparação com a solução analítica.
O resultado da solução numérica com o método explícito está mostrado na figura 13.6: ele é impressionantemente bom, embora seja computacionalmente muito caro. A escolha judiciosa de ∆ t e ∆ x para obeder ao critério (13.38) foi fundamental para a obtenção de um bom resultado “de primeira”, sem a necessidade dolorosa de fi- car tentando diversas combinações até que o esquema se estabilize e produza bons resultados.
Esquemas implícitos Embora o esquema explícito que nós utilizamos acima seja acurado, ele é lento — se você programou e rodou difusao1d-exp.py, deve ter no- tado alguma demora para o programa rodar. Embora nossos computadores estejam ficando a cada dia mais rápidos, isso não é desculpa para utilizar mal nossos recursos computacionais (é claro que, ao utilizarmos uma linguagem interpretada — Python — para programar, nós já estamos utilizando muito mal nossos recursos; no entanto, nosso argumento é didático: com uma linguagem mais simples, podemos aprender mais rápido e errar menos. Além disso, todos os ganhos relativos que obtivermos se manterão em qualquer outra linguagem) Vamos portanto fazer uma mudança fundamental nos nossos esquemas de dife- renças finitas: vamos calcular a derivada espacial no instante n + 1:
un i +^1 − uni ∆ t
= D
un i ++ 11 − 2 un i +^1 + un i −+ 11 ∆ x^2
un i +^1 − uni = Fo( un i ++ 11 − 2 un i +^1 + un i −+ 11 ), −Fo un i −+ 11 + ( 1 + 2Fo) un i +^1 − Fo un i ++ 11 = uni. (13.39)
Reveja a discretização (13.5)–(13.9): para i = 1,... , Nx − 1, (13.39) acopla 3 valores das incógnitas un +^1 no instante n + 1. Quando i = 0, e quando i = Nx , não
Matemática Aplicada 252
podemos utilizar (13.39), porque não existem os índices i = −1, e i = Nx +1. Quando i = 1 e i = Nx − 1, (13.39) precisa ser modificada, para a introdução das condições de contorno : como u 0 n = 0 e unNx = 0 para qualquer n , teremos
( 1 + 2Fo) un 1 + 1 − Fo un 2 + 1 = un 1 , (13.40) −Fo un N + x −^12 + ( 1 + 2Fo) un N + x −^11 = unNx − 1. (13.41)
Em resumo, nossas incógnitas são un 1 +^1 , un 2 + 1 ,... un N + x −^11 ( Nx −1 incógnitas), e seu cálculo envolve a solução do sistema de equações
1 + 2Fo −Fo 0... 0 0 −Fo 1 + 2Fo Fo 0... 0 .. .
0... 0 −Fo 1 + 2Fo −Fo 0 0... 0 −Fo 1 + 2Fo
un 1 +^1 un 2 +^1 .. . un N + x −^12 un N + x −^11
un 1 un 2 .. . unNx − 2 unNx − 1
A análise de estabilidade de von Neumann procede agora da maneira usual:
εn i +^1 = εni + Fo( εn i ++ 11 − 2 εn i +^1 + εn i −+ 11 ) ξl ea ( tn +∆ t )^ e i kl^ i ∆ x^ = ξl eatn^ e i kl^ i ∆ x
Ä
ξl ea ( tn +∆ t )^ e i kl^ ( i +^1 )∆ x^ − 2 ξl ea ( tn + ∆ t ) e i kl^ i ∆ x
- ξl ea ( tn +∆ t )^ e i kl^ ( i −^1 )∆ x^
ä , ea ∆ t^ = 1 + ea ∆ t^ Fo
Ä
e i kl^ ∆ x^ − 2 + e −i kl^ ∆ x^
ä , ea ∆ t^ = 1 + ea ∆ t^ 2Fo
cos( kl ∆ x ) − 1
ea ∆ t^ = 1 − ea ∆ t^ 4Fo sin^2
kl ∆ x 2
ea ∆ t
1 + 4Fo sin^2
kl ∆ x 2
| ea ∆ t^ | =
1 + 4Fo sin^2
Ä (^) k l ∆ x 2
ä ≤ 1 sempre. (13.43)
Portanto, o esquema implícito (13.39) é incondicionalmente estável, e temos con- fiança de que o programa correspondente não se instabilizará. Existem várias coisas atraentes para um programador em (13.42). Em primeiro lugar, a matriz do sistema é uma matriz banda tridiagonal; sistemas lineares com este tipo de matriz são particularmente simples de resolver, e estão disponíveis na lite- ratura (por exemplo: Press et al., 1992, seção 2.4, subrotina tridag). Em segundo lugar, a matriz do sistema é constante : ela só precisa ser montada uma vez no pro- grama, o que torna a solução numérica potencialmente muito rápida. Nós vamos começar, então, construindo um pequeno módulo, convenientemente denominado alglin.py, que exporta a função tridag, que resolve um sistema tridi- agonal, mostrado na listagem 13.9. Em seguida, o programa difusao1d-imp.py resolve o problema com o método implícito. Ele está mostrado na listagem 13.10. A principal novidade está nas linhas 42–46, e depois novamente na linha 56. Em Python e Numpy, é possível especi- ficar sub-listas, e sub-arrays, com um dispositivo denominado slicing , que torna a
