Solving a matrix equation with fixed point iteration method

I want to solve the following equation for

where

I believe Equ.1 can be solved using fixed iteration method. The twist here is that the term itself have a self-consistent equation:

So, to solve Equ.1, I have to solve Equ.3 for and then put the value of in Equ.2 and calculate and finally put in Equ.1 to solve it.

In Equ.1, And H and are 2-by-2 matrices given in below code. is identity matrix and E is a scalar parameter.

I'm working on this code. The loop for is converging well, but the loop for is very slow for certain u values, and it never converges for others. Any tips to speed up convergence or alternative solution methods?

clear; clc;

% parameters of equations:
E = 1;
n = 0.1;
u = 0.2;


% parameters of this script:
Nk = 1000; % number of points for integrating over kx and ky
max_iter = 2000; % # of maximum iterations
convergence_threshold = 1e-6;




% k-points and limits
xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;

kxs = linspace(xmin,xmax,Nk);
dkx = kxs(2) - kxs(1);
kys = linspace(ymin,ymax,Nk);
dky = kys(2) - kys(1);




%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_0 %%%%%%%%%%%%%%%%%%%%%%%%

Sigma_0 = (0.1 + 0.1i)*eye(2); % initial guess

for iter = 1:max_iter

    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );

    integral_term_Equ3 = zeros(2);
    parfor ikx = 1:Nk
        qx = kxs(ikx);
        for iky = 1:Nk
            qy = kys(iky);
            integral_term_Equ3 = integral_term_Equ3 + G_0(qx,qy) * dkx * dky;
        end
    end

    %new value of Sigma_0 (Equ.3):
    new_Sigma_0 =  n * u * inv( eye(2) - 1/(4*pi^2) * integral_term_Equ3 * u); 

    diff = norm(new_Sigma_0 - Sigma_0); %difference

    fprintf('G_0 Iteration: %d, Difference: %0.9f\n', iter, diff);

    if diff < convergence_threshold
        fprintf('G_0 converged after %d iterations\n', iter);
        break;
    end

    Sigma_0 = new_Sigma_0; % update Solution

end

G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 ); %the chosen G_0 


%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_x %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_x = Sigma_0;   %taking Sigma_0 as initial guess

for iter = 1:max_iter

    % Calculating the integration in Sigma (Equ.1) via sum:

    integral_term_Equ1 = zeros(2);
    parfor ikx = 1:Nk
        qx = kxs(ikx);
        for iky = 1:Nk
            qy = kys(iky);
            integrant = G_0(qx,qy) * (Sigma_x - 1i*E*v_x(qx,qy)) * G_0(qx,qy)' + 1i*E/2 * (G_0(qx,qy)*v_x(qx,qy)*G_0(qx,qy) + G_0(qx,qy)'*v_x(qx,qy)*G_0(qx,qy)');
            
            integral_term_Equ1 = integral_term_Equ1 + integrant * dkx * dky;
        end
    end

    %new value of Sigma_x (Equ.1):
    new_Sigma_x =  -1/(4*pi^2*n) * Sigma_0 * integral_term_Equ1 * Sigma_0'; 

    diff = norm(new_Sigma_x - Sigma_x); %difference

    fprintf('G Iteration: %d, Difference: %0.9f\n', iter, diff);

    if diff < convergence_threshold
        fprintf('G converged after %d iterations\n', iter);
        break;
    end

    Sigma_x = new_Sigma_x; % update Solution

end








%%%%%%%%%%%%%%%%%%%%%%%% H and v_x functions %%%%%%%%%%%%%%%%%%%%%%%%
function H = H(kx,ky)
J = 1;
D = 0.5;
S = 1;

a1 = [0,-1]';
a2 = [sqrt(3)/2,1/2]';
a3 = [-sqrt(3)/2,1/2]';

b1 = [sqrt(3)/2,-3/2]';
b2 = [sqrt(3)/2,3/2]';
b3 = [-sqrt(3),0]';

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


k = [kx,ky];
h0 = 3*J*S;
hx = -J*S*(  cos(k*a1) + cos(k*a2) + cos(k*a3)  );
hy = -J*S*(  sin(k*a1) + sin(k*a2) + sin(k*a3)  );
hz = 2*D*S*(  sin(k*b1) + sin(k*b2) + sin(k*b3)  );

H = s0*h0 + sx*hx + sy*hy + sz*hz;

end

function v_x = v_x(kx,ky)
J = 1;
D = 0.5;
S = 1;

sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];

dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);

v_x = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;

end

The result of above code are:

