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

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.

Answer 1

[Mike Hosea]

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]

Mike, thanks a lot. it probably works, but i will edit my code to my real problem, and add some function that is the limit of an integral. i try do edit by myself, but doenst work. if you can help me, i think it is simple

[Mike Hosea]

I doubt anything is simple anymore. :)

[Mike Hosea]

I'm not sure what you're trying to accomplish with

f1 =  ...;
rmin = fzero(f1,1.1);
k=1;
n_iter = length(rmin); 
for i=1:n_iter
    vec = zeros(1,n_iter);
    if(isreal(rmin(i))==1)
        vec(k)=rmin(i);
        k=k+1;
     end
end

The problem is that FZERO is only capable of returning one element, and it's always going to be real, so all of that code just returns the result of the FZERO call. That said, since f1 is a function of not one but three variables, we need to call it in a context where the other two variables can be dynamically bound to whatever they happen to be when we evaluate the f1 function. It turns out that we're reaching the limits of what makes sense with anonymous functions. We need to define a sub-function to do the work, and to do that we need to make the whole thing a function. For example:

function I_value = doit
ss=0.6;
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);

[Lucas Pollito]

Hi Mike. Again, thanks a lot. It works. God Bless you ! you are a great man and a great numeric analyst. So, if i could ask for something more, take a look at this

when ss=0.01 and T = 0.1, like this

then it DONT WORK, and the error is like this :

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.

and the answer is NAN

Why occurs that ? Should i use gausslegendre for integral and newtonraphson for fzero ?

[Lucas Pollito]

Hi Mike. Again, thanks a lot. It works. God Bless you ! you are a great man and a great numeric analyst. So, if i could ask for something more, take a look at this

when ss=0.01 and T = 0.1, like this

then it DONT WORK, and the error is like this :

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.

and the answer is NAN

Why occurs that ? Should i use gausslegendre for integral and newtonraphson for fzero ?

[Lucas Pollito]

and mike, this

f1 =  ...;
rmin = fzero(f1,1.1);
k=1;
n_iter = length(rmin); 
for i=1:n_iter
    vec = zeros(1,n_iter);
    if(isreal(rmin(i))==1)
        vec(k)=rmin(i);
        k=k+1;
     end
end

is because I have many roots of the f1 equation. fzero find a lot of roots, and i want the max value of a real root. I see that you take it away from equation. Can I have it back ? rs

[Lucas Pollito]

and mike, this

f1 =  ...;
rmin = fzero(f1,1.1);
k=1;
n_iter = length(rmin); 
for i=1:n_iter
    vec = zeros(1,n_iter);
    if(isreal(rmin(i))==1)
        vec(k)=rmin(i);
        k=k+1;
     end
end

is because I have many roots of the f1 equation. fzero find a lot of roots, and i want the max value of a real root. I see that you take it away from equation. Can I have it back ? rs

[Mike Hosea]

Yes, you can easily have something like that back, but that particular section of code never did what you say/think it did. FZERO finds exactly one root or it errors. It will never find more than one root at a time. To do that, you will need to call it iteratively and do something clever to give it starting intervals that bracket each root. But if you only need the largest one, perhaps a search followed by a call to FZERO could work, but searching for a bracket on the last sign change is not a sure thing.

[Lucas Pollito]

@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 ??