Faster method for numerical integration in 3D

I have a very complex expression which includes three integrations:

Please bear with me. This looks very complicated but indeed it is very simple. Allow me to explain.

here ℏ is a constant, Tr means trace of product of matrices which depend upon . Also, and . And β is also a constant.

I want to evaluate these three integration numerically, but the method I am using is very slow and also it is not giving me expected results. I will be very thankful to you if you could have a look at my method and suggest any faster method for this 3D integration. My codes are given below:

main file:

clear; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a = 1;
JN = 1;
Dp = 0.75*JN;
T = 600;
eta = 0.01*JN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Limits
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xmin = -2*pi/(3*a); 
xmax = 4*pi/(3*a);
ymin = -2*pi/(a*sqrt(3)); 
ymax = 2*pi/(a*sqrt(3));
Emin = -10000; %The function goes to approximately zero beyond this range
Emax = 10000;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Numerical Integration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic

EBZsum = integral(@(E) integral(@(x) integral(@(y) (FUN_Exy(E,x,y,JN,Dp,T,eta)), ymin, ymax, 'ArrayValued', 1) , xmin, xmax, 'ArrayValued', 1), Emin, Emax, 'ArrayValued', 1);

tim = toc;

EBZsum = EBZsum*hbar; 
fprintf('time = %f mins, sigma = %f q2/hbar\n',tim/60,EBZsum)

The helping function:

% FUN_Exy.m
function out = FUN_Exy(E,x,y,JN,Dp,T,eta)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constants & Helping Expressions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a = 1;
s = 1;
t = pi;

kB = 0.0863;
hbar = 6.582e-13;
beta = 1/(kB*T);
ny = cos(t);


H = [                                                   4*JN*s, -(s*cos((a*x)/4 - (3^(1/2)*a*y)/4)*(4*JN - Dp*ny*4i))/2,                    -(s*cos((a*x)/2)*(4*JN + Dp*ny*4i))/2
    -(s*cos((a*x)/4 - (3^(1/2)*a*y)/4)*(4*JN + Dp*ny*4i))/2,                                                    4*JN*s, -(s*cos((a*x)/4 + (3^(1/2)*a*y)/4)*(4*JN - Dp*ny*4i))/2
    -(s*cos((a*x)/2)*(4*JN - Dp*ny*4i))/2, -(s*cos((a*x)/4 + (3^(1/2)*a*y)/4)*(4*JN + Dp*ny*4i))/2,                                                    4*JN*s];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Expressions of vx, vy, GA, GR, fE, dfdE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vx = 1/hbar * [                                                         0, (a*s*sin((a*x)/4 - (3^(1/2)*a*y)/4)*(4*JN - Dp*ny*4i))/8,                    (a*s*sin((a*x)/2)*(4*JN + Dp*ny*4i))/4
    (a*s*sin((a*x)/4 - (3^(1/2)*a*y)/4)*(4*JN + Dp*ny*4i))/8,                                                          0, (a*s*sin((a*x)/4 + (3^(1/2)*a*y)/4)*(4*JN - Dp*ny*4i))/8
    (a*s*sin((a*x)/2)*(4*JN - Dp*ny*4i))/4, (a*s*sin((a*x)/4 + (3^(1/2)*a*y)/4)*(4*JN + Dp*ny*4i))/8,                                                          0];

vy = 1/hbar * [                                                                  0, -(3^(1/2)*a*s*sin((a*x)/4 - (3^(1/2)*a*y)/4)*(4*JN - Dp*ny*4i))/8,                                                                  0
    -(3^(1/2)*a*s*sin((a*x)/4 - (3^(1/2)*a*y)/4)*(4*JN + Dp*ny*4i))/8,                                                                   0, (3^(1/2)*a*s*sin((a*x)/4 + (3^(1/2)*a*y)/4)*(4*JN - Dp*ny*4i))/8
    0,  (3^(1/2)*a*s*sin((a*x)/4 + (3^(1/2)*a*y)/4)*(4*JN + Dp*ny*4i))/8,                                                                  0];


GR = inv((E+1i*eta)*eye(size(H))-H);
GA = GR';
dGR = -GR*GR; %this is derivative of GR w.r.t. E
dGA = -GA*GA; %this is derivative of GA w.r.t. E

fE = (exp(beta*E) - 1)^(-1);
dfdE = -beta*exp(beta*E)*(exp(beta*E) - 1)^(-2);

TR1 = trace( vx*(GR*vy*GR + GA*vy*GA - 2*GR*vy*GA) )*dfdE; %first term in the given expression
TR2 = trace( vx*(GA*vy*dGA - dGA*vy*GA - GR*vy*dGR + dGR*vy*GR) )*fE; %second term in the given expression

out =  real(hbar/(2*(2*pi)^3) * (TR1-TR2) ); 
end

[John D'Errico]

Why not just use integral3?

[Luqman Saleem]

@John D'Errico integral3 was giving me an error because of matrices (I don't know how to work around that error).

[Torsten]

I don't know if it's faster, but you might want to try

EBZsum = integral3(@(E,x,y)arrayfun(@(E,x,y)FUN_Exy(E,x,y,JN,Dp,T,eta),E,x,y),Emin, Emax,xmin, xmax,ymin, ymax)

instead of

EBZsum = integral(@(E) integral(@(x) integral(@(y) (FUN_Exy(E,x,y,JN,Dp,T,eta)), ymin, ymax, 'ArrayValued', 1) , xmin, xmax, 'ArrayValued', 1), Emin, Emax, 'ArrayValued', 1);

[John D'Errico]

I suspect integral3 should be faster. But I've never truly tested it to verify that prediction. That it was giving you an error just means you either should tell us what the error was, or do as I suspect might work, to use the code Torsten suggests.

[Luqman Saleem]

@Torsten thank you, it worked. I have run the code with integral3 and it indeed is faster.

@John D'Errico yes, integral3 is faster.

[Torsten]

thank you, it worked. I have run the code with integral3 and it indeed is faster.

But the unsatisfactory results shouldn't have changed, have they ?