Getting NaN while using the ratio of power and factorial

조회 수: 9 (최근 30일)
Yash Khandelwal
Yash Khandelwal 2022년 7월 13일
댓글: James Tursa 2022년 7월 13일
1
With the following piece of code, I get NaN
l=0;
lp=2;
d=200;
a2= ((-sqrt(d))^lp);
a3 =((-sqrt(d))^l);
sum1=0
for m=0:500
a1= ((2^m)*(d^m))/(factorial(m))
mul=a1*a2*a3;
sum1=sum1+mul;
end
Is there any way to modify the code to get rid of the issue?
  댓글 수: 1
James Tursa
James Tursa 2022년 7월 13일
Cross posted on stackoverflow and answered there over 4 hours ago. If you had simply checked the answer there you could have saved yourself the trouble of posting here and getting essentially the same answer.

댓글을 달려면 로그인하십시오.

답변 (4개)

Torsten
Torsten 2022년 7월 13일
편집: Torsten 2022년 7월 13일
exp(400)*200 is quite large ...
format long
l=0;
lp=2;
d=200;
a2= ((-sqrt(d))^lp);
a3 =((-sqrt(d))^l);
a1 = 1.0;
sum1 = a1;
for m=0:500
a1 = a1* 2*d/(m+1);
sum1=sum1+a1;
end
sum1 = sum1*a2*a3
sum1 =
1.044293406962028e+176
exp(400)*200
ans =
1.044293937952829e+176

Steven Lord
Steven Lord 2022년 7월 13일
m = 500;
d = 200;
The numerator and denominator of your expression (for sufficiently large m) both overflow to inf. Dividing infinity by infinity results in a NaN.
numerator = 2^m*d^m
numerator = Inf
denominator = factorial(m)
denominator = Inf
x = numerator / denominator
x = NaN
One potential approach to avoid this is to avoid explicitly computing 2^m, d^m, or factorial(m).
numeratorVector = repmat(2*d, 1, m); % prod(numeratorVector) would effectively give numerator
denominatorVector = 1:m; % prod(denominatorVector) would effectively give denominator
xVector = numeratorVector./denominatorVector;
format longg
x2 = prod(xVector)
x2 =
8.78187252741829e+166
Let's check symbolically.
numeratorSymbolic = sym(2*d)^m;
vpa(numeratorSymbolic) % Pretty big
ans = 
1.07150860718626732094842504906e+1301
denominatorSymbolic = factorial(sym(m));
vpa(denominatorSymbolic) % Also pretty big
ans = 
1.220136825991110068701238785423e+1134
x3 = vpa(numeratorSymbolic/denominatorSymbolic) % Big but not quite as big as above
x3 = 
8.7818725274182843596325298923868e+166
You could be a little more sophisticated / clever if you wanted (preemptively cancelling out factors of 2 in numeratorVector by dividing even values in denominatorVector by 2.) Or you could keep track of x for each value of m then figure out what you need to multiply it by to get x for the next value of m.

Benjamin Thompson
Benjamin Thompson 2022년 7월 13일
The factorial function output increases very fast as input increases. See "doc factorial" for details. The output is "inf" for input of 171 or larger.
  댓글 수: 1
Steven Lord
Steven Lord 2022년 7월 13일
But 400^m overflows for m > 118 so this whole expression becomes infinite or undefined beyond that point.
d = 200;
(2*d)^118
ans = 1.1043e+307
(2*d)^119
ans = Inf

댓글을 달려면 로그인하십시오.


Vijeta Singh Yadav
Vijeta Singh Yadav 2022년 7월 13일
Data overflow and underflow problem exists in the code
visit this link for avoid data overflow and underflow problems.

카테고리

Help CenterFile Exchange에서 Debugging and Analysis에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by