Da eliminação de Gauss à fatoração LU

Veja esta apresentação em www.bordoni.info.

Transcript

1. Da eliminação de Gauss à fatoração LU Prof. Paulo R. G. Bordoni LNCC UFRJ
2. Jovens, precisamos reduzir a velocidade para evitar acidentes!
3. Apontarei claramente, agora, a diferença entre um curso de Cálculo Numérico e algo tipo: “Me engana que eu gosto, você finge que ensina e eu que aprendo”.
4. Vamos examinar um caso limite nos dados do programa visto no final da aula passada. Efetuaremos uma modificação mínima no código: apenas na forma de testar se o pivot é não-nulo.
5. A B Assinalei a única diferença.
6. A modificação é realmente mínima. Apenas arredondamos o valor do pivot para 14 casas decimais após a vírgula. Entretanto, Surfistas, observem a diferença no retorno da execução de cada código, na próxima transparência!
7. Para os dados marcados em azul, o código em A oferece a solução abaixo! A
8. Para os mesmos dados o código B oferece a resposta correta: o sistema é impossível. B
9. Observe, Surfista, que a 3ª linha da matriz M, na ETAPA 2, é nula. Em particular o pivot, ali na posição (2,2), aparece como zero - apenas, e tão somente, porque o código arredonda seus elementos para 2 casas decimais. A
10. A Na verdade, para os dados em azul, o pivot é minúsculo: menor que . . .5*10**(-14) Repetindo: Entretanto é diferente de zero, o que é testado e aceito como True no if da linha 29 do código! pivot < 0.000000000000005
11. E produzir, em seguida, o erro marcado aos meus pés. B Erro evitado pelas linhas marcadas do código B.
12. Mestre, eu nunca teria posto esse detalhe em meu programa! Fiquei assustada. Como eu faço? Shht! Acabei de perder totalmente a confiança no meu taco!
13. Pois é jovens. O momento de aprender é este. Não permitam que os enganem! Imaginem se a consequência do erro tivesse sido esta! Há mais de 40 anos o Rio entrava em choque com o desabamento do Elevado Paulo de Frontin, no Rio Comprido, Zona Norte do Rio. A tragédia causou a morte de 28 pessoas, deixou 30 feridas, e destruiu 17 carros. Foto: CPDoc JB
14. O que causou o desastre no Paulo de Frontin eu não sei, mas o erro apresentado no código A deve-se ao desconhecimento da representação de números no computador, estabelecida pelo padrão IEEE 754 de ponto flutuante, internacionalmente adotado para computação numérica.
15. Na verdade, a forma mais indicada para evitar esse tipo de desastre é usar a função isclose() da NumPy.
16. Fui buscar informações sobre ela em Logic Functions, na refgerência da NumPy.
17. B Nossa proposta é substituir as linhas 29 e 30 do código B, que estão corretas, por algo mais pythônico, usando isclose().
18. Nele compararemos uma variável nomeada pivot com zero, via: np.isclose(pivot,0.) Faremos um programa para exemplificar a utilização de isclose() .
19. • Iniciar atribuindo um valor ao pivot, • Definir a medida de proximidade atribuindo um valor à prox, • Depois reduzir pivot sucessiva/. por um fator de 1/10 e em cada passo: • testar se |pivot-0.|< prox através da isclose(pivot,0.), • devolver uma mensagem com o resultado da comparação. Uma descrição (algoritmo) do programa:
20. Nesta execução definimos os valores 0.002 e 3.0e-8 para pivot e prox.
21. Este é o código do programa. Note que prox envolve tanto a tolerância absoluta, atol, como a relativa, rtol.
22. Agora vamos repetir o programa do início destas transparências. Nesta 1ª parte apenas as mensagens mudam.
23. A única diferença do programa do início é a utilização de isclose.
24. A execução com o mesmo sistema linear que gerou problemas: Perfeito, Mestres!
25. Chega de andar a 40 km/h Agora poderemos retomar nossa velocidade normal!
26. Vamos retornar ao método da eliminação de Gauss e focar todos os holofotes no pivot. Pivot ? O que é pivot? Até agora eles são os divisores o divisor da kª etapa.
27. 1 /11 31 /11 Multiplicadores 21 /11 Assumam que o valor do pivot 11 é 11 = 0,001 = ൗ 1 1000 Portanto ao dividirmos elementos 21 , 31 , … , 1 por 11 estaremos, de fato, multiplicando-os por 1.000.
28. Que horror Mestra. Isto faz propagar muito o erro absoluto. Nossos estudos sobre número de condicionamento mostraram isso! Sim minha Loirinha querida! Faça como eu, despreze pivots pequenos, eles estragam o jogo, no final. Sempre vá em busca dos maiores possíveis.
29. Uma forma de minimizar a propagação do erro absoluto consiste em localizar o 1 de maior valor absoluto e trocar a linha i onde ele está com a linha 1 (na 1ª etapa). Depois repetir a ideia para as etapas seguintes nas quais o sistema restante sempre fica com 1 equação e 1 incógnita a menos.
30. É verdade Mestra, pois então os multiplicadores ficam na faixa de -1 à +1 e não ampliam os erros! Exatamente Loirinha querida!
31. O valor máximo de um array pode ser obtido com a amax:
32. É verdade Manuel, eu mesmo já usei a np.amax() várias vezes. Mas como fazemos para determinar a posição do valor máximo no array?
33. Usem esta função para buscar a posição no array do seu valor máximo
34. Existe um método da classe matrix que faz praticamente a mesma coisa. A diferença está marcada lá embaixo.
35. Vejam elas funcionando:
36. Mestre, eis como trocar duas linhas de uma matriz:
37. Agora apresentaremos o código do método de eliminação de Gauss com pivotamento parcial por coluna. Precisamos localizar a posição do pivot na coluna, em seguida realizar a troca das linhas, para depois efetuar a eliminação.
38. O código só muda onde marquei. Em: 1. -> Busco o elemento de módulo máximo da coluna. 2. -> Faço a troca das linhas. 1 2
39. Este é um 1º exemplo de execução. Surfista, faça outros, semelhantes aos exemplos anteriores. 1 2 3
40. Mestre, como posso saber quais foram as trocas de linhas efetuadas? Muito simples, Loirinha. Basta criar um vetor de inteiros, trocas = [0,1, … , − 1] e anotar as trocas nele. Veja o código. O exemplo está na próxima transparência.
41. Valeu, Mestre! 1 3 2
42. Voltaremos, oportunamente, à resolução de sistemas lineares.
43. Nas próximas 31 transparências mostraremos que o método de eliminação de Gauss pode ser descrito com operações matriciais. Teremos então a fatoração PLU.
44. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = | = 11 12 ⋯ 1 ⋮ 1 21 22 ⋯ 2 ⋮ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ ⋮ ⋮ ⋮ Lembrem-se, escrevemos todos os dados do sistema linear na matriz aumentada:
45. | = 11 12 ⋯ 1 ⋮ 1 22 ⋯ 2 ⋮ 2 ⋱ ⋮ ⋮ ⋮ ⋮ | = 11 12 ⋯ 1 ⋮ 1 21 22 ⋯ 2 ⋮ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ ⋮ ⋮ ⋮ Depois efetuamos operações elementares sobre as linhas da matriz aumentada | até transformar A numa matriz triangular superior U.
46. 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.
47. • P é uma matriz de permutações (p/. trocas de linhas) • L é uma matriz triangular inferior (L de lower) • U é a matriz triangular superior ao final da eliminação (U de upper) A fatoração PLU é o resultado final dessas multiplicações. Como o Mestre anunciou todas essas três operações elementares podem ser descritas através de multiplicações matriciais!
48. Começaremos pelas trocas, Loirinha!
49. 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.
50. 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 ⟷ .
51. = 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.
52. Este programa faz o que marquei aqui no início.
53. Ele permite conferir que o produto de permutações é comutativo. Quero lembrar que, em geral, o produto matricial não é comutativo.
54. Ele permite mostrar, também, que a inversa da permutação ↔ , é a permutação ↔ .
55. Ele permite, ainda, mostrar que as permutações são idempotentes: 2 = . Observe, Surfista, que neste exemplo escolhi = .
56. Surfista, você já reparou que as matrizes de permutação são matrizes ortogonais? É verdade, Mestre! As colunas de uma permutação P são a base canônica de ℝ 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
57. 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.
58. As operações elementares no processo de eliminação foram realizadas 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
59. 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:
60. Um programa que mostra o produto da 1 = − 1 por A.
61. 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 é forte! Eu tenho que conferir!
62. O programa que fiz para conferir a afirmação da Mestra!
63. = 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 0 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.
64. −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
65. Eu fiz este programa para efetuar a eliminação, passo a passo!
66. Vejam o resultado! A matriz A é atualizada em cada etapa.
67. Na realidade, ao programar, NÃO efetuamos multiplicações matriciais. É o que mostraremos nos dois programas a seguir. Então por que motivo mostraram tudo isto? Medite um pouco sobre tua pergunta, jovem Surfista imediatista!
68. Primeiro um programa para a fatoração LU sem pivotamento.
69. Ele mostra, passo a passo, as matrizes Mult, M e b. Ao final a solução x e o produto LU. Em cada etapa as matrizes Mult e M e o vetor b são atualizados!
70. As mudanças no código da eliminação de Gauss com pivotamento, para mostrar a fatoração PLU. Comparem com o código anterior.
71. Observe, Surfista, que o produto PLU já traz as trocas de linhas incluídas. Em cada etapa as matrizes M e L e o vetor b são atualizados!
72. Quando temos a fatoração = , resolver um sistema linear = é equivalente a resolver o sistema = , usando apenas a fatoração LU. Então Mestra, basta efetuar as trocas das linhas no termo independente b e usar a fatoração LU sem pivotamento na matriz A já com as trocas?
73. Sim, ao invés de resolver 1 1 −2 0 1 −5 2 −1 3 1 3 1 1 3 −2 2 0 1 2 3 = 3 −2 0 0 Resolvemos 3 3 1 −2 1 −5 2 −1 1 1 1 1 3 −2 2 0 0 1 2 3 = 0 −2 0 3
74. Foi o que eu fiz aqui. Usei a matriz M já com as linhas trocadas e troquei, também, as linhas do termo independente b. Depois resolvi o sistema sem pivotamento.
75. Tchau, até a próxima aula!