Método diretos para problemas inversos e a NumPy

Transcript

1. UFRJ Métodos diretos para problemas inversos e a NumPy Prof. Paulo R. G. Bordoni
2. A resolução de sistemas de equações lineares é um problema inverso extremamente importante em Engenharia.
3. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Surfista, confira que o sistema de equações lineares acima pode ser escrito como o produto = , abaixo e vice-versa. 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋯ 1 2 ⋮ = 1 2 ⋮ =
4. 11 12 ⋯ 1 1 2 ⋮ = 11 1 + 12 2 + ⋯ + 1 Sim, e isso vale para todas as outras linhas, gerando todas as equações do sistema ! Claro Mestra, ao multiplicar a 1ª linha da matriz pelo vetor , obtenho exatamente o lado esquerdo da 1ª equação:
5. Mostre-nos, Manuel. Sim, e na NumPy, é o produto interno de pelo vetor : np. dot (, )
6. Aqui está ela:
7. Mestre, por quê a resolução de sistemas lineares é um problema tão importante em engenharia? Porque, conforme veremos, ela é usada para obter a solução aproximada de muitos outros problemas.
8. Para começar, Surfista, mostre como você resolveria uma equação do 1º grau + b = 0, O sistema linear mais simples de todos, apenas uma equação = − com uma incógnita .
9. = = + Τ − E no caso de = 0 a reta é paralela ao eixo-x, cortando o eixo-y no ponto = , ou coincidindo com o eixo-x, caso = 0. Geometricamente a solução, = − Τ , é o ponto onde a reta de equação = + corta o eixo-x.
10. Facílimo, Mestra ! Fiz o programa a seguir, usando float64 para as entradas , .
11. E tirou um zero redondo, Surfista. Seu programa admite a divisão por zero e pipoca em seguida:
12. Modifiquei seu programa para evitar esse erro, Surfista. Veja:
13. Uma forma bem mais elaborada é incluir uma rotina de tratamento de erro no programa:
14. Ela não só fornece a solução como trata corretamente as duas possibilidades de divisão por zero (impossível e indeterminado), Loirinha:
15. Inclusive “entende” alguns outros tipos de entrada:
16. Além disso, ela aponta outros erros de como de distração, p/ex.
17. “Forcei a barra”, Mestra e a rotina apontou!
18. Detalhes da seterr()
19. A continuação:
20. Por exemplo para ≈ 10−7, o que pode ocorrer ? Quando é muito pequeno, realmente próximo de zero, a reta definida por = + é praticamente paralela ao eixo-x, isto trará algum problema ?
21. Fui buscar os limites da IEEE 754/2008 para testarmos o programa da Mestra, que usou float32.
22. Para testar esses valores extremos, precisamos acrescentar apenas uma linha ao programa da Mestra:
23. Temos situações em que o valor de é enorme ou de overflow.
24. E também situações onde o valor de é realmente minúsculo ou de underflow.
25. Vamos explorar um pouco do que há na scipy.linalg para resolução de sistemas de equações lineares via métodos diretos. Depois veremos métodos iterativos.
26. Uma observação importante !
27. Os sistemas lineares mais simples de resolver são aqueles cuja matriz é diagonal. 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 3 0 0 0 4 0 0 0 5 1 2 3 4 5 = 1 2 3 4 5 Sim, a solução do sistema acima é: 1 2 3 4 5 = 1 /1 2 /2 3 /3 4 /4 5 /5
28. = 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 3 0 0 0 4 0 0 0 5 E o determinante de uma matriz diagonal é o produto dos termos da diagonal: det = 1 × 2 × ⋯ × Claro, é só pensar no desenvolvimento por linhas, ou por colunas !
29. Ele vai detonar uma bomba atômica ! O Torquinho vai armar uma bruta confusão com determinantes e sistemas lineares !
30. Nada disso, apenas um programinha inocente:
31. É intriga do Mestre, não vejo nada de mais no programa !
32. Executei com esses dados. Tudo Ok !
33. E funciona até quando o sistema não tem solução !
34. Se A é uma matriz ordem n, então det = det(). Um sistema linear = possui uma única solução quando, e apenas, quando det() ≠ 0. Os dois teoremas abaixo, são fatos sobejamente conhecidos sobre sistemas lineares e determinantes. Torno a avisar, o Torquinho é traiçoeiro !
35. 0 ⋮ 0 ⋯ 0 ⋮ ⋱ 0 0 0 0 1 2 ⋮ = 1 2 ⋮ Para = 0.0001 e 1 = 2 = ⋯ = = 0.001 a solução do sistema linear acima é, obviamente, 1 = 2 = ⋯ = = 10.
36. Entrei com os dados escolhidos pelo Sherlock !
37. Eu avisei !
38. Afinal, o determinante é nulo ou não ? E o sistema, tem ou não tem solução ?
39. Observem que, de fato, det = 0.000112 = 10−48 entretanto, no programa det _32 = 0.0 e det _64() = 1. −48 ≪ tiny32 Trata-se de um problema de precisão, jovens !
40. Vou explodir o mundo !
41. Cuidado Surfista, ele é traiçoeiro ! Vai nada, você só aumentou o tamanho do sistema !
42. Para a NumPy agora o determinante é zero mesmo ! Entretanto o sistema tem solução. Ele cumpriu a ameaça !
43. Observem que det = 0.000184 = 10−336 ≪ tiny64 e portanto, no programa, det_ 64 = 0.0 Mas, ainda assim, det > 0 .
44. Não há o que fazer. Teremos que aprender a lidar com a finitude do computador! Quem ele detonou foi a NumPy!
45. Agora apresentaremos mais um problema gerado pela finitude do computador e como contorná-lo !
46. Mestres, tentei melhorar o programa “sistema diagonal”, compararando o vetor diag com o vetor nulo e obtive esse erro : Como deveremos proceder ?
47. A mensagem de erro já traz sua resposta Loirinha: use any( ) or all( ) . Mestra, analise a situação.
48. A B Em: • Em (A) criamos um aviso com o valor False; caso algum elemento de diag for nulo, mudamos aviso para True. • Em (B) testamos aviso e mostramos o resultado segundo seu valor. No programa “sistema diagonal” procedemos como na linguagem C:
49. Com a NumPy trocamos (A) e (B) pelo código abaixo usando a função any( ), que é vetorizada ! Conferimos ainda a igualdade com outra função vetorizada, a equal( ).
50. A rotina any() da numpy
51. Cuidado com a tradução de any Surfista. Trata-se de um adverbio:
52. A all( ) é irmã da any( ).
53. A rotina equal( ) da numpy e suas irmãs, que permitem comparar dois arrays elemento a elemento.
54. Usei o novo código do Mestre e repeti os dados da “explosão do mundo” !
55. A conferência com equal( ) bateu com a conferência visual !
56. Entretanto, comparar números de ponto flutuante visualmente ou com igualdade frequentemente conduz a erros! Mestre, você está procurando cabelos na testa de um ovo!
57. Preste atenção neste exemplo, Surfista ! Cacilda ! A conferência com equal( ) indicou erro !
58. Quando estudarmos “Algoritmos e estabilidade”, Surfista ! Veremos em breve que o computador é uma luneta rapidíssima, mas defeituosa!
59. E como saímos dessa Mestres? Sim, expliquem tudo. Tim-tim por tim-tim !
60. Vocês nunca deverão comparar números de ponto flutuante IEEE 754/2008 visualmente ou via igualdade, mas sim por proximidade!
61. Mas e com os inteiros ? Com números inteiros você pode usar testes de igualdade sem nenhum problema.
62. Mostre as cartas escondidas na sua manga, Manuel. Seu pedido é uma ordem Galileu!
63. A rotina isclose():
64. E a rotina allclose():
65. Atenção, nas duas, a forma de medir a proximidade entre e é dada por:
66. Jamais esqueçam a fala do Mestre: Nunca comparem números de ponto flutuante IEEE 754/2008 via igualdade! Vamos refazer o exemplo anterior usando as isclose() e allclose().
67. O código do programa “sistema linear diagonal” só muda onde marquei:
68. Novos dados para calar a boca do Surfista !
69. Eis Escher tentando nos enganar visualmente, a verdade logo em seguida, e a comparação por proximidade recuperando o “quase igual”.
70. Repetindo Surfista: os Mestres compararam dois vetores y e b componente a componente, perguntando se elas diferiam por um erro minúsculo ! Exatamente, da ordem de 32 . E isto foi feito com:
71. Voltaremos ao tribar de Penrose e às ilusões visuais maravilhosas de Escher oportunamente !
72. Tchau, até a próxima aula!