Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Caderno 10

Objetivos de aprendizagem

  • Explicar o conceito de diferenças divididas;

  • Calcular diferenças divididas a partir de um conjunto de pontos e montar o polinômio interpolador de Newton correspondente;

  • Implementar algoritmo para gerar as funções de base e o polinômio interpolado de Newton

  • Aplicar interpolação polinomial a dados originados de várias áreas do conhecimento;

Em Busca do Desconhecido: Interpolação Polinomial

Interpolação Linear

Excluindo-se o caso de função constante, a forma mais simples de interpolação é ligar dois pontos dados com uma reta. Usando semelhança de triângulos entre nós e valores de função, obtemos

f1(x)f(x0)xx0=f(x1)f(x0)x1x0\dfrac{f_1(x) − f(x_0)}{x − x_0} = \dfrac{f(x_1) − f(x_0)}{x_1 − x_0}

a qual pode ser reorganizada para fornecer

f1(x)=f(x0)+f(x1)f(x0)x1x0(xx0)(1)f_1(x) = f(x_0) + \dfrac{f(x_1) − f(x_0)}{x_1 − x_0} (x − x_0) \qquad (1)

A notação f1(x)f_1(x) indica que esse é um polinômio interpolador de primeiro grau. Observe que, além de representar a inclinação da reta ligando os pontos, o termo [f(x1)f(x0)]/(x1x0)[f(x_1) − f(x_0)]/(x_1 − x_0) é uma aproximação por diferenças divididas da primeira derivada.

Interpolação Quadrática

Com três pontos, a interpolação quadrática é obtida a partir de

f2(x)=b0+b1(xx0)+b2(xx0)(xx1)(2)f_2(x) = b_0 + b_1(x − x_0) + b_2(x − x_0)(x − x_1) \qquad (2)

Um procedimento simples pode ser usado para determinar os valores dos coeficientes.

Para b0b_0, a Equação (2) com x=x0x = x_0 pode ser usada para calcular

b0=f(x0)(3)b_0 = f(x_0) \qquad (3)

A Equação (3) pode ser substituída na Equação (2), a qual pode ser calculada em x=x1x = x_1 para

b1=f(x1)f(x0)x1x0(4)b_1 = \dfrac{f(x_1) − f(x_0)}{x_1 − x_0} \qquad (4)

Finalmente, as Equações (3) e (4) podem ser substituídas na Equação (2), a qual pode ser calculada em x=x2x = x_2 e resolvida (depois de algumas manipulações algébricas) por

b2=f(x2)f(x1)x2x1f(x1)f(x0)x1x0x2x0(5)b_2 = \dfrac{\dfrac{f(x_2) − f(x_1)}{x_2 − x_1} − \dfrac{f(x_1) − f(x_0)}{x_1 − x_0}}{x_2 − x_0} \qquad (5)

Forma Geral do Interpolante de Newton

A análise anterior pode ser generalizada para ajustar um polinômio de grau nn a n+1n + 1 pontos dados. O polinômio de grau nn é

fn(x)=b0+b1(xx0)++bn(xx0)(xx1)(xxn1)(6)f_n(x) = b_0 + b_1(x − x_0) + \dots + b_n(x − x_0)(x − x_1) \dots (x − x_{n − 1}) \qquad (6)

Como foi feito anteriormente com as interpolações linear e quadrática, os pontos dados podem ser usados para calcular os coeficientes b0b_0, b1b_1, \dots , bnb_n. Para um polinômio de grau nn, n+1n + 1 pontos dados são necessários: (x0,f(x0)),(x1,f(x1)),,(xn,f(xn))(x_0, f(x_0)), (x_1, f(x_1)), \dots , (x_n, f(x_n)). Usamos esses pontos dados e as seguintes equações para calcular os coeficientes

b0=f(x0)(7)b1=f[x1,x0](8)b2=f[x2,x1,x0](9)bn=f[xn,xn1,,x1,x0](10)b_0 = f(x_0) \qquad (7) \\ b_1 = f[x_1, x_0] \qquad (8) \\ b_2 = f[x_2, x_1, x_0] \qquad (9) \\ \vdots \\ b_n = f[x_n, x_{n - 1}, \dots , x_1, x_0] \qquad (10)

onde as expressões entre colchetes correspondem a diferenças divididas.

Diferenças divididas

A primeira diferença dividida finita é representada em geral por

f[xi,xj]=f(xi)f(xj)xixjf[x_i, x_j] = \dfrac{f(x_i) − f(x_j)}{x_i − x_j}

A segunda diferença dividida finita, que representa a diferença das duas primeiras diferenças divididas, é expressa em geral por

