cap07

Prévia do material em texto

FIS045
Métodos Computacionais em física
Cap 7: Equações diferenciais parciais
Prof. Gustavo Guerrero (sala 4120)
e-mail para enviar as listas de exercicios:
�s045.ufmg@gmail.com
Universidade Federal de Minas Gerais
2017
Física Computacional
Introdução
Em ciências naturais comunmente encontramos problemas com muitas
variáveis limitadas por condições de iniciais e de contorno. A grande
maioria desses problemas pode ser modelado por equações diferenciais
parciais (EDP). Um caso familiar de EDP é a equação de onda que em 1
dimensão tem a forma:
∂2u
∂x2
= A
∂2u
∂t2
(1)
Onde A é uma constante. A solução u depende de variáveis espaciais e
temporais, i.e., u = u(x , t)
Física Computacional
Eq. de onda em duas dimensões, (equação hiperbólica)
Em duas dimensões temos u = u(x , y , t) (u por simplicidade). Essa
equação aparece em e.g., ondas em uma mola, ondas de pressão, ondas
de super�cie, ondas eletromagnéticas, ondas de som.
∂2u
∂x2
+
∂2u
∂y2
= A
∂2u
∂t2
,
no caso de ondas eletromagnéticas, A = c , onde c é a velocidade da luz.
Física Computacional
Equação de difusão (equação parabólica)
A equação de difuão em uma dimensão é:
∂2u
∂x2
= A
∂u
∂t
(2)
Nesse caso A é chamada de constante de difusão. Pode ser usada para
um grande número de processos difusivos, desde a difusão de moléculas,
calor ou campo magnético.
As equações hiperbólica e parabólica são chamados problemas de valor
inicial (ou problemas de Cauchy). O valor de u é conhecido em um
tempo inicial e a equação descreve como u(x , t) se propaga no espaço e
no tempo.
Física Computacional
Equação de Laplace (equação elíptica)
A equação de Laplace aparece em electrostática, é similar a equação de
onda, exeto que A = 0:
∂2u
∂x2
+
∂2u
∂y2
= 0, (3)
Ou, em caso de existur uma carga representada por uma densidade de
carga ρ(x), temos a equação de Poisson:
∂2u
∂x2
+
∂2u
∂y2
= −4πρ(x). (4)
Na equação elíptica queremos encontrar uma solução estática u(x , y) que
satisfaz a equação numa região (x , y) de interesse. Esse tipo de
problemas são chamado problemas de contorno, com:
I Condições de Dirichlet: especi�ca valores nas fronteiras como
função do tempo.
I Condições de Neumann: espe�cica os valores dos gradientes
normais.
I Condições abertas: permitem o �uxo das quantidades através das
fronteiras.
Física Computacional
Equação de Schroedinger
Em mecanica quantica temos a equação de Schroedinger que em duas
dimensões tem a forma:
−∂
2u
∂x2
− ∂
2u
∂y2
+ f (x , y)u = ı
∂u
∂t
.
Física Computacional
Equações de Maxwell
Um importante sistema de equações diferenciais parciais lineares são as
equações de Maxwell, que no caso de carga nula são:
∂E
∂t
= curlB,
−curlE = B,
divE = divB = 0.
Física Computacional
Equações de Euler
Um conhecido sistema de EDP não lineares são as equações de Eular, as
quais para um caso de um �uido incompresivel é não viscoso são:
∂u
∂t
+ (u · ∇)u = −Dp; divu = 0,
onde p é a pressão, e
∇ = ∂
∂x
êx +
∂
∂y
êy ,
é o operador gradiente. No caso de um �uido viscoso temos a equação de
Navier-Stokes:
∂u
∂t
+ u∇u− ν∆u = −Dp; divu = 0.
Física Computacional
Equações de �uxo-conservativo, problemas de valor inicial
Um grande número de problemas de valor inicial são as equações de �uxo
conservativo:
∂u
∂t
= −∂F(u)
∂x
= 0,
onde o vetor F é chamado de �uxo conservado. Como exemplo, considere
uma equação hiperbólica padrão, i.e., a equação de onda 1D com
velocidade de propagação v :
∂2u
∂t2
= v2
∂2u
∂x2
,
que pode ser escrita como o sistema de equações:
∂r
∂t
= v
∂s
∂x
,
∂s
∂t
= v
∂r
∂x
(eq. de Maxwell em 1D para propagação de uma onda EM???? ) com
r ≡ v ∂u
∂x
, s ≡ ∂u
∂t
Física Computacional
Neste caso r e s são as duas componentes de u e o �uxo é dado por
F =
[
0 −v
−v 0
]
· u
Logo, uma possível solução numérica para uma equação hiperbólica
consiste em resolver a equação de �uxo conservativo.
∂u
∂t
= −v ∂u
∂x
A partir da solução analítica, sabemos a solução de essa equação é uma
onda que se propaga na direção positiva do eixo x : u = f (x − vt).
Vamos estudar inicialmente como várias estrategias para resolver esse
problema numericamente.
Física Computacional
Diferenças �nitas
Usamos diferenças �nitas para discretizar a equação acima com pontos
igualmente espaçados nos eixos t e x :
xi = x0 + j∆x , j = 0, 1, ..., J
tn = t0 + n∆t, n = 0, 1, ...,N
Temos varias alternativas para representar a derivada temporal, a mais
obvia a diferenciação de eular para frente.
∂u
∂t
=
un+1j − unj
∆t
+O(∆t)
Calculamos um valor desconhecido, n + 1 em termos dos valores
desconhecidos, n. Para a derivada espacial podemos usar a rep. de
segunda ordem forward time centered space (FTCS):
∂u
∂x
=
unj+1 − unj−1
2∆x
+O(∆x2)
Física Computacional
O que resulta em
un+1j − unj
∆t
= −v
unj+1 − unj−1
2∆x
(5)
É um algoritmo explícito, simple, ocupa pouca memoria e é rapido.
Infelizmente não funciona.
Analise de estabilidade de von Neumann: Imagine os coe�cientes das
equações variar tão suavemente que podem ser considerados constantes.
Nesse caso, soluções independentes da equação tem a forma:
unj = ζ
ne ikj∆x (6)
Com k : número de onda espacial e ζ = ζ(k) é um número complexo. A
dependência temporal dessa solução para um automodo simples, é uma
suscesão de numeros complexos que dependem de k . Assim, se
|ζ(k)| > 1 implica modos que crescem exponencialmente (uma
instabilidade numérica).
Substituindo (6) em (5) temos e dividindo por ζn temos (ver
demonstração no quadro):
Física Computacional
ζ(k) = 1− i v∆t
∆x
sin(k∆x)
com modulo > 1 para todo k . Assim o esquema é incondicionalmente
instável.
Método de Lax: O método FTCS pode "ser curado"fazendo uma
mudança simples introduzida por Lax. Remplazamos unj na derivada
temporal por sua média:
unj →
1
2
(unj+1 + u
n
j−1),
o que resulta em:
un+1j =
1
2
(unj+1 + u
n
j−1)− v∆t
unj+1 − unj−1
2∆x
(7)
Substituindo (6) em (7) encontramos:
ζ(k) = cos(k∆x)− i v∆t
∆x
sin(k∆x),
Física Computacional
que para representar um sistema estável requer:
|v |∆t
∆x
≤ 1.
Esse o famoso criterio de estabilidade de Courant-Friedrichs-Lewy, mais
conhecido como a condição de Courant.
Mas como uma mudança tão simples estabiliza o código? (resolver EDP
é uma arte mais do quem uma ciência.) Podemos desmiti�car essa arte
reescrevendo a equação (7) na forma:
un+1j − unj
∆t
= −v
unj+1 − unj−1
2∆x
+
1
2
unj+1 − 2unj + unj−1
∆t
, (8)
que é a representação exata em FTCS da equação:
∂u
∂t
= −v ∂u
∂x
+
(∆x)2
2∆t
∇2u
i.e., adicionamos um termo difusivo a equação de �uxo conservativo.
O esquema de Lax têm difusão numérica ou viscosidade numérica.
Física Computacional
Segunda ordem de precisao em tempo
Metodo de Lax-Wendro�: é um método de segunda ordem em tempo
que elimina grande dissipação numérica. Devemos de�nir passos
intermediarios, uj+1/2 em passos de tempo intermediarios, tn+1/2 e no
ponto intermediario da grade, xj+1/2. Esses pontos são calculados pelo
método de Lax:
u
n+1/2
j+1/2 =
1
2
(unj+1 + u
n
j )−
∆t
2∆x
(F nj+1 − F nj ) (9)
Usando essas variáveis podem se calcular, F n+1/2j±1/2 , depois o valor de u
n+1
j
pode ser calculado usando diferenças centradas:
un+1j = u
n
j −
∆t
∆x
(F
n+1/2
j+1/2 − F
n+1/2
j−1/2 ) (10)
Substituindo (9) em (10), com F = vu, e α = v∆t/∆x , temos:
un+1j = u
n
j −
α
2
[
1
2
(unj+1 + u
n
j−1)−
α
2
(unj+1 − 2unj − unj−1)
]
Física Computacional
Figura: Representação do método de 2 passos de Lax-Wendro�. Os
pontos com a × são calculados com o método de Lax. Esses pontos,
mais os 3 pontos originais, permitem calcular o ponto un+1j . O método
tem segunda ordem de precisão em tempo. O salto no ponto un+1/2j é
conhecido como leapfrog.
Infelizmente o método de Lax-Wendro� pode levar a outra clase de erros
conhecida como erros de fase, gerando oscilações em comprimentos de
onda pequenos.
Física Computacional
Podemos aprofundar mais no problema da advecção: Nos métodos
de Euler e Lax, assumimos que dentro de cada celda o �uxo, vu, se
mantémconstante. Mas isso é só uma suposição. Podemos supor
também que o estado dentro de cada celula é uma função linear da
posição, ou um picewise linear subgrid scale model.
Dentro de cada celda, o estado no começo de um passo temporal é:
u(x , t = tn) = u
n
j + σ(x − xj) , para xj−1/2 < x < xj+1/2
Onde σnj é uma inclinação (slope), e xj = (xj−1/2 − xj+1/2)/2. Na
interface, o �uxo varia no tempo entre tn e tn+1:
Fj−1/2(t) = vu(x = xj−1/2, t)
= vunj−1 + vσ
n
j−1(xj−1/2 − xj−1 − v(t − tn))
= vunj−1 + vσ
n
j−1
(
∆x
2
− v(t − tn)
)
Física Computacional
De forma similar para Fj+1/2(t). O �uxo médio durante o passo temporal
∆t é:
< Fj−1/2(t) >
tn+1
tn =
1
∆t
∫ tn+1
tn
vu(x = xj−1/2, t)
= vunj−1
1
2
vσnj−1(∆x − v∆t)
Com um resultado similar para < Fj+1/2(t) >
tn+1
tn . A diferença é:
< Fj+1/2(t) >
tn+1
tn − < Fj−1/2(t) >
tn+1
tn = v(u
n
j − unj−1)
+
1
2
v(σnj − σni−1)(∆x − v∆t)
Física Computacional
Finalmente usando a equação (10) acima:
un+1j = u
n
j −
∆t
∆x
(F
n+1/2
j+1/2 − F
n+1/2
j−1/2 )
temos:
un+1j = u
n
j −
v∆t
∆x
(unj − unj−1)−
v∆t
∆x
1
2
(σnj − σni−1)(∆x − v∆t)
A questão é agora, como escolher a inclinação apropriada?
centered slope : σj =
unj+1−u
n
j−1
2∆x Fromm
′s method
upwind slope : σj =
unj −u
n
j−1
2∆x Beam−Warming method
downwind slope : σj =
unj+1−u
n
j
2∆x Lax−Wendro� method
Física Computacional
Limitadores de inclinação
Métodos de alta ordem como Beam-Warming ou Lax-Wendro� temdem
a produzir oscillações como consequência de erros de fase ou super
estimação da inclinação.
Para corrigir esse problema se utilizam os limitadores de inclinação. A
inclinação, σj é modi�cada unicamente quando necessário.
Para medir oscilações em uma solução numérica de�nimos a variação
total (TV) do estado u:
TV (u)−
N∑
i=1
|qi − qi−1|,
Física Computacional
Se a quantidade u varia monotonicamente, TV se mantém constante, se
u desenvolve máximos ou minimos, então TV aumenta. Um esquema é
chamado de redutor da variação total (TVD), se:
TV (un+1) ≤ TV (un)
Alguns limitadores de inclinação são:
I inclinação minmod
σn = minmod
(
unj − unj−1
∆x
,
unj+1 − unj
∆x
)
onde
minmod(a, b) =
{ a sim |a| < |b| e ab > 0
b sim |a| > |b| e ab > 0
0 sim ab < 0
}
Física Computacional
I inclinação superbee
σn = maxmod
(
σ
(1)
j , σ
(2)
j
)
Com:
σ
(1)
j = minmod
(
unj+1 − unj
∆x
, 2
unj − unj−1
∆x
)
σ
(2)
j = minmod
(
2
unj+1 − unj
∆x
,
unj − unj−1
∆x
)
Física Computacional
Equação de difusão
Procuramos a seguir a solução numérica da equação de difusão:
∂u
∂t
= −D ∂
2u
∂x2
(11)
Onde D é o coe�ciente de difusividade. Podemos resolver a equação
usando diferenças �nitas de forma explicita:
un+1j − unj
∆t
= D
unj+1 − 2unj + unj−1
∆x2
A analise de estabilidade de von-Neumann indica que esse método é
estavel sob a condição:
D∆t
∆x2
<
1
2
A qual limita os modelos a passos temporais muito curtos.
Física Computacional
Outra possibilidade é expresar (11) em forma implícita:
un+1j − unj
∆t
= D
un+1j+1 − 2u
n+1
j + u
n+1
j−1
∆x2
assim:
un+1j − ηu
n+1
j+1 + 2ηu
n+1
j − ηu
n+1
j−1 = u
n
j
−ηun+1j+1 + (1 + 2η)u
n+1
j − ηu
n+1
j−1 = u
n
j
com η = D∆t/∆x2. Pode se ver o sistema tridiagonal:
ÂUj+1 = Uj
onde
A =

