Тот же код работает в Mathematica, но не в Matlab

Question

У меня проблема с кодом, который я публикую ниже.Из этого кода я знаю, чего ожидать: мне нужно иметь только одно (положительное) собственное значение, равное единице, и значение, которое я называю qfi, должно быть равно 4, независимо от значения a (поскольку значение a должно увеличиваться, я долженувеличьте количество режимов, чтобы получить правильный результат, но это ожидаемо. И это не так, но, например, если у меня есть = 10 с 150 режимами, получим qfi 3.99986, что нормально).

Этот код отлично работает в Mathematica.Дело в том, что мне нужно написать что-то более сложное, чем это, и я не знаю, как правильно использовать mathematica, поэтому я хочу использовать matlab.Но в Matlab это не работает вообще, даже в тех случаях, когда все четко определено (потому что, если я увеличу количество мод в Matlab, я получу NaN в матрицах).Например, для a = 2 для mode = 100 (что является хорошим числом режимов) я получаю 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) )