3410.sagews

from pylab import *

#---------------------------------#
#           Funcoes               #
#---------------------------------#
#Funcao a ser interpolada.
def f(x):
   return 1.0/(1+x**2)
   #return sin(x)

#Algoritmo de Horner para avaliar
#o polinomio. 
def pol(x,d,n,z):
   r = d[n]
   for k in range(n-1,-1,-1):
      r = r*(z-x[k])+d[k]
   return r

#Algoritmo para calcular os coeficientes 
#do polinomio interpolador na forma de Newton.     
def newton(x,y,n):
   d = y.copy()
   for k in range(1,n+1):
      for i in range(n,k-1,-1):
         d[i] = (d[i]-d[i-1])/(x[i]-x[i-k])
   return d

#---------------------------------#
#           Inicio                #
#---------------------------------#
n = 6                     #grau do polinomio interpolador
a = -5                    #inicio do intervalo de pontos
b = 5                     #fim do intervalo de pontos
x = linspace(a,b,num=n+1) #gera pontos entre [a,b]
y = f(x)                  #calcula a funcao sobre os pontos
d = newton(x,y,n)         #calcula os coeficientes do polinomio

#faz os graficos da funcao
clf()                     #limpa a tela
xp = linspace(a,b,1000)   #gera mil pontos entre [a,b] para fazer grafico
plot(xp,f(xp))            #faz grafico da funcao sobre esses pontos
plot(xp,pol(x,d,n,xp))    #faz grafico do polinomio sobre os pontos
plot(x,y,'.')             #faz grafico dos pontos usados.
savefig('grafico.png')    #mostra a figura