matlab 数值积分 奇点,一个关于电磁场的包含奇点、振荡特性的函数的数值积分...

h=0.0007;

c=3e8;%光速

er=2.33;%相对介电常数

e0=8.8542*10^(-12);%真空介电常数

f=2e9;%frequency

k0=2*pi*f/c;%自由空间波数

probex=x1;probey=y1;probez=z1;%监测点坐标

dipolex=x2;dipoley=y2;dipolez=z2;%等效偶极子坐标

for m=1:882;

for n=1:128;

%求柱坐标系

fi(m,n)=atan((probey(m)-dipoley(n))/(probex(m)-dipolex(n)));

r(m,n)=sqrt((probex(m)-dipolex(n))^2+(probey(m)-dipoley(n))^2);

z(m,n)=probez(m)-dipolez(n);

%转移系数矩阵

if mod(m,2)==1

if mod(n,2)==1

fun1=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*(besselj(0,x*r(m,n))-(besselj(1,x*r(m,n))./(x*r(m,n))));

integral11(m,n)=quadgk(fun1,0,1.0001*k0);

integral21(m,n)=integral11(m,n);

integral12(m,n)=quadgk(fun1,1.000145*k0,inf);

integral22(m,n)=integral12(m,n);

a(m,n)=(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral11(m,n)+integral12(m,n));

else

fun21=@(x)besselj(0,x*r(m,n)).*x.*sqrt(x.^2-k0^2).*exp(-z(m,n)*sqrt(x.^2-k0^2))./(sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2)));

integral11(m,n)=quadgk(fun21,0,1.0001*k0);

integral12(m,n)=quadgk(fun21,1.000145*k0,inf);

fun22=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*((sin(fi(m,n)))^2*besselj(0,x*r(m,n))+cos(2*fi(m,n))*(besselj(1,x*r(m,n))./(x*r(m,n))));

integral21(m,n)=quadgk(fun22,0,1.0001*k0);

integral22(m,n)=quadgk(fun22,1.000145*k0,inf);

a(m,n)=(integral11(m,n)+integral12(m,n))/(2*pi)-(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral21(m,n)+integral22(m,n));

end

else

if mod(n,2)==1

fun31=@(x)besselj(0,x*r(m,n)).*x.*sqrt(x.^2-k0^2).*exp(-z(m,n)*sqrt(x.^2-k0^2))./(sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2)));

integral11(m,n)=quadgk(fun31,0,1.0001*k0);

integral12(m,n)=quadgk(fun31,1.000145*k0,inf);

fun32=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*((cos(fi(m,n)))^2*besselj(0,x*r(m,n))-cos(2*fi(m,n))*(besselj(1,x*r(m,n))./(x*r(m,n))));

integral21(m,n)=quadgk(fun32,0,1.0001*k0);

integral22(m,n)=quadgk(fun32,1.000145*k0,inf);

a(m,n)=(integral11(m,n)+integral12(m,n))/(2*pi)-(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral21(m,n)+integral22(m,n));

else

fun4=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*(besselj(0,x*r(m,n))-(besselj(1,x*r(m,n))./(x*r(m,n))));

integral11(m,n)=quadgk(fun4,0,1.0001*k0);

integral21(m,n)=integral11(m,n);

integral12(m,n)=quadgk(fun4,1.000145*k0,inf);

integral22(m,n)=integral12(m,n);

a(m,n)=(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral11(m,n)+integral12(m,n));

end

end

end

end