Solving self-consistent equations that involve matrices

I want to calculate for a given value of and E. The self-consistent equation is

And

here is identity matrix of order 3, is 3-by-3 matrix, and are constants. The other variables are given in my code below. To solve this equation, I have wrote the following code (with help of ChatGPT) but this code does not coverge. Can someone please help me.

clear; clc;

% System Parameters
kx = pi/3;
ky = pi/5;
E = 8.5;
n_imp = 0.2;
u = 5;
JN = 1;
s = 1;
a = 1;
Dz = 0.75;


% Convergence Parameters
max_iterations = 1000;
convergence_threshold = 1e-6;


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


H = @(kx, ky) [4*JN*s, -(s*cos((a*kx)/4 + (3^(1/2)*a*ky)/4)*(Dz*4i + 4*JN))/2, (s*cos((a*kx)/4 - (3^(1/2)*a*ky)/4)*(Dz*4i - 4*JN))/2
    (s*cos((a*kx)/4 + (3^(1/2)*a*ky)/4)*(Dz*4i - 4*JN))/2, 4*JN*s, -(s*cos((a*kx)/2)*(Dz*4i + 4*JN))/2
    -(s*cos((a*kx)/4 - (3^(1/2)*a*ky)/4)*(Dz*4i + 4*JN))/2, (s*cos((a*kx)/2)*(Dz*4i - 4*JN))/2, 4*JN*s];


% Initial guess for G^R
G_R_fun = @(kx, ky) inv(E * eye(3) - H(kx, ky)); % Initial G^R function
G_R = inv(E * eye(3) - H(kx, ky)); %value at given kx, ky and E
 


figure();
xlim([1 max_iterations])
diffs = [];

for iter = 1:max_iterations

    % Calculate the denominator integration involving the GR
    integral_term = zeros(3);
    xs = xmin:0.01:xmax;
    Nx = length(xs);
    parfor ix = 1:Nx
        qx = xs(ix);
        for qy = ymin:0.01:ymax
            integral_term = integral_term + 1/(4*pi^2) * G_R_fun(qx, qy) * u * 0.01 * 0.01;
        end
    end

    %Evaluate Sigma^R:
    SigmaR = n_imp * u * inv(eye(3) - integral_term); 

    % Calculate the updated GR using the self-consistent equation
    G0_R_new_fun = @(kx, ky) inv(E * eye(3) - H(kx, ky) - SigmaR); %function
    G0_R_new = G0_R_new_fun(kx, ky); %value at given kx, ky and E
    
    % Check for convergence
    diff = norm(G0_R_new - G_R, 'fro');
    fprintf('Iteration: %d, Difference: %0.8f\n', iter, diff);
    
    if diff < convergence_threshold
        fprintf('Converged after %d iterations\n', iter);
        break;
    end
    
    % Update Green's function for the next iteration
    G_R = G0_R_new;
    G_R_fun = G0_R_new_fun;

    % Plotting differences:
    diffs(iter) = diff;
    plot(1:iter,diffs,'Marker','*','LineWidth',2,'Color','Black');
    drawnow;

end


fprintf('Final Green''s function:\n');
disp(G_R);

[Torsten]

You cannot "calculate G^R for a given value of k_y, k_y and E". As written, it requires G^R to be known for all xmin <= q_x <= x_max and ymin <= q_y <= ymax.

[Luqman Saleem]

@Torsten wait. wait wait. Are you are saying I can't solve this equation self-consistently? I have an initial guess of G^R for all xmin <= q_x <= x_max and ymin <= q_y <= ymax as given by

G_R_fun = @(kx, ky) inv(E * eye(3) - H(kx, ky)); % Initial G^R function

and then I update G^R for all xmin <= q_x <= x_max and ymin <= q_y <= ymax as given by

G_R = G0_R_new;

Is it not how self-consistent equations are solved?

[Torsten]

I wrote this because in your question, you write that you want to "calculate G^R for a given value of k_y, k_y and E".

"a" would mean: one.

I didn't know whether it's obvious for you that you have to compute G^R for all values xmin <= q_x <= x_max and ymin <= q_y <= ymax in order to be able to evaluate it in one specified point.

I don't know what you exactly mean with "self-consistent" in this context.

[Bruno Luong]

I'm confused, in your equation the denominator of is a matrix, and ChatGPT obviously understood as inverse of the matrix. Same for .

Honestly I never seen in the math literature the inverse is written as denominator of the expression.

I have a doubt how to correctly interpret your equation.

[Luqman Saleem]

@Bruno Luong It is the inverse. In Many-Body Physics's Green's function formalism, it is common practice to denote inverse as 1/something, that's why I wrote it like that.

[Luqman Saleem]

@Torsten you are right, it is for all k_x, k_y at a given E. And it is self-consistent in a sense that one needs to guess G^R first then put it in given equation and then check if the guess is the actual solution, if not, update the guess and keep going. is not it?

[Bruno Luong]

How physics peoples interpret ? Is it or ?

And why you write () at the end?

so the is not relevant in your equation. And why not pull out u outside the double integral, put it over outside ?

[Torsten]

And it is self-consistent in a sense that one needs to guess G^R first then put it in given equation and then check if the guess is the actual solution, if not, update the guess and keep going. is not it?

In mathematics, the method is called "fixed-point iteration" and it converges if you have a contraction mapping for G^R:

https://en.wikipedia.org/wiki/Contraction_mapping

