Cálculo Numérico, trabalho final

O trabalho final, para encerrar o curso.

Transcript

1. LNCC UFRJ Trabalho final do curso Prof. Paulo R. G. Bordoni
2. Este trabalho envolve a resolução numérica de EDP’s utilizando a resolução de sistemas lineares. EDP’s Equações Diferenciais Parciais
3. Serão dois problemas, propostos no texto na mesma cor que esta daqui. Para devolver até as 12:00 h do dia 15/07, da mesma forma que no trabalho anterior.
4. + ℎ () ( + ℎ) ′() ′ = → + − () Lembrem-se, nesta definição de derivada, o h pode ser positivo ou negativo, nunca nulo. Iniciaremos analisando algumas aproximações para a derivada 1ª de uma função : (, ) → ℝ.
5. + ℎ () ( + ℎ) ′() ′ = → + − () 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. Para ℎ > 0, a aproximação por diferença avançada para a derivada 1ª de f é definida por: + ℎ ′() − ( − ℎ) ℎ ′ ≅ + − ()
7. − ℎ − ( − ℎ) ′() ℎ Novamente para ℎ > 0, a aproximação por diferença atrasada para a derivada 1ª de f é definida por: ′ ≅ − ( − )
8. − ℎ + ℎ − ( − ℎ) ′() ℎ ℎ + ℎ A aproximação por diferença centrada para a derivada 1ª de f é definida por (novamente ℎ > 0): ′ ≅ + − ( − )
9. Mestra, pelas imagens percebo que a aproximação por diferença centrada é melhor. Há garantias matemáticas disto? Sim Filhota, é algo que passaremos a examinar graficamente.
10. A reta secante, das diferenças avançadas (em vermelho), para ℎ = 0.5, 0.2 0.05.
11. O programa que traça os gráficos exibidos. Surfista, faça programas para as outras duas aprox.
12. Neste programa mostro, graficamente, a convergência das aproximações por diferenças avançadas, atrasadas e centradas para ′ , = /4 e = sen().
13. O programa permite escolher o valor de x, a quantidade n de pontos e a escala-y. Observem que as diferenças centradas convergem bem mais rapidamente. ′(/4 ) = 0.70711
14. Usando o fato que ′′ = (′ )′ com diferenças avançadas, obtemos: ′′() ≅ ′ + ℎ − ′() ℎ Então, aproximando tanto ′( + ℎ) como ′() por diferenças atrasadas, obtemos ′′ ≅ + ℎ − 2 + ( − ℎ) ℎ2
15. Adoro festa de São João!
16. É 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?
17. 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.
18. Vamos começar modelando 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
19. 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
20. 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 ∆ ∙ ∆ T alta T baixa (, ) ( + ∆, ) ∆
21. Portanto a equação de balanço de energia térmica fica: , − + ∆, ∆ = = ∆ ∙ ∆(, ) T alta T baixa (, ) ( + ∆, ) ∆
22. 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 (, ) ( + ∆, ) ∆
23. Utilizando a Lei de Fourier, obtemos a equação condução de calor: (, ) = 2(, ) 2 (, )
24. Mais condições no contorno: • (0, ) = ℎ() • , = , t > 0 Uma condição inicial: • , 0 = , ∈ [0. ] = 2 2 , (, ) ∈ 0, × (0, ∞) x = 0 x = L O problema de evolução da temperatura T (x, t) numa barra Um problema de valor inicial e de contorno.
25. 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) os pontos internos ao intervalo [, ] são dados por = + ∆, com = 1, ⋯ , . Então, para = ∆, = 1, 2, ⋯, teremos: , = 2 2 , , = 1, ⋯ , , = = 1,2, ⋯
26. Notação simplificadora: , = ( , ) 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 = (2 , 3 )
27. Lembrem-se Mestres, 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
28. Tais esquemas são descritos pelas moléculas de discretização: t i+1 t i −1 +1 Explícito t i+1 t i −1 +1 Cranck-Nicholson t i+1 t i −1 +1 Implícito
29. 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 Em cada ponto , do interior do retângulo, usaremos as aproximações: , ≅ , − ,−1 ∆ 2 2 , ≅ +1, − 2, + , − −1, (∆)2
30. 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 Assim, em cada ponto , do interior do retângulo, vale a igualdade: , − ,−1 ∆ = +1, − 2, + −1, (∆)2 Para = 0 e = + 1 a temperatura é prescrita e também para = 0.
31. 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 Então, para = 1,2, ⋯ , e = 1,2, ⋯, vale a igualdade: −−1, + 2 + , − +1, = ,−1 . Além disso, para: = 1,2, ⋯ , → 0, = ℎ = ℎ( ) = 1,2, ⋯ , → +1, = = ( ) = 1,2, ⋯ , , → ,0 = = ( ) = (∆)2 ∆
32. x 0 = a x 3 x 31 = b x 1 x 2 t = 0 t = t 1 x 29 x 30 3,1 Portanto, os valores da temperatura no nível de tempo 1 são dados pela solução do sistema linear −−1,1 + 2 + ,1 − +1,1 = , = 1,2, ⋯ , com 0,1 = ℎ1 e +1,1 = 1 = (∆)2 ∆
33. Na forma matricial, ele fica = onde: = 1 2 ⋮ , = ,1 = −1 −1 −1 −1 ⋱ ⋱ ⋱ −1 −1 = 1 2 ⋮ −1 + ℎ1 0 ⋮ 0 1 = 2 + , = (∆)2 ∆
34. Para os outros níveis de tempo 2 , 3 , ⋯ basta substituir o termo independente: 1 = 1 2 ⋮ −1 + ℎ1 0 ⋮ 0 1 2 = 1 2 ⋮ −1 + ℎ2 0 ⋮ 0 2
35. • Escolha um intervalo 0, e uma função inicial : 0, → ℝ; • Faça uma discretização espacial 0 = 0, 1 , … , = com pelo menos 30 sub-intervalos igualmente espaçados de 0, ; • Escolha um instante de tempo final e faça uma discretização temporal 0 = 0, 1 , … , = compatível com a espacial, para o intervalo temporal [0, ]; • Construa a solução numérica do problema de valor inicial e de contorno, com condições de contorno homogêneas, conforme descrito nas últimas transparências; • Apresente a solução (, ) graficamente, no domínio 0, × [0, ]. Eis o 1º problema:
36. T = 0 ºC T = 0 ºC T = 10 ºC T = 20 ºC Surfista, imagine uma placa de metal, quadrada, cujas bordas são mantidas às temperaturas indicadas. Como será a distribuição de temperatura no interior da placa em estado estacionário?
37. = 2 2 + 2 2 + Poderia existir uma fonte de calor Trata-se de um problema de propagação de calor em duas dimensões: Insisto. Eu quero o estado estacionário!
38. Depois que você estaciona um carro, ele não se move mais, certo? Então, quando a temperatura não varia mais com o tempo, ela está em estado estacionário. Ih, Mestre, o que é estado estacionário?
39. T = 0 ºC T = 0 ºC T = 100 ºC T = 200 ºC 2 2 + 2 2 = 0 no interior da placa Mais as condições no contorno Se não há variação de temperatura então = 0 e as equações que modelam o fenômeno físico são:
40. Essas equações são conhecidas como: • equação de Laplace (qdo. = 0) • e equação de Helmholtz (qdo. ≠ 0). Uma justa homenagem a esses cientistas!
41. Construiremos uma discretização do quadrado 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
42. 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.
43. 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
44. Aplicando o esquema em cruz para os pontos internos 1, 2, ... , 9, ficaremos com o sistema de 9 equações a 9 incógnitas T 1 , T 2 , ..., T 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 ...
45. Simplificando obtemos o sistema linear abaixo, que é penta-diagonal.
46. É 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( ).
47. Um programa para obter a solução x desse sistema linear.
48. Agora Surfista, imagine construirmos uma malha com 50 pontos internos na horizontal por outros 50 pontos internos na vertical. Ao todo teremos 2.500 pontos internos e, consequentemente, um sistema linear de 2.500 equações à 2.500 incógnitas! Impressionante Mestra! Mas continuará sendo um sistema penta- diagonal. A única diferença é que a 2ª sobre diagonal iniciará na coluna 51 e a 2ª sub diagonal na linha 51.
49. Pois é Surfista para um sistema esparso como esse, os métodos apropriados para resolver o sistema linear são os métodos iterativos.
50. No conjunto de transparências sobre métodos iterativos estacionários o Mestre apresentou uma forma de descrever a solução dos métodos de Jacobi e Gauss-Seidel através de equações. Inclusive o Sherlock explicou porquê e o Manuel apresentou onde buscar informações na Scipy.
51. • Escolha um retângulo = [0, ] × [0, ] e uma fonte de calor descrita por uma função : → ℝ; • Faça uma discretização espacial 0 = 0, 1 , … , = com pelo menos 25 sub-intervalos igualmente espaçados de 0, e outra com pelo menos 30 sub-intervalos de [0, ]; • Imponha condições de contorno com fluxo nulo em duas laterais de e de temperatura prescrita nas outras duas laterais; • Construa a solução numérica do problema desse problema de contorno; • Apresente a solução (, ) graficamente, no domínio 0, × [0, ]. • Use o método dia_matrix( ) da classe scipy.sparse para guardar a matriz do problema. Eis o 1º problema:
52. Mãs à obra, é realmente fácil!
53. O exemplo ajuda. As matrizes do problema possuem estrutura com diagonais não nulas.
54. Tchau, até nossa happy- hour!