Self-consistent solution of integral equations using fsolve
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)). Is it effective way to approach the problem? Please suggest any improvement to solve above self-consistent problem effectively? Thanks in advance. % 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
No answer yet