Interpolalaçao 1D por partes e splines

Por efetuar.

Transcript

1. Interpolação unidimensional por partes e splines 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. A ideia é muito simples: 1. Dividimos o intervalo [, ] em vários subintervalos. 2. Em cada sub-intervalo usamos uma polinomial de grau baixo. Dividir para conquistar é um mantra da computação.
5. Troquei o gigante por 5 anões. Melhorou muito nos extremos, mas piorou um pouco no centro. Poucos anões!
6. Sim, quanto mais anões, melhor:
7. Foi uma das primeiras lições que aprendemos na MatPlotLib. Lembro-me perfeitamente!
8. Com mais de 50 pontos as curvas ficam suaves.
9. Passaremos a analisar uma classe de espaços vetoriais importantíssima em teoria de aproximação.
10. Um desses espaços é o ℒ( , ). ℒ , é o subespaço de , constituído pelas funções que são lineares em cada sub-intervalo de uma partição [, ] = { = 0 < 1 < ⋯ < = } do intervalo [, ].
11. Veja a cara de uma função do espaço ℒ(9 0,2 ) no qual 9 0,2 = {0. , 0.25, 0.5, 0.75, 1. , 1.25, 1.5, 1.75, 2. } Os correspondentes foram sorteados.
12. O programa que produziu esse gráfico. Em vermelho a função da scipy.interpolate que gera a interpoladora. Já, já vamos estudá-la.
13. Todo espaco vetorial tem, ao menos, uma base. Mostre-nos uma base do ℒ( , ), Mestra. Ah, minha filha, uma delas é uma coleção de chapéus de bruxa!
14. Sim Loirinha. A base canônica de ℒ( , ) é a coleção ℋ = {ℎ0 , ℎ1 , … , ℎ } constituída pelas N+1 “hat functions” que o Filósofo vai te mostrar nas próxima transparências. Chapeus de bruxa, Mestra? Amei!
15. As funções chapéu, “hat functions”, formam uma base para esse espaço.
16. O programa que desenha as “hat functions”
17. A continuação do programa.
18. Aqui estão todos os elementos da base. hat(0) hat(3) hat(2) hat(1) hat(4) hat(5)
19. 3 1 4 3 4 Nesta figura mostro uma combinação linear de duas delas: = 3 ℎ3 + 4 ℎ4
20. A função () = = 0.4ℎ0 () + 0.7ℎ1 () + 1.2ℎ2 () + 1.3ℎ3 () + 0.9ℎ4 () + 0.6ℎ5 ()
21. 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.
22. Um elemento finito é uma terna (, , Σ) na qual: • K é um subconjunto fechado de Ω; • é um espaço de funções polinomiais para as funções de forma; • Σ é o conjunto de graus de liberdade local.
23. 0 1 e 1 1 0 1 e 1 0 Na construção dos elementos finitos lineares usamos as duas funções da base de Lagrange linear definidas no elemento = [0 , 1 ]: 0 = 1 ℎ (1 − ) 1 = 1 ℎ ( − 0 ) No método dos elementos finitos elas recebem o nome de coordenadas baricêntricas, e são usadas as letras gregas 0 , 1 no lugar de 0 , 1 .
24. 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: • Dλ0 x = −1/h • Dλ1 x = 1/h
25. ℎ5 ℎ6 5 6 0 1 5 6 A correspondência entre as funções básicas: ℎ5 = 0 () ℎ6 () = 1 () para ∈ = [5 , 6 ]. ou ℎ5 | = 0 , ℎ6 | = 1 Base global X Base local
26. Você sabe porque o chamei aqui, Mr. Charlie Chaplin? Sei sim, Mestre. Para mostrar meu chapéu coco!
27. 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?
28. 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!
29. Vou acordar este mexicano que está dormindo. Preciso do chapéu dele! Mestra, não faça isto, é o Mestre!
30. É que além dos chapéu-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.
31. Veja o gráfico da “chapéu-coco” nº 3
32. Zzz..., “es lo sombrero” mexicano nº4
33. Todos os 5 chapéus mexicanos
34. E todas as 4 chapéu-coco:
35. O programa que mostra as funções básicas Quad(no).
36. Continuação.
37. O final do programa.
38. Faremos interpolações quadráticas por partes em subintervalos sucessivos (= elementos finitos) definidos por ternas de pontos (= nós) de uma partição (malha) [, ] = { = 0 < 1 < ⋯ < = } de , .
39. O programa que faz o gráfico das interpoladoras quadráticas. Estamos usando a KroghInterpolator( ).
40. , é o subespaço de , constituído pelas funções que são quadráticas em cada 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. Observem a grande melhoria na aproximação por elementos finitos quadráticos com relação à brincadeira de trocar um gigante por vários anões. Incrível, é a mesma função do leão!
42. O programa da interpolação quadrática por partes. Ele é uma pequena modificação do anterior.
43. A continuação do programa.
44. Vamos começar com a interpolação polinomial por partes de Lagrange e de Hermite.
45. Mestres, vamos ver o que há mais na scipy.interpolate!
46. Assinalei outros dois métodos disponibilizados pela scipy.interpolate:
47. Um deles.
48. E o outro: a classe PiecewisePolynomial.
49. 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 ).
50. 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...
51. 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.
52. Para efetuar a interpolação de com o método ( ) 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.
53. A interpolação de Hermite. O programa a seguir sorteia 7 pontos e 7 direções tangentes e constrói a interpoladora.
54. 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.
55. O programa que traçou o gráfico da transparência anterior .
56. = = 3 O default = é uma polinomial de grau 3 em cada sub-intervalo.
57. 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.
58. O programa que traça o gráfico da função e da interpoladora, com = 3.
59. É só usar = 0. A interpoladora constante por partes usando a função ( ) do pacote interpolate da SciPy.
60. Agora usando = 1 temos a interpolação linear por partes.
61. Com = 2 temos interpolação quadrática em cada sub-intervalo. Neste caso a tangência só ocorre no extremo direito de cada subintervalo.
62. A interpolação de Hermite Especificando o parâmetro = 3 ou deixando = .
63. 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
64. 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ó.
65. Observe que agora são 3n condições para n pontos. Neste exemplo 12.
66. Agora vamos examinar outro tipo de interpolação polinomial por partes, conhecido como “splines”.
67. p() q() 2 2 3 1 Observem que, no extremo comum, (2 , 2 ), a “colagem” entre p e q não é suave. Matematicamente temos ′(2 ) ≠ ′(2 ). Na interpolação polinomial por partes, de Lagrange, a única exigência é a continuidade na "colagem” das interpoladoras p e q. No gráfico: 2 = 2 = (2 ).
68. 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. Já vimos que a interpolação 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.
69. p() q() 2 2 3 1 A imposição da condição ′′ 2 = ′′ 2 obriga a mesma curvatura local para as duas curvas. Nessa interpolação por partes de Hermite, são necessárias 3n informações para interpolar n pontos.
70. A interpolação polinomial spline, é obtida “forçando a barra”. Uma régua “spline”, muito usada antigamente (antes da computação gráfica) por desenhistas, projetistas e arquitetos para traçar curvas suaves.
71. Já explico. Suponha que pretendemos interpolar uma função ∶ , → ℝ usando uma spline () constituída por polinomiais quadráticas em cada subintervalo em cada um dos 4 subintervalos definidos pela partição = 0 < 1 < 2 < 3 < 4 = () 1 2 0 3 4 Como assim?
72. 1 2 0 3 4 2 1 3 4 Assim, nossa spline () será constituída por 4 polinomiais de grau 2: 1 = 0 + 1 + 2 2, 2 = 0 + 1 + 2 2, 3 = 0 + 1 + 2 2, 4 = 0 + 1 + 2 2. Ao todo 4 × 3 = 12 coeficientes a determinar.
73. () 1 2 0 3 4 1 2 0 3 4 Para aliviar a notação, vou escrever (0 ) = 0 , (1 ) = 1 , (2 ) = 2 , (3 ) = 3 , (4 ) = 4 , como mostro na figura.
74. Como queremos que spline () interpole a função f, teremos 8 condições (8 equações lineares) para as 12 incógnitas. Nelas já estão embutidas 3 condições de continuidade da spline () nos nós internos 1 2 0 3 4 1 2 3 4 1 2 0 3 4 1 0 = 0 1 1 = 1 2 1 = 1 2 2 = 2 , 3 2 = 2 , 3 3 = 3 , 4 3 = 3 , 4 4 = 4 . 8 equações
75. 1 2 3 4 E a continuidade para as derivadas da interpoladora nos nós internos fornece outras três. 1 2 0 3 4 1 ′ (1 ) = 2 ′ (1 ) 2 ′ (2 ) = 3 ′ (2 ) 3 ′ (3 ) = 4 ′ (3 )
76. Ao todo temos 11 condições para determinar as 12 incógnitas 0 , 1 , 2 , 0 , 1 , 2 , 0 , 1 , 2 , 0 , 1 , 2 . Falta uma condição, que pode ser escolhida impondo uma condição sobre o valor da derivada da última polinomial 4 () no ponto 4 = . Por exemplo, 4 ′′ 4 = ou 4 ′ 4 = , onde e são valores escolhidos.
77. Então é só resolver o sistema linear para obter a spline interpoladora.
78. Seja () uma função definida num intervalo [, ]. Uma função () é um spline de grau p, com nós nos pontos = 0 , < 1 , < ⋯ < = , quando satisfaz as condições: • Em cada subintervalo , +1 é uma polinomial de grau p, • () é de classe −1 , . () é um spline interpolante quando, além dessas duas condições, satisfaz também: • = , = 0,1, ⋯ , . A definição geral de “spline”. Acabamos de detalhar a construção de um spline quadrático.:
79. Uma “spline” cúbica interpolante: É uma função : [, ] → ℝ, de classe 2[, ] que, restrita a cada subintervalo , +1 , = 1,2, ⋯ , é uma polinomial de grau 3. De forma mais explícita, ela satisfaz: 1) = () para ∈ . +1 , k = 1,2, ⋯ , 2) = , = 1,2, ⋯ , 3) = +1 , = 1,2, ⋯ , − 1 4) ′ = ′+1 , = 1,2, ⋯ , − 1 5) ′′ = ′′+1 , = 1,2, ⋯ , − 1
80. O classe interp1d ...
81. O parâmetro kind permite escolhermos interpolação constante, linear, quadrática e cúbica por partes.
82. O resto dos parâmetros. Observem a indicação de outro método, mais recente.
83. Duas polinomiais do 2º grau (5 pontos).
84. Note que com 3 polinomiais do 2º grau (7 pontos), já temos uma aproximação razoável.
85. Com 5 polinomiais do 2º grau (9 pontos), a aproximação já é muito boa.
86. Esse é o código do programa. Marquei em vermelho as linhas de código envolvidas na definição da polinomial quadrática por partes.
87. Usando 2 polinomiais cúbicas, verificamos que o erro já é muito pequeno. O código é muito parecido.
88. Usando 3 polinomiais cúbicas, já não distinguimos, visualmente, a função da interpoladora por partes!
89. Agora vamos usar a classe dos splines da scipy.
90. Mais detalhes, que usaremos agora!
91. Esta função retorna a representação tck da spline cúbica (default).
92. Esta função avalia a spline representada por tck no array x.
93. Os gráficos da função = 2 e da spline cúbica passando por pontos escolhidos.
94. O programa. Marquei em vermelho as linhas de código envolvidas na definição e na avaliação da spline.
95. O parâmetro der da função splev( ) permite mostrar a derivada 1ª ou 2ª da spline interpolante.
96. O código do programa. Observe as chamadas a splev( )
97. A spline cúbica paramétrica, para curvas definidas parametricamente no plano x-y, p/ex.
98. A interpolação da borboleta com uma spline cúbica paramétrica.
99. Vejam, são necessários mais pontos para melhorar a precisão:
100. Este é o programa que interpola a borboleta.
101. Para interpolar a espiral de Arquimedes, basta mudar a definição da curva paramétrica!
102. Com 20 pontos a interpolação com a spline paramétrica já está boa:
103. Tchau! Até a próxima.