Diferenças finitas e EDP's

...

Transcript

1. LNCC UFRJ Diferenças finitas e EDP’s Prof. Paulo R. G. Bordoni
2. Vou iniciar apresentando aproximações por diferenças divididas finitas (DDF) para as derivadas de uma função : (, ) → ℝ. Assumirei que as funções envolvidas são deriváveis tantas vezes quantas necessário nos intervalos em questão.
3. + ℎ () ( + ℎ) ′() ′ = lim ℎ→0 + ℎ − () ℎ Lembrem-se, a derivada num ponto ∈ (, ) de uma função : , → ℝ é definida por:
4. + ℎ () ( + ℎ) ′() ′ = lim ℎ→0 + ℎ − () ℎ A intuição por trás dessa imagem é que: “se a curva (azul) definida pelo gráfico da função f é suficientemente suave, então a reta secante (vermelha) tenderá à reta tangente (azul) quando ℎ → 0”.
5. As diferenças divididas finitas (DDF) são definidas pelos quocientes de Newton (h > 0). ′() ≅ − ( − ℎ) ℎ diferença atrasada ′() ≅ + ℎ − ( − ℎ) 2ℎ diferença centrada ′() ≅ + ℎ − () ℎ diferença avançada
6. + ℎ ′() + ℎ − () ℎ diferença avançada − ℎ ′() − ( − ℎ) ℎ diferença atrasada − ℎ ′() + ℎ − ( − ℎ) ℎ + ℎ ℎ diferença centrada Elas correspondem aos gráficos:
7. Quando h diminui a secante “vai para” a tangente! A reta secante, das aproximações por diferenças avançadas (em vermelho), para ℎ = 0.5, 0.2 0.05.
8. Vamos traçar, numa mesma figura, os gráficos das aproximações da derivada () de uma função num ponto . Plotaremos os gráficos das sequências das diferenças avançadas, atrasadas e centradas, para ℎ = 2− e = 0,1, … , 50, com e fixos.
9. Eis o gráfico prometido. Note Surfista, que a DDF centrada aproxima mais rapidamente o valor da derivada que as DDF avançada e retardada. Mestra, você pediu com 50 pontos e não 40.
10. Nesse caso, os erros são causados pelo cancelamento catastrófico e ampliado ao dividirmos por números minúsculos Quando usamos 50 pontos, aparece a propagação de erro da IEEE 754/2008.
11. Veja quando ampliamos a escala vertical:
12. Observem como o erro vai decaindo:
13. Mas depois ele aumenta até atingir um valor máximo:
14. Estou mostrando o gráfico dos erros, em escala logarítmica para que possamos entender melhor o que ocorre.
15. 1. Os erros decaem exponencialmente (retas na escala logarítmica). 2. Os erros por DDF avançadas e retardadas é igual e eles decaem na mesma velocidade (estão superpostos na descida). 3. As DDF centradas convergem mais rapidamente que as outras duas (vermelho abaixo de azul e verde e com inclinação maior). 4. O erro mínimo nas aproximações por DDF centrada é da ordem de 10−11, muito menor que nas outras duas, entre 10−8 e 10−7. 5. Após o mínimo todos os erros crescem, estacionando num valor máximo (horizontal à direita). Observem que:
16. Usando o fato que ′′ = (′ )′ aplicando diferenças avançadas na função ′(), obtemos: ′′() ≅ ′ + ℎ − ′() ℎ Então, aproximando tanto ′( + ℎ) como ′() por diferenças atrasadas, obtemos ′′ ≅ + ℎ − 2 + ( − ℎ) ℎ2
17. Adoro festa de São João!
18. É uma forma de energia, a térmica. Calor É a transferência de energia térmica de uma região onde a temperatura é maior para uma região onde ela é menor. Condução de calor Como modelar a propagação do calor?
19. Existem duas grandezas físicas envolvidas: • A temperatura = (, ) e • O fluxo de calor = (, ) A lei de Fourier informa que o fluxo de calor é proporcional ao gradiente da temperatura.
20. Vamos modelar a propagação de calor numa barra metálica, envolta num isolante, para não haver troca de calor com o meio externo através de sua superfície cilíndrica. T alta (, ) T baixa Mais adiante vamos retirar essa restrição.
21. Indicando por φ o fluxo de calor na direção x a Lei de Fourier fica: , = − (, ) Nessa fórmula: • k é o coeficiente de difusividade térmica, • ρ é a densidade do material, • C é a capacidade calórica do material. T alta (, ) T baixa
22. Considerando um volume elementar ∆ = ∙ ∆ temos: • ∙ é o fluxo de calor entrando • ∙ + ∆ é o fluxo de calor saindo • O calor acumulado nesse volume elementar, na ausência de fontes ou sumidouros, é proporcional a ∆ ∙ ∆ u alta u baixa (, ) ( + ∆, ) ∆
23. Portanto a equação de balanço de energia térmica fica: , − + ∆, ∆ = = ∆ ∙ ∆(, ) u alta u baixa (, ) ( + ∆, ) ∆
24. Portanto, passando ao limite quando ∆ → 0 e ∆ → 0, obtemos a equação de balanço de energia térmica (, ) = (, ) Mas então − + ∆ − ∆ = ∆ ∆ T alta T baixa (, ) ( + ∆, ) ∆
25. Utilizando a Lei de Fourier, obtemos a equação condução de calor: (, ) = 2(, ) 2 (, )
26. Assumindo que a troca de calor com o meio, através da superfície lateral da barra, é proporcional à diferença de temperatura entre a barra e o meio no ponto e no instante , podemos escrever (, ) = 2(, ) 2 + + +
27. Um dos problemas que se coloca é descobrir a distribuição de temperatura na barra num instante final conhecendo a distribuição inicial de temperatura num instante 0 . 0 (, 0 = 0) = 0 ()
28. Para assegurar a unicidade de solução, precisamos prescrever a temperatura da barra nos seus extremos ou o fluxo de calor através deles, durante o tempo transcorrido.
29. Mais condições de contorno: 0, + (0, ) = ℎ() , + (0, ) = , t > 0 Uma distribuição inicial de temperatura: , 0 = 0 , ∈ [0. ] A lei que rege o fenômeno: (, ) = 2(, ) 2 + + + para (, ) ∈ 0, × (0, ∞) Temos então um problema de valor inicial e de contorno: Uma fonte ou sumidouro f de calor: , , ∈ 0, , ≥ 0
30. • No problema de Dirichlet, quando ℎ = = 0 para ≥ 0, temos o problema homogêneo. • No problema de Neumann, quando, ℎ = = 0 para ≥ 0, dizemos que as fronteiras são reflexivas. • Quando = = 0, temos o problema de Dirichlet. • Quando = = 0, temos o problema de Neumann. • Caso , , e sejam diferentes de zero temos o problema de Robin.
31. Vamos resolver o problema de Dirichlet homogêneo. O mais simples de todos: quando = = 0. Para obter uma solução aproximada do problema de valor inicial e de contorno, precisaremos efetuar discretizações espaciais e temporais.
32. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t 3 t 1 t 2 t 18 t 19 t 20 = t F x 29 x 30 Assumindo ∆ = − (+1) construímos uma malha espacial cujos pontos internos ao [, ] são dados por = + ∆, com = 1, ⋯ , . Então, para ∆ dado e = ∆, = 0, 1, 2, ⋯, teremos uma malha temporal num intervalo finito [0 , ], onde = ∆.
33. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t 3 t 1 t 2 t 18 t 19 t 20 = t F x 29 x 30 2,3 = (2 , 3 ) Adotaremos a notação simplificadora: , = ( , ) Com essa notação, a equação , = 2 2 , + ( , ) pode ser reescrita como , = 2, 2 + , , uma igualdade válida para cada ponto interno da malha. 3,1 = (3 , 1 )
34. Mestre, você abusou da notação, não? Sim, na realidade deveria ter escrito , = 2 2 , + , Ah, Loirinha, é mais fácil entender com , = 2, 2 + ,
35. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t i+1 t i t F x 29 x 30 A ideia para resolver numericamente o problema é: conhecidos os valores da temperatura num nível de tempo calcular os valores da temperatura no nível de tempo seguinte +1 usando as aproximações em diferenças finitas. +1, ,
36. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t i+1 t i t F x 29 x 30 +1, , • Aprox. espacial no nível de tempo : 2 , 2 ≅ +1,−2,+−1, (∆)2 • Aprox. espacial no nível de tempo + 1: 2,+1 2 ≅ +1,+1−2,+1+−1,+1 (∆)2 • A passagem entre níveis de tempo (mesmo ): • diferenças avançadas: , ≅ ,+1−, ∆ • diferenças retardadas: , ≅ ,−,−1 ∆
37. Para a equação de calor = 2 2 +, existem três possibilidades de discretização temporal: Discretização temporal Algoritmo Ordem Estabilidade Diferenças avançadas Explícito ( + ℎ2) 2 ℎ2 ≤ 1/2 Diferenças atrasadas Implícito ( + ℎ2) incondicional Cranck-Nicholson Média (2 + ℎ2) incondicional
38. Tais esquemas são descritos pelas moléculas de discretização: Cranck-Nicholson t i+1 t i −1 +1 Implícito t i+1 t i −1 +1 Explícito t i+1 t i −1 +1 Nos 3 esquemas: indica valor conhecido indica valor a calcular
39. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t i+1 t i t F x 29 x 30 No método explícito, você calcula diretamente o valor de +1, em função dos valores conhecidos de ,−1 , , , ,+1 . +1, ,−1 , , , ,+1
40. Como a igualdade , = 2, 2 + , é válida para cada ponto interno da malha, usando diferenças avançadas no tempo, teremos: ,+1 − , ∆ = +1, − 2, + −1, (∆)2 + , Método explícito
41. Isolando o termo +1, envolvendo o tempo +1 na posição , obtemos o método (algoritmo) explícito: ,+1 = + 1 − 2 , + +1, + −1, +, para = 0,1, … , e = 1,2, … , , e ,0 = (0 ) para = 1, … , onde = _ ∗ ∆/∆2 = ∆/(_ ∗ ℎ)
42. Adotando a notação: = 0 , 1, 2, ⋯ , , e considerando a matriz = 1 − 2 1 − 2 1 − 2 Então a equação da transparência anterior fica +1 = + para = 0, … , com 0 dado.
43. Notem que então podemos calcular 1 , 2 , 3 , … pela fórmula +1 = + para = 0, … , com 0 dado.
52. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t i+1 t i t F x 29 x 30 No método implícito, você precisará resolver um sistema linear envolvendo os valores de +1,+1 , +1, , +1,+1 em função dos valores conhecidos de , . +1,−1 , +1, , +1,+1 ,
53. Como a igualdade , = 2, 2 + , é válida para cada ponto interno da malha, usando diferenças avançadas no tempo, teremos: ,+1 − , ∆ = +1,+1 − 2,+1 + −1,+1 (∆)2 + , Método implícito
54. ,+1 − , ∆ = 2(∆)2 ൛ (+1,+1 −2,+1 + −1,+1 ) ൟ + +1, − 2, + −1, Fazendo a média da parte espacial, obtemos o método de Crank-Nicholson:
55. x 0 = a x 3 x 31 = b x 1 x 2 t 0 = 0 t i+1 t i t F x 29 x 30 No método de Crank-Nicholson, você precisará resolver um sistema linear envolvendo uma média entre os valores do método implícito e explícito. +1,−1 , +1, , +1,+1 ,−1 , , , ,+1
56. Atenção com as condições iniciais e as de contorno. Sim, quando = 0, a condição inicial fornece , 0 = ,0 = ( ), para = 0, … , + 1. Com as condições homogêneas, para = 0, … , e: • = 0 teremos , = 0, = 0, • = + 1 teremos , = +1 , = 0.
57. • −+1 ,+ + [2 − ൗ (∆)2 (∆)] ,+1 − −1 , +1 = ൗ [ (∆)2 ∆ ] , , para k = 1, … , • ,0 = , para = 0,1, … , + 1 • ,0 = 0, para = 0, … , • , +1 = 0, para = 0, … , Resolvendo a equação do método implícito para o nível de tempo, + 1 recairemos num sistema linear de equações à incógnitas:
58. Trata-se de um sistema tridiagonal, que pode ser resolvido usando-se métodos iterativos, para valores muito grandes de .
59. 2 − −1 −1 2 − −1 ⋱ −1 2 − −1 −1 2 − 1 , 1 2 , 1 ⋮ −1 , 1 , 1 = 1 2 ⋮ −1 A solução do problema de valor inicial e de contorno homogêneo de Dirichlet, no 1º nível de tempo, é a solução do sistema linear:
60. O sistema linear que o Surfista achou pode ser reescrito na forma matricial como ∙ 1 = 0 . Ele fornece a solução no 1º nível de tempo. Nos demais, é só acharmos a solução de ∙ +1 = .
61. Agora vamos aproximar numericamente a solução de uma equação do tipo = 2 2 + 2 2 +
62. • = 2 2 + 2 2 + em Ω × (0, 0 ) • , , 0 = 0 (, ) em Ω • , , prescrita na borda Ω × (0, 0 ) Uma fonte de calor na placa Ω. Trata-se de um problema de propagação de calor numa placa de metal retangular Ω = (0, 0 ) × (0, 0 ).
63. = = = Surfista, imagine que na placa de metal, as bordas são mantidas às temperaturas indicadas (em graus centígrados). =
64. Construiremos uma discretização do retângulo 0, 0 × [0, 0 ], através de uma malha uniforme envolvendo n subintervalos na direção x e outros m subintervalos na direção y. No exemplo = = 4. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
65. Nos pontos internos da malha (em vermelho) aproximaremos as derivadas parciais por diferenças finitas. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j O valor da temperatura U nos pontos da borda da placa (em amarelo, azul e verde) é conhecido.
66. A aproximação em diferenças finitas para as derivadas 2ªs da temperatura U em pontos internos da malha obedece o esquema em cruz indicado abaixo para o ponto 3. 2 3 2 ≃ − 23 + 6 ℎ2 2 3 2 ≃ 2 − 23 + ℎ2 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
67. Aplicando o esquema em cruz para os pontos internos 1, 2, … , 9 ficaremos com o sistema de 9 equações a 9 incógnitas 1 , 2 , … , 9 indicado abaixo: 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j − 21 + 2 ℎ2 + − 21 + 4 ℎ2 = 0 1 − 22 + 3 ℎ2 + − 22 + 5 ℎ2 = 0 2 − 23 + ℎ2 + − 23 + 6 ℎ2 = 0 ... 8 − 29 + ℎ2 + 6 − 29 + ℎ2 = 0
68. Assumindo = = = = = = 0 , = = = 100 e = ℎ = = 200 graus centígrados e simplificando obtemos o sistema linear abaixo, que é penta-diagonal. Troque T por U
69. Ele foi feito apenas para mostrar a parte espacial das aproximações. Para a parte temporal use um esquema de Crank-Nicholson do problema de 1 dimensão espacial (da barra). Muitíssima atenção que esse sistema não está levando em consideração a variação temporal da temperatura.
70. É claro que para resolver esse sistema linear usaremos o pacote scipy.linalg. Como a matriz desse sistema linear, com sinal trocado, é semi-definida- positiva e de banda, poderemos utilizar cholesky_ banded( ) e cho_solve_banded( ).
71. Um programa para obter a solução x desse sistema linear.
72. Pois é Surfista para um sistema esparso como esse, os métodos apropriados para resolver o sistema linear são os métodos iterativos. Entretanto, como as matrizes envolvidas serão penta-diagonais e definidas-positiva (por quê?) fica mais fácil usar cholesky_banded() e cho_solve_banded().
73. Dizem que, “Simpatia é quase amor”. Meu bloco preferido!
74. Condição inicial + Variação temporal + Condições de contorno + (, ) : Meus blocos até parecem “confeti” !
75. Tchau, até a próxima aula!