# Enter data points to interpolate
# xpts = [2,-1,5,8]
# ypts = [3,4,7,-1]
xpts = [2,-1,5,8,-2,10,13]
ypts = [3,4,7,-1,0,0,6]
deg = len(xpts)-1
degp1 = deg + 1
degp2 = deg + 2
pts = []
for k in range(degp1):
pts.append((xpts[k],ypts[k]))
end
G = point(pts,size=50,color="purple")
show(G)
|
The Lagrange Polynomail determined using linear equations. The MATRIX APPROACH.
A = zero_matrix(degp1,degp2)
for row in range(degp1):
for col in range(degp1):
A[row,col] = xpts[row]^col
# print xpts[row]^col
end
A[row,degp1] = ypts[row]
end
show(A)
|
|
coefs = (A[:,0:degp1]).inverse()*A[:,degp1]
show(coefs)
|
|
def p(x):
y = 0
for k in range(degp1):
y += coefs[k]*x^k
end
return (x,y[0])
end
form = str(coefs[0])
for k in range(1,degp1):
form += ("$+%s$"%str(coefs[k])+"$x^%s$"%str(k))
end
m = min(xpts)
M = max(xpts)
n = 200
delx = (M-m)/n
poly = []
for k in range(n+1):
x = m+delx*k
poly.append( p(x=x) )
end
P = point(poly)
pretty_print(html(form))
show(G+P)
|
(36481/12474)
(36481/12474)
|
# Same thing but now evaluating using a Horner's Approach
def p(x):
y = coefs[deg]
for k in range(1,degp1):
y = y*x+coefs[deg-k]
end
return (x,y[0])
end
form = str(coefs[0])
for k in range(1,degp1):
form += ("$+%s$"%str(coefs[k])+"$x^%s$"%str(k))
end
m = min(xpts)
M = max(xpts)
n = 200
delx = (M-m)/n
poly = []
for k in range(n+1):
x = m+delx*k
poly.append( p(x=x) )
end
H = point(poly)
pretty_print(html(form))
show(G+H)
|
(36481/12474)
(36481/12474)
|
Now, use the basis polynomial approach
var('x')
p = 0
for k in range(degp1):
L=ypts[k]
for j in range(degp1):
if k!=j:
L = L*(x-xpts[j])/(xpts[k]-xpts[j])
end
end
p += L
show(L)
end
show(p)
H = plot(p,(x,m,M))
show(xpts)
show(G+H)
|
print xpts
print ypts
[2, -1, 5, 8, -2, 10, 13] [3, 4, 7, -1, 0, 0, 6] [2, -1, 5, 8, -2, 10, 13] [3, 4, 7, -1, 0, 0, 6] |
# Now, divided differences
dd = [ypts[0]]
for k in range(1, degp1):
term = ypts[k]
for j in range(k):
term = (term - dd[j]) / (xpts[k] - xpts[j])
dd.append(term)
end
print dd
[3, -1/3, 5/18, -17/162, -97/3240, 1289/285120, -1009/1995840] [3, -1/3, 5/18, -17/162, -97/3240, 1289/285120, -1009/1995840] |
# Using these divided difference coefficients and evaluating with a Horner's approach.
PolyD = dd[deg]
for k in range(0,degp1):
PolyD = PolyD*(x-xpts[deg-k])+dd[deg-k]
end
H = plot(PolyD,(x,m,M))
show(G+H)
show(xpts)
show(PolyD)
|
|
|