Decimal and very small values returning zeros

32 views (last 30 days)
I am running a integral calculation which is returning all zeros.
So i am running an integral over a meshgrid for space and time which i believe is working correctly.
My time scale is 0 to 1e-11 but when i use these values I get all zeros. If I reduce the it to 0 to 1e-9 I get numerical answers.
Is Matlab having rounding errors? If so how do I tell matlab to keep the sig figs?
t=0:2e-13:9e-11;
x=0:1:5;
[time space]=meshgrid(t,x);
fun=@(c) (pi.*time).^(-1/2).*exp(-(space-c).^2./(time)).*exp(-abs(space)/1);
answer=integral(fun,-inf,inf,'ArrayValued',true)
  1 Comment
James Tursa
James Tursa on 19 Jul 2022
You will need to show us your code. There is very little meaningful advice we can give you until you do this.

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 19 Jul 2022
Edited: Bruno Luong on 19 Jul 2022
Is is more subtle than I though.
You should not use 'ArrayValued'==true to evaluate a series of integral. The use case is to integrate a vector function depends on the domain, therefore MATLAB stop when the integrate vector converges.
  1. You should do just regular integral inside a loop, vectorize it is a very BAD idea;
  2. As your function is has the support so narrow (order of time which is tiny) and you integrate over the entire real axis, matlab indeed canot detect and miss where your function is ~= 0, and integral returns 0 for most of the time. You should help integral by narrowing the integration interval
t=0:2e-13:9e-11;
x=0:1:5;
[time, space]=meshgrid(t,x);
answer = zeros(size(time));
for i=1:numel(time)
ti = time(i);
si = space(i);
fun=@(c) (pi.*ti).^(-1/2).*exp(-(si-c).^2./(ti)).*exp(-abs(si)/1);
d = abs(log(realmin)); % 708.3964, theoretically exp(-x) for x >= d is zero numerically
lo = si-sqrt(d*ti);
hi = si+sqrt(d*ti); % fun(x) for x < lo or x > hi must be tiny
answer(i)=integral(fun,lo,hi);
end
Warning: Inf or NaN value encountered.
Warning: Inf or NaN value encountered.
Warning: Inf or NaN value encountered.
Warning: Inf or NaN value encountered.
Warning: Inf or NaN value encountered.
Warning: Inf or NaN value encountered.
answer
answer = 6×451
NaN 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 NaN 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 0.3679 NaN 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 0.1353 NaN 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 0.0498 NaN 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 0.0183 NaN 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067 0.0067
max(answer, [], 'all')
ans = 1.0000
min(answer, [], 'all') % not 0 anymore
ans = 0.0067
  7 Comments
Bruno Luong
Bruno Luong on 20 Jul 2022
Edited: Bruno Luong on 20 Jul 2022
I explain already, I repeat it here
  • As your function is has the support so narrow (order of time which is tiny) and you integrate over the entire real axis, matlab indeed cannot detect and miss where your function is ~= 0, and integral returns 0 for most of the time. You should help integral by narrowing the integration interval
If I take the support of your function in meter, 1e-13 is about the size of a nucleus of an atom, and you want to integrate on -Inf to Inf. So you tell MATLAB to look from the otherside of the universe and find where is this tiny pic and integrate it.
Try to peak a random c and evaluate your function, you'll see that most of the time it returns 0.
The lo and hi helps MATLAB to look into this smaller interval on which it integrates your function.
Bruce Griffin
Bruce Griffin on 20 Jul 2022
thank you very much I understand that alot better and that makes alot of sense.

Sign in to comment.

More Answers (2)

Jan
Jan on 19 Jul 2022
Edited: Jan on 19 Jul 2022
No, Matlab uses the standard IEEE754 conventions and has no rounding errors.
The question is funny: A huge number of scientists use Matlab for decades to get reliable results. You can be sure, that such problem would have been found before.
How do you check, if the output is zero? Remember that the standard output to the command window has a limited size of decimal figures. So maybe all you need is:
format long g
See: format
x = [25 56.31156 255.52675 9876899999];
format short
x % Some zeros are shown, which are not zeros numerically:
x = 1×4
1.0e+09 * 0.0000 0.0000 0.0000 9.8769
format long g
x
x = 1×4
1.0e+00 * 25 56.31156 255.52675 9876899999

John D'Errico
John D'Errico on 19 Jul 2022
Edited: John D'Errico on 19 Jul 2022
Um, why are you doing this numerically at all??????? With the variable of integration as c, your kernel is:
(pi*t)^(-1/2).*exp(-(x-c).^2./t).*exp(-abs(x))
with limits from -inf to inf. c enters in there only in one place. And this will be a known integral. I'll do it effectively by hand though.
syms t real positive
syms x real
syms c
K = (pi*t)^(-1/2).*exp(-(x-c).^2./t).*exp(-abs(x))
K = 
First, I'll transform the problem to make it a little more clear.
u = (x-c)/sqrt(t)
Then we will have
du = -dc/sqrt(t)
Since the limits of integration were -inf to inf, the only impact is to swap the limits, since we have x-c, but then we have an extra factor of -1 in there in the du term. So the net result is we will still integrate over the limits [-inf,inf].
Ku = 1/sqrt(pi) * exp(-u^2)*exp(-abs(x))
But what is the integral of exp(-u^2) over -inf,inf]? (Yes, I know I can do this on paper, but this is a MATLAB forum.)
syms u
int(exp(-u^2),u,[-inf,inf])
ans = 
Do you see that cancels the 1/sqrt(pi)? And since x is not a function of c, it just comes out of the integral.
The result is, EVERYTHING COLLAPSES. The value of the integral you want to compute is just
exp(-abs(x))
It is not even apparently a function of t. Using a numerical integration to solve the problem is just wrong, since you then have serious numerical problems solving it. Even if the kernel were a little more complicated, I might at worst suggest the use of a Gauss-Hermite integration. As a check, see that Bruno computed a nmerical solution.
format long g
x = (0:5)';
[x,exp(-abs(x))]
ans = 6×2
0 1 1 0.367879441171442 2 0.135335283236613 3 0.0497870683678639 4 0.0183156388887342 5 0.00673794699908547
Indeed, that is the same as what Bruno found. There is no dependence on t at all.
  1 Comment
Bruce Griffin
Bruce Griffin on 20 Jul 2022
I apologize I coded the wrong integral. Thank you for catching this for me. The real integral is
fun=@(c) (pi.*ti).^(-1/2).*exp(-(si-c).^2./(ti)).*exp(-abs(c)/1);

Sign in to comment.

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!