Interpolação por partes e splines

Estudo de interpolação polinomial por partes e splines com utilização do pacote scipy.interpolate.

Transcript

1. Interpolação por partes e splines Prof. Paulo R. G. Bordoni UFRJ
2. Recomeçamos com o exemplo terrível no qual o erro cresceu assustadoramente com o grau. Ele disparou nos extremos!
3. Dividir para conquistar! É um mantra em computação e uma ideia muito conhecida dos generais e políticos. E aí Mestres, o que faremos para conseguir mais precisão?
4. A ideia é muito simples: Começamos dividindo o intervalo em vários subintervalos. Depois, no lugar de uma polinomial de grau alto, no intervalo grande, usamos uma polinomial de grau baixo em cada subintervalo. Saboreiem a conquista na próxima transparência!
5. É claro que, quanto mais subdivisões, melhor a aproximação:
6. Esse é programa. Assinalei em vermelho as linhas de código envolvidas na construção da polinomial interpoladora.
7. Aliás é assim que a MatPlotLib engana nossos olhos. Vejam só:
8. No meu programa só uso a MatPlotLib para desenhar o gráfico de () em [0, 2] .
9. Passaremos a analisar uma classe de espaços vetoriais importantíssima em teoria de aproximação.
10. São os espaços das funções contínuas e lineares por partes. Veja a cara de uma função desse espaço:
11. O programa que gera esse gráfico.
12. No exemplo, usamos a uma polinomial interpoladora de grau 1 para cada subintervalo −1 , , = 1, ⋯ , . Dito de outra forma, em cada subintervalo −1 , , = 1, ⋯ , , temos uma retinha = + −1 −1 () onde = −−1 −−1 e −1 = − −1− .
13. 3 1 4 3 4 ℎ3 3 1 As funções chapéu, “hat functions”, formam uma base para esse espaço. Nesta figura mostro = 3 ℎ3 + 4 ℎ4
14. A cada dois subintervalos consecutivos [−1 , ], [ , +1 ], usamos a uma polinomial interpoladora de grau 2. Em particular, o número total de pontos deverá ser ímpar. Outra possibilidade é efetuar interpolação quadrática por partes.
15. Vamos descobrir o conteúdo do Tutorial da SciPy sobre interpolação: Já exploramos parte da scipy.linalg e da numpy.polinomials Agora vamos partir para a scipy.interpolate.
16. O início do Tutorial da SciPy sobre interpolação:
17. Eis a scipy.interpolate na Reference.
18. Há muita coisa aqui!
19. O Mestre já detalhou isto! Experimentem usar estas rotinas da scipy.interpolate.
20. Vamos pular esta parte.
21. Veremos splines mais adiante.
22. Usaremos parte disto mais adiante.
23. Vamos pular esta parte também.
24. Vamos começar com a interpolação polinomial por partes de Lagrange e de Hermite.
25. Interpolação de com a PiecewisePolynomial.
26. 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 ).
27. 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...
28. 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.
29. 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.
30. A interpolação de Hermite. O programa a seguir sorteia 7 pontos e 7 direções tangentes e constrói a interpoladora.
31. 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.
32. O programa que traçou o gráfico da transparência anterior .
33. = = 3 O default = é uma polinomial de grau 3 em cada sub-intervalo.
34. 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.
35. O programa que traça o gráfico da função e da interpoladora, com = 3.
36. É só usar = 0. A interpoladora constante por partes usando a função ( ) do pacote interpolate da SciPy.
37. Agora usando = 1 temos a interpolação linear por partes.
38. Com = 2 temos interpolação quadrática em cada sub-intervalo. Neste caso a tangência só ocorre no extremo direito de cada subintervalo.
39. A interpolação de Hermite Especificando o parâmetro = 3 ou deixando = .
40. 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
41. 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ó.
42. Observe que agora são 3n condições para n pontos. Neste exemplo 12.
43. Agora vamos examinar outro tipo de interpolação polinomial por partes, conhecido como “splines”.
44. 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 ).
45. 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.
46. 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.
47. 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.
48. 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?
49. 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 incógnitas .
50. () 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.
51. 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 incluídas 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 = 2 2 2 = 2 , 3 2 = 2 , 3 3 = 3 , 4 3 = 3 , 4 4 = 4 .
52. 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 )
53. 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 = 0.
54. Então é só resolver o sistema linear para obter a spline interpoladora.
55. 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.:
56. 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
57. O classe interp1d ...
58. O parâmetro kind permite escolhermos interpolação constante, linear, quadrática e cúbica por partes.
59. O resto dos parâmetros. Observem a indicação de outro método, mais recente.
60. Duas polinomiais do 2º grau (5 pontos).
61. Note que com 3 polinomiais do 2º grau (7 pontos), já temos uma aproximação razoável.
62. Com 5 polinomiais do 2º grau (9 pontos), a aproximação já é muito boa.
63. Esse é o código do programa. Marquei em vermelho as linhas de código envolvidas na definição da polinomial quadrática por partes.
64. Usando 2 polinomiais cúbicas, verificamos que o erro já é muito pequeno. O código é muito parecido.
65. Usando 3 polinomiais cúbicas, já não distinguimos, visualmente, a função da interpoladora por partes!
66. Na Scipy encontramos:
67. Esta função retorna a representação tck da spline cúbica (default).
68. Esta função avalia a spline representada por tck no array x.
69. Os gráficos da função = 2 e da spline cúbica passando por pontos escolhidos.
70. O programa. Marquei em vermelho as linhas de código envolvidas na definição e na avaliação da spline.
71. O parâmetro der da função splev( ) permite mostrar a derivada 1ª ou 2ª da spline interpolante.
72. O código do programa. Observe as chamadas a splev( )
73. A spline cúbica paramétrica, para curvas definidas parametricamente no plano x-y, p/ex.
74. A interpolação da borboleta com uma spline cúbica paramétrica.
75. Vejam, são necessários mais pontos para melhorar a precisão:
76. Este é o programa que interpola a borboleta.
77. Para interpolar a espiral de Arquimedes, basta mudar a definição da curva paramétrica!
78. Com 20 pontos a interpolação com a spline paramétrica já está boa:
79. Tchau! Até a próxima.