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

Question

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

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

所以，我post把这两个代码，大家有错的可以找找看。我想这与精度有关，但我在网上找到的关于提高精度的内容提到了符号工具箱。还有什么我可以做的吗？

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

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

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

The same code works in Mathematica but not in Matlab

matlab

wolfram-mathematica