G_0 Iteration: 1, Difference: 0.127690191
G_0 Iteration: 2, Difference: 0.001871435
G_0 Iteration: 3, Difference: 0.000038569
G_0 Iteration: 4, Difference: 0.000000739
G_0 converged after 4 iterations
G Iteration: 1, Difference: 0.046150821
G Iteration: 2, Difference: 0.050105083
G Iteration: 3, Difference: 0.050493792
G Iteration: 4, Difference: 0.050542706
G Iteration: 5, Difference: 0.050547627
G Iteration: 6, Difference: 0.050543314
G Iteration: 7, Difference: 0.050536343
G Iteration: 8, Difference: 0.050528515
G Iteration: 9, Difference: 0.050520402
G Iteration: 10, Difference: 0.050512195
G Iteration: 11, Difference: 0.050503958
G Iteration: 12, Difference: 0.050495711
G Iteration: 13, Difference: 0.050487461
G Iteration: 14, Difference: 0.050479212
G Iteration: 15, Difference: 0.050470964
G Iteration: 16, Difference: 0.050462717
G Iteration: 17, Difference: 0.050454471
G Iteration: 18, Difference: 0.050446226
G Iteration: 19, Difference: 0.050437983
G Iteration: 20, Difference: 0.050429741
G Iteration: 21, Difference: 0.050421501
G Iteration: 22, Difference: 0.050413262
G Iteration: 23, Difference: 0.050405024
G Iteration: 24, Difference: 0.050396788
G Iteration: 25, Difference: 0.050388553
G Iteration: 26, Difference: 0.050380319
G Iteration: 27, Difference: 0.050372087
G Iteration: 28, Difference: 0.050363856
G Iteration: 29, Difference: 0.050355626
G Iteration: 30, Difference: 0.050347398
G Iteration: 31, Difference: 0.050339171
G Iteration: 32, Difference: 0.050330945
G Iteration: 33, Difference: 0.050322721
G Iteration: 34, Difference: 0.050314498
G Iteration: 35, Difference: 0.050306276
G Iteration: 36, Difference: 0.050298056
G Iteration: 37, Difference: 0.050289837
G Iteration: 38, Difference: 0.050281619
G Iteration: 39, Difference: 0.050273403
G Iteration: 40, Difference: 0.050265188
G Iteration: 41, Difference: 0.050256975
G Iteration: 42, Difference: 0.050248762
G Iteration: 43, Difference: 0.050240552
G Iteration: 44, Difference: 0.050232342
G Iteration: 45, Difference: 0.050224134
G Iteration: 46, Difference: 0.050215927
G Iteration: 47, Difference: 0.050207721
G Iteration: 48, Difference: 0.050199517
G Iteration: 49, Difference: 0.050191315
G Iteration: 50, Difference: 0.050183113
G Iteration: 51, Difference: 0.050174913
G Iteration: 52, Difference: 0.050166714
G Iteration: 53, Difference: 0.050158517
G Iteration: 54, Difference: 0.050150321
G Iteration: 55, Difference: 0.050142126
G Iteration: 56, Difference: 0.050133932
G Iteration: 57, Difference: 0.050125740
G Iteration: 58, Difference: 0.050117549
G Iteration: 59, Difference: 0.050109360
G Iteration: 60, Difference: 0.050101172
G Iteration: 61, Difference: 0.050092985
G Iteration: 62, Difference: 0.050084800

Answer 1

[Torsten]

main()

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Sigma_0 = 
  0.020586935783072 + 0.002472080743781i  0.000584326409329 + 0.001608817728934i
  0.000584326409329 + 0.001608817728934i  0.020586736770647 + 0.002650780039849i

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Sigma_x = 
 -0.000000000000000 + 0.024841208312721i  0.000769082661170 + 0.019717976813217i
 -0.000769082661170 + 0.019717976813218i  0.000000000000000 + 0.019770775411212i

