Newton’s Divided-difference and Lagrange interpolating polynomials Python implementation
This blog was inspired when I’m studying for my ENGG 414 (Numerical Methods) exam.
Maybe you are doing an experiment or you have datasets that want to explore, but upon reviewing your datasets you noticed that it missing some datapoints that might be important for whatever you are doing. So, what should you do? Should you panic? Absolutely no, because Newton and Lagrange got your back.
All you need to do is to interpolate.
Interpolation is a process of estimating intermediate values between precise data points. The most common method used for this purpose is polynomial interpolation.
Using interpolation, you can now estimate datapoints that might be missing from your data. But how exactly are you going to do that?
Newton’s Divided-difference and Lagrange interpolating polynomials provide a simple and easy algorithm that can be implemented in a computer.
Newton’s Divided-difference interpolating polynomial
Newton polynomial can be written as
enough of the crazy part let’s go to implementation!
Note: The implementation was done in Python 3 so we could also visualize some graphs using python libraries like matplotlib, pandas and numpy. Jupyter notebook was used to write these codes.
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
This is the algorithm for Newton’s polynomial. It accepts an array of x and y values, and the value (xi) you want to interpolate. This function returned a multidimensional array of the interpolated values for different orders and their respective relative errors that would be used for visualization. (But you can modify the code to return an array of interpolated values in different orders if you don’t need to visualize)
def newton_interpolation(x, y, xi):
#length/number of datapoints
n = len(x)
#divided difference initialization
fdd = [[None for x in range(n)] for x in range(n)]
#f(X) values at different degrees
yint = [None for x in range(n)]
#error value
ea = [None for x in range(n)]
#finding divided difference
for i in range(n):
fdd[i][0] = y[i]
for j in range(1,n):
for i in range(n-j):
fdd[i][j] = (fdd[i+1][j-1] - fdd[i][j-1])/(x[i+j]-x[i])
#just printing dd here
fdd_table = pd.DataFrame(fdd)
print(fdd_table)
#interpolating xi
xterm = 1
yint[0] = fdd[0][0]
for order in range(1, n):
xterm = xterm * (xi - x[order-1])
yint2 = yint[order-1] + fdd[0][order]*xterm
ea[order-1] = yint2 - yint[order-1]
yint[order] = yint2
#returning a map for pandas dataframe
return map(lambda yy, ee : [yy, ee], yint, ea)
For our sample problem, we will try to estimate the value of ln(2) using the given datapoints:
x = [1,4,6, 5, 3, 1.5, 2.5, 3.5]
y = [0, 1.3862944, 1.7917595, 1.6094379, 1.0986123 , 0.4054641, 0.9162907, 1.2527630]
We can use our newton_interpolation function to solve for ln(2)
a = newton_interpolation(x, y, 2)
df = pd.DataFrame(a, columns=['f(x)','error'])
and it will show the following data frame:
And now we have estimated the value of ln(2) using our eight datapoints. The highest order will yield the nearest approximation which is ln(2) = 0.69348.
We can visualize these polynomial functions of different order and compare it to the actual function of f(x) = ln(x).
It’s amusing, isn’t it?
Now let’s hear more about Lagrange and his approach to the interpolation world.
Lagrange interpolating polynomial
Given a set of k + 1 data points
where no two xj are the same, the interpolation polynomial in the Lagrange form is a linear combination
of Lagrange basis polynomials
Geez, what was that? Okay, enough of foreplay let’s get to the more exciting part.
Below is the implementation of that crazy math expression of Lagrange polynomial. The function takes an array of x and y values and the value you want to interpolate, xx. Compared to the Newton’s polynomial function, this is set up to compute a single nth-order prediction, where n + 1 is the number of data points.
def lagrange_interpolation(x,y,xx):
n = len(x)
sum = 0
for i in range(n):
product = y[i]
for j in range(n):
if i != j:
product = product*(xx - x[j])/(x[i]-x[j])
sum = sum + product
return sum
Using the same datapoints from the latter approach, we could also solve for the estimated value of ln(2) using Lagrange polynomial and graph each value of different order.
plt.scatter(x,y)
xx = np.linspace(0, 8, 100)
yy = [lagrange_interpolation(x,y,i) for i in xx]xo= np.linspace(0.01,8,200)
plt.plot(xo, np.log(xo), label="f(x)=lnx", color="black")plt.grid(True, which='both')
plt.axvline(x=0, color='k')
plt.axhline(y=0, color='k')
plt.plot(xx,yy)fig = plt.figure()for i in range(1,len(x)):
g1 = fig.add_subplot(len(x)/2,2,i)
g1.plot(xx, lagrange_interpolation(x[i-1:],y[i-1:],xx))
g1.scatter(x,y)
plt.grid(True, which='both')
plt.axvline(x=0, color='k')
plt.axhline(y=0, color='k')
The above code will give us these beautiful graphs.
I was on the edge of my seat seeing these graphs!
Now,
Which of the two should you use?
For that question, let me quote Mark from math.stackexchange.com
If you want to have an easy formula for the remainder of the interpolation then it is much better to work with Newton’s method. Another advantage is that if you found the interpolation polynomial in the points x0, x1,…,xn and then you want to add the point xn+1 then using Newton’s method you can use the polynomial in the points x0,…,xn to easily find the polynomial for the points x0,…,xn+1. You don’t have to find all the coefficients all over again.
On the other hand, if you want to do Numerical differentiation or Numerical integration then it is much easier to work with Lagrange’s method. So having different ways to calculate the interpolation polynomial gives a lot more options. It all depends on what exactly you need to do.
In terms of the algorithm complexity, both algorithms run on O(n²) time but Newton’s method need extra O(n²) for space.
Have fun interpolating!
References
Chapra, S. C., & Canale, R. P. (2015). Numerical Methods for Engineers (7th ed.). New York: McGraw-Hill Education.
Lagrange polynomial. (2020, May). Retrieved from Wikipedia: https://en.wikipedia.org/wiki/Lagrange_polynomial
Mathematics. (2018). Retrieved from StackExchange: https://math.stackexchange.com/questions/2884888/what-is-difference-between-newton-interpolation-and-lagrange-interpolation
Newton polynomial. (2020, May). Retrieved from Wikipedia: https://en.wikipedia.org/wiki/Newton_polynomial