Faster method for numerical integration in 3D

9 次查看（过去 30 天）

显示更早的评论

Luqman Saleem 2023-3-24

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1934749-faster-method-for-numerical-integration-in-3d

评论： Torsten 2023-3-24

在 MATLAB Online 中打开

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