Mathematica数值积分结果为零？特征向量求解报错排查

在Mathematica中对依赖矩阵特征值与特征向量的函数执行数值积分时，调用NIntegrate得到结果为0，同时出现报错：Eigenvectors::eivec0 无法找到全部特征向量，且该报错后续输出将被抑制。相关代码如下：

\[Gamma]1 = 0.381;
u = 100 10^-3;
\[Phi] = \[Pi]/4;
v = 3 10^-3;
n = 1 10^-3;

H[k_, q_, \[Theta]_] := {
    {u/2, v (k Cos[\[Theta]] + q Cos[\[Phi]] - I (k Sin[\[Theta]] + q Sin[\[Phi]])), 0, 0},
    {v (k Cos[\[Theta]] + q Cos[\[Phi]] + I (k Sin[\[Theta]] + q Sin[\[Phi]])), u/2, \[Gamma]1, 0},
    {0, \[Gamma]1, -(u/2), v (k Cos[\[Theta]] + q Cos[\[Phi]] - I (k Sin[\[Theta]] + q Sin[\[Phi]]))},
    {0, 0, v (k Cos[\[Theta]] + q Cos[\[Phi]] + I (k Sin[\[Theta]] + q Sin[\[Phi]])), -(u/2)}
};

sortedEigenvaluesAndVectors[kk_, \[Theta]_, q_] := Module[{eigenvalues, eigenvectors, eigenpairs},
    {eigenvalues, eigenvectors} = Eigensystem[H[kk, q, \[Theta]]];
    eigenpairs = Transpose[{eigenvalues, eigenvectors}];
    eigenpairs = SortBy[eigenpairs, First];
    {eigenpairs[[All, 1]], eigenpairs[[All, 2]]}
];

Xintegrand[kk_, \[Theta]_, qq_, w_] := (
    4 kk)/(2 \[Pi])^2 (Abs[
        ConjugateTranspose[sortedEigenvaluesAndVectors[kk, \[Theta], 0][[2]][[1]]] .
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[2]][[3]]
    ]^2 / (w +
        sortedEigenvaluesAndVectors[kk, \[Theta], 0][[1]][[1]] -
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[1]][[3]] + I n) +
    Abs[
        ConjugateTranspose[sortedEigenvaluesAndVectors[kk, \[Theta], 0][[2]][[1]]] .
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[2]][[4]]
    ]^2 / (w +
        sortedEigenvaluesAndVectors[kk, \[Theta], 0][[1]][[1]] -
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[1]][[4]] + I n) +
    Abs[
        ConjugateTranspose[sortedEigenvaluesAndVectors[kk, \[Theta], 0][[2]][[2]]] .
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[2]][[3]]
    ]^2 / (w +
        sortedEigenvaluesAndVectors[kk, \[Theta], 0][[1]][[2]] -
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[1]][[3]] + I n) +
    Abs[
        ConjugateTranspose[sortedEigenvaluesAndVectors[kk, \[Theta], 0][[2]][[2]]] .
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[2]][[4]]
    ]^2 / (w +
        sortedEigenvaluesAndVectors[kk, \[Theta], 0][[1]][[2]] -
        sortedEigenvaluesAndVectors[kk, \[Theta], qq][[1]][[4]] + I n)
);

NIntegrate[Xintegrand[kk, \[Theta], 1, 0.1], {kk, 0, 200}, {\[Theta], 0, 2 \[Pi]}]

问题原因分析 #

1. 特征向量求解失效：矩阵H在部分参数点（如特定kk、\[Theta]值）会出现重特征值或数值奇异，导致Eigensystem无法生成完整特征向量组，后续内积计算出现数值异常，最终积分结果被错误归零。
2. 重复计算的数值波动：原代码多次重复调用sortedEigenvaluesAndVectors，不仅效率低下，还可能因参数变化时特征值排序的不连续性引入数值跳变，干扰积分收敛。
3. 复数积分的抵消问题：积分式中的分式为复数，NIntegrate对复数函数积分时，易因数值误差导致正负分量相互抵消，得到零结果。

