In [33]:
import numpy as np
In [34]:
x = np.array([1,2,6,11])
y = np.array([-3,-29,-33,-90])
In [35]:
# construct the Lagrange basis polynomials
nplus1 = len(x)
xx = np.linspace(0,12,400)  # evaluation points

l = []
for k in range(nplus1):
    # construct l_k evaluated at xx
    lk = np.ones_like(xx)
    for j in range(nplus1):
        if j!=k:
            lk *= (xx-x[j])/(x[k]-x[j])
    l.append(lk)
In [36]:
type(l),len(l)
Out[36]:
(list, 4)
In [39]:
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
plt.axhline(0,color='k',alpha=.25)
#plt.plot(x,y,'o')
plt.plot(x,np.zeros_like(x),'ko')
for xk in x:
    plt.axvline(xk,alpha=0.2)
for lk in l:
    plt.plot(xx,100*lk)
plt.savefig('lagrange_polys.pdf')
In [40]:
# now construct the interpolating polynomial

p = np.zeros_like(xx)
for k in range(nplus1):
    p += y[k]*l[k]
    
# and plot it
plt.plot(x,y,'o')
#plt.plot(x,np.zeros_like(x),'ko')
for xk in x:
    plt.axvline(xk,alpha=0.2)
plt.plot(xx,p)
plt.savefig('lagrange_polys.pdf')

monomial basis instead¶

In [28]:
A = np.empty((nplus1,nplus1))
for k in range(nplus1):
    # write the kth column
    A[:,k] = x**k   # RHS is b_k
np.array(A,dtype=int)
Out[28]:
array([[   1,    1,    1,    1],
       [   1,    2,    4,    8],
       [   1,    6,   36,  216],
       [   1,   11,  121, 1331]])
In [29]:
c = np.linalg.solve(A,y)
c
Out[29]:
array([ 40.38666667, -53.31111111,  10.54      ,  -0.61555556])
In [31]:
# construct p
p = np.zeros_like(xx)
for k in range(nplus1):
    p += c[k]*xx**k
# and plot it
plt.plot(x,y,'o')
#plt.plot(x,np.zeros_like(x),'ko')
for xk in x:
    plt.axvline(xk,alpha=0.2)
plt.plot(xx,p)
plt.savefig('interpolant_using_monomial_basis.pdf')
In [ ]: