# Self-consistent solution of integral equations using fsolve

15 Ansichten (letzte 30 Tage)
SACHIN VERMA am 11 Apr. 2023
Bearbeitet: SACHIN VERMA am 23 Apr. 2023
Hello all
I tried to solve the the self-consistent problem using numerical data integration. The matlab code (attached below) shows finite output which changes randomly as i increased number of data points for numerical integration and final results "G" diverges (or shows large error) for small "T" (T<10^(-2)).
% numerical self-consistent calculation
clc
clear variables
warning('off')
fileID = fopen('data.txt','w');
KbT = logspace(-4,2,21);
for i = 1:length(KbT)
Dd = 10.0;
tn = 0.2;
ed = -0.5;
delta = 10^(-6);
a = KbT(i)./tn;
syms g Gr11 Ga11 G11r G11a
x = (1:10000000);
x(1)=0.4; %initial guess
n = 1;
while true
m = 100; % number of data points for integration wrt to 'z' using trapezoidal rule
z=linspace(-10.0,10.0,m);
y = zeros(1,100);
for k=1:m
Gr = fun(z(k),Dd,ed,tn,KbT(i),delta,x(n));
%y = (-1./pi).*((1./2).*(1-tanh(z(k)./(2.*KbT(i))))).*imag(Gr11);
y = (-1./pi).*(1./(exp(z(k)./KbT(i))+1)).*imag(Gr);
Wanted_sol(k) = double(y);
end
%x(n+1) = integral(@(t) interp1(z,Wanted_sol,t,'linear','extrap'), z(1), z(end),'ArrayValued',true);
x(n+1) = quadgk(@(t) interp1(z,Wanted_sol,t,'linear','extrap'), z(1), z(end),'RelTol',0,'AbsTol',1e-9);
if (abs(x(n+1)-x(n))<0.001)
ndown = x(n+1);
nup = ndown;
m = 10; % number of data points for integration wrt to 'z' using simpson rule
omega=linspace(-5.*KbT(i),5.*KbT(i),m);
f1 = zeros(1,10);
f2 = zeros(1,10);
f3 = zeros(1,10);
for j = 1:length(omega)
Gr = fun(omega(j),Dd,ed,tn,KbT(i),delta,nup);
Ga = conj(Gr);
T1 = (tn.^2).*(Gr.*Ga);
fdd = (1./KbT(i)).*(exp(omega(j)./KbT(i))./(exp(omega(j)./KbT(i))+1).^2);
f1(j) = fdd.*T1;
end
% Use quadgk to integrate the data
L01 = 2.*quadgk(@(t) interp1(omega,double(f1),t,'linear','extrap'), omega(1), omega(end),'RelTol',0,'AbsTol',1e-9);
% Use Gaussian quadrature to integrate the data
%L01 = 2.*integral(@(t) interp1(omega,f1,t,'linear','extrap'), omega(1), omega(end),'ArrayValued',true);
G = L01;
fprintf(fileID,'%8.6e %8.6e %8.6e\n',[a,G,nup]');
fprintf('%8.6e %8.6e %8.6e\n',[a,G,nup]')
break
end
x(n)=x(n+1);
n=n+1;
end
end
fclose(fileID);
%setting the Matlab figure
a = d(:,1);
b = d(:,2);
c = d(:,3);
semilogx(a,b,'o-K','Linewidth', 2.0,'Markersize',4.0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Gr = fun(z,Dd,ed,tn,KbT,delta,x)
options = optimoptions('fsolve','Display','off','TolFun',1e-9,'TolX',1e-9);
% Integration limit
lower_lmt = -10.0;
upper_lmt = 10.0;
y_0 = [0.1; 0.1];
% self-consistent equations
F = @(y) double([
./(z-z1+1i.*delta.*heaviside(z)-1i.*delta.*heaviside(-z))), lower_lmt, upper_lmt, 'RelTol',0,'AbsTol',1e-9);
./(z-z1+1i.*delta.*heaviside(z)-1i.*delta.*heaviside(-z))), lower_lmt, upper_lmt, 'RelTol',0,'AbsTol',1e-9)...
]);
sol = fsolve(F,y_0,options);
eng_1 = vpa(sol(1));
eng_2 = vpa(sol(2));
g0 = z-ed+(tn./pi).*log(abs((Dd-z)./(Dd+z)))+1i.*tn;
Gr = ((1-x-eng_1)./(g0-eng_2-1i.*tn.*eng_1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Is it effective way to approach the problem? Please suggest any improvement to solve above self-consistent problem effectively?
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (1)

Gokul Nath S J am 21 Apr. 2023
Bearbeitet: Gokul Nath S J am 21 Apr. 2023
Hi Sachin,
It seems that you are trying to solve an self consistent solution of an integral equation. As an alternative approach, I can recommend you to write the equation in terms of an differential equation and apply ODE45 which can replace many redundant section.
Thanks,
Gokul Nath S J
##### 1 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden
SACHIN VERMA am 22 Apr. 2023
Bearbeitet: SACHIN VERMA am 23 Apr. 2023
Thanks sir,
I will also try to implement the alternative approach (i.e. by using "ODE45") to solve the problem. I think the previous method (using fsolve and numerical integration) is also a correct approach if one choose the integration limit, and other parameters correctly.
Thank you for the suggestion.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Numerical Integration and Differentiation finden Sie in Help Center und File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by