numerical curl vs symbolic curl

I am tring to calculate curl of an expression both numerically and symbolically. When I compare them both, they are very different from each other. I get expected results using symbolic formalism but it gives me physically wrong results and a annoying kink at point. I have given both codes and results below.

Is there any way to get correct curl (like symbolic formalism)?

Symbolic curl:

% ============= Symbolic curl =============

clear; clc;

syms kx ky kz

% ============ system constants ============
aR = 2.95;
muBx = 0.1;
me = 9.1e-31/(1.6e-19*(1e10)^2);
m = 0.32*me;
hbar = 6.5911e-16;
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];


% ============ things required to build J ============
k = sqrt(kx^2+ky^2);

HR = hbar^2*k^2/(2*m)*eye(2) + aR * (kx*sy -ky*sx);
Hz = muBx*sx;
H = HR+Hz;

[Uk,~] = eig(H);

vx = 1/hbar*diff(H,kx);
vy = 1/hbar*diff(H,ky);

Sx = Uk'*hbar*sx*Uk;

X = 0.5*(Sx*vx + vx*Sx);
Y = 0.5*(Sx*vy + vy*Sx);

% ============ J = (Jx,Jy,0) ============
Jx = Uk(:,1)'*X*Uk(:,1); %x component of J
Jy = Uk(:,1)'*Y*Uk(:,1); %y component of J


vec = [kx,ky,kz]; %vector w.r.t. to I want to take curl
J = [Jx,Jy,0]; %J

sym_curl = curl(J,vec); %calculatig curl of J

curl_fun = matlabFunction( sym_curl(end) );

% ============ kx and ky points ============
NBZ = 300;
kx1 = -pi; kx2 = pi;
ky1 = kx1; ky2 = kx2;
kx_points = linspace(kx1,kx2,NBZ);
ky_points = linspace(ky1,ky2,NBZ);

% === calculating curl at each kx and ky ====

z_curl = zeros(NBZ,NBZ);
for ikx = 1:NBZ
    kx = kx_points(ikx);
    for iky = 1:NBZ
        ky = ky_points(iky);
        z_curl(ikx,iky) = curl_fun(kx,ky);
    end
end
[X,Y] = meshgrid(kx_points,ky_points);
%%
figure
surf(X,Y,real(z_curl),'linestyle','none');

xlabel('k_x');
ylabel('k_y');
zlabel('(\nabla \times J)_{symbolic}');
axis([kx1 kx2 ky1 ky2])
set(gca,'fontname','times','fontsize',14);

Numerical curl:

% ============= Numerical curl =============
clear; clc;

NBZ = 300;

% ============ constants ============
aR = 2.95;
muBx = 0.1;

me = 9.1e-31/(1.6e-19*(1e10)^2);
m = 0.32*me;
hbar = 6.5911e-16; %eV.s
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];

% ============ kx and ky points ============
kx1 = -pi; kx2 = pi;
ky1 = kx1; ky2 = kx2;
kxs = linspace(kx1,kx2,NBZ);
kys = linspace(ky1,ky2,NBZ);


Jx_comps = zeros(NBZ,NBZ);
Jy_comps = Jx_comps;
parfor ikx = 1:NBZ
    kx = kxs(ikx);
    for iky = 1:NBZ
        ky = kys(iky);
        
        % ============ things required to build J ============
        k = sqrt(kx^2+ky^2);
        
        HR = hbar^2*k^2/(2*m)*eye(2) + aR * (kx*sy -ky*sx);
        Hz = muBx*sx;
        H = HR+Hz;
        
        [Uk,~] = eig(H);
        
        vx = 1/hbar*[ (hbar^2*kx)/m,        -aR*1i; aR*1i, (hbar^2*kx)/m];
        vy = 1/hbar*[ (hbar^2*ky)/m,           -aR; -aR, (hbar^2*ky)/m];
        
        Sx = Uk'*hbar*sx*Uk;
        
        X = 0.5*(Sx*vx + vx*Sx);
        Y = 0.5*(Sx*vy + vy*Sx);
        
        % ============ J = (Jx,Jy,0) ============
        Jx = Uk(:,1)'*X*Uk(:,1); %x component of J
        Jy = Uk(:,1)'*Y*Uk(:,1); %y component of J
        
        
        Jx_comps(ikx,iky) = real(Jx);
        Jy_comps(ikx,iky) = real(Jy);
    end
end


[X,Y]=meshgrid(kxs,kys);
% === calculating numerical curl at each kx and ky ====
[z_curl,~] = curl(X,Y,Jx_comps,Jy_comps);

%%
figure

surf(X,Y,real(z_curl),'linestyle','none');

xlabel('k_x');
ylabel('k_y');
zlabel('(\nabla \times J)_{numerical}');
axis([kx1 kx2 ky1 ky2])
set(gca,'fontname','times','fontsize',14);

[David Goodmanson]

Hi Luqman,

what's eig_sort ?

[Luqman Saleem]

Hey, it basically calculates sorted eigenvalues and eigenvectors. I sorry I forget to change it to "eig()". I have updated it now. It should be working fine now.

[Paul]

I ran the code and did not get the same plot as shown above for the numeric case. Was the code edited after posting the figures?

Also, it looks like meshgrid() isn't being used correctly, which would affect the surf plots. When doing something like:

[X,Y]=meshgrid(x,y)
Z = f(X,Y)

Variation in x goes across the columns of Z and the variation in y goes down the rows of Z. The documentation of curl() says that the X,Y,U,V must all be in meshgrid form, so I think the call to curl should be

[z_curl,~] = curl(X,Y,Jx_comps.',Jy_comps.');

