Triple integral with absolute value limits

I would like to solve the following integral numerically:

I wrote this code, but the result is not coming. Can anyone tell me what my mistake is?

lambda1=25/(1000)^2;
lambda1=100/(1000)^2

lambda1 = 1.0000e-04

fun =integral3(@(l0,l2,l1) (l1.^-3).*exp(-pi*(lambda2.*l2.^(2)+lambda1.*l0.^(2)))./(sqrt(1-((l0.^(2)+l2.^(2)-l1.^(2))/(2*l0.*l2)).^2)),0,Inf,0,@(l0,l2) abs(l0-l2),@(l0,l2) (l0+l2))

Error using integral3
Invalid argument list. Function requires 1 more input(s).

Answer 1

[Walter Roberson]

lambda1=25/(1000)^2;
lambda2=100/(1000)^2

lambda2 = 1.0000e-04

fun =integral3(@(l0,l2,l1) (l1.^-3).*exp(-pi*(lambda2.*l2.^(2)+lambda1.*l0.^(2)))./(sqrt(1-((l0.^(2)+l2.^(2)-l1.^(2))/(2*l0.*l2)).^2)),0,Inf,0,Inf,@(l0,l2) abs(l0-l2),@(l0,l2) (l0+l2))

Warning: Inf or NaN value encountered.

Warning: The integration was unsuccessful.

fun = NaN

[Walter Roberson]

With l0 and l2 both ranging from 0 to +inf, then when they are both 0, 4*l0^2*l2^2 would be 0 and you would have a division by 0. The numerator in that situation would be (-l1^2)^2 which would be l1^4

Is it possible that in the limit case, the overall numerator compensates? Well, in this situation l0^2*lambda1 + l2^2*lambda2 would be 0, so exp(-pi*0) which would be 1, so the overall numerator would be l1^-3. What is l1 ? Well, l1 ranges from abs(l0-l2) to l0+l2 and both of those are 0 so l1 is 0

Stepping back a second... that means that in the other part, the l1^4 numerator is going to 0^4 and the denominator is going to plain 0, so in the limit I guess that would be an overall 0, leading to sqrt(1-0) and so in the limit case the denominator might be defined after all... in the limit case.

But then when we apply that logic to the overall numerator and see that we dealing with 0^(-3) that is a "faster" division by 0 than problems in the denominator, so we get problems all over again.

But l0 == l2 an infinite number of times with those limits, so l1 would have a lower bound of 0 an infinite number of times, and that gives you the problems noted above even when l0 and l2 are not 0.

[mohammed saleh]

Thank you for your response.

In fact, the value of l2 should not be greater than 10, so l2 < l0.

Since we are dealing with an area, the terms of integration must be positive |l0-l2| .

Values of l0 can be any value from 0 to inf.

Here I am trying to calculate pdf of l1 Where the value of l1 depends on the values of l0&l2

Should the code be in this way with changing limits depending on the above data?

fun =integral3(@(l0,l2,l1) (l1.^-3).*exp(-pi*(lambda2.*l2.^(2)+lambda1.*l0.^(2)))./(sqrt(1-((l0.^(2)+l2.^(2)-l1.^(2))/(2*l0.*l2)).^2)),0,Inf,0,@(l0) l0,@(l0,l2) (l0-l2),@(l0,l2) (l0+l2))

[Walter Roberson]

@(l0) l0 would not prevent l2 from being exactly equal to l0 . You would need something like @(l0) lo*(1-eps)

[mohammed saleh]

How can I simplify this integration? If I say l1<(10+12)

Because integration is now very slow

[Walter Roberson]

@mohammed saleh

Could you confirm that you want l1 < (10+12) which is l1 < 22 -- or do you want l1 < (l0+l2) ?

The current call

  fun =integral3(@(l0,l2,l1) (l1.^-3).*exp(-pi*(lambda2.*l2.^(2)+lambda1.*l0.^(2)))./(sqrt(1-((l0.^(2)+l2.^(2)-l1.^(2))/(2*l0.*l2)).^2)),0,Inf,0,@(l0) l0,@(l0,l2) (l0-l2),@(l0,l2) (l0+l2))

already restricts l1 to be between l0-l2 (inclusive) and l0+l2 (inclusive)

[Walter Roberson]

If the problem is that the bounds for integral3() are ">=" for the lowerbound and "<=" for the upper bound, and you need ">" and "<" then:

• multiply any negative lowerbound by (1-eps) to get a result that is slightly less negative
• multiply any positive upperbound by (1-eps) to get a result that is slightly less positive
• multiply any positive lowerbound by (1+eps) to get a result that is slightly more positive
• multiply any negative upperbound by (1+eps) to get a result that is slightly more negative

This assumes that the bounds have not been reversed, that lower bound is < upperbound

