Método diretos para problemas inversos

...

Transcript

1. UFRJ Métodos diretos para problemas inversos e scipy.linalg Prof. Paulo R. G. Bordoni
2. A resolução de problemas inversos é a parte mais desafiadora e interessante da Análise Numérica.
3. Uma forma de resolver problemas inversos é utilizando soluções de problemas diretos. Mestra, sinceramente, não entendi o trocadilho!
4. Não é trocadilho, Cabelos de Fogo. O que você faz para dividir um inteiro por outro inteiro ? Sim, como você calcula Τ 20 5 ?
5. “Tá brincando”, Mestra, o quociente é 4. “E fiz de cabeça”, nem pensei ! 20 5 4 0
6. Estou “vendo” seu pensamento, Surfista. Você resolveu um problema direto 5 x 4 = 20 e conferiu o resto. 20 5 4 -5x4 0
7. A bem da verdade, várias vezes. Usou a tabuada da multiplicação, até obter a igualdade. 5 x 1 = 5 5 x 2 = 10 5 x 3 = 15 5 x 4 = 20 É !
8. É, mas e no caso de 23/5 em que a divisão não dá exata – o resto não é zero ? 23 5 4 -20 3 Sim, e usou a definição de múltiplo inteiro para fazer a divisão: Τ 20 5 = 4 ↔ 5 × 4 = 20
9. 23 5 4 -20 3 Um quociente 4 e resto 3. Portanto o resultado é 3,6. 30 5 6 0 30 = 3 x 10 e
10. Vocês recorreram a três problemas diretos e ao conhecimento da representação decimal dos números: Não Mestra, foi feita uma subtração, que é um problema inverso !
11. Será sempre assim, Mestres ? É verdade, Cabelos de Fogo, mas 3 = 23 − 20 ↔ 23 = 20 + 3
12. Seja mais explícita na sua pergunta Loirinha ! Eu quero saber se sempre resolveremos problemas inversos por meio de uma quantidade finita de problemas diretos.
13. Não coleguinha, as dízimas periódicas respondem tua pergunta ! É verdade, Cabelos de Fogo! Uma dízima periódica é um número decimal com um grupo de dígitos após a vírgula se repetindo, sem parar !
14. Claro: Τ 7 3 = 2,33 … = 2, ത 3 Τ 13 7 = 1,451328 451328 … = 1, 451328 Você repete as operações até obter uma aproximação satisfatória !
15. Este é um exemplo de método iterativo. A solução aproximada é obtida por uma repetição de métodos diretos. Entra em jogo a noção de convergência de uma sequência.
16. A resolução de sistemas de equações lineares é um problema inverso extremamente importante em Engenharia.
17. 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 ⋮ =
18. 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:
19. Mostre-nos, Manuel. Sim, e na NumPy, é o produto interno de pelo vetor : np. dot (, )
20. Aqui está ela:
21. 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.
22. 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 .
23. = = + Τ − 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.
24. Facílimo, Mestra ! Fiz o programa a seguir, usando float64 para as entradas , .
25. E tirou um zero redondo, Surfista. Seu programa admite a divisão por zero e pipoca em seguida:
26. Modifiquei seu programa para evitar esse erro, Surfista. Veja:
27. Uma forma bem mais elaborada é incluir uma rotina de tratamento de erro no programa:
28. Ela não só fornece a solução como trata corretamente as duas possibilidades de divisão por zero (impossível e indeterminado), Loirinha:
29. Inclusive “entende” alguns outros tipos de entrada:
30. Além disso, ela aponta outros erros de como de distração, p/ex.
31. “Forcei a barra”, Mestra e a rotina apontou!
32. Detalhes da seterr()
33. A continuação:
34. 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 ?
35. Fui buscar os limites da IEEE 754/2008 para testarmos o programa da Mestra, que usou float32.
36. Para testar esses valores extremos, precisamos acrescentar apenas uma linha ao programa da Mestra:
37. Temos situações em que o valor de é enorme ou de overflow.
38. E também situações onde o valor de é realmente minúsculo ou de underflow.
39. 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.
40. Uma observação importante !
41. 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
42. = 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 !
43. Ele vai detonar uma bomba atômica ! O Torquinho vai armar uma bruta confusão com determinantes e sistemas lineares !
44. Nada disso, apenas um programinha inocente:
45. É intriga do Mestre, não vejo nada de mais no programa !
46. Executei com esses dados. Tudo Ok !
47. E funciona até quando o sistema não tem solução !
48. 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 !
49. 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.
50. Entrei com os dados escolhidos pelo Sherlock !
51. Eu avisei !
52. Afinal, o determinante é nulo ou não ? E o sistema, tem ou não tem solução ?
53. 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 !
54. Vou explodir o mundo !
55. Cuidado Surfista, ele é traiçoeiro ! Vai nada, você só aumentou o tamanho do sistema !
56. Para a NumPy agora o determinante é zero mesmo ! Entretanto o sistema tem solução. Ele cumpriu a ameaça !
57. Observem que det = 0.000184 = 10−336 ≪ tiny64 e portanto, no programa, det_ 64 = 0.0 Mas, ainda assim, det > 0 .
58. Não há o que fazer. Teremos que aprender a lidar com a finitude do computador! Quem ele detonou foi a NumPy!
59. Agora apresentaremos mais um problema gerado pela finitude do computador e como contorná-lo !
60. Mestres, tentei melhorar o programa “sistema diagonal”, compararando o vetor diag com o vetor nulo e obtive esse erro : Como deveremos proceder ?
61. A mensagem de erro já traz sua resposta Loirinha: use any( ) or all( ) . Mestra, analise a situação.
62. 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:
63. 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( ).
64. A rotina any() da numpy
65. Cuidado com a tradução de any Surfista. Trata-se de um advérbio:
66. A all( ) é irmã da any( ).
67. A rotina equal( ) da numpy e suas irmãs, que permitem comparar dois arrays elemento a elemento.
68. Usei o novo código do Mestre e repeti os dados da “explosão do mundo” !
69. A conferência com equal( ) bateu com a conferência visual !
70. 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!
71. Preste atenção neste exemplo, Surfista ! Cacilda ! A conferência com equal( ) indicou erro !
72. Já advertimos sobre esse tipo de problema em “Algoritmos e estabilidade”, Surfista ! Lembrem-se o computador é uma luneta rapidíssima, mas defeituosa!
73. E como saímos dessa Mestres? Sim, expliquem tudo. Tim-tim por tim-tim !
74. Vocês nunca deverão comparar números de ponto flutuante IEEE 754/2008 visualmente ou via igualdade, mas sim por proximidade!
75. Mostre as cartas escondidas na sua manga, Manuel. Seu pedido é uma ordem Galileu!
76. A rotina isclose():
77. E a rotina allclose():
78. Atenção, nas duas, a forma de medir a proximidade entre e é dada por:
79. 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().
80. O código do programa “sistema linear diagonal” só muda onde marquei:
81. Novos dados para calar a boca do Surfista !
82. Eis Escher tentando nos enganar visualmente, a verdade logo em seguida, e a comparação por proximidade recuperando o “quase igual”.
83. 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:
84. Voltaremos ao tribar de Penrose e às ilusões visuais maravilhosas de Escher oportunamente !
85. Tchau, até a próxima aula!
86. Curiosidade computacional: como um computador (processador, máquina virtual) executa uma subtração de inteiros.
87. O problema inverso mais simples é computar o oposto − de um número inteiro ∈ ℤ. Lembrem-se, o oposto de é o número ∈ ℤ que resolve a equação + = 0.
88. Sim para obter o oposto basta trocar o sinal de z ! Mestres, trata-se de um problema trivial !
89. A solução computacional é outra. Lembrem-se os únicos símbolos computacionais primários são o 0 e o 1. Todo o resto é feito com cadeias de 0’s e 1’s de tamanho pré-fixado. Como?
90. Ainda assim é trivial Sherlock! Se o computador opera com cadeias de 8 bits, ele reserva 7 para a magnitude do número e o 8º bit para o sinal, pondo 0 para + e 1 para −. É isso aí Surfista! Mestres, meu colega comprovou que trata-se, realmente, de um problema trivial !
91. Surfista, você apresentou um algoritmo para solucionar para o problema. Como o Sherlock disse, o algoritmo computacional é outro. Sim, a solução utilizada pelos cientistas de computação é de uma engenhosidade magnífica!
92. Para facilitar o entendimento descreveremos a solução nas bases 10 e 2. Começamos definindo o complemento à 9, 9 (), de um número decimal de um dígito. Na base 2 falamos em complemento à 1, 1 .
93. O 9 () é o número que somado à dá 9: + 9 = 9. Veja na tabela: 9 () 0 9 1 8 2 7 3 6 4 5 5 4 6 3 7 2 8 1 9 0 Na base 2 é mais simples ainda: 1 0 = 1 e 1 1 = 0.
94. Por exemplo 9 (1.036) = 8.963 e 1 1.001.101 = 0.110.010. Quando é um inteiro com vários dígitos o 9 () é obtido substituindo cada dígito de pelo seu 9 (). Na base 2 obtemos o 1 () “virando bits”: 0 → 1 e 1 → 0.
95. O correspondente na base 2 é o complemento a 2: 2 = 1 + 1. O complemento a 10 () de um número com vários dígitos é, simplesmente: 10 () = 9 () + 1.
96. Fácil, fácil: 10 2.751 = 7.248 + 1. Na base 2 é mais fácil ainda: 10 1.101.001 = 0.010.110 + 1.
97. Fácil fazer, mas não entendo aonde vocês dois querem chegar! Na próxima transparência farei mágica. Vou tirar um coelho da cartola!
98. Nossa, dá zero! O complemento a 10 de 4.532 é 5.468. Agora some esses dois números de 4 dígitos e despreze o 5º dígito do resultado. 4.532 + 5.468 1 0.000 lixo
99. Em outras palavras, para um número com dígitos, o complemento a 10 de é o “oposto” de (basta desprezar o dígito + 1 da soma). É por esse motivo que um número inteiro negativo é armazenado no computador, como 2 ().
100. Sim, quando você digita o número − 1.001.101 no teclado, o computador guarda 0.110.010 + 1. Para o sistema operacional, bastam duas operações triviais: • “Virar os bits” do número dado; • E depois somar 1.
101. É importante lembrar que, por padrão, os inteiros no computador já são armazenados em cadeias com 8, 16 ou 32 bits. Assim, se você digita 1 1001, o computador armazena 0001 1001 (com 8 bits). Logo −1 1001 é armazenado como 1110 0110 + 1.
102. 0001 1001 + 1110 0111 1 0000 0000 Gostei da mágica: Surfista, em Matemática não há mágica. Mestres expliquem isso! lixo
103. Loirinha, observe que o complemento a 10 de = 3.741 é 10 () = 6.258 + 1 e portanto: 1 + 6.258 = = 1 + 9 − 3 ∙ 1.000 + 9 − 7 ∙ 100 + + 9 − 4 ∙ 10 + 9 − 1 ∙ 1 = = 1 + 9 ∙ 1000 + 9 ∙ 100 + 9 ∙ 10 + 9 ∙ 1 − − 3 ∙ 1000 + 7 ∙ 100 + 4 ∙ 10 + 1 ∙ 1 = = 1 + 9.000 + 900 + 90 + 9 −(3.000 + 700 + 40 + 1) = 1 + 9.999 − 3741 = 10.000 − 3.741
104. Resumindo: 1 + 6.258 + 3.741 = 10.000 Ou seja, 10 + = 10.000 e desprezando o 5º dígito: 10 + = 0.000.
105. Resumindo tudo: O computador não faz subtrações, mas sim uma soma com o oposto do minuendo: Moral da estória: O problema inverso foi resolvido através de uma sequência de dois problemas diretos.
106. Agora sim, até a próxima aula!