Integration output is NaN

4 views (last 30 days)
Amit  Patel
Amit Patel on 20 Sep 2019
Commented: Walter Roberson on 20 Sep 2019
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;

Answers (1)

Walter Roberson
Walter Roberson on 20 Sep 2019
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.
  2 Comments
Amit  Patel
Amit Patel on 20 Sep 2019
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
Walter Roberson on 20 Sep 2019
Sorry, it is time for me to head to bed.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!