Lagrange interpolation is a useful polynomial interpolation formula that lets us fit to an arbitrary function. While we might expect that increasing the number of points would decrease the average error of our polynomial approximation, this intuition fails because of Runge's phenomenon. If we pick our points carefully, though, we can overcome this phenomenon and still achieve a great approximation.

Specifically, in this post, we're going to approximate the function \(\frac{1}{1 + x^6}\) over \([0, 5]\) by using two different sets of points: a set of equidistant points across the interval vs. Chebyshev nodes. We'll compare the average error (MSE) and the maximum error over the interval.

## Helper functions

We'll add a function that returns the Lagrange interpolating polynomial for a set of points:

from operator import mul from functools import reduce def lagrange(pts, f): """ Return a Lagrange interpolating polynomial for f at the points pts. """ def basis(others): # w_{i}(x) = (x - x_{1}) * ... * (x - x_{i-1}) * (x - x_{i+1}) * ... return lambda x : reduce(mul, [x - o for o in others]) # Get each basis polynomial and coefficients bases = [basis(pts[:i] + pts[i + 1:]) for i in range(len(pts))] # a_{i} = f(x_{i}) / w_{i}(x_{i}) weights = [f(pt) / bases[i](pt) for i, pt in enumerate(pts)] return lambda x : sum([weights[i] * bases[i](x) for i in range(len(pts))])

This is a basic implementation of Lagrange interpolation - nothing special here.

Next, we'll add a couple of helper functions to compute the average and maximum error:

def average_error(p, f): # Get the MSE of p over 1000 points on the interval [0, 5] N = 1000 error = 0 for pt in [5 * i / N for i in range(N)]: error += (p(pt) - f(pt))**2 return error / N def max_error(p, f): error = float("-inf") for pt in [5 * i / 1000 for i in range(1000)]: curr_error = abs(p(pt) - f(pt)) if curr_error > error: error = curr_error return curr_error

## Computing for equidistant vs. Chebyshev points

Finally, we'll add our function to approximate:

def g(x): return 1 / (1 + x**6)

And some functions to compute and evaluate the interpolating polynomial for N points, whether equidistant or Chebyshev:

def evaluate(pts, f, tag): print("Interpolating for {0}".format(tag)) poly = lagrange(pts, f) print("Average error: {0}".format(average_error(poly, f))) print("Max error: {0}".format(max_error(poly, f))) def equidistant(n): # Interpolate for n points from [0, 5] pts = [5 * i / n for i in range(n)] evaluate(pts, g, "{0} equidistant points".format(n)) import math def chebyshev(n): # Interpolate for n points from [0, 5] pts = [(5 + 5 * math.cos((2 * i + 1) * math.pi / (2 * n + 2))) / 2 for i in range(n)] evaluate(pts, g, "{0} Chebyshev points".format(n)) n_values = [3, 5, 10, 20, 30] for n in n_values: equidistant(n) print("") chebyshev(n) print("")

## Evaluation

Here's what the results show:

Mean Squared Error:

Number of points | Equidistant | Chebyshev |

3 | 0.083449016830 | 0.123504372195 |

5 | 0.082739268475 | 0.08960932097 |

10 | 0.64207073475 | 0.0127952042780 |

20 | 26.4328328320 | 0.0005982536404 |

30 | 435.529425177 | 9.3148785757e-06 |

Maximum Error:

Number of points | Equidistant | Chebyshev |

3 | 0.86442582784 | 0.0030308889852 |

5 | 1.34512664557 | 0.0063971895249 |

10 | 6.2268209092 | 0.000121486335635 |

20 | 61.205287745 | 2.1722048813e-05 |

30 | 312.690373601 | 4.1559874016e-06 |

I initially wanted to graph these results, but as you can see, the error for the Chebyshev nodes is so low that it would hardly be visible.

All in all, clearly it's preferable to use Chebyshev nodes. This is especially clear as the number of points increases. I hope this was helpful!