Newton’s polynomial interpolation in three iteration lines using Python

Newton’s interpolation or divided differences method commonly refers to the algorithm that obtains the interpolating polynomial function of a set of points of the form (x, f (x)). Interpolation in this context means adjusting a curve given by a polynomial of degree n, so that it contemplates the points based on the degree required for a better approximation of the set. The divided difference operator handles, as its name suggests, the ratio between two values. In the numerator is the difference between the values ​​of f calculated for a point and the previous point, and the difference of the values ​​of x in the denominator are associated with the same two points, respectively. In this way we can write the following equation that relates what was said

Therefore, from the table given below, it appears that the divided differences take the following form

Importance of divided differences

From them we can obtain the factors that will be included in the construction of a polynomial p(x) that interpolates the set of points, and only one item — first in each column — in the hierarchy of divided differences is used in the calculation, as follows

In light green, the factors used in the construction of the interpolator polynomial are highlighted

If we think about this problem, we will see that for a set with n points the scalability of the calculations can become quite exhaustive without the computational approach. Thus, when analyzing the behavior of the divided differences, it is noted that for each column that we advance the number of items in the next column is reduced by one unit, while the values ​​of y and x follow a logic defined by what was discussed earlier.

The idea is to create a list that stores each column of divided differences into sublists. In other words, a list of sublists will be obtained, where the first sublist includes n-1 differences divided for the first level — column of f(x) — and the last one presents only one element, relative to the iterations that resulted from the combination of only two divided differences.

The problem ends when we obtain these sublists and are able to access the first element of each one for the construction of p(x). After a lot of scribbling and analysis of indexes, we have that the following loop solves this step

where deltas is the list with the previously mentioned sublists and coord is the list with ordered pairs (x, y). It’s worth mentioning that the value of k associated with the second for in the scope is the least trivial to be positioned since at the same time that it’s necessary to iterate for each list when adding the calculated items, it’s also necessary to go through empty lists of deltas and fill them taking into account the change in length of the sublists. This is aimed at avoiding errors such as

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: list index out of range

Generic approach for the same set of points

If we substitute the numerical values ​​of the previous table, we will notice that each coordinate assumes the following indices for x and y, already the uppercase deltas represent in its index the position of the sublist within the list that stores the delta columns— deltas, while the second index value of lowercase deltas refers to the item’s position within the deltas sublist for the same column.

so we get the following interpretation for the program input and the output given by printing the deltas list

coord = [(x0, y0),(x1, y1), (x2, y2), (x3, y3)]
deltas = [delta0, delta1, delta2]
# Before iterations, deltas is still a list with empty sublists

For the table presented with the numeric values ​​replaced, if we apply this to the right side of p(x), being

we’ll have the following

simplifying

To improve visualization, with the help of Mathematica, if we type the following lines into any .nb file

plot = Plot[x^3-3*x^2+2*x,{x,0,5.5}];
points = Graphics[{PointSize[Large],Red,Point[{{1,0},{3,6},{4,24}, {5,60}}]}];
Show[plot,points]

we arrive that p(x) coherently approaches the four points studied

Conclusion 🤔

Curve fitting is an efficient way to describe the behavior of certain studied phenomena. From the moment we make an analysis that requires graphic completion by adjusting a set of points, it is important to know the artifice presented in order to proceed to more advanced stages.

If you want to access the interpolation algorithms following the methods of Newton and Lagrange access my repository on GitHub.

There you will see that the algorithms require the use of two libraries that need to be installed via pip (SymPy and Matplotlib), but some instructions have been left to facilitate the execution of the project.