How to calculate the numerical integration that contains singular points?

11 views (last 30 days)
T(x,y,afa) is a generated integrand, and the codes are as following.When I calculate M=arrayfun(@(D) integral2(@(x,y) T(x, y, D), 0,pi/2,-pi/6,pi/6,'reltol', 1e-6), afa) with varying afa=0:0.005:pi/6, the curve of output is not smooth and seems like noise. This is because the integrand has singular points. How to solve this problem? Many thanks!
function U=T(x,y,afa)
d1=1.34e-9;
d2=1.34e-9;
mu=5.5;
vh=1;
HBAR=1.05457266e-34;
ME=9.1093897e-31;
ELEC=1.60217733e-19;
Kh=2.95e10;
kc=sqrt(2.*ME.*ELEC./HBAR.^2);
k=kc.*sqrt(mu);
kh=sqrt(k.^2-(Kh-k.*sin(x).*cos(y)).^2-k.^2.*sin(x).^2.*sin(y).^2);
khg=sqrt(k.^2-(2.*Kh.*sin(afa./2).*sin(afa./2)-k.*sin(x).*cos(y)).^2-(2.*Kh.*sin(afa./2).*cos(afa./2)+k.*sin(x).*sin(y)).^2);
khpl=sqrt(k.^2-(Kh-k.*sin(x).*cos(y)).^2-k.^2.*sin(x).^2.*sin(y).^2+kc.^2.*vh);
khplpl=sqrt(k.^2-(Kh-k.*sin(x).*cos(y)).^2-k.^2.*sin(x).^2.*sin(y).^2+2.*kc.^2.*vh);
khgplpl=sqrt(k.^2-(2.*Kh.*sin(afa./2).*sin(afa./2)-k.*sin(x).*cos(y)).^2-(2.*Kh.*sin(afa./2).*cos(afa./2)+k.*sin(x).*sin(y)).^2+2.*kc.^2.*vh);
A2=exp(i.*khpl.*d1)./(exp(i.*(kh+khgplpl-khg).*d1)+exp(i.*khplpl.*d1));
U=abs(A2).^2;
end

Answers (1)

Raynier Suresh
Raynier Suresh on 24 Mar 2020
The quadgk function can handle singularity if the singularity is present at the boundary. In case if your singularity is not at the boundary you can split the integration domain to place the singularity at the boundary. Refer to the below links for more information,

Categories

Find more on QSP, PKPD, and Systems Biology in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!