使用Numpy求解振子链本征向量结果异常的问题排查

问题描述 #

我尝试用Python的Numpy库求解并绘制振子链（弹簧）的本征向量，但结果和预期不符：运行代码得到的第一本征向量存在负值，而正确的第一本征向量（对应整体平移模式）应该全为非负值。

原代码 #

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()

问题排查与修正 #

1. 矩阵构造完全颠倒 #

振子链的动力学矩阵（描述每个振子受力的矩阵）是三对角矩阵，正确的赋值逻辑为：

• 对角元：每个振子受左右两个弹簧的拉力，值为 2*spring_tension/mass（即你的变量k）
• 次对角元（上下相邻位置）：每个振子受相邻振子的拉力，值为 -spring_tension/mass（即你的变量h）

原代码把对角元和次对角元的赋值完全搞反了，修正后的matrix_creation函数：

def matrix_creation(dimension,h,k,l):
    matrix=np.zeros((dimension,dimension))
    # 第一行：仅右邻居，对角元为k，次对角元为h
    matrix[0][0]=k
    matrix[0][1]=h
    # 最后一行：仅左邻居，对角元为k，次对角元为h
    matrix[dimension-1][dimension-1]=k
    matrix[dimension-1][dimension-2]=h
    # 中间行：左右均有邻居
    for j in range(1,dimension-1):
        matrix[j][j-1]=h
        matrix[j][j]=k
        matrix[j][j+1]=h
    return matrix

2. 本征值与本征向量未排序 #

np.linalg.eig返回的结果是无序的，默认不会按本征值从小到大排列。振子链的最低本征值（对应整体平移模式）才会对应全同符号的本征向量，需要先对本征值排序，再重新排列本征向量：

# 获取本征值排序后的索引
sorted_indices = np.argsort(eigenvalues)
# 按索引重新排列本征值和本征向量
sorted_eigenvalues = eigenvalues[sorted_indices]
sorted_eigenvectors = eigenvectors[:, sorted_indices]

3. 本征向量符号处理 #

本征向量的符号是任意的（若v是本征向量，则-v也是同一本征值的本征向量）。如果排序后的本征向量是全负的，只需乘以-1即可匹配预期的非负模式：

# 确保每个本征向量的首元素为正
for i in range(sorted_eigenvectors.shape[1]):
    if sorted_eigenvectors[0, i] < 0:
        sorted_eigenvectors[:, i] *= -1

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

def matrix_creation(dimension,h,k,l):
    matrix=np.zeros((dimension,dimension))
    # 第一行
    matrix[0][0]=k
    matrix[0][1]=h
    # 最后一行
    matrix[dimension-1][dimension-1]=k
    matrix[dimension-1][dimension-2]=h
    # 中间行
    for j in range(1,dimension-1):
        matrix[j][j-1]=h
        matrix[j][j]=k
        matrix[j][j+1]=h
    return matrix


def plot_eigenvectors(rest_position, eigenvectors, N):
    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)
    eigenvalues,eigenvectors=np.linalg.eig(matrix)
    
    # 排序本征值与本征向量
    sorted_indices = np.argsort(eigenvalues)
    sorted_eigenvectors = eigenvectors[:, sorted_indices]
    
    # 修正本征向量符号
    for i in range(sorted_eigenvectors.shape[1]):
        if sorted_eigenvectors[0, i] < 0:
            sorted_eigenvectors[:, i] *= -1
    
    rest_position=[a for a in range(1,N+1)]
    plot_eigenvectors(rest_position,sorted_eigenvectors,N)

if __name__=="__main__":
    main()

内容的提问来源于stack exchange，提问作者tristan ledet ledet