function main
clear; clc;
format long
% parameters of equations:
E = 1;
n = 0.1;
u = 0.2;
% parameters of this script:
Nk = 300; % number of points for integrating over kx and ky
% k-points and limits
xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;
kxs = linspace(xmin,xmax,Nk);
dkx = kxs(2) - kxs(1);
kys = linspace(ymin,ymax,Nk);
dky = kys(2) - kys(1);
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_0 %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_0 = (0.1 + 0.1i)*eye(2); % initial guess
sigma0 = [Sigma_0(:,1);Sigma_0(:,2)];
sigma0 = fsolve(@fun_Sigma0,sigma0,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_0 = [sigma0(1:2),sigma0(3:4)]
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_x %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_x = Sigma_0;
sigmax = [Sigma_x(:,1);Sigma_x(:,2)];
sigmax = fsolve(@fun_Sigmax,sigmax,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_x = [sigmax(1:2),sigmax(3:4)]

function res = fun_Sigma0(sigma0)

    Sigma_0 = [sigma0(1:2),sigma0(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );
    integral_term_Equ3 = zeros(2);
    for ikx = 1:Nk
        qx = kxs(ikx);
        for iky = 1:Nk
            qy = kys(iky);
            integral_term_Equ3 = integral_term_Equ3 + G_0(qx,qy) * dkx * dky;
        end
    end
    Res = Sigma_0 - n * u * inv( eye(2) - 1/(4*pi^2) * integral_term_Equ3 * u);
    res = [Res(:,1);Res(:,2)];
    
end
function res = fun_Sigmax(sigmax)

    Sigma_x = [sigmax(1:2),sigmax(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );

    integral_term_Equ1 = zeros(2);
    for ikx = 1:Nk
        qx = kxs(ikx);
        for iky = 1:Nk
            qy = kys(iky);
            G_0_num = G_0(qx,qy);
            v_x_num = v_x(qx,qy);
            integrant = G_0_num * (Sigma_x - 1i*E*v_x_num) * G_0_num' + 1i*E/2 * (G_0_num*v_x_num*G_0_num + G_0_num'*v_x_num*G_0_num');
            integral_term_Equ1 = integral_term_Equ1 + integrant * dkx * dky;
        end
    end

    Res = Sigma_x - (-1/(4*pi^2*n) * Sigma_0 * integral_term_Equ1 * Sigma_0');
    res = [Res(:,1);Res(:,2)];
    
end
%%%%%%%%%%%%%%%%%%%%%%%% H and v_x functions %%%%%%%%%%%%%%%%%%%%%%%%
function H = H(kx,ky)
J = 1;
D = 0.5;
S = 1;
a1 = [0,-1]';
a2 = [sqrt(3)/2,1/2]';
a3 = [-sqrt(3)/2,1/2]';
b1 = [sqrt(3)/2,-3/2]';
b2 = [sqrt(3)/2,3/2]';
b3 = [-sqrt(3),0]';
s0 = eye(2,2);
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
k = [kx,ky];
h0 = 3*J*S;
hx = -J*S*(  cos(k*a1) + cos(k*a2) + cos(k*a3)  );
hy = -J*S*(  sin(k*a1) + sin(k*a2) + sin(k*a3)  );
hz = 2*D*S*(  sin(k*b1) + sin(k*b2) + sin(k*b3)  );
H = s0*h0 + sx*hx + sy*hy + sz*hz;
end
function v_x = v_x(kx,ky)
J = 1;
D = 0.5;
S = 1;
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);
v_x = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;
end
end

[Luqman Saleem]

@Torsten OMG...!! how did you do this. I spent a lot of time on stupid fixed piont iteration, and you solved it within a minute. wow. just wow. Thank you very much.

[Torsten]

I think there is a slight mistake in the approximation of the double integrals in your original code. This should work out correct:

main()

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Sigma_0 = 
  0.020579735521490 + 0.002458782459241i  0.000589986563542 + 0.001606680960751i
  0.000589986563542 + 0.001606680960751i  0.020584658399803 + 0.002637599312309i

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Sigma_x = 
  0.000000000000000 + 0.028225712037855i  0.000969520609637 + 0.018676280051416i
 -0.000969520609637 + 0.018676280051416i  0.000000000000000 + 0.020095681609324i

function main
clear; clc;
format long
% parameters of equations:
E = 1;
n = 0.1;
u = 0.2;
% parameters of this script:
Nk = 300; % number of points for integrating over kx and ky
% k-points and limits
xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;
kxs = linspace(xmin,xmax,Nk);
dkx = kxs(2) - kxs(1);
kys = linspace(ymin,ymax,Nk);
dky = kys(2) - kys(1);
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_0 %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_0 = (0.1 + 0.1i)*eye(2); % initial guess
sigma0 = [Sigma_0(:,1);Sigma_0(:,2)];
sigma0 = fsolve(@fun_Sigma0,sigma0,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_0 = [sigma0(1:2),sigma0(3:4)]
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_x %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_x = Sigma_0;
sigmax = [Sigma_x(:,1);Sigma_x(:,2)];
sigmax = fsolve(@fun_Sigmax,sigmax,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_x = [sigmax(1:2),sigmax(3:4)]

function res = fun_Sigma0(sigma0)

    Sigma_0 = [sigma0(1:2),sigma0(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );
    
    G_0_num = zeros(Nk,Nk,2,2);
    for ikx = 1:Nk
      qx = kxs(ikx);  
      for iky = 1:Nk
        qy = kys(iky);
        G_0_num(ikx,iky,:,:) = G_0(qx,qy);
      end
    end
    integral_term_Equ3 = zeros(2);
    %integral_term_Equ3(1,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,1)),1));
    %integral_term_Equ3(2,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,1)),1));
    %integral_term_Equ3(1,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,2)),1));
    %integral_term_Equ3(2,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,2)),1));
    integral_term_Equ3(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),2));
    integral_term_Equ3(2,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,1)),2));
    integral_term_Equ3(1,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,2)),2));
    integral_term_Equ3(2,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,2)),2));
    %integral_term_Equ3 = zeros(2);
    %for ikx = 1:Nk
    %    qx = kxs(ikx);
    %    for iky = 1:Nk
    %        qy = kys(iky);
    %        integral_term_Equ3 = integral_term_Equ3 + G_0(qx,qy) * dkx * dky;
    %    end
    %end
    Res = Sigma_0 - n * u * inv( eye(2) - 1/(4*pi^2) * integral_term_Equ3 * u);
    res = [Res(:,1);Res(:,2)];
    
