Interpolação Polinomial

O problema da interpolação

Dados \(n+1\) pontos distintos \((x_i, y_i)\), busca-se uma função que passe exatamente por eles.

• Interpolação polinomial - busca-se \(P_n(x)\) de grau \(\leq n\):

\[P_n(x_i) = y_i, \qquad i = 0, 1, \ldots, n\]

Aplicações: reconstrução de tabelas (propriedades de materiais, dados experimentais), integração e diferenciação numéricas, projeto de trajetórias e curvas em CAD/CAM.

Aproximação polinomial: sistema de Vandermonde

Escrevendo \(P_n(x) = a_0 + a_1 x + \cdots + a_n x^n\) e impondo \(P_n(x_i) = y_i\) para cada nó, obtém-se um sistema linear nos coeficientes \(a_j\) — o sistema de Vandermonde:

\[\begin{bmatrix} 1 & x_0 & \cdots & x_0^n \\ 1 & x_1 & \cdots & x_1^n \\ \vdots & & & \vdots \\ 1 & x_n & \cdots & x_n^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix}\]

O determinante de Vandermonde é \(\prod_{i<j}(x_j - x_i)\): não nulo sempre que os nós são distintos — daí a existência e unicidade de \(P_n\).

Vandermonde em Python

import numpy as np

x = np.array([0., 1., 2., 4.])
y = np.array([1., 3., 2., 5.])

V = np.vander(x, increasing=True)   # matriz de Vandermonde
a = np.linalg.solve(V, y)           # coeficientes a0..an

print(a)                            # [a0, a1, a2, a3]
print(np.polynomial.polynomial.polyval(1.5, a))   # avalia P(1.5)

print(np.linalg.cond(V))            # número de condição — cresce rápido com n

Forma de Lagrange

Evita resolver o sistema linear: constrói \(P_n\) diretamente por combinação dos valores \(y_i\).

\[P_n(x) = \sum_{i=0}^{n} y_i\, L_i(x), \qquad L_i(x) = \prod_{\substack{j=0 \\ j\neq i}}^{n} \frac{x - x_j}{x_i - x_j}\]

Cada base \(L_i\) vale \(1\) em \(x_i\) e \(0\) nos demais nós — por isso \(P_n(x_i) = y_i\) é satisfeito automaticamente.

Lagrange em Python

import numpy as np
from scipy.interpolate import lagrange

x = np.array([0., 1., 2., 4.])
y = np.array([1., 3., 2., 5.])

P = lagrange(x, y)                  # objeto poly1d
print(P(1.5))

# implementação manual da base L_i(x)
def base_lagrange(x, i, t):
    termos = [(t - x[j])/(x[i] - x[j]) for j in range(len(x)) if j != i]
    return np.prod(termos, axis=0)

t = 1.5
Pt = sum(y[i]*base_lagrange(x, i, t) for i in range(len(x)))
print(Pt)

Erro de interpolação e o fenômeno de Runge

Se \(f\in C^{n+1}[a,b]\), o erro de \(P_n\) interpolando \(f\) é:

\[f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^{n}(x - x_i), \qquad \xi\in(a,b)\]

O produtório cresce próximo às extremidades do intervalo, especialmente com nós equiespaçados e grau alto — origem do fenômeno de Runge.

\[f(x) = \frac{1}{1 + 25x^2}, \qquad x\in[-1,1] \quad \text{(função de Runge)}\]

Runge: nós equiespaçados × Chebyshev

import numpy as np
import matplotlib.pyplot as plt

f = lambda x: 1 / (1 + 25*x**2)       # função de Runge
xx = np.linspace(-1, 1, 400)

# 11 nós equiespaçados
xe = np.linspace(-1, 1, 11)
ce = np.polyfit(xe, f(xe), 10)

# 11 nós de Chebyshev — suprimem a oscilação
k  = np.arange(11)
xc = np.cos((2*k + 1)*np.pi / 22)
cc = np.polyfit(xc, f(xc), 10)

plt.plot(xx, f(xx), 'k-', label='f(x)')
plt.plot(xx, np.polyval(ce, xx), label='equiespaçados')
plt.plot(xx, np.polyval(cc, xx), label='Chebyshev')
plt.plot(xe, f(xe), 'o'); plt.legend()

Interpolação por splines

Em vez de um polinômio de grau \(n\) para todos os pontos, usa-se um polinômio de grau baixo por subintervalo, unidos com continuidade nos nós.

Para o spline cúbico \(S(x)\), com \(n+1\) nós (\(n\) subintervalos), cada trecho \(S_i\) é um cúbico e exige-se, em cada nó interno:

\[S_i(x_i) = S_{i+1}(x_i), \qquad S_i'(x_i) = S_{i+1}'(x_i), \qquad S_i''(x_i) = S_{i+1}''(x_i)\]

Continuidade até a 2ª derivada \(\Rightarrow\) curva suave (\(C^2\)), sem os saltos de curvatura de uma interpolação linear por partes.

Splines em Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

f = lambda x: 1 / (1 + 25*x**2)      # função de Runge
x = np.linspace(-1, 1, 11)
y = f(x)

S = CubicSpline(x, y, bc_type='natural')

xx = np.linspace(-1, 1, 400)
plt.plot(xx, f(xx), 'k-', label='f(x)')
plt.plot(xx, S(xx), label='spline cúbico')     # sem oscilação de bordas
plt.plot(x, y, 'o'); plt.legend()

Aplicação em Engenharia: perfil de came

Projeto de cames (motores, automação industrial): define-se o deslocamento do seguidor em pontos de projeto ao longo da rotação, e a curva entre eles precisa ser suave.

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

# ângulo de came (graus) e deslocamento do seguidor (mm) — pontos de projeto
theta = np.array([0, 60, 120, 180, 240, 300, 360])
desloc = np.array([0, 8, 15, 15, 8, 0, 0])

S = CubicSpline(theta, desloc, bc_type='periodic')  # came fecha o ciclo em 360°

ang = np.linspace(0, 360, 400)
plt.plot(ang, S(ang), label='deslocamento')
plt.plot(ang, S(ang, 1), label="velocidade  S'")
plt.plot(ang, S(ang, 2), label='aceleração  S\'\'')
plt.plot(theta, desloc, 'o'); plt.legend()

bc_type='periodic' garante que a came se fecha suavemente após uma volta completa.

Exercícios

1. Construa \(P_3(x)\) para 4 pontos pelas formas de Vandermonde e de Lagrange e verifique que os polinômios coincidem.
2. Estime o número de condição da matriz de Vandermonde para \(n = 5, 10, 15\) nós equiespaçados e comente o crescimento.
3. Reproduza o fenômeno de Runge e compare nós equiespaçados, nós de Chebyshev e spline cúbico (gráfico Matplotlib único com as três curvas).
4. Para o exemplo do perfil de came, compare bc_type='natural' e 'periodic': qual evita descontinuidade de velocidade em \(\theta=0°/360°\)?
5. Interpole \(\sin x\) em \([0,\pi]\) com 5 nós igualmente espaçados; estime o erro máximo pela fórmula do erro e compare ao erro real.