because Jx_comps and Jy_comps are defined with kx and ky varying down the rows and across the columns respectively.

In the symbolic code, it would be easier to just do:

[X,Y] = meshgrid(kx_points,ky_points);
z_curl = curl_fun(X,Y);

Answer 1

[Paul]

In addition to fixing the meshgrid per the comment above ...

The eigenvectors Uk are normalized differently between the symbolic and numeric solutions. I normalized the numeric Uk to be the same as in the symbolic case. The figures now are a prett close match, but I have no idea if this is the correct result you seek.

clear; clc;
syms kx ky kz
% ============ system constants ============
aR = 2.95;
muBx = 0.1;
me = 9.1e-31/(1.6e-19*(1e10)^2);
m = 0.32*me;
hbar = 6.5911e-16;
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];
% ============ things required to build J ============
k = sqrt(kx^2+ky^2);
HR = hbar^2*k^2/(2*m)*eye(2) + aR * (kx*sy -ky*sx);
Hz = muBx*sx;
H = HR+Hz;
[Uk,~] = eig(H);
vx = 1/hbar*diff(H,kx);
vy = 1/hbar*diff(H,ky);
Sx = Uk'*hbar*sx*Uk;
X = 0.5*(Sx*vx + vx*Sx);
Y = 0.5*(Sx*vy + vy*Sx);
% ============ J = (Jx,Jy,0) ============
Jx = Uk(:,1)'*X*Uk(:,1); %x component of J
Jy = Uk(:,1)'*Y*Uk(:,1); %y component of J
vec = [kx,ky,kz]; %vector w.r.t. to I want to take curl
J = [Jx,Jy,0]; %J
sym_curl = curl(J,vec); %calculatig curl of J
curl_fun = matlabFunction( sym_curl(end) );
% ============ kx and ky points ============
NBZ = 300;
kx1 = -pi; kx2 = pi;
ky1 = kx1; ky2 = kx2;
kx_points = linspace(kx1,kx2,NBZ);
ky_points = linspace(ky1,ky2,NBZ);
% % === calculating curl at each kx and ky ====
% z_curl = zeros(NBZ,NBZ);
% for ikx = 1:NBZ
%     kx = kx_points(ikx);
%     for iky = 1:NBZ
%         ky = ky_points(iky);
%         z_curl(ikx,iky) = curl_fun(kx,ky); % I think this indexing is backward
%     end
% end
[X,Y] = meshgrid(kx_points,ky_points);
z_curl = curl_fun(X,Y);
%%
figure
surf(X,Y,real(z_curl),'linestyle','none');
xlabel('k_x');
ylabel('k_y');
zlabel('(\nabla \times J)_{symbolic}');
axis([kx1 kx2 ky1 ky2])
view(3);
set(gca,'fontname','times','fontsize',14);

% numerical code
%%
% ============= Numerical curl =============
clear; clc;
NBZ = 300;
% ============ constants ============
aR = 2.95;
muBx = 0.1;
me = 9.1e-31/(1.6e-19*(1e10)^2);
m = 0.32*me;
hbar = 6.5911e-16; %eV.s
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];
% ============ kx and ky points ============
kx1 = -pi; kx2 = pi;
ky1 = kx1; ky2 = kx2;
kxs = linspace(kx1,kx2,NBZ);
kys = linspace(ky1,ky2,NBZ);
Jx_comps = zeros(NBZ,NBZ);
Jy_comps = Jx_comps;
for ikx = 1:NBZ
    kx = kxs(ikx);
    for iky = 1:NBZ
        ky = kys(iky);
        
        % ============ things required to build J ============
        k = sqrt(kx^2+ky^2);
        
        HR = hbar^2*k^2/(2*m)*eye(2) + aR * (kx*sy -ky*sx);
        Hz = muBx*sx;
        H = HR+Hz;
        
        [Uk,~] = eig(H);
        
        % normalize Uk to match the normalization of the symbolic eigenvectors
        Uk(:,1) = Uk(:,1)/Uk(2,1);
        Uk(:,2) = Uk(:,2)/Uk(2,2);
        
        vx = 1/hbar*[ (hbar^2*kx)/m,        -aR*1i; aR*1i, (hbar^2*kx)/m];
        vy = 1/hbar*[ (hbar^2*ky)/m,           -aR; -aR, (hbar^2*ky)/m];
        
        Sx = Uk'*hbar*sx*Uk;
        
        X = 0.5*(Sx*vx + vx*Sx);
        Y = 0.5*(Sx*vy + vy*Sx);
        
        % ============ J = (Jx,Jy,0) ============
        Jx = Uk(:,1)'*X*Uk(:,1); %x component of J
        Jy = Uk(:,1)'*Y*Uk(:,1); %y component of J
        
%         didn't take the real part in the symbolic approach        
%         Jx_comps(ikx,iky) = real(Jx); 
%         Jy_comps(ikx,iky) = real(Jy);
        Jx_comps(ikx,iky) = Jx;
        Jy_comps(ikx,iky) = Jy;
    end
end
[X,Y]=meshgrid(kxs,kys);
% === calculating numerical curl at each kx and ky ====
[z_curl,~] = curl(X,Y,Jx_comps.',Jy_comps.');  % have to transpose to match meshgrid definition
%%
figure
surf(X,Y,real(z_curl),'linestyle','none');  % real part only to match symbolic
xlabel('k_x');
ylabel('k_y');
zlabel('(\nabla \times J)_{numerical}');
axis([kx1 kx2 ky1 ky2])
view(3)
set(gca,'fontname','times','fontsize',14);