numerical curl vs symbolic curl

1 次查看（过去 30 天）
Luqman Saleem
Luqman Saleem 2021-5-25
回答： Paul 2021-5-28
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);
  3 个评论
显示 1更早的评论
Luqman Saleem
Luqman Saleem 2021-5-26
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
Paul 2021-5-28
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);

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Paul
Paul 2021-5-28
Ran in:
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);

类别

在 Help Center 和 File Exchange 中查找有关 Conversion 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by