As I said, I've never heard of self-consistent in this context.

The task reminds me of solving integral equations for all 9 components of your matrix G^R.

Answer 1

[Bruno Luong]

Use fixpoint to find good startting point for lsqnonlin: this code converges and find a solution:

% System Parameters
E = 8.5;
n_imp = 0.2;
u = 5;
JN = 1;
s = 1;
a = 1;
Dz = 0.75;

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

% Integration step and discretization vectors of the domain 
dx = 0.01;
qx = xmin:dx:xmax;
dy = 0.01;
qy =  ymin:dy:ymax;

% Number of iterations for fix points
fixpointniter = 10;

params = struct('E', E,...
    'n_imp', n_imp,...
    'u', u, ...
    'JN', JN, ...
    's', s, ...
    'a', a, ...
    'Dz', Dz, ...
    'dx', dx, ...
    'dy', dy, ...' ...
    'qx', qx, ...
    'qy', qy, ...
    'fixpointniter', fixpointniter);

[SigmaROut, Gxy] = EstimateSigmaR(params);

fix point iteration: 1/10
fix point iteration: 2/10
fix point iteration: 3/10
fix point iteration: 4/10
fix point iteration: 5/10
fix point iteration: 6/10
fix point iteration: 7/10
fix point iteration: 8/10
fix point iteration: 9/10
fix point iteration: 10/10

                                            Norm of      First-order 
 Iteration  Func-count      Resnorm            step       optimality
     0         10           153.904                               32
     1         20           73.8192         9.77849         1.17e+06      
     2         30           22.2866         3.02885             1.41      
     3         40           22.2866         11.2746             1.41      
     4         50           11.4495             2.5             1.72      
     5         60           11.4495         4.73589             1.72      
     6         70           5.31462         1.18397             1.14      
     7         80          0.973527         2.36795            0.744      
     8         90         0.0099212         0.60406           0.0641      
     9        100       1.81649e-05       0.0907756          0.00139      
    10        110       2.25454e-09      0.00990778         1.07e-05      
    11        120       7.90605e-17     0.000113883         5.79e-09      

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

relative error = 0.000000

SigmaROut

SigmaROut = 
   0.1045 - 1.4188i   0.9007 - 0.3963i  -0.2477 - 0.4365i
  -0.2478 - 0.4370i   0.0740 - 1.3894i   0.8528 - 0.4035i
   0.9003 - 0.3964i  -0.2184 - 0.4167i   0.0736 - 1.3892i

%%
function [SigmaROut, Gxy] = EstimateSigmaR(params)

% Use fixed point to find the first guess of SigmaR
SigmaR = zeros(3,3);
fixpointniter = params.fixpointniter;
for i = 1:fixpointniter
    fprintf('fix point iteration: %d/%d\n', i, fixpointniter)
    [~, SigmaR] = SigmaRError(SigmaR, params);
end

% Use lsqnonlin to refine the solutuon
x0 = SigmaR(:);
opt = struct('Display', 'iter', 'StepTolerance', 1e-10);
unconstrainedarg = cell(1,7);
x = lsqnonlin(@(x) SigmaRError(x, params), x0, unconstrainedarg{:}, opt);
SigmaROut = reshape(x,3,3);

[errorSigma, ~, Gxy] = SigmaRError(x, params);
relerr = norm(errorSigma) / norm(x);
fprintf("relative error = %f\n", relerr);

end

%%
function H = Hfun(params)
kx = params.qx;
ky = params.qy;
JN = params.JN;
s = params.s;
a = params.a;
Dz = params.Dz;

[kx, ky] = ndgrid(kx, ky);
kx = reshape(kx, [1 1 size(kx)]);
ky = reshape(ky, [1 1 size(ky)]);
z = zeros(size(kx));

% you should factorize this code, many terms computed over and over again,
% such as cos(a*kx)
H = [z+4*JN*s, -(s*cos((a*kx)/4 + (3^(1/2)*a*ky)/4)*(Dz*4i + 4*JN))/2, (s*cos((a*kx)/4 - (3^(1/2)*a*ky)/4)*(Dz*4i - 4*JN))/2;
    (s*cos((a*kx)/4 + (3^(1/2)*a*ky)/4)*(Dz*4i - 4*JN))/2, z+4*JN*s, -(s*cos((a*kx)/2)*(Dz*4i + 4*JN))/2;
    -(s*cos((a*kx)/4 - (3^(1/2)*a*ky)/4)*(Dz*4i + 4*JN))/2, (s*cos((a*kx)/2)*(Dz*4i - 4*JN))/2, z+4*JN*s];
end

%%
function [dx, SigmaR, Gxy] = SigmaRError(x, params)
dx = params.dx;
dy = params.dy;
E = params.E;
u = params.u;
n_imp = params.n_imp;
SigmaR = reshape(x, 3, 3);
% Calculate the updated GR using the self-consistent equation
Hxy = Hfun(params);
I3 = eye(3);
Gxy = pageinv(E*I3 - Hxy - SigmaR);
integral_term = (dx*dy*u)/(4*pi^2) * sum(Gxy, [3 4]);
SigmaR = (n_imp * u) * inv(I3 - integral_term); %#ok shutup stupid mlint on INV
dx = x(:)-SigmaR(:); % when it's consistent, dx is close to 0
end

[Luqman Saleem]

@Bruno Luong Thank you very much.