381 - A7 - Interpolating Polynomials

218 days ago by Professor381

# 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)