Why do I get an error when I ran the integral?

Hi everyone, I got an error when I run following the code, I think because of the symbolic piecewise functions or the nested integrals, or maybe something else. Could anyone help here, please?

clear all; clear; clc;


hU = 100;
NU = 10;
x0 = 150;
a1 = 11.95;
b1 = 0.135;
drm = sqrt(hU^2 + x0^2);
DML = 0.57 * drm^1.6;

% AL = zeros(size(50:10:300));
j = 1;
Rs = 50:10:300;

for i = 50:10:300
   
    xl = sqrt(i^2 - x0^2 + hU^2);
    xh = sqrt(i^2 + x0^2 + hU^2);
     syms x r
    
    f(x) = piecewise(x >= hU & x <= xl, (2 * x) / i^2, x > xl & x <= xh, (2 * x) / (pi * i^2) * acos(((x^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(x^2 - hU^2))));
    f(r) = piecewise(r >= hU & r <= xl, (2 * r) / i^2, r > xl & r <= xh, (2 * r) / (pi * i^2) * acos(((r^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(r^2 - hU^2))));
    
  

    DLM(r) = r^1.6;
    
    PL(r) = 1 / (1 + a1 * exp(-b1 * ((180/pi)*asin(hU / r) - a1)));
    PN(r) = 1 - PL(r);
    
    PL(x) = 1 / (1 + a1 * exp(-b1 * ((180/pi)*asin(hU / x) - a1)));
    PN(x) = 1 - PL(x);
    
    
    AL(j) = NU * vpa(int(PL(r) * f(r) * (int(PL(x) * f(x), r, xh) + int(PN(x) * f(x), DLM(r), xh))^((NU)-1), hU, DML), 6);
    
    j = j + 1;
  
end


plot(Rs, AL, 'b--o')
grid on
xlabel('Avg. Cell radius in meters')
ylabel('Association Probability')
legend('LOS Association Prob.')

[Walter Roberson]

clear all; clear; clc;


hU = 100;
NU = 10;
x0 = 150;
a1 = 11.95;
b1 = 0.135;
drm = sqrt(hU^2 + x0^2);
DML = 0.57 * drm^1.6;

% AL = zeros(size(50:10:300));
j = 1;
%Rs = 50:10:300;
Rs = 50:10:100;

for i = Rs
   
    xl = sqrt(i^2 - x0^2 + hU^2);
    xh = sqrt(i^2 + x0^2 + hU^2);
     syms x r
    
    f(x) = piecewise(x >= hU & x <= xl, (2 * x) / i^2, x > xl & x <= xh, (2 * x) / (pi * i^2) * acos(((x^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(x^2 - hU^2))));
    f(r) = piecewise(r >= hU & r <= xl, (2 * r) / i^2, r > xl & r <= xh, (2 * r) / (pi * i^2) * acos(((r^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(r^2 - hU^2))));
    
  

    DLM(r) = r^1.6;
    
    PL(r) = 1 / (1 + a1 * exp(-b1 * ((180/pi)*asin(hU / r) - a1)));
    PN(r) = 1 - PL(r);
    
    PL(x) = 1 / (1 + a1 * exp(-b1 * ((180/pi)*asin(hU / x) - a1)));
    PN(x) = 1 - PL(x);
    
    
    AL(j) = NU * vpa(int(PL(r) * f(r) * (int(PL(x) * f(x), r, xl, xh) + int(PN(x) * f(x), DLM(r), xh))^((NU)-1), hU, DML), 6);
    AL(j)
    j = j + 1;
  
end

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

ans = 

NaN

ans = 

NaN

ans = 

NaN

ans = 

NaN

ans = 

NaN

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

ans = 

NaN

dAL = double(AL);

plot(Rs, dAL, 'b--o')
grid on
xlabel('Avg. Cell radius in meters')
ylabel('Association Probability')
legend('LOS Association Prob.')

[Pooja Kumari]

What is the error you are getting?

[Hadeel Obaid]

Hi @Walter Roberson. Rs (the radius) should be at least 300m or 500m.

[Hadeel Obaid]

Below is the error as @Walter Roberson posted:

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

[Walter Roberson]

Several of your int() do not indicate which variable to integrate over, and are guessing incorrectly.

Answer 1

[Walter Roberson]

    f(x) = piecewise(x >= hU & x <= xl, (2 * x) / i^2, x > xl & x <= xh, (2 * x) / (pi * i^2) * acos(((x^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(x^2 - hU^2))));

There is no clause for x > xh or x < hU

[Hadeel Obaid]

Thanks @Walter Roberson, but what do you mean?

[Walter Roberson]

piecewise() returns NaN if none of the conditions are met.

You do not have conditions covering the cases x < hU or x > xh

[Hadeel Obaid]

the same problem though.

[Wafaa]

Warning: Unable to check whether the integrand

exists everywhere on the integration interval.

> In symengine

In sym/int (line 162)

In pp (line 31)

Warning: Unable to check whether the integrand

exists everywhere on the integration interval.

R(x,y)=piecewise(y <= x & x <= 1 & 0 <= y, y + 1, x < y & 0 <= x & y <= 1, x + 1);

%define xi
xi=zeros(1,m);
for i=1:m
    xi(i)=(i-1)/(m-1);   
end
%define psi(i)=L(R(x,y)) where y=xi
Lpsi=sym(zeros(1,m));
psi=sym(zeros(1,m));
for i=1:m
    % psi(i)=ff(xi(i))-(R(x,xi(i))^3)*(x/2)+diff(R(x,xi(i))^3)*((x/3)-(x^2/2))+diff(R(x,xi(i))^3,2)*((x/4)-(2*x^2/4)+(x^3/3));%int(K(x,y)*(R(x,y))^3,y,0,xi(i))+int(K(x,y)*(R(x,y))^3,y,xi(i),1);  
    psi(i)=ff(xi(i))-int(K(x,y)*(R(x,y))^3,y,0,1);
    Lpsi(i)=ff(x)-int(K(x,y)*subs(psi(i),x,y)^3,y,0,1);
    % Lpsi(i)=psi(i)^3*(x/2)+diff(psi(i)^3)*((x/3)-(x^2/2))+diff(psi(i)^3,2)*((x/4)-(2*x^2/4)+(x^3/3));
end

please help me

[Walter Roberson]

your xl is mostly complex-valued, which is outside the valid piecewise() range, so you mostly get NaN

hU = 100;
NU = 10;
x0 = 150;
a1 = 11.95;
b1 = 0.135;
drm = sqrt(hU^2 + x0^2);
DML = 0.57 * drm^1.6;

% AL = zeros(size(50:10:300));
j = 1;
%Rs = 50:10:300;
Rs = 50:10:100;

for i = Rs
   
    xl = sqrt(i^2 - x0^2 + hU^2)
    xh = sqrt(i^2 + x0^2 + hU^2)
     syms x r
    
    f(x) = piecewise(x >= hU & x <= xl, (2 * x) / i^2, x > xl & x <= xh, (2 * x) / (pi * i^2) * acos(((x^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(x^2 - hU^2))))
    %f(r) = piecewise(r >= hU & r <= xl, (2 * r) / i^2, r > xl & r <= xh, (2 * r) / (pi * i^2) * acos(((r^2) - hU^2 + x0^2 - i^2) / (2 * x0 * sqrt(r^2 - hU^2))));
    
  

    DLM(r) = r^1.6;
    
    %PL(r) = 1 / (1 + a1 * exp(-b1 * ((180/pi)*asin(hU / r) - a1)));
    %PN(r) = 1 - PL(r);
    
    PL(x) = 1 / (1 + a1 * exp(-b1 * ((180/pi)*asin(hU / x) - a1)));
    PN(x) = 1 - PL(x);
    inner3 = PL(x) * f(x)
    inner1 = int(inner3, x, xl, xh)
    inner4 = PN(x) * f(x)
    inner2 = int(inner4, x, DLM(r), xh)
    outer = int(PL(r) * f(r) * (inner1 + inner2)^((NU)-1), r, hU, DML)
    AL(j) = NU * vpa(outer, 6);
    AL(j)
    j = j + 1;
  
end

xl = 0.0000e+00 + 1.0000e+02i

xh = 187.0829

f(x) = 

inner3 = 

inner1 = 

NaN

inner4 = 

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

inner2 = 

outer = 

NaN

ans = 

NaN

xl = 0.0000 +94.3398i

xh = 190

f(x) = 

NaN

inner3 = 

NaN

inner1 = 

NaN

inner4 = 

NaN

inner2 = 

NaN

outer = 

NaN

ans = 

NaN

xl = 0.0000 +87.1780i

xh = 193.3908

f(x) = 

NaN

inner3 = 

NaN

inner1 = 

NaN

inner4 = 

NaN

inner2 = 

NaN

outer = 

NaN

ans = 

NaN

xl = 0.0000 +78.1025i

xh = 197.2308

f(x) = 

NaN

inner3 = 

NaN

inner1 = 

NaN

inner4 = 

NaN

inner2 = 

NaN

outer = 

NaN

ans = 

NaN

xl = 0.0000 +66.3325i

xh = 201.4944

f(x) = 

NaN

inner3 = 

NaN

inner1 = 

NaN

inner4 = 

NaN

inner2 = 

NaN

outer = 

NaN

ans = 

NaN

xl = 0.0000 +50.0000i

xh = 206.1553

f(x) = 

inner3 = 

inner1 = 

NaN

inner4 = 

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

inner2 = 

outer = 

NaN

ans = 

NaN

dAL = double(AL);

plot(Rs, dAL, 'b--o')
grid on
xlabel('Avg. Cell radius in meters')
ylabel('Association Probability')
legend('LOS Association Prob.')