Barycentric interpolation, Runge's phenomenon and Gibb's phenomenon¶

We experiment with polynomial interpolants.

One can rewrite the Lagrange interpolating formula in an astonishing form:

$$ p_n(x) = \frac{\sum_{k=0}^n \frac{w_k y_k}{(x-x_k)}}{\sum_{k=0}^{n} \frac{w_k}{(x - x_k)}} $$ for $x \neq x_i$, where $\{x_i, y_i\}_{i=0}^{n}$ interpolation of points and

$$ w_k = \frac{1}{\prod_{i=0 \dots n, i\neq k} (x_k - x_i)} $$

are called the interpolating weights. (Yes, $p_n(x)$ is a polynomial! and it exactly equal to formula you derived in HW2).

The uniform and Chebysev points on $[-1, 1]$ are defined as :

$$ x_i = -1 + 2 \frac{i}{n} $$

$$ x_i = \cos(\frac{i\pi}{n}) $$

for $i = 0, \dots n$

Implement a polynomial interpolant of the following functions using Barycentric interpolation. Check that the point to evaluate is not one of the nodes $x_i$, or (that the distance to one is less than $10^{-15}$). If it is, return the corresponding value $y_i$.

Test it on uniform points and the function $f_1(x) = \sin(2\pi x)$ below.

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
def compute_barycentric_weights(x):
    # Compute barycenter weights 
    w = np.ones_like(x)

    return w


def eval_polynomial_at_one_z(w,x,y,z):
    return z


def eval_polynomial_at_many_z(w,x,y,z):
    # Evaluate polynomial interpolant p_n(x) at many points
    pnz = np.zeros_like(z)
    return pnz

def estimate_max_error(f1, f2):
    # Approximate the l_infty difference using 1e4 points.
    dense_grid = np.linspace(-1,1,10000)
    return np.max(np.abs(f1(dense_grid) - f2(dense_grid)))


n = 10
f1 = lambda x: np.sin(2*np.pi *  x)
uniform_grid = np.linspace(-1,1,n)
f1x = f1(uniform_grid)

# Plot the function and its interpolant on the fine grid
w = compute_barycentric_weights(uniform_grid)
dense_grid = np.linspace(-1,1,10000)
pn1x = eval_polynomial_at_many_z(w,uniform_grid,f1x,dense_grid)

# Plot the function and its interpolant
plt.plot(dense_grid,pn1x)
plt.plot(dense_grid,f1(dense_grid))

pn_x = lambda dense_grid: eval_polynomial_at_many_z(w,uniform_grid,f1x, dense_grid)
max_error = estimate_max_error(f1,pn_x)
# Should be around 0.08 for n = 10
print("max error:", max_error)

Test your code on:

  • $ f_1(x) = \cos( 10 \pi x) $
  • Runge's function $ f_2(x) = \frac{1}{ 1 + 25 x^2 } $

At uniform, Chebyshev and random point. You can get random points on $[-1,1]$ by running np.random.rand(n)*2 -1 .

Visualize the functions are their interpolant for n = 5, 10, 30.

Then, make a convergence plot of the max error as n increases from n = 1 to n = 50. To do this, approximate and plot the error between the interpolant and the function $\max_x |f(x) - p_n(x)| $. You can approximate the max by taking the max over a fine uniform grid with $n = 10^4$ points (see code above).

Do the interpolant converge for $f_1(x)$ ? What about $f_2(x)$?

Dos the interpolant at random points behave like the uniform interpolant or like the Chebyshev interpolant?

In [ ]:
n = 30
# TODO: Change the function, , see how result changes
f = lambda x: np.sin(2*np.pi *  x)
# f = lambda x: 1 / ( 1 + 25 * x**2)

# TODO: Change the grid, see how result changes
grid = np.linspace(-1,1,n)


fx = f(grid)

# Plot the function and its interpolant on the fine grid
w = compute_barycentric_weights(grid)
dense_grid = np.linspace(-1,1,10000)
pn1x = eval_polynomial_at_many_z(w,grid,fx,dense_grid)

# Plot the function and its interpolant
plt.plot(dense_grid,pn1x)
plt.plot(dense_grid,f(dense_grid))

pn_x = lambda dense_grid: eval_polynomial_at_many_z(w,grid,fx, dense_grid)
max_error = estimate_max_error(f1,pn_x)
# Should be around 0.08 for n = 10
print("max error:", max_error)
In [ ]:
errors = []
ns = np.arange(5,50)
f = lambda x: np.sin(2*np.pi *  x)
for n in ns:
    # Make convergence plot


plt.semilogy(ns,errors)

Test it on:

  • $ f_3(x) = \begin{cases} 1 & \text{ if }x < 1 \\ 1 & \text{ else} \end{cases} $

Plot the the Chebyshev interpolant of $f_3$ for different values of n. How does the oscillation near $0$ change with $n$? Make a convergence plot for n = 5 to 500 in increment of 10.

This tendency of polynomial interpolant to oscillate near discontinuities is called Gibb's phenomenon.

In [ ]:
n = 20
# TODO: Change the function, , see how result changes
f = lambda x: np.sign(x)
In [ ]:
# scaling plot
errors = []
ns = np.arange(5,500,10)

plt.semilogy(ns,errors)