Stability of polynomial interpolation

def compute_divided_diff_coef(x, y):
    n = x.size
    coef = y.copy()
    for i in range(1,n):
        coef[i:] = (coef[i:] - coef[i-1:-1])/(x[i:] - x[:-i])
    return coef

def newton_interp(xi, yi, x):
    coef = compute_divided_diff_coef(xi, yi)
    n = xi.size
    p = np.zeros(x.size)
    for i in range(n-1,0,-1):
        p = (coef[i]+p) * (x-xi[i-1])
    p = p +  coef[0]
    return p

def improved_lag_interp(xk, yk, x):
    n = xk.size
    l = np.ones(x.size)
    p = np.zeros(x.size)
    
    wk = np.zeros(n) 
    for k in range(n):
        wk[k] = 1. / (np.prod( xk[k] - xk[np.arange(n)!=k] ))
    
    mask = np.isin(x, xk[:])
    p[np.where(mask)] = yk[np.where(np.isin(xk, x))]
    mask = np.invert(mask)    

    for k in range(n):
        l[mask] *= x[mask] - xk[k]
        p[mask] += (wk[k]/(x[mask] - xk[k])) * yk[k]
        
    return l*p

def cheb_points(xmin, xmax, n):
    x = np.zeros(n+1)
    for i in range(n+1):
        x[i] = (xmin+xmax)/2 + ((xmax-xmin)/2) * np.cos(((2*i+1)*np.pi)/(2*n + 2))
    return x

def noise(eps, x): 
    return eps * (2*np.random.rand(x.size)-1) 

def lagrange(k, xk, x):
    lag = np.ones(x.size)
    for j in range(xk.size):
        if (j!=k): lag *= (x-xk[j])/(xk[k]-xk[j])
    return lag

def compute_cond(xk, yk, pn, x):
    cond = np.zeros(x.size)
    for k in range(xk.size):
        cond += np.abs(lagrange(k, xk, x) * yk[k]) / np.abs(pn)
    return cond

Consider the function the polynomial interpolation of $f(x) = \sin(\pi x)$ on $[-1,1]$ .

def f(x):
    return np.sin(np.pi*x)

xmin = -1.; xmax = 1.

Newton interpolation¶

We interpolate the function $f(x)$ using Chebyshev points.

# x use to display f and pn
xplot = np.linspace(xmin, xmax, 500)

# maximal degree of Newton interpolating polynomial
n_max = 71

# array of degree
n = np.arange(19, n_max+1, 2)

# compute for each degree xi, yi and pn
xk = []; yk = []; pn = []
for i, ni in enumerate(n):
    xk.append(cheb_points(xmin, xmax, ni))
    yk.append(f(xk[i]))
    pn.append(newton_interp(xk[i], yk[i], xplot))
    
# Create figure
fig = go.Figure(layout_title="Polynomial interpolation using Chebyshev points")

# add f(x) plot
fig.add_trace(go.Scatter(visible=True, x=xplot, y=f(xplot), name="f(x)"))

# add yi and pn(x) invisible plots
for i, ni in enumerate(n):
    fig.add_trace(go.Scatter(visible=False, x=xplot, y=pn[i], name=f"p{ni}(x)"))
    fig.add_trace(go.Scatter(visible=False, x=xk[i], y=yk[i], mode='markers', name="interpolation pts"))

# Make plot visible for n[0]
fig.data[1].visible = True
fig.data[2].visible = True

# Create and add slider
steps = []
for i, ni in enumerate(n):
    step = dict(method="update", label = f" {ni+1}",
                args=[{"visible": [el==0 or el==2*i+1 or el==2*i+2 for el in range(len(fig.data))]}])
    steps.append(step)
        
sliders = [dict(currentvalue={"prefix": "nb points: "}, steps=steps)]
legend = dict(orientation="h", y=1.1)

fig.update_layout(sliders=sliders, legend=legend, height=500, modebar=dict(orientation='v'))
fig.update_yaxes(range=[-1.1, 1.1])

fig.show()

For this point distribution, the problem is well-conditioned formulation. The errors observed starting at 68 points are due to the stability of the interpolation method used.

First form of the barycentric interpolation formula¶

We interpolate the function $f(x)$ using Chebyshev points.

# maximal degree of Newton interpolating polynomial
n_max = 520

# array of degree
n = np.arange(19, n_max+1, 40)

# x use to display f and pn
xplot = np.linspace(xmin, xmax, 2000)

# compute for each degree xi, yi and pn
xk = []; yk = []; pn = []
for i, ni in enumerate(n):
    xk.append(cheb_points(xmin, xmax, ni))
    yk.append(f(xk[i]))
    pn.append(improved_lag_interp(xk[i], yk[i], xplot))
    
# Create figure
fig = go.Figure(layout_title="Polynomial interpolation using Chebyshev points")

# add f(x) plot
fig.add_trace(go.Scatter(visible=True, x=xplot, y=f(xplot), name="f(x)"))

# add yi and pn(x) invisible plots
for i, ni in enumerate(n):
    fig.add_trace(go.Scatter(visible=False, x=xplot, y=pn[i], name=f"p{ni}(x)"))
    fig.add_trace(go.Scatter(visible=False, x=xk[i], y=yk[i], mode='markers', name="interpolation pts"))

# Make plot visible for n[0]
fig.data[1].visible = True
fig.data[2].visible = True

# Create and add slider
steps = []
for i, ni in enumerate(n):
    step = dict(method="update", label = f" {ni+1}",
                args=[{"visible": [el==0 or el==2*i+1 or el==2*i+2 for el in range(len(fig.data))]}])
    steps.append(step)
        
sliders = [dict(currentvalue={"prefix": "nb points: "}, steps=steps)]
legend = dict(orientation="h", y=1.1)

fig.update_layout(sliders=sliders, legend=legend, height=500, modebar=dict(orientation='v'))
fig.update_yaxes(range=[-1.1, 1.1])

fig.show()

In this case, the number of points can be increased without introducing errors. This is the expected behavior when using a stable algorithm for a problem with good local conditioning.