Speeding up code to refine grid points

I want to calculate following 2D integration as sums.

E and η are constants, is identity matrix of order 2, and is a 2 by 2 matrix whose elements are functions of . Tr means trace of matrix.

For small values of η, both real and imaginary parts of the integrand are not smooth (that's why integral2() fails), as shown below (real part). Fortunately, this non-smoothness appears only in a small portion of the space. To achieve a better approximation of the integration, I have decided to increase the number of grid points only around the non-smooth parts. While I feel like my code accomplishes what I want, it is quite messy. I believe that the same task could be performed much more efficiently. Could anyone please help me optimize my code? (I have tried to add comments to each line, and I welcome any questions.)

Code:

clear; clc;

% parameters
E = 4.4;
Dz = 0.5;
eta = 1e-3;
refined_by = 20; %factor by which coarse dkx will be divided


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


% Coarse grid points
dkx = 0.01;
dky = dkx;

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

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


%--------------------------------------------------------------------------
% calculating sums with coarse grid
%--------------------------------------------------------------------------
coarse_values = zeros(NBZx,NBZy);
parfor ikx = 1:NBZx
    kx = kxs(ikx);
    for iky = 1:NBZy
        ky = kys(iky); %#ok
        coarse_values(ikx,iky) = fun(kx,ky,E,Dz,eta);
    end
end

I_coarse = sum(coarse_values(:))*dkx*dky; %answer with coarse grid

%--------------------------------------------------------------------------
% calculating a threshold
%--------------------------------------------------------------------------

Real_coarse_values = real(coarse_values); %taking only the real part

% Deciding a threshold number. A dense grid will be taken 
% around the grid points for which integrant's absolute values is above
% threshold:
threshold = ( max(Real_coarse_values(:)) - min(Real_coarse_values(:)) )/100; 

% the values that abe within range of threshold:
okay_coarse_values = coarse_values;
okay_coarse_values(abs(Real_coarse_values) > threshold) = 0;

% contribution of these values to the total integration:
I_okay_coarse = sum(okay_coarse_values(:))*dkx*dky;


%--------------------------------------------------------------------------
% finding location of grid points which needs refinement
%--------------------------------------------------------------------------
[xx,yy] = find(abs(Real_coarse_values) > threshold);
List = [xx,yy];
List = sortrows(List,1);
%this is list of all points (of coarse grid) that with be refined.

LList = length(List);


%--------------------------------------------------------------------------
% refined grid points and calculating sums with refined grid
%--------------------------------------------------------------------------

new_dkx = dkx/refined_by; 
new_dky = new_dkx;

new_kxs = xmin:new_dkx:xmax;
new_kys = ymin:new_dky:ymax;

new_NBZx = length(new_kxs);
new_NBZy = length(new_kys);


% I plan to got to each k-point in List one and one and then calculating
% the corresponding refined k-points as:

refine_values = zeros(new_NBZx,new_NBZy);
for iList = 1:LList
    fprintf('iList/LList: %d/%d.\n',iList,LList);

    xList = List(iList,1); %kx coarse value
    yList = List(iList,2); %ky coarse value

    % correspoding kx and ky points in new refined grid are:
    xs = refined_by*xList - 2*refined_by + 1 : refined_by*xList + 1;
    ys = refined_by*yList - 2*refined_by + 1 : refined_by*yList + 1;
    
    %discarding the points which lie outside the new_kxs and new_kys arrays
    xs(xs<1 | xs>new_NBZx ) = []; 
    ys(ys<1 | ys>new_NBZy ) = [];

    %evaluating fun at refined points:
    for ix = xs
        kx = new_kxs(ix);
        for iy = ys
            ky = new_kys(iy); 
            refine_values(ix,iy) = fun(kx,ky,E,Dz,eta);
        end
    end
end
I_refined = sum(refine_values(:))*new_dkx*new_dky;

total_refined = I_okay_coarse + I_refined;

fprintf('Coarse: %0.4f, Refined: total_refined: %0.4f \n',I_coarse,total_refined);










%--------------------------------------------------------------------------
% function
%--------------------------------------------------------------------------

function out = fun(kx,ky,E,Dz,eta)
J = 1;
S = 1;
h0 = eye(2)*(3 * J * S);
hx = [0, 1; 1, 0]*(-J * S * (cos(ky / 2 - (3^(1/2) * kx) / 2) + cos(ky / 2 + (3^(1/2) * kx) / 2) + cos(ky)));
hy = [0, -1i; 1i, 0]*(-J * S * (sin(ky / 2 - (3^(1/2) * kx) / 2) + sin(ky / 2 + (3^(1/2) * kx) / 2) - sin(ky)));
hz = [1, 0; 0, -1]*(-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 = h0 + hx + hy + hz;

out = trace( inv( (E+1i*eta)*eye(2) - H ) );
end

[Torsten]

For small values of η, both real and imaginary parts of the integrand are not smooth (that's why integral2() fails)

If det(E+1i*eta)*eye(2) - H) is different from 0 in the region of interest for x,y and z, your function "out" as a function of x,y and z looks analytic to me. What makes you think it is not "smooth" ? Just by visual inspection ?

Did you try dblquad or quad2d for integration ? But looking at your plot, I think every integration method will fail in this case.

And the integration rule is not completely correct as stated: Either you have to evaluate your function in the cell centers - then all function values can be multiplied by deltakx*deltaky. Or - if you evaluate the function on the cell boundaries as you do - the contribution of the borders of the rectangle must be halfed.

[Luqman Saleem]

@Torsten I need to calculate this integral at various E values. For small eta and some values of E, the denominator becomes very huge at some (kx,ky) values which makes the integrand not smooth.

I have tried integral2, quad2d, and dblquad. All of them fail.

Could you please elaborate your last point.. about cell centers and cell boundaries?

[Torsten]

I need to calculate this integral at various E values. For small eta and some values of E, the denominator becomes very huge at some (kx,ky) values which makes the integrand not smooth.

In a mathematical sense, your integrand is smooth as smooth can be if det((E+1i*eta)*eye(2) - H) <> 0 in your region of integration. Maybe your function is not easy to integrate, but it is infinitly often differentiable, even analytic.

Could you please elaborate your last point.. about cell centers and cell boundaries?

I already corrected your previous ChatGPT solution in this respect, but here again for illustration:

f = @(x) x.^2;
x = 0:0.1:1;
dx = x(2) - x(1);
x_half = (x(1:end-1)+x(2:end))/2;
integral_centers = sum(f(x_half))*dx   % Your method, but with f applied to cell centers

integral_centers = 0.3325

integral_boundaries = trapz(x,f(x))    % trapz method which evaluates f at cell boundaries, but only halfed for x(1) and x(end)

integral_boundaries = 0.3350

integral_wrong = sum(f(x))*dx          % Your wrong integration with full weight at x(1) and x(end)

integral_wrong = 0.3850

1/3*x(end)^3-1/3*x(1)^3

ans = 0.3333

[Luqman Saleem]

@Torsten thank you so much clarification. I have just one more question. In you example, if we define "x" as x = linspace(0,1,10), and dx = x(2)-x(1). Then we don't need to to calculate x_half. Am I correct?

[Torsten]

x_half are the cell centers of the x-grid. You need them to evaluate your function at these values.

Here is an explanation how to use the so called "midpoint rule" for 2d-integration:

https://www.kristakingmath.com/blog/midpoint-rule-for-double-integrals