Eigenvalues of a Matrix for a mesh of values and 3D plot.

Problem: Here is what i want to do. I need to construct an Hamiltonian matrix of dimension 2 by 2, then diagonalise it and find the eigenvalues for set of input values. I can do the problem for input values on a line, but i want a 3D plot and hence want to use a meshgrid.

The following working code plots for values on a line.

clear all
clc

a0 = 1.42; alpha0=1;
%-------------------------------------Define k range----------------------------------
kx = -2*pi/(3*a0):2*pi/(3*49.5*a0):2*pi/(3*a0);
ky = -2*pi/(3*a0):2*pi/(3*49.5*a0):2*pi/(3*a0);
a= 3*a0/2; b= sqrt(3)*a0/2;
N=2;
% % % % Hamiltonian generation
H=zeros(N,N);
for c = 1:100
    hk(c)=exp(1i*ky(c)*a0)+exp(1i*(a*kx(c)- b*ky(c)))+exp(-1i*(a*kx(c)+ b*ky(c)));%h(k) for kx, ky ne 0.
    
    %*****Hamiltoninan generation starts here**************
 for i=1:N
     for j=1:N
         if (j-i)==1
             H(i,j)= -alpha0*hk(c);
         elseif (j-i)==-1
             H(i,j)= -alpha0*conj(hk(c));
         end 
     end 
 end 
 E(:,c) = eig(H); %Eigenvalues
 
end

plot(kx,E);

Now for the above eigenvalues, i can't use surf to 3D plot and hence i tried the following code to generate eigenvalues in a meshgrid.

This code has an error.

a0 = 1.42; alpha0=1;
%-------------------------------------Define k range----------------------------------
kx = -2*pi/(3*a0):2*pi/(3*49.5*a0):2*pi/(3*a0);
ky = -2*pi/(3*a0):2*pi/(3*49.5*a0):2*pi/(3*a0);

[x,y]=meshgrid(kx,ky);
a= 3*a0/2; b= sqrt(3)*a0/2;
N=2;
H=zeros(N,N);
for l=1:100
    for m=1:100
   hk(l,m)=exp(1i*y(l,m)*a0)+exp(1i*(a*x(l,m)- b*y(l,m)))+exp(-1i*(a*x(l,m)+ b*y(l,m)));%h(k) for kx, ky ne 0.
    
   %*****Hamiltoninan**************
 for i=1:N
     for j=1:N
         if (j-i)==1
             H(i,j)= -alpha0*hk(l,m);
         elseif (j-i)==-1
             H(i,j)= -alpha0*conj(hk(l,m));
         end 
     end 
 end 
 E(l,m) = eig(H); %Here is the problem because it can't store two values in one grid
H=zeros(N,N);
    end 
 
end
 figure
 meshgrid(kx,ky,100);
 surf(kx,ky,E);

So the core of the problem is that, i want to solve the Hamiltonian matrix for the values scattered in a mesh and obtain the eigenvalues and plot those to obtain a 3D plot.

Any help or suggestion will be helpful.

Thank you.

Answer 1

[J. Alex Lee]

It seems you already understand the problem

E(l,m) = eig(H); %Here is the problem because it can't store two values in one grid

and have solved it in your 1D case

E(:,c) = eig(H); %Eigenvalues

If you want to extend your solution, just augment E by one more dimension

E(l,m,:) = eig(H);

As for plotting, I don't think surf/mesh/contour can handle a multi-dimensional Z variable, so you'll probably have to split

figure(1)
meshgrid(kx,ky);
surf(kx,ky,E(:,:,1));

So a full solution is (after some simplification of your code0

clear;
close all;
clc

a0 = 1.42;
alpha0 = 1;

% Define k range
k = -2*pi/(3*a0):2*pi/(3*49.5*a0):2*pi/(3*a0);
% it was redundant to do this twice

% don't hard-code the length of k
M = length(k);

[x,y] = meshgrid(k);

a = 3*a0/2;
b = sqrt(3)*a0/2;

% compute hk for all pairs of kx and ky
hk = exp(1i*y*a0) + exp(1i*(a*x- b*y)) + exp(-1i*(a*x+ b*y));

% by the way, the need for a0 is unclear, as it seems to divide out

% pre-allocate your results matrix
E = nan(M,M,2);

for i = 1:M
	for j = 1:M
		
		% your H computation appears needlessly complex
		H = -alpha0*[0,hk(i,j);conj(hk(i,j)),0];
		
		E(i,j,:) = eig(H);
		
	end
end

figure; hold on;
surf(x,y,E(:,:,1));
surf(x,y,E(:,:,2));

[J. Alex Lee]

And by the way I guess the solution is pretty trivial by paper and pencil...

The eigenvalues are

lam1 = alpha0*sqrt(hk.*conj(hk))
lam2 = -lam2

[Vira Roy]

Thank you so much ,J. Alex Lee . It worked like a charm. And thanks for all the inputs in the code. This was actually a prototype code for a different problem where i need that complex part of the code to form the matrix of higher dimensions.

And yeah you are right the example is trivial for this case. And once again Thank you so very much.

I managed to get the plot i wanted using surf and for anyone else who might benefit from this thread, here is the plotting command

figure; 
surf(x,y,E(:,:,1));
surf(x,y,E(:,:,2));

hold on
surf(x,y,-E(:,:,1));
surf(x,y,-E(:,:,2));

view(30,30)