Decimal and very small values returning zeros

23 次查看（过去 30 天）
Bruce Griffin
Bruce Griffin 2022-7-19
编辑： Bruno Luong 2022-7-20
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 个评论
James Tursa
James Tursa 2022-7-19
You will need to show us your code. There is very little meaningful advice we can give you until you do this.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2022-7-19
编辑：Bruno Luong 2022-7-19
Ran in:
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 个评论
显示 5更早的评论
Bruno Luong
Bruno Luong 2022-7-20
编辑：Bruno Luong 2022-7-20
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 2022-7-20
thank you very much I understand that alot better and that makes alot of sense.

请先登录，再进行评论。

更多回答（2 个）

Jan
Jan 2022-7-19
编辑：Jan 2022-7-19
Ran in:
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 2022-7-19
编辑：John D'Errico 2022-7-19
Ran in:
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 个评论
Bruce Griffin
Bruce Griffin 2022-7-20
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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by