end
function res = fun_Sigmax(sigmax)

    Sigma_x = [sigmax(1:2),sigmax(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );

    G_0_num = zeros(Nk,Nk,2,2);
    for ikx = 1:Nk
      qx = kxs(ikx);  
      for iky = 1:Nk
        qy = kys(iky);
        g_0_num = G_0(qx,qy);
        v_x_num = v_x(qx,qy);
        G_0_num(ikx,iky,:,:) = g_0_num * (Sigma_x - 1i*E*v_x_num) * g_0_num' + 1i*E/2 * (g_0_num*v_x_num*g_0_num + g_0_num'*v_x_num*g_0_num');
      end
    end
    integral_term_Equ1 = zeros(2);
    %integral_term_Equ1(1,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,1)),1));
    %integral_term_Equ1(2,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,1)),1));
    %integral_term_Equ1(1,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,2)),1));
    %integral_term_Equ1(2,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,2)),1));
    integral_term_Equ1(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),2));
    integral_term_Equ1(2,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,1)),2));
    integral_term_Equ1(1,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,2)),2));
    integral_term_Equ1(2,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,2)),2));
    
    %integral_term_Equ1 = zeros(2);
    %for ikx = 1:Nk
    %    qx = kxs(ikx);
    %    for iky = 1:Nk
    %        qy = kys(iky);
    %        G_0_num = G_0(qx,qy);
    %        v_x_num = v_x(qx,qy);
    %        integrant = G_0_num * (Sigma_x - 1i*E*v_x_num) * G_0_num' + 1i*E/2 * (G_0_num*v_x_num*G_0_num + G_0_num'*v_x_num*G_0_num');
    %        integral_term_Equ1 = integral_term_Equ1 + integrant * dkx * dky;
    %    end
    %end

    Res = Sigma_x - (-1/(4*pi^2*n) * Sigma_0 * integral_term_Equ1 * Sigma_0');
    res = [Res(:,1);Res(:,2)];
    
end
%%%%%%%%%%%%%%%%%%%%%%%% H and v_x functions %%%%%%%%%%%%%%%%%%%%%%%%
function H = H(kx,ky)
J = 1;
D = 0.5;
S = 1;
a1 = [0,-1]';
a2 = [sqrt(3)/2,1/2]';
a3 = [-sqrt(3)/2,1/2]';
b1 = [sqrt(3)/2,-3/2]';
b2 = [sqrt(3)/2,3/2]';
b3 = [-sqrt(3),0]';
s0 = eye(2,2);
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
k = [kx,ky];
h0 = 3*J*S;
hx = -J*S*(  cos(k*a1) + cos(k*a2) + cos(k*a3)  );
hy = -J*S*(  sin(k*a1) + sin(k*a2) + sin(k*a3)  );
hz = 2*D*S*(  sin(k*b1) + sin(k*b2) + sin(k*b3)  );
H = s0*h0 + sx*hx + sy*hy + sz*hz;
end
function v_x = v_x(kx,ky)
J = 1;
D = 0.5;
S = 1;
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);
v_x = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;
end
end

[Luqman Saleem]

Thanks for the update on integration. But Should not the direction of trapz() over "kxs" be "1", i.e., instead of the following

