Newton's method for polynomial interpolation

Question

I've seen that in Newton's method for interpolating polynomials, the coefficients can be found algorithmically using (in Python-ish):

a = Y_DataPoints.copy()
m = length(X_DataPoints)
for k in range(1,m):
       a[k: m] = (a[k:m] - a[k-1]) / (X_Data[k:m] - X_Data[k-1]) 

But I don't really understand the model is subtracting all points past $k_{i}$ by $k_{i-1}$. It seems like you would only subtract $k_{i}$ by $k_{i-1}$, not the entire vector.

Can someone shed some light on this?

Answer 1

This is clever. I looked in Numerical Analysis by Kinkaid and Cheney, and did not see this algorithm in the section of divided differences. Here is how it works. We look for coefficients in the expansion $$f(x) = a_0+a_1(x-x_0) + a_2(x-x_0)(x-x_1)+a_3(x-x_0)(x-x_1)(x-x_2)+\dots \tag{1} $$ Clearly, $a_0=f(x_0)$. This is already recorded in our a array, so the zeroth entry there stays unchanged. To find $a_1$, subtract $a_0$ from both sides of (1) and divide the result by $x-x_0$: $$\frac{f(x)-a_0}{x-x_0} = a_1 + a_2 (x-x_1)+a_3 (x-x_1)(x-x_2)+\dots \tag{2}$$ But wait a moment. Relation (2) is the same as (1), only for a different function and with indices shifted. Great! The step $a_k\leftarrow \dfrac{a_k-a_0}{x_k-x_0}$ is the algorithmic version of replacing $f$ with the left side of (2), which we could call $f_1$. After that, the process repeats: $a_1$ is already found, and we move on to $$\frac{f_1(x)-a_1}{x-x_1} = a_2 +a_3 (x-x_2)+\dots \tag{3}$$ The step $a_k\leftarrow \dfrac{a_k-a_1}{x_k-x_1}$ amounts to replacing $f_1$ with the left side of (3). And so on.