How to calculate the numerical integration that contains singular points?

조회 수: 11 (최근 30일)
Henan Fang
Henan Fang 2019년 12월 19일
답변: Raynier Suresh 2020년 3월 24일
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에서 QSP, PKPD, and Systems Biology에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by