integral3 reached limit error
显示 更早的评论
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.
2 个评论
Walter Roberson
2013-8-26
It seems to me that if the bound on error is 2.1e+138 then you probably have something that is approaching infinity.
Walter Roberson
2013-8-26
In all three of your functions, teta drops out when the expression is simplified, so you are doing a 3 dimensional integral of an expression with two free variables. That is certainly allowed, but suggests that possibly you might have an error in your equations.
采纳的回答
Walter Roberson
2013-8-27
0 个投票
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. )
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