f[xi,xj,xk]=f[xi,xj]f[xj,xk]xixkf[x_i, x_j, x_k] = \dfrac{f[x_i , x_j] − f[x_j, x_k]}{x_i − x_k}

Analogamente, a nn-ésima diferença dividida é dada por

f[xn,xn1,,x1,x0]=f[xn,xn1,,x1]f[xn1,xn2,,x0]xnx0f[x_n, x_{n − 1}, \dots , x_1, x_0] = \dfrac{f[x_n, x_{n − 1}, \dots , x_1] − f[x_{n − 1}, x_{n − 2}, \dots , x_0]}{x_n − x_0}

Funções de Base de Newton

Código gerador de funções de base de Newton de grau nn por computação simbólica.

from sympy import Symbol

def symbolic_vector(n,var):
    """Cria uma lista com n variáveis simbólicas.
    
        entrada:
            n: numero de pontos
            var: uma string (ex. 'x')
        saida:
            V: ['var0','var1',...,'varn-1']
    """
    if not isinstance(var,str):
        raise TypeError("{0} must be a string.".format(var))                 
        
    V = [Symbol(var  + str(i)) for i in range(0,n)]
        
    return V

def N_nj(X,j):
    """ Calcula a função de base de Newton N_{n,j}(x). 
        
        entrada:
            X: um vetor contendo variáveis simbólicas
    """
    # pega a variavel base passada e converte para simbólica
    x = X[1]
    x = str(x)
    x = Symbol(x[0:-1])
    N = x/x;
    if j > 0: 
        for k in range(0,j):
            N *= (x - X[k])
            
    return N
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# número de nós de interpolação: interpolação de (n-1)-ésimo grau
n = 5

# domínio de interpolação
x0,x1 = -1,1

# constroi vetor simbolico em x
X = symbolic_vector(n,'x')

# constroi pontos numericos 
xp = np.linspace(x0,x1,num=n,endpoint=True)

# cria malha numérica
xv = np.linspace(x0,x1)

# matriz das funções
Y = np.zeros((n,len(xv)))
for i in range(0,n):
    Y[i,] = np.zeros(np.shape(xv))

# montagem de dict para substituição: [xk,x0,x1,x2,...]
k = [str(i) for i in X]
k.insert(0,'x')

# preenche matriz
for i in range(0,Y.shape[0]):
    for j in range(0,np.size(xv)):
        v = list(np.concatenate([np.asarray([xv[j]]),xp]))
        d = dict(zip(k,v))
        Y[i,j] = N_nj(X,i).subs(d)

# plotagem das funções

# nós 
plt.scatter(xp,np.zeros(xp.shape),c= 'k')

leg = []
for i in range(0,Y.shape[0]):
    plt.plot(xv,Y[i,])
    s = '$N_{' + str(i) + '}(x)$'
    leg.append(s)
        
plt.grid()        
plt.legend(leg,loc='best')
plt.title('Funções de base de Newton até ordem ' + str(n-1) + ' em ['+str(x0)+','+str(x1)+']');
<Figure size 640x400 with 1 Axes>

Exemplo:

Encontre o polinômio interpolador de Newton de ordem 1 P1(x)P_1(x) para a tabela abaixo

xy
-14
01

Compute o valor de P1(0.35)P_1(0.35).

import numpy as np
import matplotlib.pyplot as plt

# interpolação linear

# coeficiente a0 = y0
# coeficiente a1 = (y1-y0)/(x1-x0)

# pontos
x0,y0 = -1,4
x1,y1 =  0,1

# ordem 0
a0 = y0

# interpolador de Newton
a1 = (y1-y0)/(x1-x0)
P1 = lambda x: a0 + a1*(x-x0)

# ponto interpolado
xp = -0.35
yp = P1(xp)
yp

# plotagem

# nós
plt.plot([x0,x1],[0,0],'ok')

# valores nodais
plt.plot([x0,x1],[y0,y1],'ok')

# interpolador
x = np.linspace(x0,x1,30,endpoint=True)
plt.plot(x,P1(x),label='$P_1(x)$')

# ponto interpolado
plt.plot(xp,0,'sr')
plt.plot(xp,yp,'sr')
plt.axvline(xp,0,yp,c='r',ls='dashed')

plt.grid()
plt.legend(loc='best');
<Figure size 640x400 with 1 Axes>
P1(xp)
2.05

Exemplo:

Encontre o polinômio interpolador de Newton de ordem 2 P2(x)P_2(x) para a tabela abaixo

xy
-14
01
2-1

Compute o valor de P2(0.35)P_2(0.35).

# interpolação quadrática

# Usando tabela DD: 
# https://vnicius.github.io/numbiosis/interpolador-newton/index.html


# par adicional 
x2,y2 = 2.,-1.