integral_term_Equ1(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),2));

we should have

integral_term_Equ1(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),1));

[Torsten]

I think you are right.

Becomes obvious if Nkx and Nky are introduced instead of Nk for both.

main()

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Sigma_0 = 
  0.020579735521490 + 0.002458782459241i  0.000589986563542 + 0.001606680960751i
  0.000589986563542 + 0.001606680960751i  0.020584658399803 + 0.002637599312309i

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Sigma_x = 
  0.000000000000000 + 0.028225712037855i  0.000969520609637 + 0.018676280051416i
 -0.000969520609637 + 0.018676280051416i  0.000000000000000 + 0.020095681609324i

function main
clear; clc;
format long
% parameters of equations:
E = 1;
n = 0.1;
u = 0.2;
% parameters of this script:
Nkx = 300; % number of points for integrating over kx 
Nky = 300; % number of points for integrating over ky
% k-points and limits
xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;
kxs = linspace(xmin,xmax,Nkx);
dkx = kxs(2) - kxs(1);
kys = linspace(ymin,ymax,Nky);
dky = kys(2) - kys(1);
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_0 %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_0 = (0.1 + 0.1i)*eye(2); % initial guess
sigma0 = [Sigma_0(:,1);Sigma_0(:,2)];
sigma0 = fsolve(@fun_Sigma0,sigma0,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_0 = [sigma0(1:2),sigma0(3:4)]
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_x %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_x = Sigma_0;
sigmax = [Sigma_x(:,1);Sigma_x(:,2)];
sigmax = fsolve(@fun_Sigmax,sigmax,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_x = [sigmax(1:2),sigmax(3:4)]

function res = fun_Sigma0(sigma0)

    Sigma_0 = [sigma0(1:2),sigma0(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );
    
    G_0_num = zeros(Nkx,Nky,2,2);
    for ikx = 1:Nkx
      qx = kxs(ikx);  
      for iky = 1:Nky
        qy = kys(iky);
        G_0_num(ikx,iky,:,:) = G_0(qx,qy);
      end
    end
    integral_term_Equ3 = zeros(2);
    %integral_term_Equ3(1,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,1)),2));
    %integral_term_Equ3(2,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,1)),2));
    %integral_term_Equ3(1,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,2)),2));
    %integral_term_Equ3(2,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,2)),2));
    integral_term_Equ3(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),1));
    integral_term_Equ3(2,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,1)),1));
    integral_term_Equ3(1,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,2)),1));
    integral_term_Equ3(2,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,2)),1));
    %integral_term_Equ3 = zeros(2);
    %for ikx = 1:Nkx
    %    qx = kxs(ikx);
    %    for iky = 1:Nky
    %        qy = kys(iky);
    %        integral_term_Equ3 = integral_term_Equ3 + G_0(qx,qy) * dkx * dky;
    %    end
    %end
    Res = Sigma_0 - n * u * inv( eye(2) - 1/(4*pi^2) * integral_term_Equ3 * u);
    res = [Res(:,1);Res(:,2)];
    
