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.

Transcript

1. LNCC Fatoração LU e aplicações Prof. Paulo R. G. Bordoni UFRJ
2. 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 + ⋯ + =
3. | = 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.
4. 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.
5. 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!
6. 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)
7. Começaremos pelas trocas, Loirinha!
8. 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.
9. 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 ⟷ .
10. = 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.
11. 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.
12. 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.
13. 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
14. 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.
15. Meu programa também permite conferir que −1 = .
16. 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?
17. 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.
18. O programa a seguir mostra isto:
19. 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
20. 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:
21. Um programa que mostra o produto da 1 = − 1 por A.
22. 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!
23. O programa que fiz para conferir a afirmação da Mestra!
24. = 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.
25. −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
26. Eu fiz este outro programa para conferir essa 2ª afirmação Mestra!
27. Vejam o resultado!
28. Agora vamos usar os superpoderes das NumPy e SciPy, para dirigir o magnífico calhambeque da capa pelo caminho das pedras.
29. Precisaremos obter fatias de ndarrays e saber operar com elas. O operador : é a ferramenta que a numpy oferece para cortar o pão!
30. 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
31. Nossa 1ª fatia será a 2ª linha da matriz 4x5: 0 1 2 3 0 1 2 3 5
32. 0 1 2 3 0 1 2 3 Agora fatiaremos uma coluna da matriz 4x4:
33. Fatiando uma coluna e multiplicando-a por um número.
34. Nosso 3º fatiamento será a sub-matriz formada pelas 3 últimas linhas da matriz: 0 1 2 3 0 1 2 3
35. 0 1 2 3 0 1 2 3 Nosso 4ª fatiamento será a sub-matriz formada pelas 2 primeiras colunas da matriz A:
36. 0 1 2 3 0 1 2 3 5 Agora fatiaremos as 2ª, 3ª e 4ª linhas de uma matriz 4x5:
37. 0 1 2 3 0 1 2 3 5 Agora fatiaremos a sub matriz indicada na figura:
38. Vou mostrar a ação NumPy de um Samurai. Amei este filme!
39. Recebendo uma matriz, um vetor e devolvendo o produto termo a termo de parte deles.
40. 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.
41. O impressionante é que tudo acontece um golpe só, como faz um samurai!
42. Repetiremos o método de eliminação, de Gauss, usando os superpoderes da NumPy e da SciPy.
43. O processo de eliminação é o mesmo. O que o mudará é a utilização de fatiamento, vetorização e difusão.
44. Vejam como tudo fica muito mais simples. Ainda assim, vou explicar cada detalhe do código.
45. O loop em k é a única estrutura de repetição utilizada. Como antes, ele é responsável pelas etapas do processo de eliminação.
46. 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.
47. Nesta linha de código, calculamos todos os multiplicadores da coluna k, abaixo do pivot. Num golpe só, como o samurai!
48. Esta é outra linha de código vetorizada. Ela atualiza, de uma vez só, toda a sub matriz que pintei de azul.
49. 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ó.
50. Na substituição retroativa, também temos vetorização e difusão.
51. Vamos examinar uma aplicação às EDP’s da resolução de sistemas lineares. EDP’s Equações Diferenciais Parciais
52. 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
53. 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:
54. 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
55. Soma Diferença Somando a fórmula com a obtemos: Avç. Atr. Já subtraindo a fórmula da obtemos: Avç. Atr.
56. 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:
57. 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:
58. Adoro festa de São João!
59. É 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?
60. 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.
61. 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
62. 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
63. 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 (, ) ( + ∆, ) ∆
64. Portanto a equação de balanço de energia térmica fica: , − + ∆, ∆ = = ∆ ∙ ∆(, ) T alta T baixa (, ) ( + ∆, ) ∆
65. 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 (, ) ( + ∆, ) ∆
66. Utilizando a Lei de Fourier, obtemos a equação condução de calor: (, ) = 2(, ) 2 (, )
67. 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.
68. 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 :
69. Agora no quadro negro!
70. 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?
71. = 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!
72. 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?
73. 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:
74. 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!
75. 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
76. 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.
77. 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
78. 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 ...
79. Simplificando obtemos o sistema linear abaixo, que é penta-diagonal.
80. É 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( ).
81. Um programa para obter a solução x desse sistema linear.
82. 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.
83. 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.
84. Dizem que, “Simpatia é quase amor”. Meu bloco preferido!
85. Meus blocos até parecem “confeti” !
86. = 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.
87. 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:
88. 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!
89. 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:
90. 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:
91. Concordo Surfista, mas assim você poderá resolver sistemas gigantes, através da solução de sistemas bem menores! Mestra, é um bocado de trabalho!
92. Então Surfista, faça isto para aquele de 2.500 equações à 2.500 incógnitas! É uma dica Mestre ?
93. Tchau, até a próxima aula!