nested integration with multiple variables

Is there away to solve this integration?

this is the latest thing I tried but it gives an error

Iout = @(y,u,v) (2*pi*lb.*y.*exp(-lb*pi.*y.^2))./(1+(t.*N.*v.^((1-ep)*a)).*y.^(ep*a).*u.^(-a))

rr = @(v) 2*pi*lb.*v.*exp(-lb*pi.*v.^2)

fun2 = integral(@(v) rr(v)*exp(integral(@(u) (-2*pi*lb.*u.*(1- (integral(@(y) Iout(y,u,v),0,v)))),0,v)),0,inf,'ArrayValued',true)

note: t,a,ep,lb,N are all arbitrary values so u can choose any value for them

error:

matrix dimensions must agree

I would appreciate if someone can spot the problem and help me with this

[Walter Roberson]

Your rr function is defined with a parameter of v, but you are calling it in fun2 with a parameter of r. Are you sure that is correct? If you do want to call it with a parameter of r then it would be constant with respect to the integration happening there so it would make more sense to move it outside of the integral.

[Hesham]

r is a typo, it should be v

[Walter Roberson]

You indicate t, a, ep, lb, N are all arbitrary values, but for our testing purposes it is important to know the correct range for them (and for r) as the choice of value will affect the convergence.

Answer 1

[Walter Roberson]

The terms v.^((1-ep)*a)) .* u.^(-a) mean that for ep > 1 and a > 0 that you end up dividing by v to a positive power times u to a positive power. As v starts at 0 and u is 0 to v, this leads to division by 0, which is making your integration impossible to calculate. You need to start at positive bound. I found that below 1E-25 that the division by 0 occurred anyhow, so you need to start above that.

Assuming that rr(r) should be rr(v):

t = 2; a = 3; ep = 5; lb =7; N = 11;
r = sin(pi/9);

LB = 1e-24;

Iout = @(y,u,v) (2*pi*lb.*y.*exp(-lb*pi.*y.^2))./(1+(t.*N.*v.^((1-ep)*a)).*y.^(ep*a).*u.^(-a));

rr = @(v) 2*pi*lb.*v.*exp(-lb*pi.*v.^2);

fun8 = @(U,v) arrayfun(@(u) integral(@(y) Iout(y,u,v),0,v), U);
fun7 = @(V) arrayfun(@(v) integral(@(u) (-2*pi*lb.*u.*(1- fun8(u,v))), LB, v), V);
fun9 = integral(@(v) rr(v).*exp(fun7(v)),LB,inf);