Speeding up code to refine grid points

Question

Luqman Saleem 2024-2-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2077901-speeding-up-code-to-refine-grid-points

编辑： Torsten 2024-2-4

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2024-2-4

编辑：Torsten 2024-2-4

在 MATLAB Online 中打开

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 2024-2-4

@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 2024-2-4

编辑：Torsten 2024-2-4

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

请先登录，再进行评论。

请先登录，再回答此问题。

Speeding up code to refine grid points

5 个评论
显示 3更早的评论隐藏 3更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Speeding up code to refine grid points

5 个评论 显示 3更早的评论隐藏 3更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论