Self-consistent solution of integral equations using fsolve

Question

SACHIN VERMA 2023-4-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1944919-self-consistent-solution-of-integral-equations-using-fsolve

编辑： SACHIN VERMA 2023-4-23

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Gokul Nath S J 2023-4-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1944919-self-consistent-solution-of-integral-equations-using-fsolve#answer_1220148

编辑：Gokul Nath S J 2023-4-21

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

SACHIN VERMA 2023-4-22

编辑：SACHIN VERMA 2023-4-23

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.

请先登录，再进行评论。

Self-consistent solution of integral equations using fsolve

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Self-consistent solution of integral equations using fsolve

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论