MATLAB: Integral3 reached limit error

numerical integration

Hi;
I have a complex numerical integration and i get the following error : "Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e+138. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy."
My Matlab codes as
function trns2 
clc
global kpp;
global teta;
global eta;
%**************************************************************

%********************P A R A M E T E R S***********************
%**************************************************************
lmd=1.55e-6;            %wavelength
%L=50;                    %distance
k=2*pi/lmd;
%******************oceanic parameters**********************
eta_s=1e-3;             %kolmogorov microscale length
x_t=1e-6;               %rate of dissipation of mean-squared temperature
eps=1e-4;               %rate of dissipation of kinetic energy per unit mass of fluid
%w=-1;                  %ratio of temperature to salinity contributions to the refractive index spectrum
A_t=1.863e-2;
A_s=1.9e-4;
A_ts=9.41e-3;
%***************************************************************


%***************************************************************
%***************************************************************
p1x=0.02;
p1y=0;
p2x=0.02;
p2y=0;
w1=-1;
w2=-2;
w3=-5;
coeff1=0.388*pi*(1e-8)*k*k*(eps^(-1/3))*x_t/(w1*w1);
coeff2=0.388*pi*(1e-8)*k*k*(eps^(-1/3))*x_t/(w2*w2);
coeff3=0.388*pi*(1e-8)*k*k*(eps^(-1/3))*x_t/(w3*w3);
j=1;
for L=10:10:50
f1=@(eta,kpp,teta)((kpp.^(-8/3)).*(1+2.35*((kpp*eta_s).^(2/3))).*...
    (w1*w1*exp(-A_t*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))+exp(-A_s*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))-2*w1*exp(-A_ts*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))).*...
    ((exp(1i*k*(p1x*p1x+p1y*p1y-p2x*p2x-p2y*p2y)./(L-eta)).*exp(1i*kpp.*cos(teta)*(p1x-p2x)).*exp(1i*kpp.*sin(teta)*(p1y-p2y)))-(exp(1i*k*(p1x*p1x+p2x*p2x+p1y*p1y+p2y*p2y)./(L-eta)).*exp(1i*kpp.*cos(teta)*(p1x-p2x)).*exp(1i*kpp.*sin(teta)*(p1y-p2y)).*exp(1i*(eta-L).*kpp.*kpp/k))));
Bx1=integral3(f1,0,L,0,inf,0,2*pi);
BxR1(j)=coeff1*(real(Bx1))
f2=@(eta,kpp,teta)((kpp.^(-8/3)).*(1+2.35*((kpp*eta_s).^(2/3))).*...
    (w2*w2*exp(-A_t*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))+exp(-A_s*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))-2*w2*exp(-A_ts*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))).*...
    ((exp(1i*k*(p1x*p1x+p1y*p1y-p2x*p2x-p2y*p2y)./(L-eta)).*exp(1i*kpp.*cos(teta)*(p1x-p2x)).*exp(1i*kpp.*sin(teta)*(p1y-p2y)))-(exp(1i*k*(p1x*p1x+p2x*p2x+p1y*p1y+p2y*p2y)./(L-eta)).*exp(1i*kpp.*cos(teta)*(p1x-p2x)).*exp(1i*kpp.*sin(teta)*(p1y-p2y)).*exp(1i*(eta-L).*kpp.*kpp/k))));
Bx2=integral3(f2,0,L,0,inf,0,2*pi);
BxR2(j)=coeff2*(real(Bx2))
f3=@(eta,kpp,teta)((kpp.^(-8/3)).*(1+2.35*((kpp*eta_s).^(2/3))).*...
    (w3*w3*exp(-A_t*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))+exp(-A_s*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))-2*w3*exp(-A_ts*(8.284*((kpp*eta_s).^(4/3))+12.978*((kpp*eta_s).^2)))).*...
    ((exp(1i*k*(p1x*p1x+p1y*p1y-p2x*p2x-p2y*p2y)./(L-eta)).*exp(1i*kpp.*cos(teta)*(p1x-p2x)).*exp(1i*kpp.*sin(teta)*(p1y-p2y)))-(exp(1i*k*(p1x*p1x+p2x*p2x+p1y*p1y+p2y*p2y)./(L-eta)).*exp(1i*kpp.*cos(teta)*(p1x-p2x)).*exp(1i*kpp.*sin(teta)*(p1y-p2y)).*exp(1i*(eta-L).*kpp.*kpp/k))));
Bx3=integral3(f3,0,L,0,inf,0,2*pi);
BxR3(j)=coeff3*(real(Bx3))
j=j+1;
end
plot(L,BxR1,L,BxR2,L,BxR3);
If you could help me to solve this problem, it will be a big pleasure for me.
Thanks a lot in advance.

Best Answer

In each of your expressions for f1, f2, f3, you have two subexpressions involving exp( 1i * T ./ (L-eta) ) where T is an expression.
In one of the two subexpressions in each case, T comes out as 0 because your p1x and p2x values are the same, and your p1y and p2y values are the same. With 0 in the numerator, the expression becomes exp(0) -- except theoretically at eta = L itself which would in theory give 0/0, but that is a removable discontinuity.
In the other of the two subexpressions in each case, the expression involving p1x, p2x, p1y, p2y gives a constant instead of 0. The numerator of the subexpression is thus non-zero except at specific values of kpp; that hints at the possibility of singularities as eta approaches L in 1/(L-eta) .
To see whether those infinity problems are real, we need to look at the rest of the subexpression within the exp() term. When we do so, we see that it is of the form (a*eta^2 + b).
Now, in the limit case, the a*eta^2 will overwhelm the b term, leading to a limit situation that is effectively a*eta^2 / eta, which would approach a*eta which would not be a problem with division by 0 singularities.
But is our eta in the range where the behavior approximates the limit? Alas, No, not unless kpp is fairly large, somewhere beyond 11000. Instead, for moderate values of kpp, the behavior of the exp() subexpression does indeed suffer from the difficulty of reaching infinity exponentially (1/x, x->0)
The subexpression is complex, and exp() of a large complex number is a complex (sin(),cos()) pair. As the complex number is reaching infinity exponentially, the (sin(), cos()) pair are going to be cycling faster and faster as eta approaches L.
... And, as you will recall, the integral() of sin() or cos() as their arguments go to infinity, is undefined.
So... your integral is undefined as eta approaches L, at least for moderate kpp. Possibly it becomes defined after a sufficiently large kpp.
(I have now confirmed that with a sufficiently large kpp, the integral over eta will go to zero within the tolerance of floating point precision. )
Related Question