How to avoid singularity while performing quadruple (4D) integration
2 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
Walter Roberson
on 6 Dec 2015
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!