Solving self-consistent equations that involve matrices
15 views (last 30 days)
Show older comments
Luqman Saleem
on 13 Aug 2023
Commented: Luqman Saleem
on 22 Aug 2023
I want to calculate
for a given value of
and E. The self-consistent equation is
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456227/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456232/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456237/image.png)
And
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456242/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456247/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456252/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456257/image.png)
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);
8 Comments
Bruno Luong
on 13 Aug 2023
Edited: Bruno Luong
on 13 Aug 2023
How physics peoples interpret
? Is it
or
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456407/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456412/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456417/image.png)
And why you write (
) at the end?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456377/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456382/image.png)
so the
is not relevant in your equation. And why not pull out u outside the double integral, put it over outside
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456387/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1456397/image.png)
Torsten
on 13 Aug 2023
Edited: Torsten
on 13 Aug 2023
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:
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.
Accepted Answer
Bruno Luong
on 14 Aug 2023
Edited: Bruno Luong
on 14 Aug 2023
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);
SigmaROut
%%
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
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!