First of all, sorry for my bad math terminology as it's not my native language and I may misuse some terms in English.
I've been tasked with writing an application which calculates the general formula of interpolating polynomial (being given a set of $X_i$ and $Y_i$), however, I have to use Neville's algorithm for that. I know Lagrange's method can be used, but I've been told to use Neville's algorithm specifically.
So, my question is – is it even possible? If so, how should I calculate the coefficients of the interpolating polynomial using this algorithm $(a_0, a_1,\dots ,a_n)$?
I feel like I'm missing a point somewhere (maybe because I'm not that good at maths), but I can't find a solution to this.
I'll be honest with it – this is homework. But I'm not asking for a solved code, just need help with understanding this.
Best Answer
I think you can calculate the coefficients of the interpolation polynomial as follows.
The algorithm given on the Wikipedia page http://en.wikipedia.org/wiki/Neville%27s_algorithm works not only for the value of the interpolation polynomial at a given point x, but also for the complete interpolation polynomial (in a mathematical sense, regarding the interpolation polynomial as an element of the polynomial ring $\mathbb R[X]$). More specifically, the initial polynomials $p_{i,i}(X) = y_i$ are constant. From there, we can calculate the intermediate polynomials $$p_{i,j}(X) = ((x_j - X) p_{i,j-1}(X) + (X - x_i) p_{i+1,j}(X))/(x_ j - x_i)$$until we reach the complete interpolation polynomial.
Now, all you have to do is to translate this recursive scheme for polynomials into a form suitable for mainstream programming languages :-). You can do this by identifying a polynomial with its coefficient vector. We know that, given n + 1 interpolation points $(x_i, y_i)$, the interpolation polynomial has degree n, and more generally, the polynomial $p_{i,j}$ from the recursive scheme with $0 \leq i < j \leq n$ has degree $j - i \leq n$. So, all the polynomials we have to consider, and their corresponding coefficient vectors are of the form $$p(X) = a_0 + a_1 X + \cdots + a_n X^n <-> (a_0, \ldots, a_n).$$ And of course, in a programming language, a coefficient vector is just an array :-). With this correspondence, we can translate the recursive scheme from above into a recursive scheme for coefficient vectors.
First, the initial polynomial $p_{i,i}(X) = y_i$ has coefficient vector $(y_i, 0, \ldots, 0)$.
Next, we have to figure out the coefficient vector of $Xp_{i,j-1}(X)$ given the coefficient vector of $p_{i,j-1}(X)$. But that's easy. Multiplying a polynomial by $X$ shifts its coefficient vector to the right. So, if we denote the coefficient vector of the polynomial $p(X)$ by $V(p(X))$ and on coefficient vectors define the shift operator $$S((a_0, \ldots, a_n)) = (0, a_0, \ldots, a_{n - 1}),$$ then we have $$V(Xp(X)) = S(V(p(X)).$$Check that :-).
Now we can rewrite the recursion for polynomials from above as a recursion for their coefficient vectors and get $$V(p_{i,j}(X)) = (x_jV(p_{i,j-1}(X)) - S(V(p_{i,j-1}(X))) + S(V(p_{i+1,j}(X))) - x_iV(p_{i+1,j}(X))) / (x_j - x_i).$$
Now you have a recursive scheme for coefficient vectors of polynomials which takes you from the coefficient vectors of the initial, constant polynomials $p_{i,i}$ all the way to the coefficient vector of the complete interpolation polynomial.