Diferenças finitas e EDP's

Transcript

1. LNCC UFRJ Diferenças finitas e EDP’s Prof. Paulo R. G. Bordoni
2. Agora vamos ter os girassóis Do fim do ano E o calor vem desumano Tudo irá se expandir Crescer com as águas Quiçá, amores nos corações E um santeiro, Milagreiro Prevê a dor De terceiros E diz que a vida É feita de ilusão ... Girassóis me lembram a Cássia Eller
3. Iniciarei com aproximações denominadas 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 forem necessárias 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 (azul) quando ℎ → 0”.
6. 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
7. + ℎ ′() + ℎ − () ℎ diferença avançada − ℎ ′() − ( − ℎ) ℎ diferença atrasada − ℎ ′() + ℎ − ( − ℎ) ℎ + ℎ ℎ diferença centrada Elas correspondem aos gráficos:
8. 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.
9. 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, … , , com e fixos e ~40 − 50.
10. 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. Para = temos ′ = e ′ 1 = , o limite da sequência.
11. Nesse caso, os erros são causados pelo cancelamento catastrófico e ampliado ao dividirmos por números minúsculos Quando usamos 50 pontos, começa a aparecer a propagação de erro da IEEE 754/2008.
12. Veja quando ampliamos a escala vertical:
13. Observem na tabela como o erro vai decaindo:
14. Mas depois ele aumenta até atingir um valor máximo:
15. Estou mostrando o gráfico dos erros, em escala logarítmica para que possamos entender melhor o que ocorre.
16. 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:
17. Usando o fato que ′′ = (′ )′ aplicando diferenças avançadas à derivada ′ da função num ponto ∈ () obtemos: ′′() ≅ ′ + ℎ − ′() ℎ Então, aproximando tanto ′( + ℎ) como ′() por diferenças atrasadas, obtemos a aproximação para a derivada 2ª: ′′ ≅ + ℎ − 2 + ( − ℎ) ℎ2
18. Adoro festa de São João!
19. É 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?
20. 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.
21. 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. u alta (, ) u baixa Mais adiante vamos retirar essa restrição.
22. 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. u alta (, ) u baixa
23. 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 (, ) ( + ∆, ) ∆
24. Portanto a equação de balanço de energia térmica fica: , − + ∆, ∆ = = ∆ ∙ ∆(, ) u alta u baixa (, ) ( + ∆, ) ∆
25. Portanto, passando ao limite quando ∆ → 0 e ∆ → 0, obtemos a equação de balanço de energia térmica (, ) = (, ) Mas então − + ∆ − ∆ = ∆ ∆ u alta u baixa (, ) ( + ∆, ) ∆
26. Utilizando a Lei de Fourier, obtemos a equação condução de calor: (, ) = 2(, ) 2 na qual resume as constantes físicas. u(, )
27. 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 + + +
28. Um dos problemas de grande interesse em engenharia (pense em vigas de metal ou em oleodutos) é descobrir a distribuição de temperatura ao longo de uma barra transcorrido um tempo = dada uma distribuição inicial de temperatura em = 0 . 0 (, = 0 ) = 0 ()
29. Para assegurar a unicidade de solução, precisamos prescrever a temperatura da barra nos seus extremos ou o fluxo de calor através de um deles, durante o tempo transcorrido.
30. 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
31. • 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 (fluxo nulo). • 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.
32. Para resolver o problema de Dirichlet não-homogêneo basta acrescentar uma solução particular (, ) à solução geral (, ) do problema homogêneo. Sim Colega! E como a derivada segunda de uma função linear na variável é nula, podemos usar a solução particular , = ℎ + Τ − ℎ ( − ) .
33. Assim, no caso do problema de Dirichlet, basta resolvermos o caso homogêneo. O mais simples de todos, aquele em que = = 0. Para obter uma solução aproximada do problema de valor inicial e de contorno, precisaremos efetuar discretizações espaciais e temporais.
34. 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 = ∆.
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 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 )
36. 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 + ,
37. 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, ,
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 +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 ∆
39. 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
40. 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
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 No método explícito, você calcula diretamente o valor de +1, em função dos valores conhecidos de ,−1 , , , ,+1 . +1, ,−1 , , , ,+1
42. 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
43. 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, … , , com = Τ ∆ ∆2 . Para iniciar, ,0 = 0 ( ) para = 1, … ,
44. 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.
45. Notem que então podemos calcular 1 , 2 , 3 , … pela fórmula +1 = + para = 0, … , com 0 dado.
49. A estabilidade exige ≪ . Confiram os valores: = 0.05 e = 0.00125.
51. Notem que a matriz é extremamente mau condicionada!
52. Ficou lindo o gráfico da solução!
53. Sim, esta é outra visualização.
54. Estourei a condição de estabilidade.
55. O início e o fim da matriz , muito bem condicionada!
56. Entretanto eis o visual da instabilidade:
57. Torno a repetir, a estabilidade exige ≪ . Observem que = Τ 40.
58. 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 O método implícito, resolve esse problema. Entretanto, em cada passo temporal, 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 ,
59. ,+1 − , ∆ = +1,+1 − 2,+1 + −1,+1 (∆)2 + , Método implícito Como a igualdade , = 2, 2 + , é válida para cada ponto interno da malha, usando diferenças avançadas no tempo, teremos (com = 1):
60. Explicitando os termos no instante + 1 em função dos valores no instante anterior obtemos: − −1, +1 + 1 + 2 , +1 − +1, +1 = , + ∆ , Método implícito = Τ ∆ ∆ 2 > 0 Essas igualdades são válidas para os pontos internos da malha, isto é, = 1, 2, … ,
61. = 1 → − 0, +1 + 1 + 2 1, +1 − 2, +1 = 1, + ∆ 1, = 2 → − 1, +1 + 1 + 2 2, +1 − 3, +1 = 2, + ∆ 2, ⋮ = 30 ⋮ → ⋮ − 29, +1 + 1 + 2 30, +1 − 31, +1 = 30, + ∆ 30, Método implícito Escrevi as igualdades ali embaixo, para = 30: Na primeira equação a parcela − 0, +1 é conhecida e na última − 31, +1 também. Elas envolvem o valor de na fronteira. = Τ ∆ ∆ 2 > 0
62. 1 + 2 − − 1 + 2 − 1 + 2 1, +1 2, +1 30, +1 = 1, 2, 30, + ∆ 1, 2, 30, + 0, +1 31, +1 Método implícito = Τ ∆ ∆ 2 > 0 Na forma matricial o sistema fica:
63. Adotando a notação: = 1, 2, ⋯ , , = 1, 2, ⋯ , , +1 0 = 1,+1 0 ⋯ 0 ,+1 , = 1 + 2 − − 1 + 2 − − − 1 + 2 Então a equação da transparência anterior fica +1 = + ∆ + α + 0 , para = 0, … , com 0 dado.
64. Trata-se de um sistema tri-diagonal, que pode ser resolvido usando-se métodos iterativos, para valores muito grandes de .
65. 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.
66. = 1 + 2 − − 1 + 2 … … − − 1 + 2 Método implícito = 1 − 2 1 − 2 … … 1 − 2 Método explícito Comparando as matrizes dos métodos fica clara a condição de estabilidade para o explícito: Sua matriz deixa de ser diagonalmente dominante, ao contrário do que acontece com o método implícito. = Τ ∆ ∆ 2 ≥ 1/2 Instável quando
67. 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. Mas isto fica para você, Surfista! +1,−1 , +1, , +1,+1 ,−1 , , , ,+1
68. ,+1 − , ∆ = 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:
69. Agora vamos aproximar numericamente a solução de uma equação do tipo = 2 2 + 2 2 +
70. • = 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 ).
71. Construiremos uma discretização do retângulo 0, 0 × [0, 0 ], através de uma malha uniforme envolvendo 4 subintervalos na direção x e outros 4 subintervalos na direção y. Utilizaremos diferenças finitas para aproximara solução. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
72. Nos pontos internos da malha (em azul) aproximaremos as derivadas parciais por diferenças finitas. O valor da temperatura U nos pontos da borda da placa (em vermelho) é conhecido. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
73. A aproximação por 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
74. 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: − 21 + 2 ℎ2 + − 21 + 4 ℎ2 = 1 1 − 22 + 3 ℎ2 + − 22 + 5 ℎ2 = 2 2 − 23 + ℎ2 + − 23 + 6 ℎ2 = 3 ... 8 − 29 + ℎ2 + 6 − 29 + ℎ2 = 9 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
75. −41 + + 2 + + 4 = ℎ21 −42 + 1 + 3 + + 5 = ℎ22 −43 + 2 + + + 6 = ℎ23 −44 + + 5 + 1 + 7 = ℎ24 −45 + 4 + 6 + 2 + 8 = ℎ25 −46 + 5 + + 3 + 9 = ℎ26 −47 + + 8 + 4 + = ℎ27 −48 + 7 + 9 + 5 + ℎ = ℎ28 −49 + 8 + + 6 + = ℎ29 Reescrevi as equações e marquei em vermelho os termos da fronteira, conhecidos e os valores da fonte em laranja. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
76. −41 + 2 + 4 = ℎ21 − ( + ) −42 + 1 + 3 + 5 = ℎ22 − −43 + 2 +6 = ℎ23 −( + ) −44 + 5 + 1 + 7 = ℎ24 − −45 + 4 + 6 + 2 + 8 = ℎ25 −46 + 5 + 3 + 9 = ℎ26 − −47 +8 + 4 = ℎ27 −( + ) −48 + 7 + 9 + 5 = ℎ28 −ℎ −49 + 8 +6 = ℎ29 −( + ) Rearranjando as equações na forma: ó = . 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
77. 1 2 3 4 5 6 U7 8 9 = ℎ2 1 2 3 4 5 6 7 8 9 − + + 0 + ℎ + −4 1 0 1 −4 1 0 1 −4 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 −4 1 0 1 −4 1 0 1 −4 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 −4 1 0 1 −4 1 0 1 −4 1 0 0 0 1 0 0 0 1 Colocando o sistema na forma matricial, obtemos um sistema = ℎ2 − onde é tridiagonal de blocos. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j A fonte de calor (, ) As condições de contorno .
78. Dizem que, “Simpatia é quase amor”. Meu bloco preferido!
79. 1 2 3 4 5 6 U7 8 9 = ℎ2 1 2 3 4 5 6 7 8 9 − + + 0 + ℎ + −4 1 0 1 −4 1 0 1 −4 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 −4 1 0 1 −4 1 0 1 −4 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 −4 1 0 1 −4 1 0 1 −4 1 0 0 0 1 0 0 0 1 Meus blocos até parecem “confeti” ! 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j
80. É claro que para resolver esse sistema linear usaremos o pacote scipy.linalg. Como a matriz desse sistema linear, com sinal trocado, é definida-positiva e de banda, poderemos utilizar cholesky_ banded( ) e cho_solve_banded( ).
81. Um programa para obter a solução x desse sistema linear pelo método de Cholesky.
82. No exemplo usaremos as seguintes condições de contorno: 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j = = = 10 ºC = = = 0. ºC = ℎ = = 0. ºC = = = 0. ºC
83. 1 4 3 5 6 7 8 9 2 a b c i h g d e f l k j = = = 100. ºC = = = 0. ºC = ℎ = = 0. ºC = = = 0. ºC A solução:
84. Vou fazer um outro exemplo para esclarecer esse ponto. Ajude aí colega! Mestra, a idéia do blocos ficou clara, só não entendi como como montá-los.
85. Considerarei agora um retângulo e farei uma discretização maior, Loirinha. Para entender como fica a matriz dos coeficientes do sistema. a b c n m r q p f g h d e l k j o i 1 3 2 4 5 6 8 7 9 10 11 13 12 14 15 16 18 17 19 20
86. a b c n m r q p f g h d e l k j o i 1 3 2 4 5 6 8 7 9 10 11 13 12 14 15 16 18 17 19 20 Aplicando o esquema em cruz para os pontos internos 1, 2, … , 20 ficaremos com o sistema de 20 equações a 20 incógnitas 1 , 2 , … , 20 . As 5 primeiras equações, indicadas abaixo, foram obtidas percorredo a 1ª linha da malha.: − 21 + 2 ℎ2 + − 21 + 6 ℎ2 = 1 1 − 22 + 3 ℎ2 + − 22 + 7 ℎ2 = 2 2 − 23 + 4 ℎ2 + − 23 + 8 ℎ2 = 3 3 − 24 + 5 ℎ2 + − 24 + 9 ℎ2 = 4 4 − 25 + ℎ2 + − 25 + 10 ℎ2 = 5
87. −41 + 2 + 6 = ℎ21 − ( + ) −42 + 1 + 3 + 7 = ℎ22 − −43 + 2 + 4 + 8 = ℎ23 − −44 + 3 + 5 + 9 = ℎ24 − −45 + 4 + 10 = ℎ25 − ( + ) − 21 + 2 ℎ2 + − 21 + 6 ℎ2 = 1 1 − 22 + 3 ℎ2 + − 22 + 7 ℎ2 = 2 2 − 23 + 4 ℎ2 + − 23 + 8 ℎ2 = 3 3 − 24 + 5 ℎ2 + − 24 + 9 ℎ2 = 4 4 − 25 + ℎ2 + − 25 + 10 ℎ2 = 5 Rearranjando essas 5 equações obtemos :
88. − 26 + 7 ℎ2 + 1 − 26 + 11 ℎ2 = 6 6 − 27 + 8 ℎ2 + 2 − 27 + 12 ℎ2 = 7 7 − 28 + 9 ℎ2 + 3 − 28 + 13 ℎ2 = 8 8 − 29 + 10 ℎ2 + 4 − 29 + 14 ℎ2 = 9 9 − 210 + ℎ2 + 5 − 210 + 15 ℎ2 = 10 −46 + 7 + 1 + 11 = ℎ26 − −47 + 6 + 8 + 2 + 12 = ℎ27 −48 + 7 + 9 + 3 + 13 = ℎ28 −49 + 8 + 10 + 4 + 14 = ℎ29 −410 + 9 + 5 + 15 = ℎ210 − A 2ª linha da malha gera estas outras 5 equações:
89. −411 + 12 + 6 + 16 = ℎ211 − −412 + 11 + 13 + 7 + 17 = ℎ212 −413 + 12 + 14 + 8 + 18 = ℎ213 −414 + 13 + 15 + 9 + 19 = ℎ214 −415 + 14 + 10 + 20 = ℎ215 − ℎ A 3ª e a 4ª linhas da malha geram as equações: −416 + 17 + 11 = ℎ216 − ( + ) −417 + 16 + 18 + 12 = ℎ217 − −418 + 17 + 19 + 13 = ℎ218 − −419 + 18 + 20 + 14 = ℎ219 − −420 + 19 + 15 = ℎ220 − ( + ) 4ª linha 3ª linha
90. −4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 −4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 −4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 −4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 1 − 4 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A matriz do sistema linear será constituída por 4 × 4 blocos de tamanho 5 × 5, como abaixo: Sim Mestra, as linhas da grade estabelecem o número de blocos (4 × 4 neste caso) e o número de pontos em cada linha a dimensão dos blocos (5 × 5 aqui).
91. O termo independente do sistema linear será constituído pela soma de dois vetores de ordem 20: ℎ2 + Um deles provém da fonte: ℎ2 = ℎ2 1 2 ⋯ 20 O outro é o vetor em frente, oriundo das condições de contorno: = − + + 0 0 0 0 0 0 ℎ + +
92. Pois é Surfista para sistemas esparsos como esse, = + , os métodos apropriados para resolver o sistema linear são os métodos iterativos lá da scipy.linalg.esparse, que veremos mais adiante. Repetindo a fala de meu colega: No caso geral de uma malha retangular com linhas e pontos internos em cada linha teremos uma matriz com blocos quadrados de tamanho , seguindo esse mesmo padrão do exemplo.
93. O programa a seguir permite que você escolha o número de linhas internas da discretização e a quantidade de pontos internos nas linhas. Retorna a solução do sistema e uma representação colorida da distribuição de calor na placa.
94. Escolhemos uma malha pequena apenas para mostrar a matriz A todinha.
95. Eis a representação colorida da distribuição de calor na placa. Experimente malhas com mais pontos, Surfista.
96. O programa: A construção da matriz A (1) e do vetor b (2). 1 2
97. A continuação, usando o método do gradiente conjugado para resolver o sistema linear esparso (3). A construção das matrizes de temperatura nos pontos da malha (4) para desenhar (5) a representação colorida das temperaturas com imshow( ). 3 4 5
98. • 2 2 + 2 2 = em Ω × (0, 0 ) • , , 0 = 0 (, ) em Ω1 • (,) = () em Ω2 • (,) = ℎ() em Ω3 • , , prescrita na borda Ω × (0, 0 ) Agora vamos resolver o problema de calor em estado estacionário numa placa Ω com o formato de um L.
99. f O método de diferenças finitas se aplica facilmente para regiões Ω como esta, constituidas por blocos retangulares contíguos com fronteira Ω = Γ = Γ ∪ Γ ∪ Γ Sim colega, é só repetir o mesmo procedimento do exemplo para cada bloco.
100. f 1 2 3 10 11 12 13 14 f1 f2 f3 f4 f5 f6 f7 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 7 8 9 4 5 6 f8 f9 f10 f11 f12 f13 f14 f15 f18 f17 f16 f22 f21 f20 f19 f29 f28 f27 f25 f24 f23 f26 Passeiem com a cruz pelos pontos azuis (os internos) para estabelecer as equações, conforme já fizemos. Ficaremos com um sistema linear constituído por: • 4 blocos 3 × 3 − 12 equações, • 3 blocos 8 × 8 − 24 equações.
101. Seu problema, Surfista é determinar a temperatura nos pontos internos, os azuis. Mais que isto, além das 36 incógnitas dos pontos internos (para as quais você tem 36 equações), acrescente outras 6 incógnitas vindas das condições de fluxo, nos pontos cinza e verdes. Mestre, como eu faço para usar a condição de fluxo prescrito na fronteira?
102. Então aproximamos a derivada usando diferença atrasada: 16 = ቉ 16 ≅ 16 − 36 ℎ , 17 = ቉ 17 ≅ 17 − 28 ℎ , 18 = ቉ 18 ≅ 18 − 20 ℎ . Bem, Loirinha, prescrever o fluxo num ponto da fronteira é prescrever o valor da derivada parcial nele. Para a fronteira cinza teremos ቃ =
103. Então teremos mais 3 equações: 16 − 36 = ℎ 16 , 17 − 18 = ℎ 17 , 16 − 20 = ℎ 16 , para as incógnitas 16 , 17 , 18 . Para a fronteira verde teremos ቃ =
104. f 1 2 3 10 11 12 13 14 f1 f2 f3 f4 f5 f6 f7 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 7 8 9 4 5 6 f8 f9 f10 f11 f12 f13 f14 f15 f18 f17 f16 f22 f21 f20 f19 f29 f28 f27 f25 f24 f23 f26 Assim, você obtém mais 3 equações para as condições de fluxo: − ቉ = ⟹ ൞ 29 − 1 = −ℎ 29 28 − 2 = −ℎ 28 27 − 3 = −ℎ 27 As direções do fluxo de calor estão indicadas por flechas, por isso o sinal + na fronteira cinza e o − na verde.
105. Atenção: • Fluxo negativo num ponto é calor entrando por esse ponto. • Fluxo positivo é calor saindo por ele. • Fluxo nulo é reflexão de calor.
106. Sobrou prá você, Surfista! Monte o sistema de equações e resolva-o.
107. Tchau, até a próxima aula!