Problem with 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)

Answer 1

[Walter Roberson]

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)

[Dharmendra]

Thanks for your explanation.I am using MATLAB 2009b. the solutions u told, its not working. it shows error.

??? Undefined function or method 'erfc' for input arguments of type 'sym'.

Error in ==> Fung at 68 fsh=(0.2821./au).*exp(-au.*au)-0.5.*erfc(au);

[Walter Roberson]

Unfortunately the symbolic toolbox in R2009b did not implement erf() or erfc().

You could try substituting the integral shown at http://www.mathworks.com/help/symbolic/erfc.html in the "more about" section.

[Dharmendra]

i have tried as u told but it didn't work. is it not possible to solve erfc(x) by any other method?

[Walter Roberson]

What difficulty did you encounter when you substituted the integral?

[Dharmendra]

Mr. Walter, i thank you for being closer to solution.

When i am replacing erfc(au(r)) by integral expression erfc(r), it is implementing (it’s not proper to replace by erfc(r) instead of erfc(au). But when replacing with integral expression of erfc(au(r)), it is showing error.

??? Undefined function or method 'erf' for input arguments of type 'double' and attributes 'full 2d complex'.

Error in ==>

Error in ==> @(phi,r)X(phi,r)

Error in ==> quad2d>tensor at 355

Z = FUN(X,Y); NFE = NFE + 1;

Error in ==> quad2d at 169

[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in ==> Fung_1_test at 76

Q(i) = quad2d(integrand,0.001,1.0,0.0,pi);

[Walter Roberson]

Please show your code with the substituted version.

Note that MATLAB's numeric erf() function is only defined over real numbers, and not over complex numbers.

[Dharmendra]

This is the code:

syms r phi t

fr=5;

er=3.5-i*0.2;

sig=0.5;

L=4.0;

k=6.283.*fr/30;

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;

   rt=sqrt(er-s2); rv=(er.*cs-rt)./(er.*cs+rt); rh=(cs-rt)./(cs+rt);

   rvh=(rv-rh)./2.0;

   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;      

   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;

   acc=exp(-2.*ks2.*cs2)./(8.*pi);

    Wmh=0; Wnh=0.0;

for n=1:4;

    for m=1:4;

    wn=n.*kL2./(n.*n + kL2.*((rx - s).^2 + ry.^2))^1.5;     

    wm=m.*kL2./(m.*m + kL2.*((rx + s).^2 + ry.^2))^1.5;     

    vhmn=((ks2.*cs2).^(n+m))./(factorial(m).*factorial(n));

    Wm=vhmn.*wn.*wm;

    Wmh=Wmh+Wm;

    end

end

Wnh=Wnh+Wmh;

VH=2.*acc.*Fvh.*Wnh.*r;

% Error function

     rss =(2.*sig)./L;                                     

     au=((q./r)./1.414)./rss;

     c=exp(-t.^2);

     er=int(c,t,au,inf);

     fsh=(0.2821./au).*exp(-au.*au)-0.5*(2./sqrt(pi)).*(er);

     sha=1./(1+fsh);

P=VH.*sha;

X=matlabFunction(P);

integrand=@(phi,r)X(phi,r);

%% Double Integration

Q(i) = quad2d(integrand,0.001,1,0.0,pi);

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.*abs(Q)');

plot(theta,sigvh,'k')

This is giving error:

Error:

??? Undefined function or method 'erf' for input arguments of type 'double' and attributes 'full 2d complex'.

Error in ==>

Error in ==> @(phi,r)X(phi,r)

Error in ==> quad2d>tensor at 355

        Z = FUN(X,Y);  NFE = NFE + 1;

Error in ==> quad2d at 169

[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in ==> Fung_1_test at 76

Q(i) = quad2d(integrand,0.001,1.0,0.0,pi);

[Walter Roberson]

The implication is that the symbolic toolbox is generating the erf() in symbolic terms, and that matlabFunction() is preserving that erf() into a numeric erf() call, but that the argument that gets passed into that numeric erf() is complex; MATLAB's numeric erf() does not handle complex values.

I have not examined to see where the complex values are coming from in your expression. Are you expecting complex values?

[Dharmendra]

i am not sure about argument whether it is complex or real. is there any other way to solve this?

[Dharmendra]

Dear Walter, i have solved this error function problem by converting complex argument into real argument.Thanks a lot for your detailed explanation.