Self-consistent solution of integral equations using fsolve
10 views (last 30 days)
Show older comments
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
d = load('data.txt');
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([
y(1)-(tn./pi).*quadgk(@(z1) ((((1./(exp(z1./KbT)+1))).*(((1-x-(y(1)))./(z1-ed+(tn./pi).*log(abs((Dd-z1)./(Dd+z1)))-1i.*tn+1i.*tn.*(y(1))-(y(2))))))...
./(z-z1+1i.*delta.*heaviside(z)-1i.*delta.*heaviside(-z))), lower_lmt, upper_lmt, 'RelTol',0,'AbsTol',1e-9);
y(2)-(tn./pi).*quadgk(@(z1) (((1./(exp(z1./KbT)+1)).*(1+1i.*tn.*(((1-x-(y(1)))./(z1-ed+(tn./pi).*log(abs((Dd-z1)./(Dd+z1)))-1i.*tn+1i.*tn.*(y(1))-(y(2)))))))...
./(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?
Thanks in advance.
0 Comments
Answers (1)
Gokul Nath S J
on 21 Apr 2023
Edited: Gokul Nath S J
on 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.
For more information, kindly refer the following links.
Thanks,
Gokul Nath S J
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!