O método de eliminação de Gauss

...

Transcript

1. Método de eliminação de Gauss Prof. Paulo R. G. Bordoni LNCC UFRJ
2. O título deste conjunto de transparências poderia ter sido “Para não dizer que não falei de flores”
3. Por quê Mestra? Ah Loirinha, porque o conteúdo desta aula é dado em Álgebra Linear!
4. Agora vamos usar os superpoderes das NumPy, para dirigir o magnífico calhambeque da capa pelo caminho das pedras.
5. Antes, precisaremos aprender a obter fatias de ndarrays e a operar com elas. Na NumPy, a faca para fatiar arrays é o operador : . É a ferramenta que ela oferece para cortar o pão!
6. 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, em NumPy, a indexação de arrays começa por zero e não por um.
7. 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]
8. Nossa 1ª fatiada será a linha 1 de uma matriz 4x5: 0 1 2 3 0 1 2 3 5 [1, : ] Esta é outra opção, Mestra!
9. Agora fatiaremos a coluna 2 de uma matriz 3x4: 0 1 2 3 0 1 2 [ ∶, 2] E agora Surfista?
10. 0 1 2 3 0 1 2 Se A é uma matriz m x n, então: • [ , : ] fornece a linha i de A, de j em diante; • [ : , ] fornece a coluna j de A, de i para baixo. [1,2: ] 0 1 2 3 0 1 2 3 [ 1: , 0]
11. Agora a fatia é uma parte da linha 1 de A: da posição 2 até o final. Faça o código para a outra figura Surfista. 0 1 2 3 0 1 2 [1,2: ]
12. 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 k-1. • Pense em j:k como o intervalo [j, k). 0 1 2 3 0 1 2 3 4 5 [1,1: 4] : k `` = ´´ [ j, k )
13. 0 1 2 3 0 1 2 3 4 5 [1,1: 4] Uma parte da linha 1, aquela correspondente aos índices [1,4) das colunas.
14. Agora vejamos como pegar submatrizes, como a abaixo: 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
15. Passei uma matriz M e peguei uma submatriz. 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
16. Agora vou passar manteiga nas fatias de pão, Mestra.
17. Multipliquei um pedaço da coluna 1 por 1.5:
18. Multipliquei uma sub-matriz por 0.5:
19. Fiz o produto direto de uma subcoluna por uma submatriz:
20. Portanto, antes, elas precisam ser convertidas, usando o modificador . A (esse A é de array): nome_da_matriz . A O produto direto de entre duas matrizes compatíveis é realizado em NumPy/SciPy através do produto delas como ndarrays.
21. Agora passarei a expor o método de eliminação de Gauss.
22. 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.
23. Fui buscar estas informações na Wikipedia!
24. As honras da Alemanha à K. F. Gauss, por muitas outras contribuições à Matemática. E as nossas também, Mestre!
25. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Vamos agora descrever como resolver um sistema linear como o abaixo, pelo método de eliminação. Descreveremos o processo, primeiramente, em linhas gerais, usando as equações.
26. 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 Nosso objetivo será atingido após várias etapas consecutivas.
27. 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 ~ .
28. 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.
29. 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.
30. O último passo consiste em resolver o sistema triangular. Algo extremamente simples. 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.
31. Sim! Como é feita, de fato, a eliminação? Mestres, vocês descreveram o processo em linhas gerais! Só papo...
32. 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:
33. 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.
34. | = 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!
35. 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
36. 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
37. 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
38. 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
39. Brilhante esta sua observação Surfista! Aliás adorei este filme.
40. Vamos proceder a contagem de operações elementares utilizadas no processo de eliminação!
41. Linha 1 Linha 1 Linha 2 Linha 2 Linha 1 21 /11 × − Linha n Linha n Linha 1 1 /11 × − Sim Sherlock, na 1ª etapa: • Efetuamos − 1 divisões, para calcular os multiplicadores, • Em cada linha − 1 multiplicações e subrações • Como são − 1 linhas: • ( − 1)2 multiplicações • ( − 1)2 subtrações
42. 1 1 + 2 2 … = {1 1 + 2 2 … } − (1 /11 ){ 11 1 + 21 2 + ⋯ } 1 = 0 Multiplicações Atenção Surfista, como já sabemos que 1 = 0 não precisamos efetuar essa conta! Discordo, são n multiplicações e n subtrações! Subtrações
43. 11 1 + 12 2 + ⋯ + 1 = 1 ෤ 22 2 + ⋯ + ෤ 2 = ෨ 2 ⋯ ෤ 2 2 + ⋯ + ෤ = ෨ Depois repetimos o processo para o subsistema sem a 1ª linha e assim por diante! Assim é só somar as operações em todas as n-1 etapas. Acabe esta conta Surfista!
44. Sim! Como é feito, de fato, o processo de eliminação? Torno a repetir vocês são bons de papo... E as contas?
45. Fiz um programa para mostrar cada etapa do processo de eliminação, seguindo a ideia do Surfista. Vejam o resultado.
46. Este é o código do programa. Marquei a observação do Surfista.
47. 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
48. 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
49. Aos meus pés uma versão sem o loop em i. Comparem com a anterior. Mais adiante, apresentaremos uma versão aperfeiçoada.
50. Neste 1º exemplo o sistema possui uma única solução.
51. Neste 2º exemplo o sistema não possui uma única solução. Na etapa 2, a 4ª equação é 0 = −10. . Portanto o sistema é impossível.
52. Neste 3º exemplo o sistema não possui uma única solução. Na etapa 1, a 3ª equação é 0. = 0. . Portanto o sistema possui uma infinidade de soluções.
53. 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.
54. Jovens, precisamos reduzir a velocidade para evitar acidentes!
55. 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”.
56. 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.
57. A B Assinalei a única diferença.
58. 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!
59. Para os dados marcados em azul, o código em A oferece a solução abaixo! A
60. Para os mesmos dados o código B oferece a resposta correta: o sistema é impossível. B
61. 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
62. 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
63. E produzir, em seguida, o erro marcado aos meus pés. B Erro evitado pelas linhas marcadas do código B.
64. 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!
65. 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
66. 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.
67. Na verdade, a forma mais indicada para evitar esse tipo de desastre é usar a função isclose() da NumPy.
68. Fui buscar informações sobre ela em Logic Functions, na referência da NumPy.
69. B Nossa proposta é substituir as linhas 29 e 30 do código B, que estão corretas, por algo mais pythônico, usando isclose().
70. Agora vamos repetir o programa do início destas transparências. Nesta 1ª parte apenas as mensagens mudam.
71. A única diferença do programa do início é a utilização de isclose.
72. A execução com o mesmo sistema linear que gerou problemas: Perfeito, Mestres!
73. Chega de andar a 40 km/h Agora poderemos retomar nossa velocidade normal!
74. 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.
75. 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.
76. 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.
77. 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.
78. É verdade Mestra, pois então os multiplicadores ficam na faixa de -1 à +1 e não ampliam os erros! Exatamente Loirinha querida!
79. O valor máximo de um array pode ser obtido com a amax:
80. É 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?
81. Usem esta função para buscar a posição no array do seu valor máximo
82. Existe um método da classe matrix que faz praticamente a mesma coisa. A diferença está marcada lá embaixo.
83. Vejam elas funcionando:
84. Mestre, eis como trocar duas linhas de uma matriz:
85. 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.
86. 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
87. Este é um 1º exemplo de execução. Surfista, faça outros, semelhantes aos exemplos anteriores. 1 2 3
88. 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.
89. Valeu, Mestre! 1 3 2
90. Voltaremos, oportunamente, à resolução de sistemas lineares.
91. 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.
92. 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:
93. | = 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.
94. 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.
95. • 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!
96. Começaremos pelas trocas, Loirinha!
97. 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.
98. 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 ⟷ .
99. = 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.
100. Este programa faz o que marquei aqui no início.
101. Ele permite conferir que o produto de permutações é comutativo. Quero lembrar que, em geral, o produto matricial não é comutativo.
102. Ele permite mostrar, também, que a inversa da permutação ↔ , é a permutação ↔ .
103. Ele permite, ainda, mostrar que as permutações são idempotentes: 2 = . Observe, Surfista, que neste exemplo escolhi = .
104. 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
105. 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.
106. 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
107. 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:
108. Um programa que mostra o produto da 1 = − 1 por A.
109. 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!
110. O programa que fiz para conferir a afirmação da Mestra!
111. = 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.
112. −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
113. Eu fiz este programa para efetuar a eliminação, passo a passo!
114. Vejam o resultado! A matriz A é atualizada em cada etapa.
115. 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!
116. Primeiro um programa para a fatoração LU sem pivotamento.
117. 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!
118. 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.
119. 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!
120. 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?
121. 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
122. 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.
123. Tchau, até a próxima aula!