I am trying to find and represent the eigenvectors of a chain of oscillators (springs) using numpy in python. The problem is that the obtained result is not exactly the expected one.
Here is my code:
import numpy as np
import matplotlib.pyplot as plt
def matrix_creation(dimension,h,k,l):
#This function creates the tridiagonal matrix
matrix=np.zeros((dimension,dimension)) #this specifies the size of the rows and columns
matrix[0][0]=h
matrix[0][1]=k
matrix[dimension-1][dimension-1]=l
matrix[dimension-1][dimension-2]=k
for j in range(1,dimension-1):
matrix[j][j-1]=h
matrix[j][j]=k
matrix[j][j+1]=l
return matrix
def plot_eigenvectors(rest_position, eigenvectors, N):
# Plot each eigenvector
for i in range(4):
plt.plot(rest_position, eigenvectors[:,i], label=f'Eigenvector {i+1}')
plt.title('Eigenvectors at time t=0')
plt.xlabel('Rest Position (x)')
plt.ylabel('Eigenvector Value (y)')
plt.legend()
plt.grid(True)
plt.show()
def main():
N=50
spring_tension=4
mass=1
h=-spring_tension/mass
k=2*spring_tension/mass
l=h
matrix=matrix_creation(N,h,k,l)
#print(np.linalg.eigvals(matrix))
#to acces al the information we can use the following comand
eigenvalues,eigenvectors=np.linalg.eig(matrix)
rest_position=[a for a in range(1,N+1)]
#plot_eigenvectors(rest_position,eigenvectors,N)
if __name__=="__main__":
main()
In the first image we can see the result of running my code.Representation of the first eigenvector using my code
It is obviously false since we should have any negative values. We can see an example of correct solution in image 2: Correct first eigenvector
This does not look like the first eigenvector of an oscillator, it looks like a higher order one …
and that is probably because np.linalg.eig
does not sort eigenvalues / eigenvectors.
I assume by "first" eigenvector, you mean sorted by increasing order of eigenvalues, in which case you might want to use
eigenvalues,eigenvectors=np.linalg.eig(matrix)
eigenvalue_order = np.argsort(eigenvalues)
sorted_eigenvalues = eigenvalues[eigenvalue_order]
sorted_eigenvectors = eigenvectors[:, eigenvalue_order]