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

...

Transcript

1. Método de eliminação de Gauss e scipy.linalg Prof. Paulo R. G. Bordoni LNCC UFRJ
2. Um facilitador para lidar com os ndarrays da NumPy é a técnica de fatiamento.
3. Na NumPy, a faca para fatiar ndarrays é o operador : . É a ferramenta que ela oferece para cortar o pão nosso de cada dia !
4. O fatiamento é simplesmente uma forma que a NumPy oferece para construir submatrizes de uma matriz dada. Inclusive para pegar linhas ou colunas da matriz.
5. 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
6. 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]
7. 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]
8. 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 )
9. Veja Surfista, como pegar uma submatriz, como a abaixo: 0 1 2 3 0 1 2 3 4 5 [1: 3, 2: 5]
10. Este é um programa que faz alguns fatiamentos:
11. Executei seu programa uma vez, Mestre.
12. O restante da execução. Repeti a matriz A para facilitar a conferência.
13. Agora Surfista, vou ensinar a passar manteiga nas fatias de pão.
14. 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.
15. Agora vamos iniciar nosso estudo do método de eliminação de Gauss, um algoritmo para resolver sistemas de equações lineares.
16. 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.
17. Fui buscar estas informações na Wikipedia!
18. As honras da Alemanha à K. F. Gauss, por muitas outras contribuições à Matemática. E as nossas também, Mestre!
19. 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.
20. 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
21. 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 !
22. 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).
23. 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.
24. Fiz até este programa Mestra, ainda mais geral:
25. Marquei a retro- substituição, com tratamento de erro.
26. Vejam uma execução do meu programa:
27. Parabens! Você excluiu a possibilidade de divisão por zero e outros erros !
28. 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.
29. Esta é rotina da SciPy para resolver sistemas triangulares.
30. O detalhamento dos seus parâmetros.
31. O código do meu programa usando a solve_triangular( )
32. Uma execução do meu programa:
33. Outras duas execuções do meu programa:
34. Agora vou explicar o processo de eliminação.
35. 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.
36. 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.
37. 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 ~ .
38. 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.
39. 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.
40. 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.
41. Sim! Como é feita, de fato, a eliminação? Mestres, vocês descreveram o processo em linhas gerais! Só papo...
42. 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:
43. 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.
44. | = 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!
45. 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
46. 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
47. 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
48. 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
49. Brilhante esta sua observação Surfista! Aliás adorei este filme.
50. 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
51. • 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
52. 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 !
53. • − 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:
54. 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.
55. • 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 .
56. 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 %
57. Agora vamos usar os superpoderes da NumPy, para simplificar o código.
58. Sim! Como é feito, de fato, o processo de eliminação? Torno a repetir vocês são bons de papo... E as contas?
59. Fiz um programa para mostrar cada etapa da eliminação. Veja Surfista, o início do processo.
60. O restante do processo. Ele acaba ao obtermos uma matriz triangular superior.
61. Nem sempre o programa “enxerga” se o sistema é impossível ou indeterminado !
62. Esta é a entrada dos dados do programa.
63. A parte final do código, onde é realizada a eliminação: Atenção para o fatiamento
64. 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
65. 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
66. 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( ).
67. Escrevi uma versão sem o laço em i, que inclui a solução do sistema. A entrada de dados é igual.
68. Observem que estou operando sobre as linhas da matriz aumentada A|b. 11 12 13 ⋯ 1 1 ã22 ã23 ⋯ ⋯ ã2 ã3 ã2 ෨ 2 ã ෨
69. Neste 1º exemplo o sistema possui uma única solução.
70. A parte final do processo de eliminação. E também a solução mais a conferência visual.
71. A 4ª linha de A é o dobro da 3ª. Neste caso o processo é interrompido.
72. Neste outro também, mas as mensagens de erro são diferentes ! A 4ª e a 5ª linha de A são iguais.
73. Para decidir se o sistema é impossível ou indeterminado precisamos examinar o termo independente. Com a matriz estendida é mais fácil !
74. Nesses dois casos o sistema é impossível 01 + 02 + 3 + 04 + 05 = −7 01 + 02 + 3 + 04 + 05 = 6
75. A 4ª linha de A é o dobro da 3ª, inclusive os termos independentes 01 + 02 + 3 + 04 + 05 = 0 O sistema é indeterminado !
76. 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.
77. Jovens, precisamos reduzir a velocidade para evitar acidentes!
78. 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
79. 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! “
80. 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: , ].
81. Claro Mestra, para impedir divisão por zero ! Mas não impede a divisão por números extremamente, próximos de zero !
82. Nesta matriz, a 3ª linha de A é a soma da 1ª com a 2ª. Portanto o sistema não terá solução única.
83. Entretanto o programa fornece uma solução. Vejam: Errada, é só conferir !!!
84. 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 !
85. O Mestre deveria ter testado a possibilidade de um pivot quase nulo !
86. Et tu Master !?!
87. 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.
88. Chega de andar a 40 km/h Agora veremos o algoritmo de eliminação com pivotamento por coluna!
89. 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.
90. 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.
91. 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.
92. 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.
93. É verdade Mestra, pois então os multiplicadores ficam na faixa de -1 à +1 e não ampliam os erros! Exatamente Loirinha querida!
94. O valor máximo de um array pode ser obtido com a amax:
95. É 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?
96. Usem esta função para buscar a posição no array do seu valor máximo
97. Usei ambas no código:
98. Eis uma execução do programa.
99. Este programinha mostra como trocar duas linhas de uma matriz. Atenção ao detalhe:
100. 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.
101. Na realidade, o código só muda onde marquei.
102. Este é um 1º exemplo de execução.
103. O exemplo do “Paulo de Frontin” e “Et tu Brute” revisitado pela 3ª vez. Agora indicando, corretamente, que a matriz é singular.
104. A repetição do exemplo onde ℎ 4 = 2 × ℎ 3
105. E a repetição do caso em que ℎ 5 = ℎ 4
106. 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 !
107. Tchau, até a próxima aula!