How to avoid singularity while performing quadruple (4D) integration

Hi every one. I am solving a 4D integral where the integrand is in a numerator/denominator form. The denominator causes a singularity when it goes to zero with in the integration limits. The complete code is as follows

Vna1= [0.1086 - 0.3148i,   0.1244 - 0.3898i,   0.1842 - 0.5670i,   0.2167 - 0.6645i,   0.2500 - 0.8014i,...
   0.2903 - 0.9211i,   0.3189 - 1.0177i,   0.3481 - 1.1504i,   0.3811 - 1.2496i,   0.4046 - 1.3420i,...
  0.4298 - 1.4713i,   0.4566 - 1.5527i,   0.4751 - 1.6427i,   0.4960 - 1.7684i,   0.5169 - 1.8319i,...
   0.5305 - 1.9213i,   0.5469 - 2.0449i,   0.5618 - 2.0878i,   0.5703 - 2.1811i,   0.5820 - 2.3081i,...
 0.5910 - 2.3275i,   0.5944 - 2.4172i,   0.6011 - 2.6462i,   0.6041 - 2.3985i,   0.6025 - 3.2024i,...
 0.6041 - 2.3985i,   0.6011 - 2.6462i,   0.5944 - 2.4172i,   0.5910 - 2.3275i,   0.5820 - 2.3081i,...
 0.5703 - 2.1811i,   0.5618 - 2.0878i,   0.5469 - 2.0449i,   0.5305 - 1.9213i,   0.5169 - 1.8319i,...
 0.4960 - 1.7684i,   0.4751 - 1.6427i,   0.4566 - 1.5527i,   0.4298 - 1.4713i,   0.4046 - 1.3420i,...
 0.3811 - 1.2496i,   0.3481 - 1.1504i,   0.3189 - 1.0177i,   0.2903 - 0.9211i,   0.2500 - 0.8014i,...
 0.2167 - 0.6645i,   0.1842 - 0.5670i,   0.1244 - 0.3898i,   0.1086 - 0.3148i];

 freq=3e9;
 lambda=(3e8)/freq;
 L=0.48*lambda;
 b=0.2*lambda;
wc=b/2;yc=b/2;
Ls1=L;
Ls2=L;
Ws1=0.04*Ls1;
Ws2=0.04*Ls2;
xmin=-0.5*Ls2; xmax=0.5*Ls2;
ymin=0; ymax=Ws2;
zmin=-0.5*Ls1; zmax=0.5*Ls1;
wmin=0; wmax=Ws1;

 spacing=0.2;
 e2=1;

 k=2*pi/lambda;
 Vna1=(1e-1)*Vna1;
 N1=49;
 del=L/(N1+1);
 X_vm= -L/2+del:del:L/2-del;
 Poly1= polyfit(X_vm,Vna1,3);
 p3=Poly1(1);p2=Poly1(2);p1=Poly1(3);p0=Poly1(4);

 %---------To calculate V1 and V2 and multiplier----------
 V1=interp1(X_vm,Vna1,0);
 V2=interp1(X_vm,Vna1,0);
 mul= 1i*freq*e2/(2*V1*conj(V2));
 %---------------------------------------------------------

 A1= @(x,y,z,w) -k^2*x.^4+4*k^2*z.*x.^3;
 B1= @(x,y,z,w) (((x-z).^2+(y-w).^2).^(1/2)).*(2*1i*k*x.^2-4*1i*z.*x);
 C1= @(x,y,z,w) (-6*k^2*(z.^2)-(y-w).^2*k^2+2).*(x.^2);
 D1= @(x,y,z,w) (4*k^2*z.^3+(2*(y-w).^2*k^2-4).*z).*x;
 E1= @(x,y,z,w) -k^2*z.^4+(2-(y-w).^2*k^2).*(z.^2)-(y-w).^2;
 T1= @(x,y,z,w) exp(-1i*k*((x-z).^2+(y-w).^2).^(1/2));
 G1= @(x,y,z,w) ((x-z).^2+(y-w).^2).^(5/2);
 for h=1:length(spacing)
     d=0.1*spacing(h);
     F1= @(x,y,z,w) ((1./(pi*sqrt(Ws1^2+(w-wc).^2))).*(p3*(z.^3)+p2*(z.^2)+p1*(z.^1)+p0*(z.^0))).*((1./(pi*sqrt(Ws2^2+(y-yc+d).^2))).*(p3*(x.^3)+p2*(x.^2)+p1*(x.^1)+p0*(x.^0))).*(1+1/k^2).*(A1(x,y,z,w)+B1(x,y,z,w)+C1(x,y,z,w)+D1(x,y,z,w)+E1(x,y,z,w)).*T1(x,y,z,w)./G1(x,y,z,w);

    %Res= integralN(F1,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax,'AbsTol',1e-3,'RelTol',1e-3);
    Res= integral2(@(x,y)arrayfun(@(x,y)integral2(@(z,w)F1(x,y,z,w),zmin,zmax,wmin,wmax),x,y),xmin,xmax,ymin,ymax);
    Res= mul*Res;
    Imp(h)=1/Res;
    Res=0;
 end

The denominator function G1(x,y,z,w) has singularity when x=y=z=w=0, or if x=z & y=w at the same time. Can any one suggest me what to do in order to avoide these singularities. I have checked the results both with IntegralN and the integral2(integral2....) command as shown in code, but the problem is not solved.

Thanks in anticipation

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/matlab/ref/integral.html

Tips:

Do not use waypoints to specify singularities. Instead, split the interval and add the results of separate integrations with the singularities at the endpoints.

[Basit]

Dear @Walter Roberson Thanks for help. The problem is how to split the interval. As far as I can see, the integral function puts arrays of values each for x,y,z and w and then solves the function for these values. As the limits for x and z are same, the array values for both will be same. Same is true for y and w. So, for each value in the array, x=z and y=w. I have tried shuffling the integrals i.e. internal integral2 command works on x & z, while external one works on y & w. Then applied the interval splitting method given in mathworks help, http://cn.mathworks.com/help/matlab/math/singularity-on-interior-of-integration-domain.html

The new command is as follows

integral2(@(w,y)arrayfun(@(w,y)integral2(@(z,x)F1(x,y,z,w),zmin,zmax,@(z)-z,xmax),w,y),wmin,wmax,@(w)-w,ymax);

However, this time I got the error

 Warning: Reached the maximum number of function evaluations (10000). The result fails the global
error test.

and Matlab held busy for more than 10 minutes, without answer.