Conditioning of Derviative Matrices¶

Construct a matrix that takes a (first) centered difference of a periodic function on $[0,1]$. Call that matrix D.

In [3]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
In [80]:
npts = 100
In [81]:
x = np.linspace(0, 1, npts, endpoint=False)

h = x[1] - x[0]
D = (
    np.roll(np.diag(np.ones(npts)), -1, axis=0)
    -
    np.roll(np.diag(np.ones(npts)), 1, axis=0)
)/(2*h)

D
Out[81]:
array([[  0.,  50.,   0., ...,   0.,   0., -50.],
       [-50.,   0.,  50., ...,   0.,   0.,   0.],
       [  0., -50.,   0., ...,   0.,   0.,   0.],
       ..., 
       [  0.,   0.,   0., ...,   0.,  50.,   0.],
       [  0.,   0.,   0., ..., -50.,   0.,  50.],
       [ 50.,   0.,   0., ...,   0., -50.,   0.]])

Or: an alternate matrix that's based on a global polynomial

In [70]:
a = np.arange(npts, dtype=np.float64)

# Chebyshev nodes
nodes = np.cos((2*(a+1)-1)/(2*npts)*np.pi)

x = nodes

Vdm = np.empty((npts, npts))
Vdm_deriv = np.zeros((npts, npts))
for i in range(npts):
    Vdm[:, i] = np.cos(i*np.arccos(x))
    Vdm_deriv[:, i] = (i*np.sin(i*np.arccos(x)))/np.sqrt(1-x**2)
        
D = Vdm_deriv @ la.inv(Vdm)

Test that the matrix actually takes derivatives:

In [82]:
alpha = 3

f = np.sin(alpha * 2*np.pi*x)
df = alpha*2*np.pi*np.cos(alpha * 2*np.pi*x)

print(la.norm(df - D@f, np.inf))

0.111424462966
In [83]:
plt.plot(x, D@f- df)
plt.plot(x, df)
Out[83]:
[<matplotlib.lines.Line2D at 0x7feddbf12278>]

Investigate the norm of this matrix.

In [84]:
la.norm(D, np.inf)
Out[84]:
100.0

What function gets amplified like that?

In [85]:
amp_func = np.ones(npts)
amp_func[2::4] = -1
amp_func[3::4] = -1
la.norm(D@amp_func, np.inf)
Out[85]:
100.0
In [86]:
plt.plot(amp_func)
plt.plot(D@amp_func)
Out[86]:
[<matplotlib.lines.Line2D at 0x7fede076bba8>]

Now, what's the conditioning of this matrix like?

In [77]:
la.cond(D)
Out[77]:
2.0564138575587231e+18

OK, this may not be completely fair. But: can look at the spectrum:

In [87]:
eigv = la.eigvals(D)

plt.plot(eigv.real, eigv.imag, "o")
plt.xlim([-2, 2])
#plt.ylim([-2, 2])
plt.grid()
In [ ]: