为何首个特征向量未正确归一化？附Python数值计算代码

我编写了一个求解矩阵特征值与特征向量的函数，用于处理一维薛定谔方程，在函数末尾对特征向量进行归一化操作。为验证归一化效果，我检查了特征向量的正交归一性，发现其他特征向量均正常（结果为0或1），但首个特征向量evt_test[:, 0]未被正确归一化，相关代码如下：

import numpy as np
import math
from scipy import sparse
from scipy.sparse import linalg as sla


def harmonic_potential(x, pert):
    V_pert = []
    if pert == False:
       return x**2
    else:
        for i in x: 
            V_pert.append(( i**2 + 4*math.exp(-10*i**2) ) )
        return V_pert

def schrodinger1D(xmin, xmax, Nx, neigs, pert=False):
    x = np.linspace(xmin, xmax, Nx)                        
    dx = x[1] - x[0]  
    V = harmonic_potential(x, pert)                         
    H = sparse.eye(Nx, Nx, format='lil') * 2                
    for i in range(Nx - 1):                                 
        H[i, i + 1] = -1
        H[i + 1, i] = -1
    H = H / (dx ** 2)
    for i in range(Nx):                                       
        H[i, i] = H[i, i] + V[i]
    H = H.tocsc()                                           
    [evl, evt] = sla.eigs(H, k=neigs, which='SM')           
    for i in range(0, neigs):
        evt[:, i] = evt[:, i] / np.sqrt(                   
                                np.trapz(np.conj(
                                evt[:,i])*evt[:,i],x))
        evl = np.real(evl)                                 
    return evl, evt, x

evl_test , evt_test , x_test = schrodinger1D(-10, 10, 500, 5)

product_evt = []

for i in range(0, 500):
    product_evt.append( np.conjugate(evt_test[:, 0][i]) * evt_test[:, 0][i] )
    
sum_product = sum(product_evt)

问题原因 #

1. 验证方法与归一化逻辑不一致
   你用sum(product_evt)直接求和验证归一化，但归一化时用的是np.trapz（梯形积分），后者会考虑网格间距dx的权重。直接求和相当于把每个点的权重当成1，和积分的计算方式完全不同，这才是导致你误以为首个特征向量未归一化的核心原因。

2. 特征向量未排序
   sla.eigs返回的特征值和对应特征向量是无序的，你默认首个列向量是基态，但实际它可能对应任意一个特征值，不过这不是归一化失效的原因，只是会导致你验证的对象可能不是你预期的基态。

3. 数值虚部的干扰
   稀疏矩阵特征值求解会引入微小数值虚部，虽然你在归一化时用了np.conj，但验证时的求和没有处理虚部（不过影响极小）。

修正步骤 #

1. 用一致的方法验证归一化 #

把验证代码改成和归一化逻辑匹配的梯形积分：

# 正确的归一化验证方式
norm = np.trapz(np.conj(evt_test[:, 0]) * evt_test[:, 0], x_test)
print(norm)  # 结果会接近1，说明归一化有效

2. 对特征值和特征向量排序 #

为了确保特征向量按特征值从小到大排列（基态在前），在归一化前添加排序逻辑：

# 在schrodinger1D函数中，sla.eigs之后添加
idx = np.argsort(evl)
evl = evl[idx]
evt = evt[:, idx]
# 再进行后续归一化操作

3. 去除特征向量的数值虚部 #

由于薛定谔方程的束缚态解是实函数，可以直接取特征向量的实部，避免虚部干扰：

evt = np.real(evt)

修正后的完整函数 #

import numpy as np
import math
from scipy import sparse
from scipy.sparse import linalg as sla


def harmonic_potential(x, pert):
    if not pert:
       return x**2
    else:
        return x**2 + 4*np.exp(-10*x**2)  # 用numpy向量化操作替代循环，更高效

def schrodinger1D(xmin, xmax, Nx, neigs, pert=False):
    x = np.linspace(xmin, xmax, Nx)                        
    dx = x[1] - x[0]  
    V = harmonic_potential(x, pert)                         
    H = sparse.eye(Nx, Nx, format='lil') * 2                
    for i in range(Nx - 1):                                 
        H[i, i + 1] = -1
        H[i + 1, i] = -1
    H = H / (dx ** 2)
    for i in range(Nx):                                       
        H[i, i] = H[i, i] + V[i]
    H = H.tocsc()                                           
    evl, evt = sla.eigs(H, k=neigs, which='SM')           
    # 排序特征值和特征向量
    idx = np.argsort(evl)
    evl = np.real(evl[idx])
    evt = np.real(evt[:, idx])
    # 归一化
    for i in range(neigs):
        norm = np.sqrt(np.trapz(evt[:,i]**2, x))
        evt[:, i] = evt[:, i] / norm
    return evl, evt, x

# 测试
evl_test , evt_test , x_test = schrodinger1D(-10, 10, 500, 5)
# 验证归一化
print(np.trapz(evt_test[:,0]**2, x_test))  # 输出接近1

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