end
function res = fun_Sigmax(sigmax)

    Sigma_x = [sigmax(1:2),sigmax(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );

    G_0_num = zeros(Nkx,Nky,2,2);
    for ikx = 1:Nkx
      qx = kxs(ikx);  
      for iky = 1:Nky
        qy = kys(iky);
        g_0_num = G_0(qx,qy);
        v_x_num = v_x(qx,qy);
        G_0_num(ikx,iky,:,:) = g_0_num * (Sigma_x - 1i*E*v_x_num) * g_0_num' + 1i*E/2 * (g_0_num*v_x_num*g_0_num + g_0_num'*v_x_num*g_0_num');
      end
    end
    integral_term_Equ1 = zeros(2);
    %integral_term_Equ1(1,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,1)),2));
    %integral_term_Equ1(2,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,1)),2));
    %integral_term_Equ1(1,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,2)),2));
    %integral_term_Equ1(2,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,2)),2));
    integral_term_Equ1(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),1));
    integral_term_Equ1(2,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,1)),1));
    integral_term_Equ1(1,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,2)),1));
    integral_term_Equ1(2,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,2)),1));
    
    %integral_term_Equ1 = zeros(2);
    %for ikx = 1:Nkx
    %    qx = kxs(ikx);
    %    for iky = 1:Nky
    %        qy = kys(iky);
    %        G_0_num = G_0(qx,qy);
    %        v_x_num = v_x(qx,qy);
    %        integrant = G_0_num * (Sigma_x - 1i*E*v_x_num) * G_0_num' + 1i*E/2 * (G_0_num*v_x_num*G_0_num + G_0_num'*v_x_num*G_0_num');
    %        integral_term_Equ1 = integral_term_Equ1 + integrant * dkx * dky;
    %    end
    %end

    Res = Sigma_x - (-1/(4*pi^2*n) * Sigma_0 * integral_term_Equ1 * Sigma_0');
    res = [Res(:,1);Res(:,2)];
    
end
%%%%%%%%%%%%%%%%%%%%%%%% H and v_x functions %%%%%%%%%%%%%%%%%%%%%%%%
function H = H(kx,ky)
J = 1;
D = 0.5;
S = 1;
a1 = [0,-1]';
a2 = [sqrt(3)/2,1/2]';
a3 = [-sqrt(3)/2,1/2]';
b1 = [sqrt(3)/2,-3/2]';
b2 = [sqrt(3)/2,3/2]';
b3 = [-sqrt(3),0]';
s0 = eye(2,2);
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
k = [kx,ky];
h0 = 3*J*S;
hx = -J*S*(  cos(k*a1) + cos(k*a2) + cos(k*a3)  );
hy = -J*S*(  sin(k*a1) + sin(k*a2) + sin(k*a3)  );
hz = 2*D*S*(  sin(k*b1) + sin(k*b2) + sin(k*b3)  );
H = s0*h0 + sx*hx + sy*hy + sz*hz;
end
function v_x = v_x(kx,ky)
J = 1;
D = 0.5;
S = 1;
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);
v_x = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;
end
end

[Luqman Saleem]

@Torsten, Thank you so much! Your assistance is truly invaluable. I'm still grappling with a problem, though. I've just realized that my method for evaluating integrals over kx and ky is inaccurate. I need to compute Sigma_0 and Sigma_x for E values ranging from 0 to 6. Unfortunately, for some E values, Sigma_0 becomes very small (around 1e-6). This tiny Sigma_0 leads to extreme oscillations in G_0 across kx and ky, making the integral evaluation very challenging. Physically, I don't expect Sigma_x components to exceed 1. However, due to this integration issue, I sometimes obtain Sigma_x values on the order of 1e6. An example is when...

n = 0.2;

u = 0.2;

E = 4.4;

Nkx = 1814;

Nky = 2095;

I obtain Sigma_0

0.0386 - 0.0012i 0.0013 + 0.0005i

0.0013 + 0.0005i 0.0386 - 0.0012i

and Sigma_x

1.0e+03 *

-0.0000 - 4.4712i -0.0051 + 1.9231i

0.0051 + 1.9231i -0.0000 - 4.4651i

This 1e3 value doesn't make sense in my calculation. I suspect the integration is the problem because the Sigma_x changes drastically when I change the values of Nkx and Nky.

I'm curious if there's a reliable numerical method to accurately compute these integrals.

[Torsten]

The usual tool for 2d-integration is "integral2", but I think it will slow down your computations drastically compared to "trapz".

You will have to call "integral2" for each matrix entry separately. Thus four calls to "integral2" in fun_Sigma0 and four calls to "integral2" in fun_Sigmax. A "matrix integration" in one call is not possible.

[Luqman Saleem]

Unfortunetly integral2 fails to compute this integral. I was calling integral2 for each element of of the matrix. It says "Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test."

[Torsten]

You could try the method in my answer given here:

https://uk.mathworks.com/matlabcentral/answers/2059374-is-there-anyway-to-increase-the-calculation-speed-of-this-sqrt-integral2-without-loosing-accuracy?s_tid=srchtitle

[Luqman Saleem]

Ahan, I will try this tomorrow. It is already 3:24AM here.

[Torsten]

Does this help ? Or do you get stuck ?

main()

