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

Fatoração LU e aplicações

Fatoração LU e aplicações

Matrizes elementares, fatoração LU, fatiamento com NUmPy e método de eliminção, de Gauss com programação vetorial. Aplicação à propagação de calor numa barra e numa placa. Fatoração LU de matrizes tridiagonais por blocos.

Paulo Bordoni

July 05, 2013
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. Ao final da aula passada Estudamos o método da eliminação,

    de Gauss, para obtenção da solução = 1 2 ⋯ de um sistema de equações lineares: 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + =
  2. | = 11 12 ⋯ 1 ⋮ 1 22 ⋯

    2 ⋮ 2 ⋱ ⋮ ⋮ ⋮ ⋮ | = 11 12 ⋯ 1 ⋮ 1 21 22 ⋯ 2 ⋮ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ ⋮ ⋮ ⋮ Efetuamos operações elementares sobre as linhas da matriz aumentada | até transformar A numa matriz triangular superior U.
  3. As operações elementares usadas foram: Somar ou subtrair linhas Multiplicar

    uma linha por um número real Trocar duas linhas de posição E elas permitiram transformar o sistema original A x = b num sistema triangular superior U x = β, com a mesma solução x.
  4. A fatoração LU é o resultado final dessas multiplicações. Todas

    essas três operações elementares podem ser descritas através de multiplicações matriciais!
  5. Assim como decompomos um número em seus fatores primos: 30

    = 2 ∙ 3 ∙ 5, 21 = 3 ∙ 7, 4 = 2 ∙ 2 Veremos que é possível decompor uma matriz A num produto = , onde: • P é uma matriz de permutações (p/. trocas de linhas) • L é uma matriz triangular inferior (L de lower) • U é uma matriz triangular superior (U de upper)
  6. Indicaremos por ↔ a matriz obtida permutando as linhas i

    e j da identidade I Observem que ↔ é uma matriz simétrica, isto é ↔ = ⟷ . É fácil observar que o produto ↔ é a matriz A com suas linhas i e j trocadas. O programa a seguir mostrar isto.
  7. Neste programa você entra com uma matriz A pelo teclado

    e com os índices i e j das linhas de A a serem trocadas. Ele retorna a matriz ⟷ e o produto ⟷ .
  8. = 0 1 0 0 0 1 0 0 0

    0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 Neste exemplo, mostro a identidade I com as linhas 1 e 2 e as 3 e 5 trocadas. Aliás, se P é uma matriz obtida da identidade I através de permutações de linhas, então o produto P A é a matriz A com as mesmas linhas trocadas.
  9. Este programa faz o que está descrito no seu início.

    Ele também mostra que o produto de permutações também é uma permutação.
  10. Mestre, corrija o erro que assinalei, não é QPA, mas

    sim PQA ! He, he, he! Eu fiz como pegadinha, só para lembrar que = , uma consequência de serem matrizes de permutação.
  11. Surfista, você já reparou que as matrizes de permutação são

    matrizes ortogonais? É verdade, Mestre! No seu exemplo, as colunas de PQ são a base canônica de ℝ5 fora de ordem. = 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
  12. Sim coleguinha. E da ortogonalidade segue que, se P é

    uma matriz de permutações então = −1. E como elas são simétricas temos −1 = = . Aliás isto é totalmente intuitivo: a inversa −1 destroca as linhas da identidade I que foram trocadas por P.
  13. Com tudo que foi provado, podemos afirmar que o conjunto

    das matrizes de permutação, munido da multiplicação matricial possui a estrutura de grupo. Será um grupo abeliano?
  14. A multiplicação da linha k de uma matriz A por

    um número α também é descrita por um produto matricial. A matriz M obtida a partir da identidade I substituindo o 1 da posição (i, i) por α faz isto.
  15. Vocês devem estar lembrados que usamos as operações elementares para

    zerar os elementos abaixo do pivot na 1ª coluna. Isto foi realizado da seguinte forma: Linha 1 Linha 1 Linha 2 Linha 2 Linha 1 21 /11 × − Linha 3 Linha 3 Linha 1 31 /11 × − Linha n Linha n Linha 1 1 /11 × − Multiplicadores
  16. 1 = − 1 = 1 0 ⋯ 0 0

    1 ⋮ ⋮ 0 ⋯ ⋱ 0 0 1 − 0 0 ⋯ 0 21 0 ⋮ ⋮ 1 ⋮ 0 ⋱ 0 ⋯ 0 = 1 0 ⋯ 0 −21 1 ⋯ 0 ⋮ −1 0 ⋱ ⋮ ⋯ 1 Pois é, isto corresponde a efetuar o produto 1 ∙ , Onde 1 = − 1 é a matriz de eliminação:
  17. E, o mais surpreendente, é que a inversa dessa matriz

    de eliminação 1 = − 1 é a matriz + 1 , obtida simplesmente trocando o sinal da 1 : 1 −1 = ( − 1 )−1 = 1 0 ⋯ 0 21 1 ⋯ 0 ⋮ 1 ⋮ 0 ⋱ ⋮ ⋯ 1 = + 1 Essa eu tenho que conferir!
  18. = 3 2 1 = 1 0 0 0 0

    1 0 0 0 0 0 0 1 0 9 1 1 0 0 0 0 1 0 ⋮ 0 0 7 4 1 0 0 1 1 0 0 0 3 1 0 ⋮ −2 5 0 0 1 0 0 1 = = 1 0 0 0 3 1 1 0 19 188 7 67 1 0 9 1 No método de eliminação usamos um produto, de matrizes eliminadoras = 3 2 1 Como no exemplo abaixo.
  19. −1 = 1 0 0 0 −3 1 0 ⋮

    2 −5 0 0 1 0 0 1 1 0 0 0 0 1 0 ⋮ 0 0 −7 −4 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 9 1 = = 1 0 0 0 −3 1 1 0 2 −5 −7 −4 1 0 −9 1 O mais surpreendente é a inversa −1 = (3 2 1 )−1= 1 −12 −1 3 −1 desse produto E. Simplesmente é a superposição das com sinal trocado
  20. Agora vamos usar os superpoderes das NumPy e SciPy, para

    dirigir o magnífico calhambeque da capa pelo caminho das pedras.
  21. Precisaremos obter fatias de ndarrays e saber operar com elas.

    O operador : é a ferramenta que a numpy oferece para cortar o pão!
  22. Representarei uma matriz por uma grade, como abaixo, porque só

    estou interessada nos índices. Isto basta para explicar o conceito de fatiamento. 0 1 2 3 0 1 2 3 5
  23. 0 1 2 3 0 1 2 3 Agora fatiaremos

    uma coluna da matriz 4x4:
  24. 0 1 2 3 0 1 2 3 Nosso 4ª

    fatiamento será a sub-matriz formada pelas 2 primeiras colunas da matriz A:
  25. 0 1 2 3 0 1 2 3 5 Agora

    fatiaremos as 2ª, 3ª e 4ª linhas de uma matriz 4x5:
  26. 0 1 2 3 0 1 2 3 5 Agora

    fatiaremos a sub matriz indicada na figura:
  27. Se M é uma matriz, o comando M . A

    é o ndarray associado à ela (o A é de Array). Esta conversão é necessária para efetuar o produto elemento a elemento.
  28. O processo de eliminação é o mesmo. O que o

    mudará é a utilização de fatiamento, vetorização e difusão.
  29. O loop em k é a única estrutura de repetição

    utilizada. Como antes, ele é responsável pelas etapas do processo de eliminação.
  30. Definimos A[k,k] como o pivot da etapa k. Se ele

    não é zero, criamos a matriz coluna Mult que guardará os multiplicadores dessa etapa.
  31. Nesta linha de código, calculamos todos os multiplicadores da coluna

    k, abaixo do pivot. Num golpe só, como o samurai!
  32. Esta é outra linha de código vetorizada. Ela atualiza, de

    uma vez só, toda a sub matriz que pintei de azul.
  33. O código da 1ª linha acima guarda os multiplicadores na

    coluna k da matriz A, abaixo do pivot. Num golpe só. A 2ª aplica as operações elementares nos elementos do vetor b, após a posição correspondente ao pivot. Também num golpe só.
  34. Vamos examinar uma aplicação às EDP’s da resolução de sistemas

    lineares. EDP’s Equações Diferenciais Parciais
  35. x x+h h h x-h h > 0 Avç. Atr.

    O desenvolvimento em série de Taylor de uma função f, x ⟼ f(x) entorno de um ponto x
  36. Esta é a aproximação por diferença avançada para a derivada

    1ª de f Da fórmula , obtemos a fórmula abaixo. Avç. Nela, passe f (x) para o 1º membro e divida por -h . Assumindo que f ” é limitada, obtemos a expressão mostrada pela Professora. O algebrismo é o seguinte:
  37. Da fórmula , obtemos a fórmula abaixo. Atr. Nela, passe

    f (x) para o 1º membro e divida por -h . Assumindo que f ” é limitada, obtemos a expressão mostrada pela Professora. O algebrismo é o seguinte: Esta é a aproximação por diferença atrasada para a derivada 1ª de f
  38. Soma Diferença Somando a fórmula com a obtemos: Avç. Atr.

    Já subtraindo a fórmula da obtemos: Avç. Atr.
  39. Divida a fórmula diferença por 2h. Assumindo que a derivada

    3ª é limitada, obtemos a expressão mostrada pelo Professor. Diferença A aproximação por diferença centrada para a derivada 1ª de f. Ela é obtida assim:
  40. Na fórmula “Soma” subtraia 2f (x) de ambos os lados

    e divida tudo por h 2 . Assumindo que a derivada 4ª é limitada, obtemos a expressão mostrada pelo Professor. Soma A aproximação por diferença centrada para a derivada 2ª de f. Ela é obtida assim:
  41. É 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?
  42. 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.
  43. 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
  44. 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 calorífica do material. T alta (, ) T baixa
  45. 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 (, ) ( + ∆, ) ∆
  46. Portanto a equação de balanço de energia térmica fica: ,

    − + ∆, ∆ = = ∆ ∙ ∆(, ) T alta T baixa (, ) ( + ∆, ) ∆
  47. 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 (, ) ( + ∆, ) ∆
  48. 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.
  49. x 0 x 3 x 4 x 5 x 1

    x 2 t 0 t 3 t 4 T t 1 t 2 Discretizando em x e t :
  50. 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?
  51. = 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!
  52. 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?
  53. 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:
  54. 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!
  55. 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
  56. 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.
  57. 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
  58. 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 ...
  59. É 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( ).
  60. 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.
  61. Pois é Surfista para um sistema esparso como esse, os

    métodos apropriados para resolver o sistema linear são os métodos iterativos. Mais adiante, no curso, estudaremos métodos iterativos para resolução de sistemas lineares. Mas agora, neste Carnaval não, vou sair no meu bloco em Ipanema.
  62. = 1 1 2 2 ⋯ 0 ⋮ ⋱ ⋮

    0 ⋯ −1 −1 = 2 ⋯ ⋮ ⋱ ⋮ ⋯ 1 1 2 ⋯ ⋮ ⋱ ⋮ ⋯ −1 = Entrando nessa folia. Vou apresentar uma fatoração T = L U por blocos para uma matriz tridiagonal por blocos T.
  63. Input: A matriz tri-diagonal por blocos T . Output: A

    matriz bi-didiagonal inferior por blocos L com I’s na diagonal e a matriz triangular superior bi-diagonal U. • Passo 1: Defina 1 = 1 • Passo 2: Para = 2, ⋯ , : • 2.1: Calcule resolvendo os sistemas −1 = . • 2.1: Calcule efetuando = − −1 . • Fim O algoritmo para efetuar a fatoração LU por blocos de T é o seguinte:
  64. Da mesma forma que aplicamos a fatoração LU sem blocos

    para calcular a solução x de = , iremos calcular a solução de = . Eu lembro, Mestra: Substituindo T por LU no sistema, ficamos com = . Então resolvemos = e, em seguida = . Mas agora com blocos!
  65. Input: A matriz triangular inferior por blocos = [ ]

    e o vetor por blocos = [ ]. Output: O vetor por blocos = [ ], solução dos sistemas lineares = . • Passo 1: Resolva 1 0 = 0. • Passo 2: Para = 1, ⋯ , repita: • = − −1 . • Fim O algoritmo para efetuar a substituição na ordem direta é o seguinte:
  66. Input: A matriz triangular superior por blocos = [ ]

    e o vetor por blocos = [ ], obtido da substituição direta. Output: O vetor por blocos = [ ], solução dos sistemas lineares = . • Passo 1: Resolva +1 = 0. • Passo 2: Para = , ⋯ , 1 repita: • = − +1 . • Fim O algoritmo para efetuar a retro-substituição é o seguinte:
  67. Concordo Surfista, mas assim você poderá resolver sistemas gigantes, através

    da solução de sistemas bem menores! Mestra, é um bocado de trabalho!
  68. Então Surfista, faça isto para aquele de 2.500 equações à

    2.500 incógnitas! É uma dica Mestre ?