calculating an integral of complex function substituting data from excel

I have an excel data whose components are all numbers ranging from 2.9 to 3(ex) 2.90234 2.90343 ~ 2.98021).

I want to import these data as variable 'e'. Also, I have other parameters and substituting those values including 'e', I want to get a series of resulting values of 'I' and eventually plot them. However, when I run the code below, I get NaN for 'I' for every data. What should I do?

x = 3.01;
k = 8.62*10^-5;
T = 160;
a = 0.2;
fun = @(x,k,T,a,b,e) exp(-b.^2./(2.*a.^2)).*heaviside(x-b-e).*exp(-(x-b-e)./(k.*T));
I = integral(@(b) fun(x,k,T,a,b,e),0,Inf)
%plot(e,I)

[dpb]

>> x = 3.01;
k = 8.62*10^-5;
T = 160;
a = 0.2;
fun = @(x,k,T,a,b,e) exp(-b.^2./(2.*a.^2)).*heaviside(x-b-e).*exp(-(x-b-e)./(k.*T));
e=2.903;fun(x,k,T,a,b,e)
Unrecognized function or variable 'b'. 
>> 

so you're missing more than just e

Use readmatrix or similar to read the spreadsheet.

Then, of course, check that your function actually produces what you expect it to for a range of cases first...

Answer 1

[Walter Roberson]

format long g
e = [2.90234, 2.90343, 2.98021];

x = 3.01;
k = 8.62*10^-5;
T = 160;
a = 0.2;
fun = @(x,k,T,a,b,e) exp(-b.^2./(2.*a.^2)).*heaviside(x-b-e).*exp(-(x-b-e)./(k.*T));
F = @(b) fun(x,k,T,a,b,e);
I = integral(F,0,Inf, 'arrayvalued', true);

Warning: Infinite or Not-a-Number value encountered.

fplot(F, [0 100])

F(9), F(10)

ans = 1×3

0 0 0

ans = 1×3

NaN NaN NaN

 b = 10,[exp(-b.^2./(2.*a.^2)), heaviside(x-b-e), exp(-(x-b-e)./(k.*T))]

b = 
    10

ans = 1×7

0 0 0 0 Inf Inf Inf

 [(x-b-e), (k.*T)], (x-b-e)./(k.*T)

ans = 1×4

-9.89234 -9.89343 -9.97021 0.013792

ans = 1×3

-717.252030162413 -717.331061484919 -722.898056844548

It appears that what is happening is that in some cases where heaviside() would be 0, that the corresponding exp() term is infinite because of the value inside the exp() becomes a large positive number. However, 0 times infinity is NaN.

I would suggest a piecewise() would be in order,

 fun = @(x,k,T,a,b,e) piecewise(x-b-e > 0, exp(-b.^2./(2.*a.^2)).*exp(-(x-b-e)./(k.*T)), 0);
 syms b
 F = arrayfun(@(E) fun(x,k,T,a,b,E), e)

F = 

 I = vpaintegral(F, b, 0, inf)

I = 

[Wonkyung Choi]

Thank you all! I did not know that zero*inf is simply NaN. I've solved the problem by just setting the range for b to be [0,1]. Again thanks a lot!