function main
clear; clc;
format long
% parameters of equations:
E = 1;
n = 0.1;
u = 0.2;
% parameters of this script:
%Nkx = 200; % number of points for integrating over kx 
%Nky = 300; % number of points for integrating over ky
% k-points and limits
xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;
%kxs = linspace(xmin,xmax,Nkx);
%dkx = kxs(2) - kxs(1);
%kys = linspace(ymin,ymax,Nky);
%dky = kys(2) - kys(1);
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_0 %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_0 = (0.1 + 0.1i)*eye(2); % initial guess
sigma0 = [Sigma_0(:,1);Sigma_0(:,2)];
%fun_Sigma0(sigma0)
sigma0 = fsolve(@fun_Sigma0,sigma0,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_0 = [sigma0(1:2),sigma0(3:4)]
%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Sigma_x %%%%%%%%%%%%%%%%%%%%%%%%
Sigma_x = Sigma_0;
sigmax = [Sigma_x(:,1);Sigma_x(:,2)];
%fun_Sigmax(sigmax)
sigmax = fsolve(@fun_Sigmax,sigmax,optimset('TolFun',1e-12,'TolX',1e-12));
Sigma_x = [sigmax(1:2),sigmax(3:4)]

function res = fun_Sigma0(sigma0)

    Sigma_0 = [sigma0(1:2),sigma0(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );
    
    [~,sol] = ode45(@(t,y)fun1_to_integrate(t,y,G_0),[xmin,xmax],zeros(4,1),odeset('RelTol',1e-10,'AbsTol',1e-10));
    integral_term_Equ3(1,1) = sol(end,1);
    integral_term_Equ3(2,1) = sol(end,2);
    integral_term_Equ3(1,2) = sol(end,3);
    integral_term_Equ3(2,2) = sol(end,4);
    
    %G_0_num = zeros(Nkx,Nky,2,2);
    %for ikx = 1:Nkx
    %  qx = kxs(ikx);  
    %  for iky = 1:Nky
    %    qy = kys(iky);
    %    G_0_num(ikx,iky,:,:) = G_0(qx,qy);
    %  end
    %end
    %integral_term_Equ3 = zeros(2);
    %%integral_term_Equ3(1,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,1)),2));
    %%integral_term_Equ3(2,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,1)),2));
    %%integral_term_Equ3(1,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,2)),2));
    %%integral_term_Equ3(2,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,2)),2));
    %integral_term_Equ3(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),1));
    %integral_term_Equ3(2,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,1)),1));
    %integral_term_Equ3(1,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,2)),1));
    %integral_term_Equ3(2,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,2)),1));
    %integral_term_Equ3 = zeros(2);
    %for ikx = 1:Nkx
    %    qx = kxs(ikx);
    %    for iky = 1:Nky
    %        qy = kys(iky);
    %        integral_term_Equ3 = integral_term_Equ3 + G_0(qx,qy) * dkx * dky;
    %    end
    %end
    Res = Sigma_0 - n * u * inv( eye(2) - 1/(4*pi^2) * integral_term_Equ3 * u);
    res = [Res(:,1);Res(:,2)];
    
