3D integral failure

2 views (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)

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by