Let's say you have a vector space (which you want to interpolate in) $\Phi$ of dimension $n$.
Then, the nodes would be $n$ points $(\xi_1, \ldots, \xi_n)$ in the domain of your basis functions which satisfy the following: for arbitrary given values $(y_1, \ldots, y_n)$, you can find a unique $f \in \Phi$ such that
$$
f(\xi_i) = y_i \qquad \forall i=1,\ldots,n.
$$
This property is called unisolvency of the points $(\xi_i)$ for the space $\Phi$. You can call any set which is unisolvent a set of interpolation nodes, since you can do interpolation in them. This is a general concept, no matter how the space looks like.
Knots are something which is particular to the way splines are constructed. For a sequence of knots, $(t_1, \ldots, t_m)$, a spline is a function which is polynomial when restricted to each nonempty knot span $(t_i, t_{i+1})$ and satisfies some additional continuity assumptions in the knots.
So knots and nodes are not the same. For a spline space, however, it is common to choose the interpolation nodes as certain averages of the knots. See for instance the definition of Gréville abscissae, which are just particular interpolation nodes, here.
However, choosing the Gréville abscissae is not the only possible choice. In fact, the Schönberg-Whitney theorem states that, if you choose the $i$-th interpolation node $\xi_i$ to lie in the interior of the support of the $i$-th B-spline function, then the resulting nodes $(x_i)$ are unisolvent.
There is no theoretical reason for the error, as that should fall like $(6\sqrt2/n)^{n+1}/(4(n+1))$. The observed error is thus due to cancellation errors and other floating point problems in the construction of the interpolation polynomial.
With a different, more compact computation, the deviation only starts to become visible at $32-35$, the error at $40+1$ nodes is $1/10$ of the error reported in the question.
q = poly(x);
dq = polyder(q)
sum=0;
for i=1:length(x)
sum = sum + y(i)*deconv(q,poly(x(i)))/polyval(dq,x(i))
end
One can get some nodes better by computing the interpolation polynomial of the error and add that to the polynomial, and then iterate that. This helps up to $41$. One might perhaps get to higher node densities using the Newton interpolating scheme, starting outside on both interval ends and moving alternatingly inwards. Or use some dedicated divide-and-conquer approach.
Best Answer
What is great with Newton's interpolation is the fact that if you add new points you don't have to re-calculate all the coefficients (see forward divided difference formula) which can be really useful ! You'll have more details at the end of this article : https://en.wikipedia.org/wiki/Newton_polynomial