Solving self-consistent equations that involve matrices

17 次查看（过去 30 天）
Luqman Saleem
Luqman Saleem 2023-8-13
评论： Luqman Saleem 2023-8-22
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);
  8 个评论
显示 6更早的评论
Bruno Luong
Bruno Luong 2023-8-13
编辑：Bruno Luong 2023-8-13
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
Torsten 2023-8-13
编辑：Torsten 2023-8-13
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.

请先登录，再进行评论。

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

采纳的回答

Bruno Luong
Bruno Luong 2023-8-14
编辑：Bruno Luong 2023-8-14
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by