Polinômios ortogonais e quadratura gaussiana

Transcript

1. LNCC UFRJ Polinômios ortogonais e quadratura Gaussiana Prof. Paulo R. G. Bordoni
2. Um conceito importante é o de ordem de precisão de uma fórmula de integração numérica.
3. O ordem de precisão de uma fórmula de integração numérica é o maior inteiro N para o qual ela integra = de forma exata, para = 0,1, ⋯ , . Mestra, é imediato comprovar que a regra trapezoidal tem ordem 1 de precisão.
4. Entretanto a regra de Simpson tem ordem 3, apesar de utilizar uma polinomial interpoladora de grau 2. Sim, pois a polinomial interpoladora utilizada para a regra trapezoidal tem grau 1.
5. Pois é. Mestre, faça um programinha para convencer o Surfista! Muito estranho Mestres, é algo totalmente inesperado!
6. Assinalei em vermelho os dados e resultados comprovadores.
7. Bem Mestra, neste caso funcionou – obtivemos o mesmo valor para o resultado por Simpson e a integral exata. E no caso geral? Respondendo pela Mestra – você tem duas possibilidades: • Conferir com várias outras polinomiais de grau ≤ 3 e deixar assim, ou • Provar a afirmação (que só é válida para pontos igualmente espaçados). Vá à luta!
8. A descrição do meu programa.
9. 1 2 3 4 Aqui eu assinalei: 1. A entrada de dados, 2. O cálculo da (), 3. O cálculo da integral exata, 4. A utilização da regra de Simpson com apenas 3 pontos.
10. Esta é a parte final do código, que gera o gráfico da () e da integral, correspondente à área sob a curva.
11. Em ℒ2(, ) temos outros conjuntos ortogonais completos além do trigonométrico, das séries clássicas de Fourier. Um deles é constituído pelos polinômios de Legendre, ℒ = 0 , 1 , 2 , … .
12. • 0 = 1, • 1 = , • +1 = 2+1 +1 − +1 −1 () para = 1,2, … As polinomiais de Legendre podem ser definidas recursivamente por:
13. • 0 = 1, • 1 = , • 2 = 3 2 2 − 1 2 , • 3 = 5 2 3 − 3 2 , • 4 = 35 8 4 − 15 4 2 + 3 8 , • 5 = 63 8 5 − 35 4 3 + 15 8 , • 6 = 231 16 6 − 315 16 4 + 105 16 2 − 5 16 . Calculei as até = 6, pela definição depois fui pegar uma onda!
14. Destaquei o módulo de Legendre no pacote polynomial da NumPy. Além do conteúdo coberto no módulo polinomial, que já estudamos, há um novo: Quadrature.
15. Mais detalhes do módulo. Recortei o início e o final, destacando o tópico Quadrature.
16. Clicando lá você vem para cá, Surfista. A classe Legendre. O domínio padrão é o intervalo [−1,1].
17. Este programa fornece a expressão e traça o gráfico da polinomial de Legendre de ordem n. Marquei os dois métodos da classe que usamos:
18. Usamos o método basis() que constrói as polinomiais da base de Legendre. Basta fornecer o grau k.
19. E o método convert() que faz a mudança de base, para a base canônica do espaço das polinomiais.
20. As duas primeiras polinomiais de Legendre:
21. As polinomiais de Legendre de grau 3 e 4:
22. As polinomiais de Legendre de graus 5 e 6:
23. Observe Loirinha, que nas polinomiais de Legendre de grau ímpar, todos os coeficientes das potências pares são nulos. É mesmo, e nas de grau par todos os coeficientes das potências ímpares são nulos!
24. Lembrem-se, uma função ↦ () é: • par se, e só se, − = , • ímpar se, e só se, − = − , sempre que −, ∈ (). Portanto as polinomiais de Legendre de grau par são funções pares e as de grau ímpar são funções ímpares.
25. Sim Mestra: • f ímpar ⟹ ׬ − = 0, • f par ⟹ ׬ − = 2 ׬ 0 , sempre que −, ⊂ (): Funções pares e ímpares desfrutam de propriedades destacadas na integração.
26. A prova dessas propriedades “é mole”, Mestre, basta usar que න − = න − 0 + න 0 mais a definição de função par ou ímpar.
27. Na NumPy e na literatura trabalha-se com as polinomiais de Legendre em ℒ2 −1, 1 . Para ℒ2 , basta efetuar a mudança de variáveis : −1,1 → [, ], definida por ↦ = + + − 2 Ou colocar = [, ] e = = , na chamada.
28. Para 0 = = 1 nem é necessário, ׬ −1 1 0 2 = 2. Logo 0 2 = 2 Fiz um programa para calcular a 2 da polinomial de Legendre de grau k.
29. No programa utilizei o método scipy.integrate.quadrature() para efetuar a integração numérica e aproximar 2 .
30. Por hora, vamos utilizá-la sem maiores explicações. Mais adiante vamos explicar o que é a quadratura Gaussiana nos mínimos detalhes.
31. Os detalhes para utilização da função quadrature()
32. As normas e a interpretação geométrica das 2 para = 1,2.
33. Idem para as 3 , 4 , 5 .
34. Prova-se que: 2 = න −1 1 ()2 = Τ 2 2 + 1 Assim podemos conferir todos os resultados calculados!
35. Vamos provar a seguir que as polinomiais de Legendre são ortogonais, isto é, que para todos os valores de , ∈ ℕ, ≠ , temos , = න −1 1 = 0
36. As polinomiais de Legendre são as autofunções do problema de Sturm-Liouville 1 − 2 () = −() correspondentes aos autovalores = ( + 1). Cruzes Galileu, de onde você tirou isso?
37. Ah coleguinha, é só pegar a expressão de uma e efetuar as operações indicadas no lado esquerdo da equação de Sturm Liouville que chegamos no lado direito. Comprovei isto para as 2 = 3 2 2 − 1 2 e 3 = 5 2 3 − 3 2 . O resultado geral, segue dessa sua ideia e da definição recursiva das , por indução finita em k.
38. Portanto satisfaz ′ = − + 1 , com = 1 − 2 . Multiplicando essa equação por e integrando em [−1, 1], obtemos න −1 1 ′ = −( + 1) න −1 1
39. Integrando por partes o lado esquerdo dessa igualdade, obtemos න −1 1 ′ = () ′ |−1 1 + න −1 1 ′ ′ Então obtemos a identidade + 1 න −1 1 = න −1 1 ′ ′ pois −1 = 1 = 0. 2 = 3 2 2 − 1 2
40. Trocando os papéis de e nessa identidade obtemos + 1 න −1 1 = න −1 1 ′ ′ Juntando as duas temos + 1 න −1 1 = + 1 න −1 1
41. Acabamos de provar que as polinomiais de Legendre são ortogonais, visto que para todos os valores de , ∈ ℕ, ≠ , temos , = න −1 1 = 0 Em outras palavras [ + 1 − ( + 1)] න −1 1 = 0
42. Modificando um pouquinho o último programa do Mestre, poderemos calcular o produto interno , = ׬ −1 1 . Em particular, para = , recalculamos a 2 . Para ≠ poderemos até confirmar a demonstração de ortogonalidade, embora isto não seja necessário!
43. A 1ª parte do novo programa. Novamente usamos a spi.quadrature().
44. E a parte gráfica.
45. Fiz 3 exemplos:
46. Claro Loirinha, definindo Τ = 2 , = 0,1,2, … temos que ∙ , ∶ ℒ2[−1,1] → ℝ é um aparelho medidor de doses de Legendre. Lembrei-me do dosador do Surfista!
47. Sim, para ∈ ℒ2[−1,1] a dose de que f possui é dada por = , = 1 2 න −1 1 E a sequência 0 , 1 , 2 , … é o espectro de Legendre de f.
48. Decidi calcular o espectro da função ∶ −1, 1 → [−1, 1] definida por = cos(). Vou te ajudar. Como f é uma função par seu espectro não terá componentes ímpares. Por exemplo: , 1 = 0.
49. Por acaso , 0 = 0., mas , 2 = −0.60793.
50. O valores espectrais decaem rapidamente, veja: , 4 = 0.12943 e , 6 = −0.01715
51. E , 8 = 0.00042, e , 10 = −0.00001.
52. Juntando tudo:
53. O espectro de uma outra função:
54. O espectro de uma polinomial de grau 7 só tem 7 elementos não-nulos.
55. O cocktail de Legendre de uma função polinomial é finito. Diferente dos cocktails para sen(x), exp(x) etc, cujos espectros são infinitos. Qualquer que seja a função polinomial p temos: • ∈ ℒ2 −1,1 , • p possui um espectro finito.
56. Em outras palavras, se é uma polinomial de grau n existem n+1 números 0 , 1 , … , tais que = ෍ =0
57. Agora exploraremos uma ideia importantíssima na integração numérica. A quadratura Gaussiana – as primeiras ideias foram apresentadas por Gauss em 1886. Gauss, C. F., "Methodus nova integralium valores per approximationem inveniendi" , Werke , 3 , K. Gesellschaft Wissenschaft. Göttingen (1886) pp. 163–196
58. Uma regra de quadratura é uma fórmula de integração aproximada de uma função () através de alguma função de ponderação (). São fórmulas do tipo: න −1 1 ≅ ෍ =1 ( )
59. A regra trapezoidal composta se enquadra na definição de regra de quadratura: න ≅ ෍ =0 , na qual = ቊ Τ ℎ 2 , = 0, ℎ , = 1,2, … − 1 e os são igualmente espaçados de h.
60. A regra de Simpson composta também se encaixa na definição de regra de quadratura; න ≅ ෍ =0 , na qual = ൞ Τ ℎ 3 , = 0, Τ 4ℎ 3 , = 1, … , − 1 í Τ 2ℎ 3 , = 1, … , − 1 e os são igualmente espaçados de h.
61. Ponderar: considerar o que é importante. Quando você vai adquirir algum bem, p/ex. uma casa, você pondera diversas coisas para efetuar a compra. O coeficiente de rendimento, CR, no seu histórico escolar é uma média ponderada das notas de cada disciplina. Os pesos são proporcionais à carga horária de cada disciplina, para valorizar as disciplinas com maior carga horária.
62. A quadradura gaussiana é uma regra de quadratura construída para produzir um resultado exato para polinômios de grau 2 − 1, ou menor, mediante uma escolha apropriada dos pontos 1 , 2 , … , e pesos 1 , 2 , … .
63. Observe que a polinomial interpoladora tem grau 1, mas os nós de interpolação não são (0. , exp 0. ) e (2. , exp(2.0)), como em Lagrange.
64. Sabemos que න 0 2 = 2 − 1 = 6.389 … e o resultado da quadratura gaussiana de grau 1 foi 6.3681..., preciso na 1ª casa decimal. Vou mostrar que uma polinomial interpoladora de Lagrange de grau 2 integra corretamente polinomiais de grau 5.
65. Conferindo: න −1 2 5 = Τ 6 6 |−1 2 = 10.5
66. Nas próximas transparências mostraremos aproximações sucessivamente melhores para න −1 1 3 usando quadratura de Gauss.
72. Marquei em vermelho as rotinas da scipy.integrate que realizam integração Gaussiana.
73. Fui la no “Tutorial” buscar informações sobre quadratura gaussiana:
74. Vamos começar examinando a fixed_quad( ).
75. O Mestre fez um programa para calcular aproximações para a integral , , = ׬ , via a rotina fixed_quad( ) de quadratura Gaussiana da scipy.integrate.
76. O programa recebe o intervalo de integração, a expressão da função a integrar e um limite superior N para a ordem de quadratura integração Gaussiana e devolve a tabela das aproximações e o gráfico.
77. Eis o programa:
78. A parte do programa que gera o gráfico.
79. Sim querida, já que 2 × 3 − 1 = 5. O programa do Mestre permite conferirmos. Mestra, uma regra de quadratura Gaussiana de ordem 3 deve, portanto, integrar corretamente uma polinomial de grau 5.
80. É mesmo, o resultado se repete a partir = 2.
81. Usando o outro programa do Mestre podemos conferir. Lá ele calcula a integral exata de polinomiais de qualquer ordem!
82. Métodos de integração numérica na scipy.integrate oriundos do QUADPACK
83. No plural porque são diversos...
84. Mais detalhes sobre os parâmetros.
85. Um programa para usar a rotina quad( ) de Integração Gaussiana da Scipy.
86. A parte gráfica:
87. Um exemplo:
88. Informações adicionais sobre a “quad” na scipy.integrate.
94. Apenas pelo tamanho das informações adicionais percebemos que seriam necessárias várias aulas para utilizar melhor essa poderosíssima rotina quad(). É o poder da abstração apresentando suas credencias.
95. Tchau, até a próxima aula!
96. Vou provar que a fórmula de integração gaussiana de grau p integra corretamente polinomiais de grau < 2. Por simplicidade vou usar = 4.
97. O Mestre vai provar que Escolhendo 0 , 1 , 2 e 3 como as raízes da polinomial de Legendre 4 (), então න −1 1 = ෍ =0 3 ( ) para toda polinomial ∈ 7 −1,1 . Nesta fórmula, = න −1 1 , onde { 0 , 1 ,2 , 3 } é a base de Lagrange de 3 ([−1,1]).
98. Trata-se de uma escolha genial! Explique porque, Filósofo. Porquê escolher logo as raízes da 4 (), Mestra?
99. O motivo fundamental é que as polinomiais de Legendre { 0 , 1 , … , 7 } formam uma base ortogonal para o espaço 7 ([−1,1]) das funções polinomiais com grau ≤ 7 em [−1,1]. Outras duas razões importantes são: • As raízes 0 , 1 , 2 e 3 de 4 estão no intervalo aberto −1,1 , • Elas são distintas entre si.
100. Você se lembra, Loirinha, como as funções básicas de Lagrange de 3 ( −1,1 ), as (), foram definidas? Claro Mestre, veja: 0 = ( − 1 )( − 2 )( − 3 ) (0 − 1 ) 0 − 2 0 − 3 1 = ( − 0 )( − 2 )( − 3 ) (1 − 0 )(1 − 2 )(1 − 3 ) 2 = ( − 0 )( − 1 )( − 3 ) (2 − 0 )(2 − 1 )(2 − 3 ) 3 = ( − 0 )( − 1 )( − 2 ) (3 − 0 )(3 − 1 )(3 − 2 )
101. Concordo contigo, Loirinha. Convençam-nos dessas afirmações, Mestres! Lembrei-me das definições, Mestre, mas acho muito estranho = න −1 1 .
102. Pois é Surfista, esta é uma das provas “mais lindas que já vi passar”, conforme nosso Vinícius de Moraes cantou a “Garota de Ipanema”. Vamos por partes. Consideremos primeiramente uma polinomial qualquer ∈ 3 ( −1,1 ).
103. Portanto, se ∈ 3 −1,1 então න −1 1 = ෍ =0 3 න −1 1 = ෍ =0 3 ( ) Bem, se ∈ 3 −1,1 então = 0 0 + 1 1 + 2 2 + 3 3 , com = , k = 0,1,2,3, conforme já vimos antes.
104. É, mas isto só funciona para ∈ 3 −1,1 . E para ∈ 7 −1,1 ? Bem Loirinha, é aí que entra a força maravilhosa das composições do Vinícius. Divida a sua polinomial pela 4 . Então teremos = 4 × + onde () ≤ 3 e também () ≤ 3.
105. Portanto න −1 1 = න −1 1 4 + න −1 1 Como () ≤ 3, por ortogonalidade temos න −1 1 4 = 0 .
106. Bem Loirinha, como ∈ 3 e { 0 , 1 , 2 , 3 } é uma base de 3 existem 0 , 1 , 2 , 3 ∈ ℝ para os quais = 0 0 + 1 1 + 2 2 + 3 3 , de forma que න −1 1 4 = ෍ =0 3 න −1 1 4 = 0 Por quê essa integral se anula ?
107. Fantástico, é aí que entrou a ortogonalidade: න −1 1 4 = 0 para ≠ 4. Bem, até agora vocês só mostraram que, න −1 1 = න −1 1 .
108. Bem Loirinha, de = 4 × + teremos = 4 × + e como as são raízes da 4 = . Agora o golpe final. Como ∈ 3 −1,1 segue, pelo primeiro resultado que mostrei, que න −1 1 = ෍ =0 3 = ෍ =0 3
109. Assim, juntando tudo num pacote só න −1 1 = න −1 1 = ෍ =0 3 , conforme prometemos.
110. As polinomiais de Legendre de ordem 2, 3 e 5, suas raízes e os pesos correspondentes.
111. Repetindo, achei fantásticas as ideias escondidas “debaixo do tapete” da ortogonalidade. Mestres, existem outras fórmulas de “integração gaussiana” tão poderosas quanto esta que usa as polinomiais de Legendre?
112. Bem, Loirinhas e Surfistas, para outros intervalos existem outras fórmulas de integração “gaussiana” apropriadas. Por exemplo para o intervalo 0, +∞ a escolha boa seriam as polinomiais de Laguerre, que formam um conjunto ortogonal completo em 2 (0, +∞).
113. Sim, mas aí é preciso utilizar o produto interno , = න 0 +∞ , onde : [0, +∞) → ℝ é uma função peso, definida pela expressão = exp(−).
114. Já para o intervalo (−∞, +∞) a escolha boa seriam as polinomiais de Hermite, que formam um conjunto ortogonal completo em 2 (−∞, +∞).
115. Sim, mas aí é preciso utilizar o produto interno , = න −∞ +∞ , onde : (−∞, +∞) → ℝ é uma outra função peso, definida pela expressão β = exp(−2).
116. Vocês precisam ir até a scipy.special para encontrar os polinômios ortogonais.
117. Vejam, além das polinomiais de Legendre, na scipy.special temos as polinomiais de Laguerre e Hermite.
118. Mais informações na scipy.special sobre Legendre, Laguerre e Hermite.
119. Como chamar a polinomial de Laguerre de grau n.
120. Como chamar a polinomial de Hermite de grau n.
121. Lembrar das equações de Legendre e Bessel