相同的代码在 Mathematica 中有效，但在 Matlab 中无效

Question

我在下面发布的代码有问题。从这段代码中我知道会发生什么：我只需要一个（正）特征值等于 1 并且我命名为 qfi 的数量必须等于 4 无论 a 的值是多少（随着 a 的值增加，我必须增加模式的数量以获得正确的结果，但这是意料之中的。不是这样，但是，例如，如果我有 150 种模式的 a=10 得到 qfi 3.99986，这没关系）。

此代码在 Mathematica 中完美运行。问题是我需要写一些比这更复杂的东西而且我不知道如何正确使用mathematica，所以我想使用matlab。但是在 Matlab 中它根本不起作用，即使在一切都定义良好的情况下（因为如果我在 Matlab 中增加模式的数量，我会在矩阵中得到一个 NaN）。例如，对于模式=100 的 a=2（这是一个很好的模式数），我得到 qfi 80.0000。

所以，我发布了这两个代码，如果你们中的任何人能找到错误。我想这与精度有关，但我在网上找到的关于提高精度的内容提到了符号工具箱。还有什么我可以做的吗？

很抱歉发布了整个代码，但我已经为此苦苦挣扎了好几天。

Mathematica 代码：（不是很好复制）

modes = 100;

a = 2;

neg = 0;

(* the (l,m) element of the matrix *)

g[l_, m_] := (a^(l + m) E^-a^2)/Sqrt[Factorial[l]*Factorial[m]]

ρ = N[Table[g[l, m], {l, 0, modes}, {m, 0, modes}]];

tr = Tr[ρ]

ρnorm = ρ*(1/tr);

(*I find the eigenvalues and vectors. We want 1 eigenvalue equal to 1 and 
  all the others 0 *)

eige = Eigenvalues[ρnorm];

eige = Chop[eige];

vec = Eigenvectors[ρnorm];

vec = Chop[vec];

For[i = 0, i <= modes, i++, x = eige[[i]]; If[x < 0, neg = neg + 1]];

neg

0

Length[eige] - Count[eige, 0]

1

Max[eige]

1.

(* A second matrix like the first one *)

deriv[k_, n_] := (a^(-1 + k + n) E^-a^2 (-2 a^2 + k + n))/ Sqrt[k! n!]

derρ = N[Table[deriv[k, n], {k, 0, modes}, {n, 0, modes}]];

derρ = Chop[derρ];

L = 0;

For[i = 0, i <= modes, i++;
 For[j = 0, j <= modes , modes, j++;
   If[eige[[i]] + eige[[j]] != 0,
    (*I multiply eigenvectors with the derρ matrix *)
    rd = vec[[i]].derρ.tra[vec[[j]]];
    rd = Chop[rd];
    rd = rd[[1]];

    (*I multiply an eigenvector with the transpose of an other one and \
         I create the L matrix *)
    L = L + rd/(eige[[i]] + eige[[j]])*tra[vec[[i]]].{vec[[j]]}
    ]]]


    L = 2*L;
    L = Chop[L];

    qfi = Tr[ρnorm.L.L]

    4.

MatLab 代码（rhodiff=derρ）：

modes=100;
acc=1e-15;
a=2;
rho=zeros(modes,modes);
rhodiff=zeros(modes,modes);

for l=1:modes
    for m=1:modes
    rho(l,m)=(a^(l+m-2))*exp(-a^2)/(sqrt(factorial(l-1)*factorial(m-1)));
    end
end

tr=trace(rho);
rhonorm=rho/tr;

[V,D]=eig(rhonorm);
D(abs(D)<acc)=0;

for l=1:modes
    for m=1:modes
    rhodiff(l,m)=(a^(l+m-1))*(-2*a^2+l+m)*exp(-a^2)/(sqrt(factorial(l-1)*factorial(m-1)));
    end
 end

 L=zeros(modes,modes);

 for i=1:modes
    for j=1:modes
        if D(i,i)+D(j,j)~=0
        rd=((V(:,i)')*rhodiff*V(:,j));   
        rd(abs(rd)<acc)=0; 

        L=L+(((rd/(D(i,i)+D(j,j))))*V(:,i)*V(:,j)');
        L(abs(L)<acc)=0;
        end
    end
end

L=2*L;
qfi=trace(rhonorm*L*L)

Answer 1

这个问题几乎肯定是由于精度。在 Matlab 中，默认为双精度。在您的代码中，有一次，您正在计算factorial(modes-1)*factorial(modes-1)。这是一个非常大的数字。

在 Matlab 中，该数字的表示将被限制为双精度。在 Mathematica 中，由于它是一个整数，我的猜测是它可以精确地表示该数字。

最好的做法是让您的计算在数值上更加稳定。做到这一点的最好方法是将表达式转换为

(a^(l+m-2))*exp(-a^2)/(sqrt(factorial(l-1)*factorial(m-1)))

到产品中，这样您就不必在内存中实际拥有一些非常大的数字。例如，上面的表达式似乎可以等效地写为

exp(-a^2) * prod( a ./ sqrt(1:l-1) ) * prod( a ./ sqrt(1:m-1) )

这个版本从不要求你计算像factorial(100) 这样的东西，它不能用双精度精确表示。

我没查过，不过和你的其他表情很像

(a^(l+m-1))*(-2*a^2+l+m)*exp(-a^2)/(sqrt(factorial(l-1)*factorial(m-1)))

可以写成

(-2*a^2+l+m)*a*exp(-a^2) * prod( a ./ sqrt(1:l-1) ) * prod( a ./ sqrt(1:m-1) )