Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt
x = np.linspace(-1, 1, 100)
pt.xlim([-1.2, 1.2])
pt.ylim([-1.2, 1.2])
for k in range(10): # crank up
pt.plot(x, np.cos(k*np.arccos(x)))
What if we interpolate random data?
n = 50 # crank up
k = n-1
i = np.arange(0, k+1)
x = np.linspace(-1, 1, 3000)
def f(x):
return np.cos(k*np.arccos(x))
nodes = np.cos(i/k*np.pi)
pt.plot(x, f(x))
pt.plot(nodes, f(nodes), "o")
[<matplotlib.lines.Line2D at 0x7ffad2e9f370>]
i = np.arange(1, n+1)
x = np.linspace(-1, 1, 3000)
def f(x):
return np.cos(n*np.arccos(x))
nodes = np.cos((2*i-1)/(2*n)*np.pi)
pt.plot(x, f(x))
pt.plot(nodes, f(nodes), "o")
[<matplotlib.lines.Line2D at 0x7ffad154d6c0>]
pt.plot(nodes, 0*nodes, "o")
[<matplotlib.lines.Line2D at 0x7ffad15ba230>]
n = 100
i = np.arange(n, dtype=np.float64)
nodes = np.cos((2*(i+1)-1)/(2*n)*np.pi)
V = np.cos(i*np.arccos(nodes.reshape(-1, 1)))
if 1:
# random data
data = np.random.randn(n)
else:
# Runge's example
data = 1/(1+25*nodes**2)
coeffs = la.solve(V, data)
x = np.linspace(-1, 1, 1000)
Vfull = np.cos(i*np.arccos(x.reshape(-1, 1)))
pt.plot(x, np.dot(Vfull, coeffs))
pt.plot(nodes, data, "o")
[<matplotlib.lines.Line2D at 0x7ffad13908b0>]
n = 10 # crank up
i = np.arange(n, dtype=np.float64)
nodes = np.cos((2*(i+1)-1)/(2*n)*np.pi)
V = np.cos(i*np.arccos(nodes.reshape(-1, 1)))
la.cond(V)
1.4142135623730963
Plot the product term from the estimate of truncation error in interpolation for the Chebyshev nodes: $$\left|\prod_{i=1}^n (x-x_i)\right| $$
def plot_err_prod(nodes, label):
eval_pts = np.linspace(-1, 1, 30000)
product = 1
for xi in nodes:
product = product*(eval_pts-xi)
pt.plot(eval_pts, np.abs(product), label=label)
n = 10 # crank up
i = np.arange(n, dtype=np.float64)
cheb_nodes = np.cos((2*(i+1)-1)/(2*n)*np.pi)
plot_err_prod(cheb_nodes, label="Chebyshev")
if 0:
nodes = np.linspace(-1, 1, n)
plot_err_prod(nodes, label="equispaced")
elif 0:
nodes = cheb_nodes.copy()
nodes[3] += 0.1
plot_err_prod(nodes, label="Perturbed")
pt.legend()
<matplotlib.legend.Legend at 0x7ffad052d6f0>