Diferenças finitas e EDP's Atualizada

As questões da P1 "online" de Cálculo Numérico 2016-2, com a questão 4 atualizada.

Transcript

1. LNCC UFRJ Diferenças finitas e EDP’s Prof. Paulo R. G. Bordoni
2. Voces terão duas semanas para devolver esta prova online.
3. 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.
4. + ℎ () ( + ℎ) ′() ′ = lim ℎ→0 + ℎ − () ℎ Lembrem-se, a derivada num ponto ∈ (, ) de uma função : , → ℝ é definida por:
5. + ℎ () ( + ℎ) ′() ′ = 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 (verde) quando ℎ → 0”.
6. A aproximação por diferença avançada, ℎ +(), para a derivada 1ª de f em ∈ (, ) é definida por: + ℎ ′() − ( − ℎ) ℎ ℎ +() = + ℎ − () ℎ , ℎ > 0
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. 1ª questão: Fazer um programa que receba uma função ↦ (), num intervalo , dado e mostre as três aproximações desenhadas na transparência anterior.
9. − ℎ − ( − ℎ) ′() ℎ A aproximação por diferença atrasada, ℎ −(), para a derivada 1ª de f em ∈ (, ) é definida por: ℎ −() = − ( − ℎ) ℎ , ℎ > 0
10. − ℎ + ℎ − ( − ℎ) ′() ℎ ℎ + ℎ A aproximação por diferença centrada, ℎ ×(), para a derivada 1ª de f em ∈ (, ) é definida por: ℎ × = + ℎ − ( − ℎ) 2ℎ , ℎ > 0
11. Mestra, pelas imagens percebo que a aproximação por diferença centrada é melhor. Há garantias matemáticas disto? Sim Filhota. É algo que, por hora, examinaremos graficamente. A prova envolve o conhecimento de série de Taylor, que veremos depois.
12. 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 ℎ +() , ℎ −() e ℎ ×() para ℎ = 2− e = 0,1, …, 50, com e fixos.
13. 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.
14. 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.
15. Veja quando ampliamos a escala vertical:
16. 2ª questão: Fazer um programa que receba uma função ↦ (), num intervalo , dado e um ponto 0 ∈ [, ] e mostre as três sequências desenhadas na transparência anterior e os gráficos em escala logarítmica mostrados a seguir.
17. Observem como o erro vai decaindo:
18. Mas depois ele aumenta até atingir um valor máximo:
19. Estou mostrando o gráfico dos erros, em escala logarítmica para que possamos entender melhor o que ocorre.
20. 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:
21. 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
22. Adoro festa de São João!
23. É 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?
24. 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.
25. 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. Talta (, ) Tbaixa
26. 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. Talta (, ) Tbaixa
27. 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 ∆ ∙ ∆ Talta Tbaixa (, ) ( + ∆, ) ∆
28. Portanto a equação de balanço de energia térmica fica: , − + ∆, ∆ = = ∆ ∙ ∆(, ) Talta Tbaixa (, ) ( + ∆, ) ∆
29. Portanto, passando ao limite quando ∆ → 0 e ∆ → 0, obtemos a equação de balanço de energia térmica (, ) = (, ) Mas então − + ∆ − ∆ = ∆ ∆ Talta Tbaixa (, ) ( + ∆, ) ∆
30. Utilizando a Lei de Fourier, obtemos a equação condução de calor: (, ) = 2(, ) 2 (, )
31. 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 + +
32. 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) = ()
33. 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.
34. Mais condições de contorno: 0, + (0, ) = ℎ() , + (0, ) = , t > 0 Uma distribuição inicial de temperatura: , 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:
35. • 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 são diferentes de zero temos o problema de Robin.
36. 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.
37. 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 ∆ = ( − )/ 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 = ∆.
38. 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.
39. Mestre, você abusou da notação, não? Sim, na realidade deveria ter escrito , = 2 2 , Ah, Loirinha, é mais fácil com , = 2, 2
40. 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, ,
41. 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 ∆
42. Para a equação de calor = 2 2 existem três possibilidades de discretização temporal: Discretização temporal Esquema 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
43. 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
44. 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, ,
45. 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, ,
46. 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, ,
47. 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 ,+1 − , ∆ = +1,+1 − 2,+1 + −1,+1 (∆)2 Método implícito
48. ,+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:
49. 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.
50. • ,+1 = ( Τ ∆ (∆)2 +1 , − 2, + −1 , , k = 1, …, − 1 • ,0 = , = 1, … , − 1 • ,0 = 0, = 0, … , • , +1 = 0, = 0, … , Resolvendo a equação do método explícito para ,+1 teremos:
51. • −+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:
52. Trata-se de um sistema tridiagonal, que pode ser resolvido usando-se métodos iterativos, para valores muito grandes de .
53. 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:
54. 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 = .
55. Fazer um programa que receba os valores numéricos de • , , , , , , e as funções: • A distribuição inicial de temperatura inicial ↦ () em [0, ], • ↦ ℎ() e ↦ () em 0, que definem as condições de contorno, e resolva o problema de valor inicial e de contorno descrito na próxima transparência pelo método de Crank- Nicholson e mostre a solução graficamente. 3ª questão:
56. (, ) = 2(, ) 2 + + para (, ) ∈ 0, × (0, ∞) • , 0 = , ∈ [0.] • 0, + (0,) = ℎ , 0 < t ≤ τ, • , + (0,) = , 0 < t ≤ τ. 3ª questão, continuação.
57. Agora vamos aproximar numericamente a solução de uma equação do tipo = 2 2 + 2 2 +
58. • = 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 ).
59. = = = Surfista, imagine que na placa de metal, as bordas são mantidas às temperaturas indicadas (em graus centígrados). =
60. 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
61. 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 T nos pontos da borda da placa (em amarelo, azul e verde) é conhecido.
62. A aproximação em diferenças finitas para as derivadas 2ªs da temperatura T 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
63. 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
64. Assumindo = = = = = = 0 , = = = 100 e = ℎ = = 200 graus centígrados e simplificando obtemos o sistema linear abaixo, que é penta-diagonal.
65. 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.
66. É 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( ).
67. Um programa para obter a solução x desse sistema linear.
68. 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().
69. Dizem que, “Simpatia é quase amor”. Meu bloco preferido!
70. Condição inicial + Variação temporal + Condições de contorno + (, ) : Meus blocos até parecem “confeti” !
71. 4ª questão: Fazer um programa que receba os dados do problema descrito pelo Mestre e apresente a solução graficamente através de uma animação colorida de curvas de nível. Usem o método de Crank-Nicholson.
72. Tchau, até a próxima aula!