I am playing around with spectral properties of differential operators. To get a feel for things I decided to start out with computing the eigenvalues and eigenvectors of the 1-D Laplacian with periodic boundary conditions
Lap =
[[-2, 1, 0, 0, ..., 1],
[ 1,-2, 1, 0, ..., 0],
[ 0, 1,-2, 1, ..., 0],
...
...
[ 0, 0, ..., 1,-2, 1],
[ 1, 0, ..., 0, 1,-2]]
So I run the following
import numpy as np
import scipy.linalg as scilin
N = 12
Lap = np.zeros((N, N))
for i in range(N):
Lap[i, i] = -2
Lap[i, (i+1)%N] = 1
Lap[i, (i-1)%N] = 1
eigvals, eigvecs = scilin.eigh(Lap)
where
> print(eigvals)
[-4.00000000e+00 -3.73205081e+00 -3.73205081e+00 -3.00000000e+00
-3.00000000e+00 -2.00000000e+00 -2.00000000e+00 -1.00000000e+00
-1.00000000e+00 -2.67949192e-01 -2.67949192e-01 9.43689571e-16]
which is what I expect. However I decide to verify that these eigenvalues and eigenvectors are correct. What I end up with is
> (Lap - eigvals[0]*np.identity(N)).dot(eigvecs[0])
array([ 0.28544445, 0.69044928, 0.83039882, 0.03466493, -0.79854101,
-0.81598463, -0.78119579, -0.7445237 , -0.769496 , -0.79741997,
-1.09625463, -0.69683007])
I expect to get the zero vector. So what is going on here?
As mentioned in the comment by @Warren, eigenvectors are columns of eigvecs. While in numpy indexing, eigvecs[0] represent first row of eigvecs. To fix it:
print((Lap-eigvals[0]*np.eye(N))@eigvecs[:,0])
[-6.66133815e-16 2.55351296e-15 -1.77635684e-15 1.11022302e-16
5.55111512e-16 -2.22044605e-16 -3.66373598e-15 -4.44089210e-16
7.77156117e-16 -1.11022302e-16 -1.66533454e-15 2.22044605e-15]
Which is basically all 0 (the numbers are there due to precision issue)