How to calculate the numerical integration that contains singular points?

4 次查看(过去 30 天)
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

回答(1 个)

Raynier Suresh
Raynier Suresh 2020-3-24
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,

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by