Complex eigenvalues for hermitian matrix

I have been trying to find the eigenenergies the Hamiltonian using the eig() function

% constants and parameters
hbar = 6.58211*10^-4;
w1 = 2000/hbar;
w2 = 2001/hbar;
syms wph;
g = 120;
pump = 10;
% the hamiltonian
h = [w1-wph pump 0; 
    pump w1-wph g; 
    0      g   w2];

This gives complex valued eigenenergies in terms of the parameter wph. This was expected because Matlab might be using an algorithm sue to which this happens. I also expected that when I plot these values against wph, they should be real but instead they come out to be complex valued with a fairly big imaginary part (code attached).

Can anyone explain why this is happening? The hamiltonian is hermitian and thus should have real eigenenergies which is not happening.

[Matt J]

What output do you see when you run the code?

[Neelesh Kumar Vij]

Here, I am attaching the plots:

Answer 1

[Matt J]

I suspect it is because you lack the float precision with which to crunch those huge integers in your expressions, like 17373258711169930298161307553886039650995152377.

Why exactly are you using symbolic (as oposed to numeric) eigenvalue analysis here?

[Neelesh Kumar Vij]

Hi Matt!

I want to plot the eigenvalues of the matrix as a function of wph. This should give me the behavior of anti-crossing, which is commonly observed in quantum mechanics. If what you are saying is true, what should be the best method of doing the above process?

[Matt J]

Why not as follows,

hbar = 6.58211*10^-4;
wph = linspace(0, 10/hbar,10000);

w1 = 2000/hbar;
w2 = 2001/hbar;

g = 120;
pump = 10;

clear eigen
for i=numel(wph):-1:1
    
    h= [w1-wph(i) pump 0; 
          pump w1-wph(i) g; 
          0      g   w2];
      
      eigen(:,i)=eig(h)*hbar;
end

eigen=sort(eigen);

or sort the eigenvalues at the end as you see fit?

[Neelesh Kumar Vij]

That does help!

Thank a lot Matt

[Matt J]

You're welcome, but please Accept-click the answer to signify that it solved your problem.