Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Diferenças finitas e EDP's

Diferenças finitas e EDP's

Paulo Bordoni

June 06, 2019
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. 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
  2. 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.
  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, … , , com e fixos e ~40 − 50.
  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. Para = temos ′ = e ′ 1 = , o limite da sequência.
  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, começa a aparecer a propagação de erro da IEEE 754/2008.
  11. Estou mostrando o gráfico dos erros, em escala logarítmica para

    que possamos entender melhor o que ocorre.
  12. 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:
  13. 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
  14. É 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?
  15. 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.
  16. 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.
  17. 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
  18. 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 (, ) ( + ∆, ) ∆
  19. Portanto a equação de balanço de energia térmica fica: ,

    − + ∆, ∆ = = ∆ ∙ ∆(, ) u alta u baixa (, ) ( + ∆, ) ∆
  20. 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 (, ) ( + ∆, ) ∆
  21. Utilizando a Lei de Fourier, obtemos a equação condução de

    calor: (, ) = 2(, ) 2 na qual resume as constantes físicas. u(, )
  22. 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 + + +
  23. 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 ()
  24. 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.
  25. 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
  26. • 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.
  27. 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 , = ℎ + Τ − ℎ ( − ) .
  28. 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.
  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 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 = ∆.
  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 = (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 )
  31. 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 + ,
  32. 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, ,
  33. 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 ∆
  34. 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
  35. 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
  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 No método explícito, você calcula diretamente o valor de +1, em função dos valores conhecidos de ,−1 , , , ,+1 . +1, ,−1 , , , ,+1
  37. 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
  38. 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, … ,
  39. 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.
  40. Notem que então podemos calcular 1 , 2 , 3

    , … pela fórmula +1 = + para = 0, … , com 0 dado.
  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 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 ,
  42. ,+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):
  43. 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, … ,
  44. = 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
  45. 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:
  46. 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.
  47. Trata-se de um sistema tri-diagonal, que pode ser resolvido usando-se

    métodos iterativos, para valores muito grandes de .
  48. 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.
  49. = 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
  50. 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
  51. ,+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:
  52. • = 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 ).
  53. 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
  54. 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
  55. 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
  56. 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
  57. −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
  58. −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
  59. 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 .
  60. 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
  61. É 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( ).
  62. 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
  63. 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:
  64. 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.
  65. 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
  66. 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
  67. −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 :
  68. − 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:
  69. −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
  70. −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).
  71. 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 ℎ + +
  72. 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.
  73. 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.
  74. Eis a representação colorida da distribuição de calor na placa.

    Experimente malhas com mais pontos, Surfista.
  75. 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
  76. • 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.
  77. 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.
  78. 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.
  79. 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?
  80. 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 ቃ =
  81. 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 ቃ =
  82. 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.
  83. Atenção: • Fluxo negativo num ponto é calor entrando por

    esse ponto. • Fluxo positivo é calor saindo por ele. • Fluxo nulo é reflexão de calor.