Solving a matrix equation with fixed point iteration method

Question

Luqman Saleem 2023-12-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2064851-solving-a-matrix-equation-with-fixed-point-iteration-method

评论： Luqman Saleem 2024-1-16

采纳的回答： Torsten

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2023-12-29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2064851-solving-a-matrix-equation-with-fixed-point-iteration-method#answer_1379636

编辑：Torsten 2023-12-29

在 MATLAB Online 中打开

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

15 个评论
显示 13更早的评论隐藏 13更早的评论

Torsten 2023-12-29

编辑：Torsten 2024-1-11

在 MATLAB Online 中打开

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

Torsten 2023-12-30

编辑：Torsten 2024-1-11

在 MATLAB Online 中打开

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 2024-1-9

@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 2024-1-11

编辑：Torsten 2024-1-11

在 MATLAB Online 中打开

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 2024-1-16

@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 2024-1-16

编辑：Torsten 2024-1-16

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 2024-1-16

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.

请先登录，再进行评论。

Solving a matrix equation with fixed point iteration method

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论
显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving a matrix equation with fixed point iteration method

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论 显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

15 个评论
显示 13更早的评论隐藏 13更早的评论