MATLAB: Problem with Error function “erfc(x)” in double integration

error function "erfc(x)" in double integration

I am looking to solve a double integration function which contain Error function (erfc(x)) in function.for clear presentation, i am giving below a code written by me.
syms r phi
er=3.5-i*0.2;
sig=0.5;
L=4.0;
k=1.0472;
theta=0.001:10:70.001;
g=length(theta);
Q=zeros(g,1);
for i=1:g;
     cs=cosd(theta(i));s=sind(theta(i));s2=s.*s;cs2=cs.*cs;ks=k.*sig;
     kL=k.*L;ks2=ks.*ks;kL2=kL.*kL;   
%(*Integration variables*)
     r2=r.*r; sf=sin(phi); csf=cos(phi); rx=r.*csf; ry=r.*sf; rx2=rx.*rx;
     ry2=ry.*ry;
%(*calculation of coefficients*)
     rt=sqrt(er-s2); rv=(er.*cs-rt)./(er.*cs+rt); rh=(cs-rt)./(cs+rt);
     rvh=(rv-rh)./2.0;
%(*calculation of field coefficients*)
rp=1.0+ rvh; rm=1.0-rvh; q=sqrt(1-r2); qt=sqrt(er-r2);
     a=rp./q; b=rm./q; c=rp./qt; d=rm./qt;      
%(* Cross pol-coeff. The formulation is in terms of B3 on p. 201*)
     B3=rx.*ry./cs; fvh1=(b-c).*(1-3.*rvh)-(b-c./er).*rp;
     fvh2=(a-d).*(1+3.*rvh)-(a-d.*er).*rm; Fvh=(abs((fvh1+fvh2).*B3)).^2;
      VH=2.*Fvh.*r ;
     rss =1.414*sig/L;
     au=q./r./1.414./rss;
     <fsh=(0.2821./au).*exp(-au.*au)-0.5.*erfc(au);>
     sha=1./(1+fsh);
P=VH.*sha;
X=matlabFunction(P);
integrand=@(phi,r)X(phi,r);
%  (Double integration)
Q(i) = quad2d(integrand,0.01,0.99,0.0,3.14);
end
ct=cotd(theta);
farg=ct./1.414./rss;
farg2=farg.*farg;
sfct=0.5*(exp(-farg2)./1.772./farg-erfc(farg));
Shdw=1./(1+sfct);
sigvh=10*log10(Shdw.*Q);
plot(theta,sigvh)
grid on
I am getting this error:
Error:
??? Error using ==> mupadmex
Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.
If the input expression contains a symbolic variable, use the VPA function instead.
Error in ==> sym.sym>sym.double at 927
Xstr = mupadmex('mllib::double', S.s, 0);
Error in ==> sym.sym>privformatscalar at 2539
    x = double(x);
Error in ==> sym.sym>privformat at 2524
s = privformatscalar(x);
Error in ==> sym.sym>sym.subsref at 1364
[inds{k},refs{k}] = privformat(inds{k});
Error in ==> Untitled_1 at 39
   fsh=(0.2821./au).*exp(-au.*au)-0.5.*erfc(au)

Best Answer

Which MATLAB version are you using?
Your computation of "au" involves the symbolic name "r", so erfc(au) is erfc() of a symbolic expression. The MuPAD symbolic engine is being invoked, but for some reason it is wanting to convert its input expression to double precision floating point, which it cannot do because of the symbolic "r".
The R2012b documentation clearly shows erfc() should operate on symbolic expressions without difficulty, so either there is a bug in your version or your version is earlier than full symbolic support was offered.
What happens if you try
sym r
erfc(r)
Related Question