3424.sagews

from pylab import *

#------------------------------------------#
#           Funcoes                        #
#------------------------------------------#

#funcao a ser interpolada
def f(x):
    return 1.0/(1+x**2)
    #return sin((x/2)**2)+exp(0.1*x)
    #return exp(x*sqrt(x+10))-exp(x**2)

#retorna valores para polinomos interpoladores 
#de grau 1 (para pontos igualmente espacados)
def p1(x,y,zs):
    r = zeros(size(zs))
    k = 0
    for z in zs:
        i = int((z-x[0])/(x[1]-x[0]))
        if(i<size(x)-1):
            d1   = (y[i+1]-y[i])/(x[i+1]-x[i])
            r[k] = y[i] + (z-x[i])*d1
        elif(i==size(x)-1):
            r[k] = y[i]
        k = k+1
    return r

#retorna valores para polinomos interpoladores 
#de grau 2 (para pontos igualmente espacados)
def p2(x,y,zs):
    r = zeros(size(zs))
    k = 0
    for z in zs:
        i = 2*int((z-x[0])/(x[2]-x[0]))
        if(i<size(x)-1):
            d11  = (y[i+1]-y[i])/(x[i+1]-x[i])
            d12  = (y[i+2]-y[i+1])/(x[i+2]-x[i+1])
            d2   = (d12-d11)/(x[i+2]-x[i])
            r[k] = y[i] + (z-x[i])*(d11 + (z-x[i+1])*d2)
        elif(i==size(x)-1):
            r[k] = y[i]
        k = k+1
    return r

#retorna valores para a spline cubica natual interpolante
#(para pontos igualmente espacados)
def spline(x,y,zs):
    g = calcg(x,y) #calculando os coefs g[i]
    r = zeros(size(zs))
    k = 0
    for z in zs:
        i = int((z-x[0])/(x[1]-x[0]))
        if(i<size(x)-1):
            a = (g[i+1]-g[i])/(6*(x[i+1]-x[i]))
            b = g[i+1]/2
            d = y[i+1]
            c = (y[i+1]-y[i])/(x[i+1]-x[i]) + (2*g[i+1]+g[i])*(x[i+1]-x[i])/6
            r[k]= d + c*(z-x[i+1]) + b*(z-x[i+1])**2 + a*(z-x[i+1])**3
        elif(i==size(x)-1):
            r[k]= y[i]        
        k = k+1
    return r

#calcula os coeficientes g[i] para o calculo da spline
def calcg(x,y):
    n = size(x)-1
    A = zeros((n+1)**2).reshape(n+1,n+1)
    b = zeros(n+1)
    A[0,0] = A[n,n] = 1
    b[0] = b[n] = 0
    for i in range(1,n):
        A[i,i-1] = (x[i]-x[i-1])
        A[i,i]   = 2*((x[i]-x[i-1])+(x[i+1]-x[i]))
        A[i,i+1] = (x[i+1]-x[i])
        b[i] = 6*((y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]))
    g = solve(A,b)
    return g        
    
#--------------------------------------------#
#        Inicio                              #
#--------------------------------------------#
    
a = -5                 #inicio do intervalo de interpolacao
b =  5                 #fim do intervalo de interpolacao
n =  6                 #numero de subintervalos (tem que ser par para p2)
x = linspace(a,b,n+1)  #gerando pontos usandos na interpolacao
y = f(x)               #gerando valores funcionais
method = spline        #metodo a ser usado (p1 ou p2 ou spline)

#plotanto as funcoes
clf()
xp = linspace(a,b,1000)
plot(xp,method(x,y,xp),"g-")
plot(x,y,"r.")
plot(xp,f(xp),"b-")
show()
savefig("interp.png")