end
function res = fun_Sigmax(sigmax)

    Sigma_x = [sigmax(1:2),sigmax(3:4)];
    % Calculation of integration in Sigma_0 (Equ.3) via sum:
    G_0 = @(kx,ky) inv( E*eye(2) - H(kx,ky) - Sigma_0 );

    [~,sol] = ode45(@(t,y)fun2_to_integrate(t,y,Sigma_x,G_0),[xmin,xmax],zeros(4,1),odeset('RelTol',1e-10,'AbsTol',1e-10));
    integral_term_Equ1(1,1) = sol(end,1);
    integral_term_Equ1(2,1) = sol(end,2);
    integral_term_Equ1(1,2) = sol(end,3);
    integral_term_Equ1(2,2) = sol(end,4);
    
    %G_0_num = zeros(Nkx,Nky,2,2);
    %for ikx = 1:Nkx
    %  qx = kxs(ikx);  
    %  for iky = 1:Nky
    %    qy = kys(iky);
    %    g_0_num = G_0(qx,qy);
    %    v_x_num = v_x(qx,qy);
    %    G_0_num(ikx,iky,:,:) = g_0_num * (Sigma_x - 1i*E*v_x_num) * g_0_num' + 1i*E/2 * (g_0_num*v_x_num*g_0_num + g_0_num'*v_x_num*g_0_num');
    %  end
    %end
    %integral_term_Equ1 = zeros(2);
    %%integral_term_Equ1(1,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,1)),2));
    %%integral_term_Equ1(2,1) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,1)),2));
    %%integral_term_Equ1(1,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,1,2)),2));
    %%integral_term_Equ1(2,2) = trapz(kxs,trapz(kys,squeeze(G_0_num(:,:,2,2)),2));
    %integral_term_Equ1(1,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,1)),1));
    %integral_term_Equ1(2,1) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,1)),1));
    %integral_term_Equ1(1,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,1,2)),1));
    %integral_term_Equ1(2,2) = trapz(kys,trapz(kxs,squeeze(G_0_num(:,:,2,2)),1));
    
    %%integral_term_Equ1 = zeros(2);
    %%for ikx = 1:Nkx
    %%    qx = kxs(ikx);
    %%    for iky = 1:Nky
    %%        qy = kys(iky);
    %%        G_0_num = G_0(qx,qy);
    %%        v_x_num = v_x(qx,qy);
    %%        integrant = G_0_num * (Sigma_x - 1i*E*v_x_num) * G_0_num' + 1i*E/2 * (G_0_num*v_x_num*G_0_num + G_0_num'*v_x_num*G_0_num');
    %%        integral_term_Equ1 = integral_term_Equ1 + integrant * dkx * dky;
    %%    end
    %%end

    Res = Sigma_x - (-1/(4*pi^2*n) * Sigma_0 * integral_term_Equ1 * Sigma_0');
    res = [Res(:,1);Res(:,2)];
    
end
%%%%%%%%%%%%%%%%%%%%%%%% H and v_x functions %%%%%%%%%%%%%%%%%%%%%%%%
function H = H(kx,ky)
J = 1;
D = 0.5;
S = 1;
a1 = [0,-1]';
a2 = [sqrt(3)/2,1/2]';
a3 = [-sqrt(3)/2,1/2]';
b1 = [sqrt(3)/2,-3/2]';
b2 = [sqrt(3)/2,3/2]';
b3 = [-sqrt(3),0]';
s0 = eye(2,2);
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
k = [kx,ky];
h0 = 3*J*S;
hx = -J*S*(  cos(k*a1) + cos(k*a2) + cos(k*a3)  );
hy = -J*S*(  sin(k*a1) + sin(k*a2) + sin(k*a3)  );
hz = 2*D*S*(  sin(k*b1) + sin(k*b2) + sin(k*b3)  );
H = s0*h0 + sx*hx + sy*hy + sz*hz;
end
function v_x = v_x(kx,ky)
J = 1;
D = 0.5;
S = 1;
sx = [0,1; 1,0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);
v_x = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;
end
function dy = fun1_to_integrate(x,~,G_0)
 fun = @(y) reshape(G_0(x,y),4,1);
 [~,dy] = ode45(@(y,~)fun(y),[ymin,ymax],zeros(4,1),odeset('RelTol',1e-10,'AbsTol',1e-10));
 dy = dy(end,:);
 dy = dy(:);
end
function dy = fun2_to_integrate(x,~,Sigma_x,G_0)
 fun = @(y) reshape(G_0(x,y)*(Sigma_x - 1i*E*v_x(x,y)) * G_0(x,y)' + 1i*E/2* (G_0(x,y)*v_x(x,y)*G_0(x,y) + G_0(x,y)'*v_x(x,y)*G_0(x,y)'),4,1);
 [~,dy] = ode45(@(y,~)fun(y),[ymin,ymax],zeros(4,1),odeset('RelTol',1e-10,'AbsTol',1e-10));
 dy = dy(end,:);
 dy = dy(:);
end
end

[Luqman Saleem]

I am running this code. It has been around 20 minutes. Let's see what is the outcome.

[Torsten]

I didn't test the code because it took too long. Maybe choosing bigger values for the error tolerances for ode45 is appropriate, but maybe if doing that you are back to "trapz".

I'd first test how long one call to "fun1_to_integrate" and "fun2_to_integrate" takes.

[Luqman Saleem]

@Torsten

Hey, the ode45 doesn't converge, and I've identified the main issue lies in the function. As becomes smaller, the function becomes increasingly challenging to integrate.

Currently, I'm focused on determining a suitable method to integrate functions like

with a small positive η. The larger η is, the easier it is to integrate over and . Fortunately, I can analytically locate points where the function becomes very large ( does not diverges on these points because we have η). Taking more points in the vicinity of points improves the integration answer. However, it's not feasible to use an infinite number of points in the entire space.

I believe that if I could write a code that takes a lot of points near these points, we can achieve accurate integration. Infact I run a test with grid of size and got quite accurate result. But it took a lot of time.

I've opened a new question to discuss it further (link: https://www.mathworks.com/matlabcentral/answers/2070611-integration-of-these-kind-of-functions)

[Torsten]

I believe that if I could write a code that takes a lot of points near these points, we can achieve accurate integration.

That's exactly what an adaptive ODE integrator like ode45 does. If it didn't succeed, I doubt you will find a way to handle this problem with existing MATLAB codes.

Are you sure that the matrix G00 has no singularities in the domain of integration ?

[Luqman Saleem]

I am sure that as long as η is non-zero there are no singularities, however, G00 can be very very very huge for some kx-ky points (I have mentioned those points in the new question that I've posted). I am aiming at case.