Integration output is NaN

I don't understand why my double integration is resulting in NaN.

%Author: Amit Patel

%Prob. of non eavesdropping event

clc; 
clear all; 
%close all;

%-----------------------------------

neta=0.75;     %Energy harvesting efficiency                          
rho=0.5;    %PSR parameter
T=1e-3;


%-----------------------------------
sigma_sq=1;                        
%-----------------------------------

Eb=(1:10:1000)*1e-3;                              
Eb_end=size(Eb);
in=Eb_end(1,2);

P=10.^(15./10);

d0=2;                              
m=3;                              
var_hsd=(d0^(-m));  

d1=3;                           
m=3;                              
var_hse=(d1^(-m));  

d2=3;                          
m=3;                             
var_hed=(d2^(-m));                
                                                                      
Rs=2;                             
SNRth=2^Rs; 

count=0;

g3=10^(-5);  %-30dB
%g3=0;
%------------------------------------

lambda_0=1/(var_hsd);
lambda_1=1/(var_hse);
lambda_2=1/(var_hed);

gammath=2^Rs;

Ravg_ana=[];
R_avg1=0;


for i=1:in
     
a=(1-rho).*neta.*rho.*P.*g3./(1-neta.*rho.*g3);
b=(1-rho).*Eb(i).*g3./((1-neta.*rho.*g3).*T) + sigma_sq;

   fun=@(x,z)  log2(1+x).*(lambda_1.*lambda_0.*(1-neta.*rho.*g3)./(P.*neta.*rho.*P.*z)).*exp(-lambda_1.*b.*x.*sigma_sq./(P.*(1-rho) -a.*x) ).*exp(-lambda_0.*x.*Eb(i).*z./((1-neta.*rho.*g3).*T.*P)).*exp(-lambda_0.*x.*sigma_sq./P).*( Eb(i).*z./((1-neta.*rho.*g3).*T)+sigma_sq + 1./( (lambda_0.*x./P) +lambda_1.*(1-neta.*rho.*g3)./(neta.*rho.*P.*z) ) )./((lambda_0.*x./P) +lambda_1.*(1-neta.*rho.*g3)./(neta.*rho.*P.*z)).*lambda_2.*exp(-lambda_2.*z);
   R_avg= integral2(fun,0,Inf,0,Inf);
    
   Ravg_ana=[Ravg_ana R_avg]
   count=count+1
   
end;

plot(Eb,Ravg_ana,'b-'); 
hold on;
xlabel('Eb');
ylabel('Average Eavesdropping rate');
grid on;

Answer 1

[Walter Roberson]

You have a sub-expression

./(neta.*rho.*P.*z)

When z is 0, that is a division by 0. When you are working numerically and z is exactly 0 (as it is because of the initial bounds) then this ends up leading to NaN.

The expression has a limit when z goes to 0, but in the form written involves a numeric division by 0.

If you have the symbolic toolbox, then

syms x z
F = simplify(expand(fun(x,z)));
Fun = matlabFunction(F, 'vars', [x, z]);

Now you should be able to use Fun with integral2: the process of expand() and then simplify() smooths out the discontinuity -- in this particular case.

[Amit Patel]

Even if I make limits of z from 0.1 to Inf then also code results in NaN. I think It has something to do with following expression:

exp(-lambda_1.*b.*x.*sigma_sq./(P.*(1-rho) -a.*x) )

When I make remove "a.*x " term from the above expression then the problem of NaN goes away:

exp(-lambda_1.*b.*x.*sigma_sq./(P.*(1-rho)) )

Any comments?

[Walter Roberson]

Sorry, it is time for me to head to bed.