I have a problem with the code I post below. From this code I know what to expect: I need to have only one (positive) eigenvalue equal to one and the quantity I name qfi must be equal to 4 no matter what the value of a is (as the value of a increases I must increase the number of modes to get the correct result, but this is expected. And it's not that but, for example if I have a=10 with 150 modes get qfi 3.99986 which is ok).
This code runs perfectly in Mathematica. The thing is that I need to write something more complicated than that and I don't know how to use mathematica properly, so I want to use matlab. But in Matlab it doesn't work at all, even in cases where everything is well defined (because if I increase the number of modes in Matlab I get a NaN in the matrices). For example, for a=2 for modes=100 (which is a good number of modes) i get qfi 80.0000.
So, I post the two codes, if any of you can find a mistake. I guess it has to do with precision, but what I found online about increasing precision mentioned the symbolic toolbox. Is there something else i could do?
I'm sorry for posting the whole code, but I'm struggling with it for days.
Mathematica code: (it's not very well copied)
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 code (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)
The issue is almost definitely due to precision. In Matlab, the default is double precision. In your code, at one point, you are computing factorial(modes-1)*factorial(modes-1). That is a very large number.
In Matlab, the representation of that number will be limited to double precision. In Mathematica, since it is an integer, my guess is that it can represent that number precisely.
The best course of action is to make your computation more numerically stable. The best way to do that is to turn expressions of the form
(a^(l+m-2))*exp(-a^2)/(sqrt(factorial(l-1)*factorial(m-1)))
into products so that you never have to actually have some really large number in memory. For instance, it seems that the above expression can be equivalently written as
exp(-a^2) * prod( a ./ sqrt(1:l-1) ) * prod( a ./ sqrt(1:m-1) )
this version never asks you to compute something like factorial(100), which cannot be represented precisely in double precision.
I haven't checked it, but it looks like your other expression
(a^(l+m-1))*(-2*a^2+l+m)*exp(-a^2)/(sqrt(factorial(l-1)*factorial(m-1)))
could be written as
(-2*a^2+l+m)*a*exp(-a^2) * prod( a ./ sqrt(1:l-1) ) * prod( a ./ sqrt(1:m-1) )