[mohammed saleh]

I want l1<(l0+l2) The problem is that it is a triple integral, and as the exponent of l1 increases, the integration becomes slow, for example, l1^3 or l1^4. I will explain to you what I want to do.

The distribution BS* and BS+ in 2D space are given as follows: The BSs* are randomly distributed in a 2D area. The user is located in the center of the area and will connect with the nearest BS*. The pdf of l0 is f(l0) = 2 * (λBS*) * π * l0 * exp(-(λBS*) * π * l0^2), where l0 is the distance between the user and the nearest BS*. The BSs+ are randomly distributed in the same area. The user will connect with the nearest BSs+. The pdf of l2 is f(l2) = 2 * (λBSs+) * π * l2 * exp(-(λBSs+) * π * l2^2), where l2 is the distance between the user and the nearest (BS+). I want to find the pdf of the distance l1 between the nearest BS* and nearest BS+ to the user, where l1=sqrt(l0^2 + l2^2-l0*l2* cos(φ)) as shown in figure below The possible positioning of the nearest BS+ within the captured annular of the inner circle with a radius of l2, centred at the origin.

l1<l0+l2

l0 0 to infinity

l2 0 to infinity

How can I simplify the integration, or is there a simpler derivation than that, as I need it to find the expectation of l1^-3 and l1^-4? so I need pdf of l1

thank you so much

[Walter Roberson]

If l1 must be strictly less than l0+l2 and both of those are positive, then set the upper bound for l1 as @(l0,l2) (l0+l2)*(1-eps)

[mohammed saleh]

its not working

[Walter Roberson]

Are you getting an error? If so then what is the error?

Or is it just not as fast as you would like? Making this change to the boundary is not expected to speed up the code measureably.

I tried the calculation symbolically, but MATLAB was not able to find an analytic solution.

[mohammed saleh]

i do it like that

fun =@(l1,l2,l0) ((l1.^(1-4)).*(exp(-pi*(lambda_2*l2.^2+lambda_1*l0.^2))./(sqrt(1-((l0.^2+l2.^2-l1.^2)./(2*l0.*l2)).^2))));

const=4*pi*lambda_2*lambda_1;

z=const*integral3(fun2,0,inf,0,inf,@(l0,l2) 1+eps*(l0-l2),@(l0,l2) 1-eps*(l0+l2));

I got it, but there are big differences with the simulation result.

if iw ant to change like this

z=@(l2) const*integral2(@(l0,l1) fun2,0,inf,@(l0,l2) 1+eps*(l0-l2),@(l0,l2) 1-eps*(l0+l2)); its possible ?

i want z in tearm of l2

[Torsten]

What is fun2 ?

Further, l0min and l0max must be functions of l1 and l2 if you define the function to be integrated as

fun =@(l1,l2,l0) ((l1.^(1-4)).*(exp(-pi*(lambda_2*l2.^2+lambda_1*l0.^2))./(sqrt(1-((l0.^2+l2.^2-l1.^2)./(2*l0.*l2)).^2))));

In your case, they are functions of l0 and l2 which makes no sense.

So redefine fun as

fun =@(l0,l2,l1) ((l1.^(1-4)).*(exp(-pi*(lambda_2*l2.^2+lambda_1*l0.^2))./(sqrt(1-((l0.^2+l2.^2-l1.^2)./(2*l0.*l2)).^2))));

[mohammed saleh]

SORRY fun=fun2

its mistake

fun =@(l1,l2,l0) ((l1.^(1-4)).*(exp(-pi*(lambda_2*l2.^2+lambda_1*l0.^2))./(sqrt(1-((l0.^2+l2.^2-l1.^2)./(2*l0.*l2)).^2))));

const=4*pi*lambda_2*lambda_1;

z=const*integral3(fun,0,inf,0,inf,@(l0,l2) 1+eps*(l0-l2),@(l0,l2) 1-eps*(l0+l2));

I got it, but there are big differences with the simulation result.

if iw ant to change like this

z=@(l2) const*integral2(@(l0,l1) fun,0,inf,@(l0,l2) 1+eps*(l0-l2),@(l0,l2) 1-eps*(l0+l2)); its possible ?

i want z as functions of l2.

[mohammed saleh]

i did like that but get error Not enough input arguments.

fun =@(l0,l2,l1) ((l1.^(1-4)).*(exp(-pi*(lambda_2*l2.^2+lambda_1*l0.^2))./(sqrt(1-((l0.^2+l2.^2-l1.^2)./(2*l0.*l2)).^2))));

z = @(l2) integral2(@(l0, l1) fun(l0,l2,l1), 0,inf, @(l0,l2) 1+eps*(l0-l2),@(l0,l2) 1-eps*(l0+l2));