Numpy eigh-function with proper signs

Question

I would like to calculate the eigenvectors of a second-derivative-matrix in python. According to math the first vector should be equal to a sin-function between 0 and pi, the second equal to a sin-function between 0 and 2*pi. Thus my code looks like

import numpy as np
from matplotlib import pyplot as plt
from scipy import sparse
import scipy.integrate as integrate
import scipy.special as special
import scipy

def create_second_deriv(size, h):
    delta_matrix_2_second_diff = (np.eye(size, dtype=np.float)*-2+np.eye(size, k=-1,dtype=np.float)*1+np.eye(size, k=1,dtype=np.float)*1)
    delta_matrix_2_second_diff /= (h*h)
    return -1*delta_matrix_2_second_diff

delta_x = 0.001
x = np.linspace(0, 1, (int)(1/delta_x))
delta_matrix = create_second_deriv(len(x), delta_x)
w, v = scipy.linalg.eigh(delta_matrix)

plt.plot(v.tolist()[0])
plt.show()
plt.plot(v.tolist()[1])
plt.show()

Now, what I get as output, is as plot for the first eigenvector, and as plot for the second eigenvector. I already know that the signs for the different values are arbitrary, but in my case they are important for later processing. Is there a way to "flip" the signs such that the resulting values are approximately equal to the expected functions? Simply using the abs()-function will not help in that case.

Answer 1

There is a gotcha with the scipy.linalg.eigh() function:

Be very careful using the eigh() routine to get eigenvalues, though. There's a "gotcha" buried in it.

The syntax is:

(eigenvalues, eigenvectors) = eigh(matrix)

This returns an array of eigenvalues and a 2D array of eigenvectors (each eigenvector consists of many components).

And there's the gotcha. Let's say you want the nth eigenvalue and eigenvector. I would write:

eigenvalues[n]
eigenvectors[n]

And I would be horribly wrong. The eigenvectors and eigenvalues do share an index, but the index on the eigenvector is SECOND column:

eigenvectors[:,n]

Thus, the last four lines of your code must be changed to:

plt.plot(v[:,0])
plt.show()
plt.plot(v[:,1])
plt.show()