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)).

% 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.

Answer 1

[Gokul Nath S J]

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.

https://www.mathworks.com/help/matlab/ref/ode45.html

https://www.mathworks.com/matlabcentral/answers/535829-solve-integral-using-ode45-by-derivative

Thanks,

Gokul Nath S J

[SACHIN VERMA]

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.