Decimal and very small values returning zeros

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)

[James Tursa]

You will need to show us your code. There is very little meaningful advice we can give you until you do this.

Answer 1

[Bruno Luong]

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

[John D'Errico]

@Bruno Luong

I note that the analytically derived result is just exp(-abs(x)). It is independent of t, as I showed.

x = 0:5;
exp(-abs(x))

ans = 1×6

1.0000 0.3679 0.1353 0.0498 0.0183 0.0067

So it looks like your computation is correct.

[Bruno Luong]

I know. The OP might use it for testing, and wonder the reasons for why MATLAB returns 0s.

[John D'Errico]

I was originally going to suggest alternatives, perhaps a Gauss-Hermite quadrature, until I looked more carefully aat the expression. But that might be another approach IF the kernel were more complicated.

[Bruno Luong]

@John D'Errico it looks like the next thing he will change is exp(-abs(si)/1);

Wonder if you have to discreatize the interval then apply Gauss quadrature and how you would do it automatically to get a good precision vs runtime.

[Bruce Griffin]

could you better explain your d lo and hi and why they are in the loop in context to -inf to inf?

Im trying to wrap my head around this but i havent dont matlab in 10+ years

[Bruno Luong]

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]

thank you very much I understand that alot better and that makes alot of sense.

Answer 2

[Jan]

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

Answer 3

[John D'Errico]

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.

[Bruce Griffin]

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);