Runge’s phenomenon - Numerical analysis

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

Consider the function $f(x) = \displaystyle \frac{1}{(1 + 25x^2)}$ .

def f(x):
    return 1/(1+25*x*x)

Interpolation using equidistant interpolation points.¶

We interpolate the function $f(x)$ using a set equidistant points with Newton’s method with Krogh interpolator from scipy : interpolate.KroghInterpolator.

xmin = -1.0
xmax = 1.0

# maximal degree of Newton interpolating polynomial
n_max = 30

# array of degree
arr_n = np.arange(1, n_max)

# compute for each degree xi, yi and pn
xi_equi = []; yi_equi = []; pn_equi = []
for i, ni in enumerate(arr_n):
    xi_equi.append(np.linspace(xmin, xmax, ni+1))
    yi_equi.append(f(xi_equi[i]))
    pn_equi.append(interpolate.KroghInterpolator(xi_equi[i], yi_equi[i]))
    
# Create figure
fig = go.Figure()

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

# 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(arr_n):
    fig.add_trace(go.Scatter(visible=False, x=xplot, y=pn_equi[i](xplot), name=f"p{ni}(x)"))
    fig.add_trace(go.Scatter(visible=False, x=xi_equi[i], y=yi_equi[i], mode='markers', name="interpolation points"))

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

# Create and add slider
steps = []
for i, ni in enumerate(arr_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)
title = "Polynomial interpolation of the Runge function using equidistant points"
fig.update_layout(sliders=sliders, legend=legend, title=title, height=500, modebar=dict(orientation='v'))
fig.update_xaxes(range=[-1.05, 1.05])
fig.update_yaxes(range=[-.25, 1.25])
fig.show()

As the number of interpolation points increases, we observe that the polynomial begins to oscillate with an ever-increasing amplitude. This is known as the Runge phenomenon.

Interpolation using Chebyshev points¶

We interpolate the function $f(x)$ using Chebyshev points points with Newton’s method.

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

# compute for each degree xi, yi and pn
xi_cheb = []; yi_cheb = []; pn_cheb = []
for i, ni in enumerate(arr_n):
    xi_cheb.append(cheb_points(xmin, xmax, ni))
    yi_cheb.append(f(xi_cheb[i]))
    pn_cheb.append(interpolate.KroghInterpolator(xi_cheb[i], yi_cheb[i]))

fig = make_subplots(rows=2, cols=1, vertical_spacing = 0.1,
                    subplot_titles=("Equidistant points ", "Chebyshev points"))

# add f(x) plot
fig.add_trace(go.Scatter(x=xplot, y=f(xplot), name="f(x)", marker=dict(color='rgb(76,114,176)')), row=1, col=1)
fig.add_trace(go.Scatter(x=xplot, y=f(xplot), name="f(x)", showlegend=False, marker=dict(color='rgb(76,114,176)')), row=2, col=1)

# add yi and pn(x) invisible plots
for i, ni in enumerate(arr_n):
    fig.add_trace(go.Scatter(visible=False, x=xplot, y=pn_equi[i](xplot), name=f"p{ni}(x)", 
                            marker=dict(color='rgb(221,132,82)')), row=1, col=1)
    fig.add_trace(go.Scatter(visible=False, x=xi_equi[i], y=yi_equi[i], mode='markers', name="interpolation points",
                             marker=dict(color='rgb(85,168,104)')), row=1, col=1)
    fig.add_trace(go.Scatter(visible=False, x=xplot, y=pn_cheb[i](xplot), showlegend=False,
                              marker=dict(color='rgb(221,132,82)')), row=2, col=1)
    fig.add_trace(go.Scatter(visible=False, x=xi_cheb[i], y=yi_cheb[i], mode='markers', showlegend=False,
                             marker=dict(color='rgb(85,168,104)')), row=2, col=1)
    
# Make plot visible for n=1
fig.data[2].visible = True
fig.data[3].visible = True
fig.data[4].visible = True
fig.data[5].visible = True

# Create and add slider
steps = []
for i, ni in enumerate(arr_n):
    step = dict(method="update", label = f" {ni+1}",
                args=[{"visible": [el==0 or el==1 or el==4*i+2 or el==4*i+3 or el==4*i+4 or el==4*i+5 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 of the Runge function ", y=0.98)
fig.update_layout(sliders=sliders, height=750, legend=legend, title=title, modebar=dict(orientation='v'))
fig.update_xaxes(range=[-1.05, 1.05])
fig.update_yaxes(range=[-.25, 1.25])
fig.show()

There is a noticeable improvement in the interpolation when using Chebyshev points.