# coeficiente a2 = f[x0,x1,x2] = ( f[x1,x2] - f[x0,x1] ) / (x2 - x0)
# a2 = ( (y2-y1)/(x2-x1) - (y1-y0)/(x1-x0) )/(x2-x0)*(xx-x0)*(xx-x1)

# interpolador de Newton
P2 = lambda xx: P1(xx) + ( ( (y2-y1)/(x2-x1) - (y1-y0)/(x1-x0) )/(x2-x0) )*(xx-x0)*(xx-x1)


# ponto interpolado
yp = P2(xp)
yp

# plotagem

# nós
plt.plot([x0,x1,x2],[0,0,0],'ok')

# valores nodais
plt.plot([x0,x1,x2],[y0,y1,y2],'ok')

# interpolador
x = np.linspace(x0,x2,30,endpoint=True)
plt.plot(x,P2(x),label='$P_2(x)$')

# ponto interpolado
plt.plot(xp,0,'sr')
plt.plot(xp,yp,'sr')
plt.axvline(xp,0,yp,c='r',ls='dashed')

plt.grid()
plt.legend(loc='best');
<Figure size 640x400 with 1 Axes>
P2(xp)
1.8983333333333332

Exemplo:

Encontre o polinômio interpolador de Newton de ordem 3 P3(x)P_3(x) para a tabela abaixo

xy
-14
01
2-1
31

Compute o valor de P3(0.35)P_3(0.35).

# interpolação quadrática

# Usando tabela DD: 
# https://vnicius.github.io/numbiosis/interpolador-newton/index.html


# par adicional 
x3,y3 = 3.,1.

# coeficiente a3 = f[x0,x1,x2,x3] 

# interpolador de Newton
P3 = lambda xxx: P2(xxx) + 1/12*(xxx-x0)*(xxx-x1)*(xxx-x2)


# ponto interpolado
yp = P3(xp)
yp

# plotagem

# nós
plt.plot([x0,x1,x2,x3],[0,0,0,0],'ok')

# valores nodais
plt.plot([x0,x1,x2,x3],[y0,y1,y2,y3],'ok')

# interpolador
x = np.linspace(x0,x3,30,endpoint=True)
plt.plot(x,P3(x),label='$P_3(x)$')

# ponto interpolado
plt.plot(xp,0,'sr')
plt.plot(xp,yp,'sr')
plt.axvline(xp,0,yp,c='r',ls='dashed')

plt.grid()
plt.legend(loc='best');
<Figure size 640x400 with 1 Axes>
P3(-0.35)
1.9428854166666665
# ponto interpolado
YP = [P1(xp),P2(xp),P3(xp)]

# plotagem
# nós
plt.plot([x0,x1,x2,x3],[0,0,0,0],'ok')

# valores nodais
plt.plot([x0,x1,x2,x3],[y0,y1,y2,y3],'ok')

# interpoladores
x = np.linspace(x0,x3,30,endpoint=True)
plt.plot(x,P1(x),'b',label='$P_1(x)$')
plt.plot(x,P2(x),'c',label='$P_2(x)$')
plt.plot(x,P3(x),'g',label='$P_3(x)$')

# ponto interpolado
plt.plot(xp,0,'sr')
plt.plot([xp,xp,xp],YP,'sr')
plt.axvline(xp,0,max(YP),c='r',ls='dashed')

plt.grid()
plt.legend(loc='best');
<Figure size 640x400 with 1 Axes>

Comparação (zoom)

# interpoladores
x = np.linspace(-0.4,-0.3,30,endpoint=True)
plt.plot(x,P1(x),'b',label='$P_1(x)$')
plt.plot(x,P2(x),'c',label='$P_2(x)$')
plt.plot(x,P3(x),'g',label='$P_3(x)$')

# ponto interpolado
plt.plot([xp,xp,xp],YP,'sr')
plt.axvline(xp,0,max(YP),c='r',ls='dashed')

plt.grid()
plt.legend(loc='best');
<Figure size 640x400 with 1 Axes>
import sympy as sym
from sympy.abc import x

P3x = 512 + 63*(x-2015) - 2.34*(x-2015)*(x-2017) - 2.45*(x-2015)*(x-2017)*(x-2020)
P3xn = sym.lambdify(x, P3x) 
x = np.array([2015,2017,2020,2022])
y = np.array([512,638,792,700])

xn = np.linspace(2015,2022,num=100)

plt.plot(x,y,'o')
plt.plot(xn,P3xn(xn))

xx = 2021
yy = P3xn(xx)
plt.plot(xx,yy,'sg')
plt.axvline(x=xx)
print(f'Censo em {xx} = {yy}')
Censo em 2021 = 775.0400000000081
<Figure size 640x400 with 1 Axes>