Points repartition and condition number

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.io as pio
pio.templates.default = "seaborn"
from scipy import interpolate
import warnings
warnings.filterwarnings('ignore')

The first form of the barycentric interpolation formula¶

Given a set of ⁠ $n + 1$ ⁠ nodes $\{ x_0, \, x_1, \, \ldots, \, x_n\}$ , the first form of the barycentric interpolation formula is defined by :

p_n(f)(x) = \left\{ \begin{aligned} &\Pi_{n+1}(x)\sum_{k=0}^n \frac{\omega_k}{x-x_k} f(x_k) \qquad{} & x\notin \{x_0,\ldots,x_n\} \\ &f(x_k) & x=x_k. \end{aligned} \right.

where the weights $w_k$ are defined by :

w_k = \frac{1}{\prod_{l\neq k} (x_k - x_l)} = \frac{1}{\Pi_{n+1}^\prime(x_k)}, \qquad{} k=0,\ldots,n,

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

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.

In this notebook, we use the first form of the barycentric interpolation formula, which has good stability properties.

Equidistant points¶

We interpolate the function $f(x)$ using a set equidistant points.

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])
    return cond / np.abs(pn)
    
# x use to display f and pn
xplot = np.linspace(xmin, xmax, 500)

# x use to display relative conditionning
eps = 0.05
xcond = np.concatenate((np.linspace(xmin+eps, -eps, 50), np.linspace(eps, xmax-eps, 50) ))

# maximal degree of Newton interpolating polynomial
n_max = 72

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

# compute for each degree xk, yk and pn
xk = []; yk = []; pn = []; cond = []
for i, ni in enumerate(n):
    xk.append(np.linspace(xmin, xmax, ni+1))
    yk.append(f(xk[i]))
    pn.append(improved_lag_interp(xk[i], yk[i], xplot))
    cond.append(compute_cond(xk[i], xk[i], improved_lag_interp(xk[i], yk[i], xcond), xcond))

# Create figure
marker=dict(size=5, symbol='x-thin', line=dict(width=1, color='rgb(76,114,176)'))
fig = make_subplots(rows=2, cols=1, vertical_spacing = 0.1,
                    subplot_titles=("Polynomial interpolation", "Relative local conditioning"))

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

# 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 points"))
    fig.add_trace(go.Scatter(visible=False, x=xcond, y=cond[i], mode='markers', showlegend=False, marker=marker), row=2, col=1)

# Make plot visible for n[0]
fig.data[1].visible = True
fig.data[2].visible = True
fig.data[3].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==3*i+1 or el==3*i+2 or el==3*i+3 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)
title = dict(text="Polynomial interpolation using equidistant points", y=0.98)

fig.update_layout(sliders=sliders, height=750, legend=legend, title=title, modebar=dict(orientation='v'))
fig.update_yaxes(range=[-1.1, 1.1], row=1)
fig.update_yaxes(type="log", exponentformat='e', row=2)
fig.update_xaxes(range=[-1.1, 1.1])

fig.show()

With this point distribution, significant errors become apparent as soon as the number of points exceeds 60. Since the algorithm used is stable, the errors observed are directly related to the local conditioning of the problem, which increases with the number of points.

Sensitivity to the perturbations

def plot_sol_pert(eps, n, xmin, xmax, xi, yi, yi_pert):

    xmin = ximin; xmax = ximax
    x = np.linspace(xmin, xmax, 1000)

    fig = go.Figure()

    fig.add_trace(go.Scatter(x=x, y=f(x), name="sin(x)"))
    fig.add_trace(go.Scatter(x=xi, y=yi, mode='markers', name="pts d'interp."))
    fig.add_trace(go.Scatter(x=x, y=improved_lag_interp(xi, yi, x), name="p(x)"))
    
    for i, epsi in enumerate(eps):
       fig.add_trace(go.Scatter(visible=False, x=xi, y=yi_pert[i], mode='markers', name=f"perturbed interp. pts"))
       fig.add_trace(go.Scatter(visible=False, x=x, y=improved_lag_interp(xi, yi_pert[i], x), name=f"p(x) with eps = {epsi}"))
   
    # Make plot visible for eps=1e-5
    fig.data[9].visible = True
    fig.data[10].visible = True
   
    # Create and add slider
    steps = []
    for i, epsi in enumerate(eps):
       step = dict(method="update", label = f" {epsi}",
                   args=[{"visible": [el==0 or el==1 or el==2 or el==2*i+3 or el==2*i+4 for el in range(len(fig.data))]}])
       steps.append(step)
           
    sliders = [dict(currentvalue={"prefix": "eps = "}, steps=steps, active=3)]
    legend = dict(orientation="h", x=0.1, y=1.14, bgcolor = 'rgba(0,0,0,0)')
    title = f"Perturbed polynomial interpolation of order {n} for sin(x)"
    
    fig.update_layout(sliders=sliders, legend=legend, title=title, height=500, modebar=dict(orientation='v'))
    
    fig.show()

eps = np.array([1e-2, 1e-3, 1e-4, 1e-5])
n = 25
ximin = -1.; ximax = 1.

xi = np.linspace(ximin, ximax, n+1)
yi = f(xi)

yi_pert = []
for i, epsi in enumerate(eps):
    yi_pert.append(f(xi))
    yi_pert[i][11] = yi_pert[i][11] + epsi

plot_sol_pert(eps, n, xmin, xmax, xi, yi, yi_pert)

Chebyshev points¶

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

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

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

# x use to display relative conditionning
eps = 0.05
xcond = np.concatenate((np.linspace(xmin+eps, -eps, 50), np.linspace(eps, xmax-eps, 50) ))
#xcond = np.linspace(xmin+0.1, xmax-0.1, 100)

# maximal degree of Newton interpolating polynomial
n_max = 520

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

# compute for each degree xk, yk and pn
xk = []; yk = []; pn = []; cond = []
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))
    cond.append(compute_cond(xk[i], xk[i], improved_lag_interp(xk[i], yk[i], xcond), xcond))

# Create figure
marker=dict(size=5, symbol='x-thin', line=dict(width=1, color='rgb(76,114,176)'))
fig = make_subplots(rows=2, cols=1, vertical_spacing = 0.1,
                    subplot_titles=("Polynomial interpolation", "Relative local conditioning"))

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

# 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 points"))
    fig.add_trace(go.Scatter(visible=False, x=xcond, y=cond[i], mode='markers', showlegend=False, marker=marker), row=2, col=1)

# Make plot visible for n[0]
fig.data[1].visible = True
fig.data[2].visible = True
fig.data[3].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==3*i+1 or el==3*i+2 or el==3*i+3 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)
title = dict(text="Polynomial interpolation using using Chebyshev points", y=0.98)
fig.update_layout(sliders=sliders, height=750, legend=legend, title=title, modebar=dict(orientation='v'))
fig.update_yaxes(range=[-1.1, 1.1], row=1)

fig.show()

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

Sensitivity to the perturbations

eps = np.array([1e-1, 1e-2, 1e-3, 1e-4, 1e-5])
n = 59
ximin = -1.; ximax = 1.

xi = cheb_points(xmin, xmax, n+1)
yi = f(xi)

yi_pert = []
for i, epsi in enumerate(eps):
    yi_pert.append(f(xi))
    yi_pert[i][49] = yi_pert[i][49] + epsi

plot_sol_pert(eps, n, xmin, xmax, xi, yi, yi_pert)