Questions about infinity in matlab
Show older comments
Hi all，
In matlab it seems that as long as a number is greater than 1e309 it will be considered inf,but my calculations tend to generate numbers larger than 1e309, but they are actually finite and not really infinite.This caused me a very distressing problem, for example, we all know that 0 * (1e310) = 0, or 0.1 * (1e310) = 1e309, but in matlab inside their results are NAN and inf, which is very unreasonable, is not it?
I know that 1e309 is undoubtedly a very large number, but this number is only the number generated by the intermediate process of my code, it is not the final result, for example, my final result is 1e309/1e300, no doubt anyone who has studied elementary mathematics knows it is 1e9, it is finite, but matlab behaves badly, it thinks it is inf.
Anyway, this behavior of matlab makes my code often appear inf or nan which makes the program crash, I don't know what to do? Does anyone have a good solution?
5 Comments
ma Jack
on 26 Oct 2022
Jan
on 26 Oct 2022
Just a side note: Your code is hard to read. A proper indentation and avoiding empty lines after each line of code is a good idea. Using too many parentheses is overdoing.
What is the purpose of these lines:
En(pp+1,Ga_2/(4*E*E))
((abs(rx)*E)^(2*pp));
Although the processing is not concerned, human readers might have less pain reading this:
function gk1d = GK1D(rx)
% These parameters are fixed and cannot be changed
H = 3; % No unit
a = 100 * 1e-9; % m(It is measured in m)
d = 15 * a; % m
NN = 5; % No unit
w = 928.5; % cm^(-1)
qy = -pi / d;
p_para = 200; % No unit
k = w * 100;
E = max([sqrt(pi/(d*d)), real(k) / (2*H)]);
ry = 0.5 * d;
gk1d = 0;
for n = -NN:NN
KY = 2 * pi * n / d - qy;
Ga_2 = KY^2 - k^2; % Gamma^2
tt = exp(1j*KY*ry) / (4*pi*d);
for pp = 0:p_para
temmpp = tt * (-1)^pp * (abs(rx)*E)^(2*pp) * ...
(En(pp + 1, Ga_2 / (4*E*E))) / factorial(pp);
gk1d = gk1d + temmpp;
% ??? En(pp+1,Ga_2/(4*E*E))
% ??? ((abs(rx)*E)^(2*pp));
end
end
end
function op = En(n,x) % nth order exponential integral
%ExpInt(1,1)
Esto = zeros(1, n);
Esto(1,1) = expint(x);
if n==1
op = Esto(1,1);
else
for jj=2:n
Esto(1,jj) = 1/(jj-1) * (exp(-x) - x * Esto(1, jj-1));
end
op=Esto(1,n);
end
end
A simplified version of En():
function Esto = En(n, x) % nth order exponential integral
%ExpInt(1,1)
Esto = expint(x);
c = exp(-x); % expensive, so compute it once only
for jj = 2:n
Esto = (c - x * Esto) / (jj - 1);
end
end
Calculating factorial(200) is a waste of time: For doubles, factorial is saturated for numbers > 171 and INF is replied and the terms xyz/factorial(pp) vanish in the sum.
ma Jack
on 26 Oct 2022
ma Jack
on 26 Oct 2022
Accepted Answer
Jan
on 25 Oct 2022
1 vote
Welcome to the world of the IEEE754 floating point format. This is not used only in Matlab, but in many other programming languages also and it is implemented in all CPUs I know (which contain a floating point unit). This format is the world wide standard.
It uses 64 bit per element and can store numbers up to 1e309, als you have observed. Remember, that for all numbers stored in "double" format, the number of valid digits is about 16. So 1e309 is a rounded number only.
You can work with arbitrary precision instead, see e.g. the symbolic toolbox and VPA . Java's BigInt can be used from Matlab also. See also:
More Answers (2)
John D'Errico
on 25 Oct 2022
Edited: John D'Errico
on 25 Oct 2022
What to do? Learn to work with large numbers like that. Sorry, but you just do. A common solution is to use logs. Don't compute those large numbers at all, but compute their logs!
For example, consider a factorial.
factorial(500)
But the log of that factorial? There are two ways to do this. First, you can just do
sum(log(1:500))
gammaln(501)
The second recognizes that the gamma function, thus gamma(n) is the same as factorial(n-1). And therefore you can compute the log of the factorial using the gammaln function, which computes the natural log of the gamma function.
Essentially, you need to learn how to work with large numbers. Similarly,
nchoosek(2000,1000)
But you can easily enough compute the log of that same binomial coefficient as:
binln = @(N,K) gammaln(N+1) - gammaln(K+1) - gammaln(N-K+1);
So is it correct for small arguments? (The warning message here comes from nchoosek, not gammaln.)
[binln(100,40),log(nchoosek(100,40))]
As you see, they agree. But binln has no problems at all for very much larger arguments, literally as large as you wish.
binln(200000,100000)
Anyway, once you have the logs of a pair of large numbers which you know lie in the numerator and denominator of a fraction, subtract them, and only then, if necessary, exponentiate to get the result you need.
Next, a common reason why numbers get big is you are using series approximations, but trying to push the convergence of that series just too far. Sorry, but this is not at all uncommon. People push things too far, without really understanding the numerical issues. For example, the simple Taylor series for say sin(x) is globally convergent in theory. Theory is great, but often useless in practice. In practice, you cannot take enough terms in that series to get convergence for say x = 100 (radians). Just look at the terms in the series. You are raising x to extremely large powers, and if x is large, things get bad long before they ever get good.
syms x
taylor(sin(x),'order',22)
You will see massive subtractive cancellation, and simply will not be able to compute the values for that series for large absolute values of x. It would be a complete waste of time to even try that in double precision arithmetic. Instead, you can do things like the use of range reduction methods, where you can now compute a solution for at least reasonable values of x. For example:
N = 1:10;
SinApprox = @(x) sum((-1).^(N-1).*x.^(2*N-1)./factorial(2*N-1));
X = 100;
Xhat = 100 - 32*pi
format long g
[SinApprox(X),SinApprox(Xhat),sin(100)]
So while I would avoid trying to evaluate SinApprox(100), SinApprox(-0.5310) is quite well behaved. And that was just a simple range reduction trick.
Another idea is to scale your variables. Too often we see people thinking that whatever units they want to work in is ok. In fact of course that is just wrong. Would you really want to measure the distance to the nearest star in millimeters? In nanometers? Of course not! Use the appropriate units so your numbers are well scaled, well behaved. Pick units that make everything near 1 in magnitude, and you will often be happier. So the distances between the planets in the solar sytem are often measured as multiples of the distance from the earth to the sun, thus an astronomical unit. The distances between nearby stars are measured in light years, not feet or meters. But wavelengths of light are often measured in nanometers. The mass of other stars themselves are measured in terms of relative solar masses, so multiples of the mass of our sun. The use of appropriate units often makes any associated computations much better behaved.
In the end, what you need to do is simply just learn the techniques to work with such problems. You may learn some of them in a numerical analysis course, though I doubt most numerical analysis courses really cover this sort of thing explicitly or in any real depth. As well, in almost all cases, going to higher precision arithmetic is not needed, though at times, it can be a fix. I usually call such things a crutch, and say they are best avoided and IMHO only used as a very last resort. (Despite the fact that I have written several versions of high/variable precision arithmetic myself for use by those who really want them.)
In matlab it seems that as long as a number is greater than 1e309 it will be considered inf,
but my calculations tend to generate numbers larger than 1e309, but they are actually finite and not really infinite.
You may consider them to be finite, but they are too large to be stored as a finite value in double precision.
This caused me a very distressing problem, for example, we all know that 0 * (1e310) = 0, or 0.1 * (1e310) = 1e309,
In variable precision arithmetic or in higher than double precision, yes. In double precision a number greater than realmax overflows to inf and by the definition of multiplication in double precision 0 times inf results in NaN.
but in matlab inside their results are NAN and inf, which is very unreasonable, is not it?
No, it is not unreasonable. It is standard double precision behavior.
I know that 1e309 is undoubtedly a very large number, but this number is only the number generated by the intermediate process of my code, it is not the final result, for example, my final result is 1e309/1e300, no doubt anyone who has studied elementary mathematics knows it is 1e9, it is finite, but matlab behaves badly, it thinks it is inf.
In exact arithmetic or higher than double precision "anyone ... knows it is 1e9". In double precision it is infinity.
One solution is to avoid computing such large intermediate values. For example, if you were to compute binomial coefficients 
using the factorial notation you could compute 
, 
, and 
but for large n and k you'd lose precision or overflow to inf.
using the factorial notation you could compute , , and but for large n and k you'd lose precision or overflow to inf.n = 200;
k = 10;
fn = factorial(n)
fk = factorial(k)
fnk = factorial(n-k)
nCk1 = fn/(fk*fnk)
Instead you could recognize that almost all the elements of 
and 
will cancel out and avoid creating such large numerators and denominators.
and will cancel out and avoid creating such large numerators and denominators.numerator = prod((n-k+1):n)
denominator = prod(1:k)
nCk2 = numerator ./ denominator
This is roughly what nchoosek does.
nCk3 = nchoosek(n, k)
Alternately you could work with higher than double precision via Symbolic Math Toolbox, though that likely will be slower than double precision.
nCk4 = nchoosek(sym(n), sym(k))
1 Comment
Jan
on 25 Oct 2022
Just a note: Even prod((n-k+1):n) and prod(1:k) can overflow, but share a lot of prime factors. There are other alogorithms, which are less fragile:
a = n - k;
b = a + 1;
for i = 2:k
b = b + (b * a) / i; % Integer values only
end
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
