matlab not a number,quadgk函数计算带有bessel函数的数值积分时遇到了警告“遇到无限值或 Not-a-Number 值”的问题。...

function F=pwd(u,rWD,RD,Sf1,Sf2,C1D,C2D,wuf11,wuf21,lanmf11,lanmf21,n11D,n21D,n12D,n22D,M112,M212,tafa11,tafa21,Pi1D,Pi2D)

f1u=(1-wuf11)/(lanmf11*n11D+(1-wuf11)*u);

f2u=(1-wuf21)/(lanmf21*n21D+(1-wuf21)*u);

a11=rWD*sqrt(f1u)*besseli(1,rWD*sqrt(f1u));

a12=-rWD*sqrt(f1u)*besselk(1,rWD*sqrt(f1u));

a21=besseli(0,RD*sqrt(f1u));

a22=besselk(0,RD*sqrt(f1u));

a23=-besselk(0,RD*sqrt(1/n12D));

a31=M112*sqrt(f1u)*besseli(1,RD*sqrt(f1u));

a32=-M112*sqrt(f1u)*besselk(1,RD*sqrt(f1u));

a33=sqrt(1/n12D)*besselk(1,RD*sqrt(1/n12D));

b11=rWD*sqrt(f2u)*besseli(1,rWD*sqrt(f2u));

b12=-rWD*sqrt(f2u)*besselk(1,rWD*sqrt(f2u));

b21=besseli(0,RD*sqrt(f2u));

b22=besselk(0,RD*sqrt(f2u));

b23=-besselk(0,RD*sqrt(1/n22D));

b31=M212*sqrt(f2u)*besseli(1,RD*sqrt(f2u));

b32=-M212*sqrt(f2u)*besselk(1,RD*sqrt(f2u));

b33=sqrt(1/n22D)*besselk(1,RD*sqrt(1/n22D));

A11=-(a22*a33-a23*a32)/(tafa11*u*(a11*a22*a33-a11*a23*a32-a12*a21*a33+a12*a23*a31));

B11=(a21*a33-a23*a31)/(tafa11*u*(a11*a22*a33-a11*a23*a32-a12*a21*a33+a12*a23*a31));

A21=-(b22*b33-b23*b32)/(tafa21*u*(b11*b22*b33-b11*b23*b32-b12*b21*b33+b12*b23*b31));

B21=(b21*b33-b23*b31)/(tafa21*u*(b11*b22*b33-b11*b23*b32-b12*b21*b33+b12*b23*b31));

X1=@(aerfa)(A11*besseli(0,sqrt(f1u*(0.732-aerfa).^2))+B11*besselk(0,sqrt(f1u*(0.732-aerfa).^2)));

X2=@(aerfa)(A21*besseli(0,sqrt(f2u*(0.732-aerfa).^2))+B21*besselk(0,sqrt(f2u*(0.732-aerfa).^2)));

Y1=quadgk(X1,-1,1,'RelTol',1e-8,'AbsTol',1e-12);

Y2=quadgk(X2,-1,1,'RelTol',1e-8,'AbsTol',1e-12);

puw1D=Y1/2+Sf1/u+Pi1D/u;

puw2D=Y2/2+Sf2/u+Pi2D/u;

pwD=((1/(puw1D*u.^2-u*Pi1D)+C1D)*Pi1D+(1/(puw2D*u.^2-u*Pi2D)+C2D)*Pi2D+1/u)/(u*(C1D+C1D)+(1/(u*puw1D-Pi1D)+1/(u*puw2D-Pi2D)));

F=pwD;

end