Interpolating $5$ points using a $5$th degree polynomial

Question

I am attempting to interpolate a set of $5$ data points $$(1, 1), (2, 3), (3, 5), (4, 7), (5, 69)$$ using a fifth degree polynomial.

So for some $$p(x)=ax^5+bx^4+cx^3+dx^2+ex^1+f$$ I can plug in the points given initially and find the matrix representation of the coefficients, which gives \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1\\\\ 32 & 16 & 8 & 4 & 2 & 1 & 3\\\\ 243 & 81 & 27 & 9 & 3 & 1 & 5\\\\ 1024 & 256 & 64 & 16 & 4 & 1 & 7\\\\ 3125 & 625 & 125 & 25 & 5 & 1 & 69 \end{bmatrix} If we rref this, we get \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \frac1{120} & \frac{59}{120}\\\\ 0 & 1 & 0 & 0 & 0 & -\frac18 & -\frac{39}8\\\\ 0 & 0 & 1 & 0 & 0 & \frac{17}{24} & \frac{403}{24}\\\\ 0 & 0 & 0 & 1 & 0 & -\frac{15}8 & -\frac{185}8\\\\ 0 & 0 & 0 & 0 & 1 & \frac{137}{60} & \frac{703}{60} \end{bmatrix} It is an system of equations with one extra variable so there is no singular solution. I attempted to solve for $p(x)$ in terms of the remaining coefficient $f=t$, which gives me $$p(x)=\left(\frac{59}{120}-\frac{t}{120}\right)x^{5}+\left(-\frac{39}{8}+\frac{t}{8}\right)x^{4}+\left(\frac{403}{24}-\frac{17t}{24}\right)x^{3}+\left(-\frac{185}{8}+\frac{15t}{8}\right)x^{2}+\left(\frac{703}{60}-\frac{173}{60}t\right)x+t$$ From my logic, this polynomial should work for all $t$, but it only works for $t=0$ (which can be seen by plugging in one of the points given initially and then solving for $t$).

Furthermore, if I plug in $t=59$, this should return the correct fourth degree polynomial to interpolate all the points, but the desmos graph says otherwise. (https://www.desmos.com/calculator/lu0rabmtoj)

What did I do wrong, or am I misunderstanding something here?

Answer 1

I recreated it in a very naive way in python and didn't have any issue, here is a graph for $t = -5, 0, 5, 10, 59$ and it looks like the theory works out perfectly.

EDIT: In your case you flipped the digits of the linear term, it should be 137 instead of 173, see here: https://www.desmos.com/calculator/rjtz8x69ho

import numpy as np
from sympy import Matrix
import matplotlib.pyplot as plt
a = np.array([1,2,3,4,5])
b = np.array([1,3,5,7,69])
exps = np.array([5,4,3,2,1,0])
B = np.array(Matrix(np.column_stack([a**k for k in exps]+[b])).rref()[0]).astype(float)
X = np.arange(0, 6, 0.001)
plt.plot(a, b, 'o')
for t in [-5, 0, 5, 10, 59]:
    poly = np.concatenate([np.squeeze(B[:, -1] - t*B[:, -2]), np.array([t])], axis=0)
    Y = ((X[:, None] ** exps) * poly) . sum(axis=1)
    plt.plot(X, Y, ':')
    plt.ylim([-5, 70])
plt.show()