Filter löschen
Filter löschen

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
clear variables
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
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);
%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
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;
% 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]')
%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));
Is it effective way to approach the problem? Please suggest any improvement to solve above self-consistent problem effectively?
Thanks in advance.

Antworten (1)

Gokul Nath S J
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.
For more information, kindly refer the following links.
Gokul Nath S J
  1 Kommentar
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.

Community Treasure Hunt

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

Start Hunting!

Translated by