Asked by Lucas Pollito
on 13 Nov 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.

Answer by Mike Hosea
on 18 Nov 2015

Accepted Answer

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 19 Nov 2015

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
on 28 Nov 2015

Lucas Pollito
on 1 Dec 2015

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.