3D integral failure

1 view (last 30 days)
Luqman Saleem
Luqman Saleem on 20 Sep 2020
I am trying to sovle the following integration but integral3 keeps failing. Is there any work around for this. The expression for integration is:
with
where Tr[...] represent trace of ... matrices multiplication, b = 1/(kB . Q), kB, hbar, Q are constants, and
and
with h(x,y) is a 3-by-3 matrix given in below code F.m. I call integral3() as shown in runner.m code.
%runner.m
clear; clc;
% parameters
Q = 34.7;
eta = 0.025;
hbar = 6.582119569e-13;
% limits of int
xmin = -2*pi/3; xmax = 4*pi/3;
ymin = -2*pi/sqrt(3); ymax = 2*pi/sqrt(3);
zmin = 3; zmax = 15;
fun = @(x,y,z) hbar/((2*pi)^3*Q) * F(x,y,z,Q,eta);
I = integral3(@(x,y,z)arrayfun(fun,x,y,z),xmin,xmax,ymin,ymax,zmin,zmax);
It gives that the result is NaN.
%F.m
function F = F(x,y,z,Q,eta)
%=============== constants ===========
J = 1;
S = 1;
Delta = 2;
DbyJ = 2/sqrt(3);
hbar = 6.582119569e-13;
kB = 8.617333262145e-2;
b = 1/(kB*Q);
% ================ h(x,y),vx,vy ============
%some parameter for h(x,y)
t1 = sqrt(3)*x;
t2 = (3*y+t1)/2;
D = 1+1i*DbyJ;
fk1 = 1+exp(1i*t1);
fk2 = 1+exp(1i*t2);
fk1k2 = 1+exp(1i*(t1-t2));
%h(x,y)
h = J*S*[4*Delta, -D*conj(fk1), -conj(D*fk2);
-conj(D)*fk1,4*Delta, -D*fk1k2;
-D*fk2, -conj(D*fk1k2),4*Delta];
%vx, vy calculated using syms expression
vx = 1/hbar*[ 0, 3^(1/2)*J*S*exp(-3^(1/2)*conj(x)*1i)*(1 + DbyJ*1i)*1i, -(3^(1/2)*J*S*exp(- (conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*1i)/2
3^(1/2)*J*S*exp(3^(1/2)*x*1i)*(- 1 + conj(DbyJ)*1i)*1i, 0, -(3^(1/2)*J*S*exp(- (y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*1i)/2
-(3^(1/2)*J*S*exp((y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*1i)/2, -(3^(1/2)*J*S*exp((conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*1i)/2, 0];
vy = 1/hbar*[ 0, 0, -(J*S*exp(- (conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*3i)/2
0, 0, (J*S*exp(- (y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*3i)/2
-(J*S*exp((y*3i)/2 + (3^(1/2)*x*1i)/2)*(1 + DbyJ*1i)*3i)/2, (J*S*exp((conj(y)*3i)/2 - (3^(1/2)*conj(x)*1i)/2)*(- 1 + conj(DbyJ)*1i)*3i)/2, 0];
% ===================== GA, GR ================
GR = inv((z+1i*eta)*eye(3)-h);
GA = GR';
G = GA-GR;
term1 = b*z/(4*(sinh(b*z/2))^2);
F = -trace((z*vy)*G*vx*G)*term1;

Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by