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

12 次查看（过去 30 天）
Basit
Basit 2015-12-5
评论： Basit 2015-12-6
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

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2015-12-6
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.
  1 个评论
Basit
Basit 2015-12-6
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by