Métodos diretos para problemas inversos

Veja esta apresentação em www.bordoni.info.

Transcript

1. LNCC Métodos diretos para problemas inversos Prof. Paulo R. G. Bordoni UFRJ
2. A resolução de problemas inversos é a parte mais desafiadora da Análise Numérica.
3. 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.
4. Sim para obter o oposto basta trocar o sinal de z ! Mestres, trata-se de um problema trivial !
5. 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?
6. 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 !
7. Surfista, você apresentou uma solução para o problema. Como o Sherlock disse, a solução computacional é outra. Sim, a solução utilizada pelos cientistas de computação é de uma engenhosidade magnífica!
8. 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 .
9. 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.
10. 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.
11. 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.
12. 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.
13. 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!
14. 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
15. 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 ().
16. 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.
17. É 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.
18. 0001 1001 + 1110 0111 1 0000 0000 Gostei da mágica: Surfista, em Matemática não há mágica. Mestres expliquem isso! lixo
19. 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
20. Portanto 1 + 6.258 + 3.741 = 10.000 Ou seja, 10 + = 10.000 e desprezando o 5º dígito: 10 + = 0.000.
21. 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.
22. Para problemas inversos a ideia a perseguir será sempre essa: buscar a solução através de problemas diretos.
23. Para equações lineares, como + = 0, , ∈ ℝ, ≠ 0, a resposta é imediata: = −/. Mas descrever a obtenção da solução via problemas diretos é algo bem mais sofisticado Mestra. O problema envolve números reais e a divisão, que é a operação inversa da multiplicação.
24. Ora Mestre é só usar uma condicional no programa para testar se ≠ 0. Não vamos entrar nessa seara, Galileu. Entretanto você falou em divisão e para calcular a solução = −/ você precisa testar antes se ≠ 0.
25. Veja só:
26. Será que você está certo Surfista? Como já falou meu conterrâneo Shakespeare: “To be or not to be, that’s the question”.
27. Este é um programa um pouco mais elaborado que o seu, Surfista:
28. Detalhes da seterr( )
29. Continuação.
30. Alguns testes com valores extremos:
31. Mais alguns valores extremos:
32. Lembre-se estamos trabalhando com aritmética finita!
33. Como você descobre, Surfista, se um número ∈ ℝ está, ou não, próximo de zero? Qual é a pergunta que não pode calar, Mestra?
34. Ou, o que dá na mesma: Como você pode testar se = ? Vou mostrar um exemplo simples. Depois chamo o Manuel para ajudar!
35. Execute este programa, Surfista! Observe as rotinas que marquei:
36. Eis como testar se = . Vou chamar meu primo Manual!
37. A rotina equal( ):
38. A rotina isclose( ):
39. Os exemplos:
40. A rotina allclose( ):
41. Para equações quadráticas, tipo 2 + + = 0, , , ∈ ℝ, ≠ 0, a solução é dada pela fórmula de Bhaskara: ± = − ± 2 − 4 2 Inclusive, dependendo do valor do discriminante ∆ = 2 − 4, podemos ter duas, uma ou nenhuma solução.
42. Meça suas palavras, Loirinha! Equações do 2º grau sempre possuirão duas raízes reais ou complexas. É verdade, esqueci que: • ∆ > 0 ⇒ − , + ∈ ℝ • ∆ = 0 ⇒ − = + ∈ ℝ • ∆ < 0 ⇒ − , + ∈ ℂ
43. Resolva uma equação do 2º grau usando a NumPy, Surfista, no caso geral! Analise também o condicionamento da solução.
44. No século XVI matemáticos italianos como Cardano e Tartaglia e outros investiram muito na obtenção de fórmulas para calcular raízes de equações de grau maior que dois.
45. Deve-se a P. Ruffini (1799) e N. H. Abel (1824), a prova do Teorema que hoje leva seus nomes: Para ≥ 5, não há uma fórmula geral, envolvendo os coeficientes da equação polinomial + ⋯ + 1 1 + 0 = 0 que permita calcular suas raízes, como no caso das equações de 1º, 2º, 3º e 4º graus.
46. Um problema inverso extremamente importante em Engenharia é a resolução de sistemas de equações lineares. Sim Galileu, e hoje para quantidades enormes de equações e incógnitas.
47. Todos esses problemas inversos citados são resolvidos diretamente, através de fórmulas. Isto consagra o nome: métodos diretos para resolver problemas inversos, em particular para sistemas lineares.
48. Cramer apresentou sua fórmula para resolver sistemas lineares em 1750. Antes disso, Cardano, no final do século XVI, trabalhou com determinantes 2x2. Leibniz usou determinantes em 1693 para resolver sistemas lineares.
49. Um determinante de 2ª ordem, é denotado e definido por: Notem que o cálculo exige 2 multiplicações e 1 subtração. = 11 12 21 22 = 11 22 − 21 12
50. 1 = 1 12 2 22 = 1 22 − 2 12 , 2 = 11 1 21 2 = 11 2 − 21 2 , ≠ 0 A regra de Cramer para resolver um sistema 2x2 como este é (prove, Surfista): ቊ 11 1 + 12 2 = 1 21 1 + 22 2 = 2 Como são 3 determinantes teremos: • 3 × 2 multiplicações, • 3 × 1 adições • e 2 divisões.
51. = 2 3 1 −1 2 1 1 3 = −5 5 = −1 = 3 1 −1 3 2 1 1 3 = 10 5 = 2 Aplicando a regra de Cramer para resolver o sistema 2x2 ቊ 2 + = 3 + 3 = −1 obtemos:
52. = 11 12 13 21 22 23 31 32 33 = 11 22 23 32 33 − 21 12 13 32 33 + 31 12 13 22 23 Um determinante de 3ª ordem, é o número definido por: Esse número é a soma/subtração dos produtos de 3 números por 3 determinantes 2 × 2. Portanto serão necessárias 3 × 2 multiplicações e 3+2 adições para calcular um determinante 3 × 3.
53. 1 = 1 12 13 2 22 23 3 32 33 , 2 = 11 1 13 21 2 23 31 3 33 , 3 = 11 12 21 22 2 31 32 3 1 = 1 , 2 = 2 , 3 = 3 , ≠ 0 ቐ 11 1 + 12 2 + 12 3 = 1 21 1 + 22 2 + 23 3 = 2 31 1 + 32 2 + 33 3 = 3 A regra de Cramer para resolver um sistema 3x3 como o acima é dada por:
54. Como são 4 determinantes, a regra de Cramer para resolver um sistema 3 × 3 exige: • 4! multiplicações, • 4 × (3 + 2) adições, • e 3 divisões. Foram necessárias: • 3! multiplicações • e 3 + 2 adições/subtrações para cada determinante 3 × 3.
55. Este é um sistema linear constituído por n equações lineares à n incógnitas. Sendo breve: um sistema n x n. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋮ 1 1 + 2 2 + ⋯ + ⋮ =
56. Aliás, Lewis Carrol conta num livro, que a notação para matrizes que usamos hoje foi proposta por Leibniz. Na forma matricial esse sistema se reescreve: 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 2 ⋯ 1 2 ⋮ = 1 2 ⋮ =
57. O determinante det de uma matriz A de ordem n é definido recursivamente, como na próxima transparência. = = 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 2 ⋯ 1
58. Para n > 1, o desenvolvimento pela − é linha do det é dado por det = 1 1 + 2 2 + ⋯ + onde = (−1)+ e é o determinante de ordem n-1 obtido a partir da matriz excluindo-se sua − é linha e − é coluna. Há um desenvolvimento semelhante a partir da − é coluna. Para n = 1, det A = a 11 .
59. det(4 × 4) α det(3 × 3) β det(3 × 3) γ det(3 × 3) δ det(3 × 3) Τ + − Τ + − Τ + − α det(2 × 2) β det(2 × 2) γ det(2 × 2) det(3 × 3) Τ + − Τ + − α det(1 × 1) β det(1 × 1) det(2 × 2) Τ + − 2 mult. 1 adiç. 3 mult. 2 adiç. 4 mult. 3 adiç. O cálculo de um determinante de ordem 4 envolve: 4 × 3 × 2 × 1 = 4! multiplicações 3 × 2 × 1 = 3! adições/subtrações
60. Se = det ≠ 0 a solução é dada por 1 = 1 , 2 = 2 , ⋯ , = onde é o determinante da matriz A com a coluna k subtituída pelo termo independente b. A regra de Cramer para resolver um sistema linear n x n como este é: 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 2 ⋯ 1 2 ⋮ = 1 2 ⋮ =
61. Portanto, para resolver um sistema linear × precisaremos efetuar o cálculo de + 1 determinantes × e divisões. Assim, o total de operações envolvidas é: • + 1 ! multiplicações • da ordem de n! adições • divisões
62. Bobagem! Para o Titan – Cray XK7 é mole! O tempo gasto, por qualquer computador, para resolver sistemas 30 x 30 é inimaginável !
63. A tabela mostra o tempo gasto pelo Cray XK7, para resolver um sistema n x n: n tempo 20 48.41 min. 25 7.37 séculos 30 1.076 idade univ. Surfista querido, para um sistema 30x30 ele demora um pouco mais que a idade do universo!
64. O programinha que fez os cálculos.
65. Jamais esqueçam de contar o número de operações envolvidas num algoritmo! Sherlock, qual sua mensagem para Surfistas e Loirinhas?
66. Definitivamente, é inviável resolver sistemas lineares pela regra de Cramer! E a sua Galileu?
67. O título deste conjunto de transparências poderia ter sido “Para não dizer que não falei de flores”
68. Por quê Mestra? Ah Loirinha, porque o conteúdo desta aula é dado em Álgebra Linear!
69. Agora vamos usar os superpoderes das NumPy, para dirigir o magnífico calhambeque da capa pelo caminho das pedras.
70. 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!
71. 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.
72. 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]
73. 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!
74. Agora fatiaremos a coluna 2 de uma matriz 3x4: 0 1 2 3 0 1 2 [ ∶, 2] E agora Surfista?
75. 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]
76. 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: ]
77. 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 )
78. 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.
79. Agora vejamos como pegar submatrizes, como a abaixo: 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
80. Passei uma matriz M e peguei uma submatriz. 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
81. Agora vou passar manteiga nas fatias de pão, Mestra.
82. Multipliquei um pedaço da coluna 1 por 1.5:
83. Multipliquei uma sub-matriz por 0.5:
84. Fiz o produto direto de uma subcoluna por uma submatriz:
85. 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.
86. Agora passarei a expor o método de eliminação de Gauss.
87. 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.
88. Fui buscar estas informações na Wikipedia!
89. As honras da Alemanha à K. F. Gauss, por muitas outras contribuições à Matemática. E as nossas também, Mestre!
90. 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.
91. 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.
92. 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 ~ .
93. 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.
94. 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.
95. 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.
96. Sim! Como é feita, de fato, a eliminação? Mestres, vocês descreveram o processo em linhas gerais! Só papo...
97. 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:
98. 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.
99. | = 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!
100. 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
101. 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
102. 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
103. 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
104. Brilhante esta sua observação Surfista! Aliás adorei este filme.
105. Vamos proceder a contagem de operações elementares utilizadas no processo de eliminação!
106. 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
107. 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
108. 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!
109. Sim! Como é feito, de fato, o processo de eliminação? Torno a repetir vocês são bons de papo... E as contas?
110. Fiz um programa para mostrar cada etapa do processo de eliminação, seguindo a ideia do Surfista. Vejam o resultado.
111. Este é o código do programa. Marquei a observação do Surfista.
112. 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
113. 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
114. Aos meus pés uma versão sem o loop em i. Comparem com a anterior. Mais adiante, apresentaremos uma versão aperfeiçoada.
115. Neste 1º exemplo o sistema possui uma única solução.
116. 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.