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

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

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

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

Paulo Bordoni

October 05, 2016
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. 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”.
  2. 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.
  3. 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!
  4. Para os dados marcados em azul, o código em A

    oferece a solução abaixo! A
  5. Para os mesmos dados o código B oferece a resposta

    correta: o sistema é impossível. B
  6. 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
  7. 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
  8. E produzir, em seguida, o erro marcado aos meus pés.

    B Erro evitado pelas linhas marcadas do código B.
  9. 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!
  10. 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
  11. 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.
  12. Na verdade, a forma mais indicada para evitar esse tipo

    de desastre é usar a função isclose() da NumPy.
  13. B Nossa proposta é substituir as linhas 29 e 30

    do código B, que estão corretas, por algo mais pythônico, usando isclose().
  14. Nele compararemos uma variável nomeada pivot com zero, via: np.isclose(pivot,0.)

    Faremos um programa para exemplificar a utilização de isclose() .
  15. • 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:
  16. Este é o código do programa. Note que prox envolve

    tanto a tolerância absoluta, atol, como a relativa, rtol.
  17. 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.
  18. 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.
  19. 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.
  20. 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.
  21. É verdade Mestra, pois então os multiplicadores ficam na faixa

    de -1 à +1 e não ampliam os erros! Exatamente Loirinha querida!
  22. É 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?
  23. Existe um método da classe matrix que faz praticamente a

    mesma coisa. A diferença está marcada lá embaixo.
  24. 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.
  25. 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
  26. Este é um 1º exemplo de execução. Surfista, faça outros,

    semelhantes aos exemplos anteriores. 1 2 3
  27. 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.
  28. 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.
  29. 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:
  30. | = 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.
  31. 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.
  32. • 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!
  33. 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.
  34. 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 ⟷ .
  35. = 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.
  36. Ele permite conferir que o produto de permutações é comutativo.

    Quero lembrar que, em geral, o produto matricial não é comutativo.
  37. Ele permite, ainda, mostrar que as permutações são idempotentes: 2

    = . Observe, Surfista, que neste exemplo escolhi = .
  38. 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
  39. 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.
  40. 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
  41. 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:
  42. 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!
  43. = 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.
  44. −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
  45. 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!
  46. 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!
  47. 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.
  48. 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!
  49. 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?
  50. 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
  51. 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.