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

1 view (last 30 days)
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>vadapt at 132
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>vadapt at 132
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.

Accepted Answer

Mike Hosea
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.
  9 Comments
Lucas Pollito
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 ??

Sign in to comment.

More Answers (0)

Categories

Tags

Community Treasure Hunt

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

Start Hunting!