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

2 views (last 30 days)

Show older comments

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>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.

##### 0 Comments

### Accepted Answer

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

Mike Hosea
on 28 Nov 2015

### More Answers (0)

### See Also

### Community Treasure Hunt

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

Start Hunting!