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

14 views (last 30 days)
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.

Accepted Answer

J. Alex Lee
J. Alex Lee on 22 Dec 2019
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));
  2 Comments
J. Alex Lee
J. Alex Lee on 22 Dec 2019
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
Vira Roy on 24 Dec 2019
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)

Sign in to comment.

More Answers (0)

Products


Release

R2016a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!