Diferenças finitas e equações diferenciais parciais

Veja a apresentação no site do Professor.

Transcript

1. LNCC UFRJ Diferenças finitas e equações diferenciais parciais Prof. Paulo R. G. Bordoni
2. Vamos iniciar estudando aproximações por diferenças divididas finitas (DDF) para as derivadas de uma função : (, ) → ℝ. Assumiremos 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 (verde) quando ℎ → 0”.
5. A aproximação por diferença avançada, ℎ +(), para a derivada 1ª de f em ∈ (, ) é definida por: + ℎ ′() − ( − ℎ) ℎ ℎ +() = + ℎ − () ℎ , ℎ > 0
6. 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.
7. O programa que traça os 3 gráficos exibidos.
8. − ℎ − ( − ℎ) ′() ℎ A aproximação por diferença atrasada, ℎ −(), para a derivada 1ª de f em ∈ (, ) é definida por: ℎ −() = − ( − ℎ) ℎ , ℎ > 0
9. − ℎ + ℎ − ( − ℎ) ′() ℎ ℎ + ℎ A aproximação por diferença centrada, ℎ ×(), para a derivada 1ª de f em ∈ (, ) é definida por: ℎ × = + ℎ − ( − ℎ) 2ℎ , ℎ > 0
10. 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.
11. 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.
12. 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.
13. 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.
14. Veja quando ampliamos a escala vertical:
15. Observem como o erro vai decaindo:
16. Mas depois ele aumenta até atingir um valor máximo:
17. Estou mostrando o gráfico dos erros, em escala logarítmica para que possamos entender melhor o que ocorre.
18. 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:
19. 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
20. Adoro festa de São João!
21. É 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?
22. 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.
23. 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
24. 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
25. 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 (, ) ( + ∆, ) ∆
26. Portanto a equação de balanço de energia térmica fica: , − + ∆, ∆ = = ∆ ∙ ∆(, ) Talta Tbaixa (, ) ( + ∆, ) ∆
27. Portanto, passando ao limite quando ∆ → 0 e ∆ → 0, obtemos a equação de balanço de energia térmica (, ) = (, ) Mas então − + ∆ − ∆ = ∆ ∆ Talta Tbaixa (, ) ( + ∆, ) ∆
28. Utilizando a Lei de Fourier, obtemos a equação condução de calor: (, ) = 2(, ) 2 (, )
29. 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 + +
30. 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) = ()
31. 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.
32. 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:
33. • 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.
34. 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.
35. 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 = ∆.
36. 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.
37. Mestre, você abusou da notação, não? Sim, na realidade deveria ter escrito , = 2 2 , Ah, Loirinha, é mais fácil com , = 2, 2
38. 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, ,
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 +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 ∆
40. 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
41. 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
42. 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, ,
43. 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, ,
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 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, ,
45. 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
46. ,+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:
47. 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.
48. • ,+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:
49. • −+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:
50. Trata-se de um sistema tridiagonal, que pode ser resolvido usando-se métodos iterativos, para valores muito grandes de .
51. 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:
52. 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 = .
53. Agora vamos aproximar numericamente a solução de uma equação do tipo = 2 2 + 2 2 +
54. = 2 2 + 2 2 + Uma fonte de calor Trata-se de um problema de propagação de calor em duas dimensões: Eu quero resolver o problema de estado estacionário!
55. 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?
56. Essa equação é conhecida como equação de Helmholtz. Para = 0 e temos a equação de Laplace. Uma justa homenagem a esses cientistas! Para o estado estacionário ( = 0 ) a equação fica 2 2 + 2 2 =
57. T = 0 ºC T = 0 ºC T = 10 ºC T = 20 ºC Surfista, imagine agora 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?
58. 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 Sem fonte de calor, a equação que modela o fenômeno físico é:
59. 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
60. 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.
61. 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
62. 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
63. Assumindo = = = = = = 0 , = = = 100 e = ℎ = = 200 graus centígrados e simplificando obtemos o sistema linear abaixo, que é penta-diagonal.
64. É 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( ).
65. Um programa para obter a solução x desse sistema linear.
66. 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.
67. Pois é Surfista para um sistema esparso como esse, os métodos apropriados para resolver o sistema linear são os métodos iterativos. Existem métodos iterativos não- estacionários para resolução de sistemas lineares desse tipo. Mas agora, neste Carnaval, vou sair no meu bloco em Ipanema.
68. Dizem que, “Simpatia é quase amor”. Meu bloco preferido!
69. Meus blocos até parecem “confeti” !
71. Então Surfista, faça isto para aquele de 2.500 equações à 2.500 incógnitas! Prometo que farei, Mestre.
72. Tchau, até a próxima aula!