Integration of these kind of functions

I want to integrate

is a 2 by 2 matrix (given below), is identity matrix of order 2, and Eand η and positive real numbers. I want to calculate this integration (for each component of the matrix) with .

Note that the integrand does not diverges as long as η is non-zero, however when η is significantly small, the integrand can be very large for some points in plane. As an example, I have plot the imaginary part of (1,2) componet of the integrand matrix for and below. (the other components also follow the same behavior):

It can be seen that when η is small the integrand is fairly smooth and can be integrated easily even with integrate2 function. But when η is small, the integrand has very sharp peaks due to which integral2 fails. To calculate the correct integration with small η by trapz, we need a lot of grid points which makes everything very slow.

Is there a way to calculate this integration with ?

I have analytically calculated the condition for which the integrand shows sharp peaks, that condition is:

where are components of matrix H. We could also find the exact points on which the condition is satisfied.

where x are the roots of equation and y is all values between -1 and +1. And

So, for each y there will be 4 roots, x, which will give eight pionts (four and four . From these eight pionts, we will take only the real ; imaginary must be neglected. Also for each y there will be two points (one and one .

Is there any way to tell MATLAB to take a lot of gird pionts near points and then calculate the integration? Or is there any way to calculate this type of integration?

Code:

clear; clc;

% parameters
E = 4.4;
eta = 0.1;5e-4; 


% kx and ky limits and points
dkx = 0.006;
dky = dkx;

xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;

kxs = xmin:dkx:xmax;
kys = ymin:dky:ymax;

NBZx = length(kxs);
NBZy = length(kys);



% H matrix:
J = 1;
Dz = 0.5;
S = 1;

s0 = eye(2);
sx = [0, 1; 1, 0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];

h0 = 3 * J * S;
hx = @(kx, ky) -J * S * (cos(ky / 2 - (3^(1/2) * kx) / 2) + cos(ky / 2 + (3^(1/2) * kx) / 2) + cos(ky));
hy = @(kx, ky) -J * S * (sin(ky / 2 - (3^(1/2) * kx) / 2) + sin(ky / 2 + (3^(1/2) * kx) / 2) - sin(ky));
hz = @(kx, ky) -2 * Dz * S * (sin(3^(1/2) * kx) + sin((3 * ky) / 2 - (3^(1/2) * kx) / 2) - sin((3 * ky) / 2 + (3^(1/2) * kx) / 2));

H = @(kx, ky) s0 * h0 + sx * hx(kx, ky) + sy * hy(kx, ky) + sz * hz(kx, ky);


%integrand:
G00 = @(kx, ky) inv(E*eye(2) - H(kx,ky) + 1i*eye(2)*eta); %integrand

[David Goodmanson]

Hi Luqman,

since G00 appears to be a 2x2 matrix, what quantity are you plotting?

[Luqman Saleem]

@David Goodmanson oops, I forget to mention that. It is the imaginary part of (1,2) componet of G00.