Polynomial fitting with equal number of data points and coefficients

I am currently experimenting with polynomial fitting using jupyter. The function below returns the least-square polynomial of degree m given the data points in xs with corresponding ys.

from numpy import *
from matplotlib import pyplot as plt

def findas(m,xs,ys):
    A = array([[0]*(m+1)]*(m+1))
    b = array([0]*(m+1))
    for k in range(m+1):
        b[k] = sum(ys*xs**k)
        for i in range(m+1):
            A[k,i] = sum(xs**(k+i))
    coefs = linalg.solve(A,b)
    print(coefs)
    def fit(x):
        return sum(coefs*(x**array(range(len(coefs)))))
    return fit

Suppose I have the following six data points and fit a polynomial of degree 5:

xs = array([1,2,3,4,5,6])
ys = array([-5.21659 ,2.53152 ,2.05687 ,14.1135 ,20.9673 ,33.5652])
ft = findas(5,xs,ys)

From my understanding, the resulting curve should pass through every single data point exactly (in fact, the Lagrange polynomial should be the result).

xdense = arange(1,6.1,0.1)
ydense = [ft(x) for x in xdense]   

plt.plot(xdense,ydense)
plt.plot(xs,ys,'rx')
plt.show()

Sample output:

However, this is not the case. The curve is quite far off! What is going on here? Does this have something to do with round-off error? Thanks in advance!

Answer

It seems there was a truncation error! The block of code

A = array([[0]*(m+1)]*(m+1))
b = array([0]*(m+1))
for k in range(m+1):
...

should read:

A = array([[0.]*(m+1)]*(m+1))
b = array([0.]*(m+1))
for k in range(m+1):
...

i.e we have to specify the zeros as float.

Moreover, round-off errors can amplify in the process of inverting a matrix. This will particularly be the case when the eigenvalues of the matrix that we want to invert differ significantly in their order of magnitude.