解决方案 #

1. 优化特征求解稳定性 #

将Eigensystem替换为数值版本，并添加Method->"Arnoldi"参数，针对奇异矩阵优化特征求解，避免无法生成完整特征向量的问题。

2. 缓存特征计算结果 #

在积分函数内一次性计算两种情况（qq=0和qq=qq）的特征值与特征向量，减少重复计算带来的数值波动。

3. 明确积分的物理分量 #

对复数项取实部（或根据需求取虚部），避免复数积分的数值抵消。

4. 提升积分鲁棒性 #

给NIntegrate添加自适应蒙特卡洛方法及精度参数，增强对奇异点的处理能力。

修改后的代码 #

\[Gamma]1 = 0.381;
u = 100 10^-3;
\[Phi] = \[Pi]/4;
v = 3 10^-3;
n = 1 10^-3;

H[k_, q_, \[Theta]_] := {
    {u/2, v (k Cos[\[Theta]] + q Cos[\[Phi]] - I (k Sin[\[Theta]] + q Sin[\[Phi]])), 0, 0},
    {v (k Cos[\[Theta]] + q Cos[\[Phi]] + I (k Sin[\[Theta]] + q Sin[\[Phi]])), u/2, \[Gamma]1, 0},
    {0, \[Gamma]1, -(u/2), v (k Cos[\[Theta]] + q Cos[\[Phi]] - I (k Sin[\[Theta]] + q Sin[\[Phi]]))},
    {0, 0, v (k Cos[\[Theta]] + q Cos[\[Phi]] + I (k Sin[\[Theta]] + q Sin[\[Phi]])), -(u/2)}
};

sortedEigenvaluesAndVectors[kk_, \[Theta]_, q_] := Module[{eigenvalues, eigenvectors, eigenpairs},
    (* 数值特征求解+Arnoldi方法提升稳定性 *)
    {eigenvalues, eigenvectors} = Eigensystem[N[H[kk, q, \[Theta]]], Method -> "Arnoldi"];
    eigenpairs = Transpose[{eigenvalues, eigenvectors}];
    (* 按特征值实部排序，避免复数排序歧义 *)
    eigenpairs = SortBy[eigenpairs, Re[First]];
    {eigenpairs[[All, 1]], eigenpairs[[All, 2]]}
];

Xintegrand[kk_, \[Theta]_, qq_, w_] := Module[{baseEvals, baseEvecs, qEvals, qEvecs, term1, term2, term3, term4},
    (* 缓存特征计算结果 *)
    {baseEvals, baseEvecs} = sortedEigenvaluesAndVectors[kk, \[Theta], 0];
    {qEvals, qEvecs} = sortedEigenvaluesAndVectors[kk, \[Theta], qq];
    
    (* 计算各项 *)
    term1 = Abs[ConjugateTranspose[baseEvecs[[1]]] . qEvecs[[3]]]^2 / (w + baseEvals[[1]] - qEvals[[3]] + I n);
    term2 = Abs[ConjugateTranspose[baseEvecs[[1]]] . qEvecs[[4]]]^2 / (w + baseEvals[[1]] - qEvals[[4]] + I n);
    term3 = Abs[ConjugateTranspose[baseEvecs[[2]]] . qEvecs[[3]]]^2 / (w + baseEvals[[2]] - qEvals[[3]] + I n);
    term4 = Abs[ConjugateTranspose[baseEvecs[[2]]] . qEvecs[[4]]]^2 / (w + baseEvals[[2]] - qEvals[[4]] + I n);
    
    (* 取实部作为积分项 *)
    (4 kk)/(2 \[Pi])^2 * Re[term1 + term2 + term3 + term4]
];

(* 自适应蒙特卡洛积分提升鲁棒性 *)
NIntegrate[Xintegrand[kk, \[Theta], 1, 0.1], {kk, 0, 200}, {\[Theta], 0, 2 \[Pi]}, 
    Method -> "AdaptiveMonteCarlo", AccuracyGoal -> 4, PrecisionGoal -> 4]

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