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

Método de eliminação de Gauss e a scipy.linalg

Método de eliminação de Gauss e a scipy.linalg

...

Paulo Bordoni

May 19, 2018
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. Na NumPy, a faca para fatiar ndarrays é o operador

    : . É a ferramenta que ela oferece para cortar o pão nosso de cada dia !
  2. O fatiamento é simplesmente uma forma que a NumPy oferece

    para construir submatrizes de uma matriz dada. Inclusive para pegar linhas ou colunas da matriz.
  3. Representarei uma matriz por uma grade, como abaixo. Por que

    só estou interessada nos índices. Isto basta para explicar o conceito de fatiamento. 0 1 2 3 0 1 2 3 5 Lembrem-se, na NumPy, a indexação de arrays começa por zero e não por um
  4. 0 1 2 3 0 1 2 3 5 Se

    A é uma matriz m x n, então: • [ , : ] fornece a linha i de A; • [ ∶ , ] fornece a coluna j de A. Esta notação é idêntica à do Matlab. [1, : ] 0 1 2 3 0 1 2 3 [ ∶, 2]
  5. 0 1 2 3 0 1 2 Se A é

    uma matriz m x n, então: • [ , : ] fornece a linha i de A, a partir de j; • [ : , ] fornece a coluna j de A, a partir de i. [1,2: ] 0 1 2 3 0 1 2 3 [ 1: , 0]
  6. No fatiamento devemos pensar em intervalos fechados à esquerda, mas

    abertos à direita: • [ , : ] fornece a linha i de A, iniciando na coluna j e acabando na coluna − 1. • Pense em : como o intervalo , . 0 1 2 3 0 1 2 3 4 5 [1,1: 4] : k `` = ´´ [ j, k )
  7. Deixe de bobagens, Mestre. Diga logo ao Surfista que podemos

    multiplicar fatias por números da mesma forma que multiplicamos matrizes e vetores por números.
  8. Agora vamos iniciar nosso estudo do método de eliminação de

    Gauss, um algoritmo para resolver sistemas de equações lineares.
  9. Carl Friedrich Gauss 30/04/1777 – 23/02/1855 Na verdade ele não

    foi criado por Gauss. Ele já aparece no cap. 8 do livro chinês “Os nove capítulos sobre a arte matemática”, composto por gerações de mestres ao longo dos séculos X a II A.C.
  10. As honras da Alemanha à K. F. Gauss, por muitas

    outras contribuições à Matemática. E as nossas também, Mestre!
  11. O algoritmo possui duas etapas. A 1ª consiste em transformar

    a matriz do sistema numa matriz triangular superior. É o processo de eliminação. E a 2ª etapa consiste em resolver esse sistema triangular superior, algo muito simples.
  12. Resolver sistemas triangulares é muito fácil. Mostre como, Surfista. Sim,

    Surfista, resolva o “sisteminha” 4x4 aos meus pés, cuja matriz é triangular inferior. 2 0 0 0 1 3 0 0 1 3 4 2 2 0 1 −1 1 3 3 4 = 2 −2 1 0
  13. 11 1 = 1 21 1 31 1 41 1

    +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 Mestra, vou pensar num caso mais geral, depois particularizo !
  14. 11 1 = 1 21 1 31 1 41 1

    +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 Mestra, da 1ª equação eu “tiro” o valor de 1 . Em seguida, substituo esse valor na 2ª equação e calculo o valor de 2 . Depois substituo esses dois valores na 3ª equação e “tiro“ 3 . Repito a ideia para a última equação (e outras, caso existam).
  15. 11 1 = 1 21 1 31 1 41 1

    +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 O processo de substituição, com fórmulas : 1ª eq. ⇒ 1 = 1 /11 2ª eq. ⇒ 2 = (2 − 21 1 )/22 3ª eq. ⇒ 3 = 3 − (31 1 + 32 2 ) /33 4ª eq. ⇒ 4 = 4 − (41 1 + 42 2 + 43 3 ) /44 claro, assumindo ≠ 0 para = 1,2,3,4.
  16. Surfista, observe no seu código que o número de: •

    adições cresce de 1 a − 1, • multiplicações/divisões cresce de 1 a n . Ambos como PAs de razão 1. Logo, o total de adições é ( − 1)/2 e o e multiplicações é ( + 1)/2. Assim total de operações é 2.
  17. 11 1 + 12 2 + ⋯ + 1 =

    1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Descreveremos o processo, primeiramente, em linhas gerais, usando as equações.
  18. Em cada etapa utilizaremos operações elementares sobre as linhas. Elas

    são três: Somar ou subtrair linhas Multiplicar uma linha por um número real Trocar duas linhas de posição As operações elementares não alteram a solução do sistema.
  19. 11 1 + 12 2 + ⋯ + 1 =

    1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Na 1ª etapa, eliminaremos os termos envolvendo 1 nas equações após a 1ª, efetuando operações elementares nas linhas. 11 1 + 12 2 + ⋯ + 1 = 1 ෤ 22 2 + ⋯ + ෤ 2 = ෨ 2 ⋯ ෤ 2 2 + ⋯ + ෤ = ෨ Os coeficientes e os termos independentes mudam após essas operações – por isto o ~ .
  20. 11 1 + 12 2 + ⋯ + 1 =

    1 ෤ 22 2 + ⋯ + ෤ 2 = ෨ 2 ⋯ ෤ 2 2 + ⋯ + ෤ = ෨ Repetimos o processo para o subsistema sem a 1ª linha: Um sistema com n-1 equações a n-1 incógnitas.
  21. Repetindo n-1 vezes o processo chegamos a um sistema triangular

    superior. 11 1 + 12 2 + ⋯ + 1 = 1 = 22 2 + ⋯ + 2 = 2 ⋯ Todos os coeficientes e termos independentes são diferentes dos originais, exceto, talvez, na 1ª linha.
  22. A parte trabalhosa do processo consiste, portanto, em eliminarmos incógnitas

    nas equações. Daí o nome: processo de eliminação. Os textos de Álgebra linear falam em escalonamento.
  23. Sim! Como é feita, de fato, a eliminação? Mestres, vocês

    descreveram o processo em linhas gerais! Só papo...
  24. 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1

    ⋮ 2 ⋱ ⋮ ⋯ 1 2 ⋮ = 1 2 ⋮ = Pois é... Agora vamos refazer tudo usando a forma matricial do sistema linear e a NumPy:
  25. Acrescentando o vetor b à última coluna da matriz A,

    obtemos uma matriz referida na literatura como “matriz aumentada”. | = 11 12 ⋯ 1 1 21 22 ⋯ 2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯ A matriz aumentada contém todos os dados do problema.
  26. | = 11 12 ⋯ 1 1 21 22 ⋯

    2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯ Como veremos, o processo de eliminação consiste em efetuar operações elementares sobre as linhas da matriz aumentada. Linha 1 Linha 1 Linha n Ahá! Por isso o fatiamento!
  27. Se 11 ≠ 0, as operações elementares são descritas por:

    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
  28. 1 1 + 2 2 … = {1 1 +

    2 2 … } − (1 /11 ){ 11 1 + 21 2 + ⋯ } 1 = 0 Nova linha i Linha 1 Antiga linha i Multiplic. Observem que os multiplicadores foram escolhidos de forma a anular todos os novos 1 abaixo do 11 . 1 /11 31 /11 Multiplicadores 21 /11 11 ≠ 0
  29. As operações elementares são operações vetoriais: a multiplicação por fator

    de escala e a adição de vetores! Linha i Linha i Linha 1 −1 /11 × + Número real Vetor Vetor Vetor 11 ≠ 0
  30. Mestre, com o fatiamento da a NumPy essa operação com

    linhas é feita num golpe só, como um Samurai! (|)[, : ] = (|) , : − (1 /11 ) (|)[1, : ] Multiplic. Linha 1 Antiga linha i Nova linha i 11 ≠ 0
  31. Vamos proceder a contagem de operações elementares utilizadas no processo

    de eliminação! Utilizaremos os resultados: ෍ =1 = + 1 2 = Τ 2 2 + Τ 2 ෍ =1 2 = + 1 2 + 1 6 = Τ 3 3 + Τ 2 2 + Τ 6
  32. • 1 divisão para calcular 1 = −1 /11 •

    multiplicações para obter 1 × ℎ 1 • − 1 adições para obter ℎ + 1 × ℎ 1 Acima temos uma linha do processo de eliminação e cada linha da matriz A é um vetor com elementos. Portanto em cada linha faremos: ℎ ℎ ℎ 1 −1 /11 × + Número real Vetor 11 ≠ 0 Vetor Vetor
  33. Entretanto, como operaremos sobre a matriz aumentada, voltamos à conta

    anterior. Detalhe: como o objetivo é zerar o 1º coeficiente da nova linha faremos uma multiplicação e uma adição a menos. Claro, já sabemos que ele é zero !
  34. • − 1 divisões para calcular os 1 = −1

    /11 • − 1 multiplicações para obter as 1 × ℎ 1 • − 1 − 1 adições para obter as ℎ + 1 × ℎ 1 Na 1ª etapa ( = 1), eliminamos (zeramos) os coeficientes da 1ª coluna de da linha 1 para baixo. Portanto na 1ª etapa repetimos esse cálculo − 1 vezes, obtendo:
  35. Na 2ª etapa ( = 2) repetimos o processo para

    o subsistema sem a 1ª linha e a 1ª coluna. Então vamos zerar os coeficientes da 2ª coluna de A, da linha 2 para baixo. Assim por diante, até transformarmos a A na matriz triangular superior U, quando = − 1.
  36. • Divisões: σ=1 −1( − 1) = 2 + ()

    • Multiplicações: σ=1 −1 − 1 = Τ (3) 3 + (2) • Adições: σ=1 −1 − 1 − 1 = Τ (3) 3 + 2 + Portanto ao final do processo de eliminação o total de operações será: Somando tudo, teremos 2 × 2 /3 mais termos com 2 e .
  37. Na realidade, estamos interessados em estimativas desses valares quando cresce.

    Confiram na tabela que o número de operações na substituição é irrelevante perante o processo de eliminação quando cresce. Substituição Eliminação Soma % da eliminação 10 100 666 766 86,95 % 100 10.000 666.666 676.666 98,52 % 1.000 106 6,6667× 108 6,6767× 108 99,85 % 10.000 108 6,6667× 1011 6,6667× 1011 99.99 %
  38. Sim! Como é feito, de fato, o processo de eliminação?

    Torno a repetir vocês são bons de papo... E as contas?
  39. Linha 1 Linha 1 A NumPy é mais poderosa ainda!

    Não precisamos do loop em i. Vejam todo o 1º passo de uma só vez! [1: , 1: ] = 1: , 1: − ([1: , 1]/11 ) × [1,1: ] Linha 2 Linha 3 Linha n − × Linha 2 Linha 3 Linha n 21 /11 31 /11 ,1 /11 Linha 1 11 ≠ 0
  40. O 11 é intitulado “pivot”. Conforme veremos adiante, ele será

    o maior da coluna. Como os pivots de basquete. Essas operações exigem que o 11 seja não-nulo: 11 ≠ 0. Caso 11 = 0 interrompemos o processo. 1 /11 31 /11 Multiplicadores 21 /11 11 ≠ 0
  41. Aos meus pés uma versão sem o laço em i.

    Comparem com a anterior. Atentem para o fatiamento e para o uso do produto exterior, outer( ).
  42. Escrevi uma versão sem o laço em i, que inclui

    a solução do sistema. A entrada de dados é igual.
  43. Observem que estou operando sobre as linhas da matriz aumentada

    A|b. 11 12 13 ⋯ 1 1 ã22 ã23 ⋯ ⋯ ã2 ã3 ã2 ෨ 2 ã ෨
  44. A parte final do processo de eliminação. E também a

    solução mais a conferência visual.
  45. A 4ª linha de A é o dobro da 3ª.

    Neste caso o processo é interrompido.
  46. Neste outro também, mas as mensagens de erro são diferentes

    ! A 4ª e a 5ª linha de A são iguais.
  47. Para decidir se o sistema é impossível ou indeterminado precisamos

    examinar o termo independente. Com a matriz estendida é mais fácil !
  48. Nesses dois casos o sistema é impossível 01 + 02

    + 3 + 04 + 05 = −7 01 + 02 + 3 + 04 + 05 = 6
  49. A 4ª linha de A é o dobro da 3ª,

    inclusive os termos independentes 01 + 02 + 3 + 04 + 05 = 0 O sistema é indeterminado !
  50. Estabelecendo uma primeira conclusão: quando possível, o final do processo

    é o mesmo, um sistema triangular superior para resolver. | = 11 12 ⋯ 1 1 22 ⋯ 2 2 ⋱ ⋮ ⋮ Senão, temos que examinar, também, como termo independente evolui no processo de eliminação.
  51. Pois é jovens. O momento de aprender é este. Estejam

    atentos ! Imaginem se a consequência de um erro (descuido) em código 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
  52. O que causou o desastre no Paulo de Frontin eu

    não sei ! Entretanto no código programa escrito pelo Mestre, há um descuido ! É um descuido envolvendo a comparação de números de ponto flutuante ! É do tipo “de grão em grão a galinha enche o papo! “
  53. Erro, Mestra ? O Mestre foi super cuidadoso usando try

    para tratamento de erros Exatamente Loirinha ! O try do código abaixo testa se pivot ≠ 0 para calcular mult[ + 1: , ].
  54. Claro Mestra, para impedir divisão por zero ! Mas não

    impede a divisão por números extremamente, próximos de zero !
  55. Nesta matriz, a 3ª linha de A é a soma

    da 1ª com a 2ª. Portanto o sistema não terá solução única.
  56. Pois é... Imaginem se fosse o programa dos mísseis Scud

    e Patriot, ou para um cálculo para a estrutura do viaduto Paulo de Frontin ! Um programa com mais de 3.000 linhas onde só aparecem os resultados !
  57. Além do erro causado pelo meu descuido na comparação por

    igualdade (que eu já havia advertido) há algo mais ! Sim, Brute, além do seu erro imperdoável, há a possibilidade do algoritmo de eliminação que você usou ser instável.
  58. Chega de andar a 40 km/h Agora veremos o algoritmo

    de eliminação com pivotamento por coluna!
  59. Vamos retornar ao método da eliminação de Gauss e focar

    todos os holofotes no pivot. Pivot ? O que é pivot? Eles são os divisores usados em cada etapa do processo.
  60. 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.
  61. Que horror Mestra. Isto faz propagar muito o erro absoluto.

    Nossos estudos sobre número de condicionamento mostraram isso! Sim minha Loirinha querida! D’agora prá frente, despreze pivots pequenos, eles estragam o jogo, no final. Sempre vá em busca dos maiores possíveis.
  62. Uma forma de minimizar a propagação do erro 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.
  63. É verdade Mestra, pois então os multiplicadores ficam na faixa

    de -1 à +1 e não ampliam os erros! Exatamente Loirinha querida!
  64. É 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?
  65. 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.
  66. O exemplo do “Paulo de Frontin” e “Et tu Brute”

    revisitado pela 3ª vez. Agora indicando, corretamente, que a matriz é singular.
  67. A repetição dos três últimos exemplos objetiva mostrar que a

    eliminação com pivotamento é um algoritmo mais estável. Pois é, o algoritmo de eliminação sem pivotamento não é estável !