1 + 2η −eta 0 0 ...
−η 1 + 2η −eta 0 ...
... ... ... ... ...
... ... ... ... −eta
0 0 ... −eta 1 + 2η

Física Computacional
O método implicito é incondicionalmente estável (ou seja, funciona para
todo ∆t e ∆x). Porem, tanto os métodos explicito quanto implícito são
de primeira ordem no tempo.
O método de Crank-Nicolson combina estabilidade com a precisão da
segunda ordem no tempo. Ele simplesmente combina os dois esquemas,
explícito e implícito. A equação discreta �caria:
un+1j − unj
∆t
=
D
2
[
(un+1j+1 − 2u
n+1
j + u
n+1
j−1 ) + (u
n
j+1 − 2unj + unj−1)
∆x2
]
Que pode ser resolvido na forma tridiagonal:
ÂUj+1 = Rj
Neste caso  é igual ao caso implicito, mas
R =
η
2
unj+1 + (1− η)unj +
η
2
unj−1
Física Computacional
Figura: Comparação entre os diferentes esquemas numéricos para resolver
a equação de difusão.
Física Computacional
Equações diferenciais em múltiples dimensões
Os métodos discutidos anteriormente podem ser geralizados a 2 ou 3
dimensões. No entanto, o tempo de computo aumenta signi�cativamente
quando os problemas são tri-dimensionais. Veremos alguns exemplos
sobre como geralizar as equações de �uxo conservativo e de difussão.
Método de Lax em duas dimensões:
∂u
∂t
= −∇F = −
(
∂Fx
∂x
+
∂Fy
∂y
)
Considerando a grade espacial:
xj = x0 + j∆x yl = x0 + l∆y
O esquema de Lax resulta em:
un+1j,l =
1
4
(
unj+1,l + u
n
j−1,l + u
n
j,l+1 + u
n
j,l−1
)
−
(
F nj+1,l − F nj−1,l + F nj,l+1 − F nj,l−1
)
Física Computacional
A analise de estabilidade leva à seguinte condição de Courant :
∆t ≤ ∆√
2(v2x + v2y )1/2
No caso N-dimensional, se |v| é a velocidade máxima de propagação no
problema, então:
∆t ≤ ∆√
N|v |
(Demonstração em Numerical Recipes, Sec. 19.3)
Física Computacional
Equação de difusão em duas dimensões:
∂u
∂t
= −D
(
∂2u
∂x2
+
∂2u
∂y2
)
O problema pode ser abordado usando diferenças �nitas centradas
(FTCS). No entanto, ja vimos que esse esquema tem limitações de
estabilidade.
Por outro lado, o método de Crank-Nicolson resulta num sistema de
equações acopladas que por sua vez não formam um sistema tri-diagonal.
Metodo alternando a direção implícita (ADI):
Divide cada passo temporal em dois passos de tamanho ∆t/2. Em cada
sub-passo, uma direção é tratada implicitamente.
u
n+1/2
j,l = u
n
j,l +
1
2
α
(
δ2xu
n+1/2
j,l + δ
2
yu
n
j,l
)
un+1j,l = u
n+1
j,l +
1
2
α
(
δ2xu
n+1/2
j,l + δ
2
yu
n+1
j,l
)
onde
δ2xu
n
j,l = u
n
j+1,l − 2unj,l + unj−1,l
e similar para δ2yu
n
j,l .
Física Computacional
What's trending in twitter this week: (source ADS-NASA)
Física Computacional
Equações de Laplace e Poisson
A equação de Laplace é:
∇2u = uxx + uyy = 0
para ser resolvida precisa de condições de contorno, tal que
u(x , y) = g(x , y) na borda do dominio, δΩ. Não há dependencia
temporal. O dominio pode ter qualquer forma e tamanho, por
simplicidade escolhemos um dominio quadrado, e uma grade igual em
cada direção:
h = ∆x = ∆y =
L
n
Onde L é o comprimento do dominio e n o número de pontos. Como
visto anteriormente, las segundas derivadas podem ser escritas de forma
discreta:
uxx '
ui+1,j − 2ui,j + ui−1,j
h
uyy '
ui,j+1 − 2ui,j + ui,j−1
h
Física Computacional
Assim, a equação de Laplace termina sendo:
ui,j =
1
4
[ui,j+1 + ui,j−1 + ui+1,j + ui−1,j ] (12)
A equação de Poisson adiciona uma complexidade na equação anterior
uxx + uyy = −ρ(x , y)
Precisamos só adicionar uma versão discreta de ρ, na equação discreta:
ui,j =
1
4
[ui,j+1 + ui,j−1 + ui+1,j + ui−1,j ] +
h2
4
ρi,j (13)
Física Computacional
As condições de contorno seriam:
ui,0 = gi,0 0 ≤ i ≤ n
ui,L = gi,L 0 ≤ i ≤ n
u0,j = g0,j 0 ≤ j ≤ n
uL,j = gL,j 0 ≤ j ≤ n
Uma "molecula"computacional para encontrar o valor de u no ponto i , i
é mostrada na �gura seguinte
Física Computacional
Para resolver a equação acima, para n = 3 (por simplicidade) podemos
reorganizar a equação acima da forma
Onde ρ̃i,j = h2ρi,j . Isolando as variáveis u11, u12, u21 e u22, i.e, os pontos
internos independentes das condições de contorno, temos:
O lado direito corresponde aos pontos conhecidos no contorno mais os
valores de ρ.
Física Computacional
O sistema pode ser resolvido usando métodos de substituição, porem, já
que a matriz A não é tridiagonal, o calculo pode ser ine�ciente.
Um método mais e�ciente é o conhecido método iterativo de Jacobi. O
método é simples,
I Começamos fazendo um chute inicial para os pontos internos (fora
da fronteira), u(0)i,j .
I Resolvemos depois as Eq. (13) ou Eq. (14), para obtermos u(1)i,j .
I Repetimos o processo iterativamente k vezes, até que u(k)i,j e u
(k−1)
i,j
sejam iguais com um determinado criterio.
Física Computacional
Como aumentar a performance de códigos em python
Porque Python é lento?
I Variáveis dinmâmicas. Contrarioa Fortran ou C++, quando uma
variável é digitada python asigna um tipo.
I Python é uma linguagem interpretada e nao compilada. Python
carece da optimização dada pelo compilador.
I Objetos em Python podem ter um accesso ine�ciente à memoria
(NUMPY melhora o accesso)
Física Computacional
I Por outro lado Fortran e C++ utilizam variáveis previamente
de�nidas, precisam de compiladores que podem aumentar a
performance do código e tém excelente accesso à memoria.
I Rotinas e funções escritas em Fortram ou C podem ser compiladas e
depois interpretadas como módulos de python. Dessa forma a
velocidade de execução de um programa aumenta
consideravelmente. Para isso utilizamos f2py.
I Dessa forma podemos desenvolver códigos rápida e facilmente em
Python e deixar as partes do código que demandam mais computo a
rotinas e�cientes escitas em Fortran (ou C).
Física Computacional
I Normalmente f2py é parte do pacote NUMPY que possue funções
vetorizadas para trabalhar com arrays.
I É importante veri�car se a rotina que queremos converter a modulo
não existe já dentro da biblioteca NUMPY.
I Para compilarmos uma rotina de Fortran usamos (usaremos como
exemplo nossa função rel_func.f90 que resolve a eq. de Laplace em
2 dimensões.
I f2py -c -m rel_func rel_func.f90
I Se tudo der certo, o arquivo rel_func.so será criado
I Para usar o módulo no python, usamos:
I import rel_func
I print rel_func.relax.__doc__
I Bibliogra�a completa de f2py pode ser encotrada em:
https://docs.scipy.org/doc/numpy-dev/f2py/
Física Computacional
Física Computacional
Exercícios
Usando os códigos desenvolvidos em sala de aula para resolver a
equação de difusão, usando os métodos, explícito, implícito e
Crank-Nicolson, resolva o problema 10.4 do livro de Hjorth-Jensen,
partes (a)-(e).
Física Computacional