How to calculate the numerical integration that contains singular points?

6 ビュー (過去 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 Center および File ExchangeNumerical Integration and Differentiation についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by