Precept 10: Barycentric interpolation, Runge's phenomenon and Gibb's phenomenon

In this precept, we experiment with polynomial interpolants.

Recall the formulas for Barycentric interpolation of points $\{x_i, y_i\}_{i=0}^{n}$:

$$ w_k = \frac{1}{\prod_{i=0 \dots n, i\neq k} (x_k - x_i)} $$$$ 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$.

The uniform and Chebysev points (of the second kind) on $[-1, 1]$ are defined as :

$$ x_i = -1 + 2 \frac{i}{n} $$$$ x_i = \cos \left( \frac{i\pi}{n} \right) $$

for $i = 0, \dots n$

Implement a polynomial interpolant of the following functions using Barycentric interpolation. Check that the point to evaluate is not of of the nodes, up to $10^{-15}$. If it is, return the corresponding value $y_i$.

• First write a function which computes the Barycentric weights in O(n^2)
• Then write a function which given the weights evaluate the polynomial at a single point in O(n)
• Then write a function which evaluates the interpolant at many points by looping through the function you just wrote.

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

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 pts = 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?

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)

errors = []
ns = np.arange(5,50)
f = lambda x: np.sin(2*np.pi *  x)
for n in ns:
    # Make convergence plot
    errors
    
plt.semilogy(ns,errors)

Test it on:

• $ f_3(x) = \begin{cases} -1 & \text{ if }x < 0 \\ 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.

n = 20
# TODO: Change n, see how result changes
f = lambda x: np.sign(x)

# scaling plot
errors = []
ns = np.arange(5,500,10)
for n in ns:
    # Make convergence plot
    errors

plt.semilogy(ns,errors)