Interpolação por partes

...

Transcript

1. Interpolação unidimensional por partes Prof. Paulo R. G. Bordoni UFRJ
2. Iniciamos revendo o exemplo do fenômeno de Runge, no qual o erro cresceu assustadoramente com o grau. O erro disparou nos extremos!
3. DIVIDIR PARA CONQUISTAR! Troquem um gigante por muitos anões! E aí Mestres, como faremos para matar esse leão?
4. Dividir para conquistar é um mantra da computação. A ideia é muito simples: 1. Dividimos o intervalo [, ] em vários subintervalos – o total de anões; 2. Em cada sub-intervalo usamos uma polinomial de grau baixo – os anões.
5. Troquei o gigante por 4 anões. Melhorou bastante nos extremos, mas piorou no centro. Não ficou bom!
6. Mestre, precisaremos de muitos anões. Vamos para Liliput com Gulliver!
7. 4 anões 14 anões 24 anões Sim, quanto mais anões, melhor:
8. Foi uma das primeiras lições que aprendemos na MatPlotLib. Lembro-me perfeitamente!
9. Do grande ao pequeno. É impossível não lembrar de M. C. Escher.
10. O pacote de interpolação da SciPy é grande. Precisaríamos de várias aulas explorá-lo por completo.
11. Risquei o que já vimos da scipy.interpolate.
12. Vamos ao 1º da lista:
13. Continuação:
14. • Recebe: 1. Um intervalo [, ], 2. Uma partição , = 0 , 1 , … , de [, ], 3. Um valor para cada ∈ , . • Retorna o gráfico de uma função linear por partes pelos pontos ( , ), utilizando a interp1d( ). Fiz um programa que:
15. Este é o programa. Note que informamos que a interpolação será linear em kind = ‘linear’
16. Os valores de ( , ) são escolhidos pelo usuário:
17. O programa: • Recebe: 1. Um intervalo [, ], 2. A expressão de uma função : [, ] → ℝ, 3. Uma partição de , . • E devolve: 1. Uma tabela dos pares dos pares , ( ) , 2. Os gráficos de e da interpoladora linear por partes de Lagrange da função sobre a partição , . Agora faremos uma interpolação linear por partes de uma função : [, ] → ℝ.
18. A única diferença é que = ( ).
19. Chamei-o com os dados à frente e veja o resultado:
20. Passaremos a analisar uma classe de espaços vetoriais importantíssima em teoria de aproximação.
21. Um desses espaços é o ℒ( , ). Ele é constituído pelas funções que são contínuas num intervalo [, ] e lineares em cada sub-intervalo de uma partição [, ] = { = 0 < 1 < ⋯ < = } de [, ].
22. A rotina interp1d( ) do pacote scipy.interpolate permite construirmos facilmente funções (vetores) desse subespaço. ℒ , é um subespaço de , .
23. Todo espaco vetorial tem uma base, ao menos uma. Mostre-nos uma base do ℒ( , ), Mestra. Ah, minha filha, uma delas é a minha coleção de chapéus de bruxa!
24. No meu tempo a Igreja não queimou apenas bruxas. Lembrem- se de Giordano Bruno! Eu mesmo escapei por pouco.
25. Chapeus de bruxa, Mestra? Amei! Sim Loirinha. A base canônica de ℒ( , ) é a coleção ℋ = {ℎ0 , ℎ1 , … , ℎ } constituída pelas N+1 “hat functions” que o Galileu vai te mostrar nas próximas transparências.
26. As funções chapéu, “hat functions”, formam uma base para esse espaço. Esta é a ℎ3 . A definição ponto-a-ponto dela é: ℎ3 = 0 ∈ [0,2) − 2 ∈ [2,3) 4 − ∈ [3,4) 0 ∈ [4,5]
27. O programa que desenha as “hat functions”
28. A continuação do programa.
29. Aqui estão todos os elementos da base. hat(0) hat(3) hat(2) hat(1) hat(4) hat(5)
30. Aspecto fundamental e importantíssimo: ℒ( , ) é um subespaço de , de dimensão finita + 1. É essa característica que permite computarmos aproximações para equações diferenciais ordinárias e parciais através do método dos elementos finitos.
31. 3 1 4 3 4 Nesta figura desenhei em azul uma combinação linear de duas delas: 3 ℎ3 + 4 ℎ4
32. A função abaixo é uma combinação linear de cinco funções-chapéu ℎ0 , ℎ1 , … , ℎ5 : = 0.3ℎ0 + 0.9ℎ1 + 0.5ℎ2 + 1.1ℎ3 + 1.3ℎ4 + 1.0ℎ5
33. Veremos adiante que essa característica é fundamental para o método dos elementos finitos. Um aspecto fundamental da “hat function” ℎ3 é que ela é nula fora dos dois subintervalos que possuem o ponto 3 . Os matemáticos dizem: uma função a suporte compacto. O suporte da ℎ3 é o intervalo fechado [2.0, 4.0]. O suporte da função ℎ3 .
34. • K é um subconjunto fechado de ഥ Ω = [a, b]; • é um espaço de funções polinomiais para as funções de forma; • Σ é o conjunto de graus de liberdade local. Para seu conhecimento: Um elemento finito é uma terna (, , Σ) na qual:
35. a b el 1 a el 1 = − − = − − No método dos elementos finitos elas recebem o nome de coordenadas baricêntricas, e são usadas as letras gregas , no lugar de , . Na construção dos elementos finitos lineares usamos as duas funções da base de Lagrange linear definidas no elemento = [, ]:
36. São propriedades básicas das coordenadas baricêntricas: • ∈ 1 (), = 1,2, = [0 , 1 ] • = = ቊ 1, = 0, ≠ , , = 1,2 • 1 + 2 = 1, ∀ ∈ As derivadas das coordenadas baricêntricas são: ቊ 0 () = − Τ 1 ℎ 1 () = Τ 1 ℎ , ∈ .
37. A correspondência global x local ℎ5 ℎ6 5 6 Base global 0 1 5 6 Base local A correspondência entre as funções básicas no intervalo [5 , 6 ] : ൜ ℎ5 = 0 () ℎ6 () = 1 () para ∈ [5 , 6 ]. ou ℎ5 | 5,6 = 0 , ℎ6 | 5,6 = 1
38. Numa próxima aula resolveremos alguns problemas de valor de contorno pelo método dos elementos finitos.
39. Agora passaremos a examinar aproximações quadráticas por partes.
40. , é o subespaço de , constituído pelas funções que são quadráticas em cada par consecutivo de sub-intervalo de uma partição [, ] = { = 0 < 1 < ⋯ < = } do intervalo [, ] – os elementos finitos quadráticos. Anotaremos o espaço de funções quadráticas por partes por ( , ).
41. Atenção Surfista: Como necessitamos de pares de subintervalos consecutivos, o número N da partição , obrigatoriamente será ímpar.
42. Para enfatizar o uso de polinomiais de grau ≤ 2 usarei um ℎ −1 , , +1 , −1 , , +1 em cada três pontos sucessivos −1 , , +1 , = 1, … , − 1 (pares de subintervalos) da partição [, ] = { = 0 < 1 < ⋯ < = } de , .
43. Pintarei cada polinomial de grau ≤ 2 com uma cor diferente, para ficar clara a ideia de “por partes”.
44. O programa que faz o gráfico das interpoladoras quadráticas. Marquei onde usamos a KroghInterpolator( ).
45. Usei o mesmo programa que chama a rotina interp1d( ) da scipy.interpolate mas agora para desenhar uma função quadrática por partes em cada subintervalo da partição. A única diferença foi no parâmetro kind.
46. Observem a grande melhoria na aproximação funções quadráticos com relação à brincadeira de trocar um gigante por vários anões. Incrível, é a mesma função do leão!
47. O programa da interpolação quadrática por partes. Ele é uma pequena modificação do anterior (o colorido).
48. A continuação do programa.
49. Você sabe porque o chamei aqui, Mr. Charlie Chaplin? Sei sim, Mestre. Para mostrar meu chapéu coco!
50. Eu explico, Surfista. Chapéus de bruxa são pontudos e retos, mas chapéus coco são arredondados – como parábolas, polinomiais de grau 2. E dá para explicar porque, Mr. Chaplin?
51. Saquei Mestre! As funções básicas agora serão parábolas com abas, como chapéus coco. 2 2 1 1 Ótimo, mas tirou só meio-certo, Surfista!
52. Vou acordar este mexicano que está dormindo. Preciso do chapéu dele! Mestra, não faça isto, é o Mestre descansando!
53. É que além dos chapéus-côco, vamos precisar também de “sombreros” mexicanos pontudos! Porque uma base do espaço das funções quadráticas por partes é constituída por chapéus-côco e “sombreros” mexicanos. Vejam nas próximas transparências.
54. Veja o gráfico da “chapéu-coco” nº 3
55. Zzz..., “es el sombrero” mexicano nº4
56. Agora vou mostrar uma base global para , . Usarei 4 pares de subintervalos. Portanto = 2 × 4 + 1 = 9, a dimensão de , . Admirem as funções básicas Loirinhas e Surfistas!
57. Todos os 5 chapéus mexicanos
58. E todos os 4 chapéus-côco:
59. Na construção dos elementos finitos quadráticos usamos as três funções da base de Lagrange quadrática (as funções de forma) definidas localmente no elemento = [, , ]. = − − ∙ − − = − − ∙ − − = − − ∙ − −
60. As funções de forma são portanto definidas no elemento = , , através do produto das coordenadas baricêntricas: = () ∙ () = () ∙ () = () ∙ ()
61. Na base exibida, a partição , foi uniforme (pontos = 0 + ∙ ℎ, igualmente espaçados) e as funções da base de Lagrange de , são: 0 = 1 2ℎ2 ( − 1 )( − 2 ) 1 = −1 ℎ2 ( − 0 )( − 2 ) 2 = 1 2ℎ2 ( − 0 )( − 1 )
62. 0 = 1 ℎ2 − Τ 1 + 2 2 1 = −1 ℎ2 − Τ 0 + 2 2 2 = 1 ℎ2 [ − Τ 0 + 1 2] As derivadas das coordenadas baricêntricas são:
63. São propriedades básicas das coordenadas baricêntricas: • ∈ 2 (), = 0,1,2, = [0 , 1 , 2 ] • = = ቊ 1, = 0, ≠ , , = 0,1,2 • 0 + 1 + 2 = 1, ∀ ∈
64. O programa que mostrou as funções básicas ().
65. Continuação.
66. O final do programa.
67. Além da interpolação polinomial por partes de Lagrange estudaremos a de Hermite.
68. p() q() 2 2 3 1 Na interpolação polinomial por partes, de Lagrange, a única exigência é a continuidade da interpolação. No gráfico: 2 = 2 = (2 ).
69. p() q() 2 2 3 1 Observem que, no extremo comum, (2 , 2 ), a “colagem” entre p e q não é suave. Matematicamente (2 ) ≠ (2 ). Ondas abruptas, mas surfistas com curvas suaves...
70. p() q() 2 2 3 1 A imposição da condição 2 = 2 suaviza a colagem, pois obriga a coincidência das retas tangentes. A interpolação polinomial por partes, de Hermite, exige, em cada ponto, a igualdade tanto do valor da função, como do valor da derivada da função. O dobro de condições.
71. Marquei a classe PchipInterpolator e sua função de conveniência na scipy,interpolate
72. Pchip é uma abreviação para “Piecewise Cubic Hermite Interpolator”.
73. Assinalei como a derivada é estimada pelo algoritmo.
74. É muito mais fácil usar a função de conveniência:
75. Usei-a neste programa:
76. Vejam dois resultados, com 5 e 13 pontos:
77. Agora a classe PiecewisePolynomial.
78. Notem que a classe PiecewisePolynomial está presente nq versão 17.1 da Scipy, mas não aparece mais na versão 19.0
79. Para efetuar a interpolação com os métodos da ( ), o parâmetro é um array unidimensional com tamanho n, o número de pontos da partição do − . Já o parâmetro é um array bidimensional, com tamanho n x k, com ≥ 2. Para interpolar uma função f, o parâmetro será constituído por n pares , . Entretanto poderiam ser dois valores quaisquer.
80. A interpolação de Hermite. O programa a seguir sorteia 7 pontos e 7 direções tangentes e constrói a interpoladora.
81. O parâmetro “orders” é mandatário com relação à condição de tangência. Sim, basta conferir os graus e os gráficos das polinomiais em cada sub-intervalo.
82. O programa que traçou o gráfico da transparência anterior .
83. = = 3 O default = é uma polinomial de grau 3 em cada sub-intervalo.
84. Repetindo, quando não especificamos o parâmetro orders, ele assume o valor 3. Em outras palavras, as polinomiais são de grau 3 em cada sub-intervalo.
85. O programa que traça o gráfico da função e da interpoladora, com = 3.
86. É só usar = 0. A interpoladora constante por partes usando a função ( ) do pacote interpolate da SciPy.
87. Agora usando = 1 temos a interpolação linear por partes.
88. Com = 2 temos interpolação quadrática em cada sub-intervalo. Neste caso a tangência só ocorre no extremo direito de cada subintervalo.
89. A interpolação de Hermite Especificando o parâmetro = 3 ou deixando = .
90. p() q() 2 2 3 1 A imposição da condição ′′ 2 = ′′ 2 obriga a mesma curvatura local para as duas curvas p e q que possuem o ponto (2 , 2 ) em comum. A curvatura local é o inverso do raio da circunferência osculadora – procure um livro de computação gráfica para mais detalhes. Observação
91. Para usar = 4 precisamos fornecer também as derivadas 2ªs. A polinomial, sua derivada e sua derivada 2ª são iguais em cada nó.
92. Observe que agora são 3n condições para n pontos. Neste exemplo 12.
93. Não deixem de estudar a classe Ppoly.
94. Continuação da classe Ppoly.
95. O final da classe Ppoly.
96. Tchau, até a próxima aula.