# Why i have this errors in my triple integral ? change numeric methods ?

1 view (last 30 days)
Lucas Pollito on 13 Nov 2015
Commented: Lucas Pollito on 1 Dec 2015
Is this triple integral, i want to get the I_value, but for ss and T below, i have many erros that the I_value is NAN. Should i use gausslegendre for integral and newton_raphson for fzero ?
function I_value = doit
ss=0.01;
T=0.1;
tic
f2 = @(r,b,g) 1./(r.^2.*sqrt(1 - (b./r).^2 - (g^-2)*((2/15)*(ss)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8* r).*(1./(r - 1).^8 - 1./(r + 1).^8)) - (ss)^3 *(1./(r -1).^3 - 1./(r + 1).^3 - 3./(2* r).* (1./(r - 1).^2 - 1./(r + 1).^2)))));
% The folloing function only works sith scalar b and g values.
X_scalar_b_scalar_g = @(b,g)real(pi - 2*b*integral(@(r)f2(r,b,g),rmin(b,g,ss),Inf,'AbsTol',1e-3,'RelTol',1e-3));
% Make X work with array inputs for b and a scalar g value.
X_scalar_g = @(b,g)arrayfun(@(b)X_scalar_b_scalar_g(b,g),b);
f3 = @(b,g) 2*(1 - cos(X_scalar_g(b,g))).*b;
qQd_scalar_g = @(g)integral(@(b)f3(b,g),0,10,'AbsTol',1e-3,'RelTol',1e-3);
% Make qQd_scalar_g work with array g inputs.
qQd = @(g)arrayfun(qQd_scalar_g,g);
f4 =@(g) g.^5.*qQd(g)./(exp(g.^2/T));
I_value = (1/T^3)*integral(f4,0,5,'AbsTol',1e-3,'RelTol',1e-3)
toc
function r = rmin(b,g,ss)
f1 = @(r) 1 - (b./r).^2 - (g^-2)*((2/15)*(ss)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8*r).*(1./(r - 1).^8 - 1./(r + 1).^8)) -(ss)^3 *(1./(r-1).^3 - 1./(r+1).^3 - 3./(2*r).*(1./(r-1).^2 - 1./(r+1).^2)));
r = fzero(f1,1.1);
and the error is
Exiting fzero: aborting search for an interval containing a sign change
because no sign change is detected during search.
Function may not have a root.
Warning: Infinite or Not-a-Number value encountered.
> In funfun\private\integralCalc>midpArea at 397
In funfun\private\integralCalc at 105
In integral at 88
In ajuda2>@(b,g)real(pi-2*b*integral(@(r)f2(r,b,g),rmin(b,g,ss),Inf,'AbsTol',1e-3,'RelTol',1e-3)) at 7
In ajuda2>@(b)X_scalar_b_scalar_g(b,g) at 9
In ajuda2>@(b,g)arrayfun(@(b)X_scalar_b_scalar_g(b,g),b) at 9
In ajuda2>@(b,g)2*(1-cos(X_scalar_g(b,g))).*b at 10
In ajuda2>@(b)f3(b,g) at 11
In funfun\private\integralCalc>iterateScalarValued at 314
In funfun\private\integralCalc at 75
In integral at 88
In ajuda2>@(g)integral(@(b)f3(b,g),0,10,'AbsTol',1e-3,'RelTol',1e-3) at 11
In ajuda2>@(g)arrayfun(qQd_scalar_g,g) at 13
In ajuda2>@(g)g.^5.*qQd(g)./(exp(g.^2/T)) at 14
In funfun\private\integralCalc>iterateScalarValued at 314
In funfun\private\integralCalc at 75
In integral at 88
In ajuda2 at 15
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 3.5e-02. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.

Mike Hosea on 18 Nov 2015
Well, theoretically like this:
ss=0.6;
T=0.1;
a = 0.0001;
f2 = @(r,b,g) 1./(r.^2.*sqrt(1 - (b./r).^2 - (g^-2)*(2/15*(ss)^9)));
% The folloing function only works sith scalar b and g values.
X_scalar_b_scalar_g = @(b,g)real(pi - 2*b*integral(@(r)f2(r,b,g),a,10,'AbsTol',1e-3,'RelTol',1e-3));
% Make X work with array inputs for b and a scalar g value.
X_scalar_g = @(b,g)arrayfun(@(b)X_scalar_b_scalar_g(b,g),b);
f3 = @(b,g) 2*(1 - cos(X_scalar_g(b,g))).*b; % somente especular
qQd_scalar_g = @(g)integral(@(b)f3(b,g),0,10,'AbsTol',1e-3,'RelTol',1e-3);
% Make qQd_scalar_g work with array g inputs.
qQd = @(g)arrayfun(qQd_scalar_g,g);
f4 =@(g) g.^5.*qQd(g)./(exp(g.^2/T));
I_value = (1/T^3)*integral(f4,0,5,'AbsTol',1e-3,'RelTol',1e-3);
But your problem has a = 0, and in that case I get non-finite numbers. I had to loosen the tolerances to get it to complete. Seems to be a difficult problem, or at least a difficult formulation of the problem. If some work can be done symbolically, it might be worth looking into. I really didn't spend any time thinking about the problem itself, just formally made the definitions you provided work.
Lucas Pollito on 1 Dec 2015
@Mike, in my second integral, qQd(g), for values of b and g small, the integral has a oscilatory behaviour. Which is the best method to integrate it ? you can do a plot 3d in qQd, g and b to see that. May i use quadgk ? Levin or Filon Rule ??