MATLAB: Difference in MATLAB function evaluation

approximationerrorgammagamma_function

According to the properties of gamma function, for any x > 0

But when I try to use it in MATLAB, it shows a very large error

gamma(150) - 149*gamma(149)

ans =

 -1.5240e+247

Can anybody explain the reason for this? I need it for an approximation algorithm I am trying to develop.

Thanks in advance.

Best Answer

For your use case, 150! is a very large number so the numerical precision issue starts to kick in. I'd suggest you to use gammaln instead, for example, you can do

gammaln(150)-(log(149)+gammaln(149))

and the result is indeed 0.

HTH.

Related Solutions

MATLAB: Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

format long g
zdata=[2370.1776112981,26665.5436328856,33999.0700566522,43581.3827551805,56913.3992480947,88844.7226527384,130991.934524133,15112.911231467,176285.501029454,223360.820798491,294455.827929450,35478.787278946,411258.878280649,422702.346106183,511111.835129998,633910.342810244,780340.236378774,619555.134715430,633756.837052416,788803.530398270,777517.401792431,773693.310713579,844333.304097923,150471.6450798029,20517.8541413645,28176.4764184613,29779.9220185300,584153.2552467975,757318.1239141541,984123.5255855076,127002.474102260,157119.325006092,211544.282909359,281369.717982640,271505.093553462,327050.316856829,344420.741097441,374147.438486281,421984.526771006,432866.381682883,554318.948982768,471859.782416242,584014.198110098,622301.344956696]
zdata = 1×44
           2370.1776112981          26665.5436328856          33999.0700566522          43581.3827551805          56913.3992480947          88844.7226527384          130991.934524133           15112.911231467          176285.501029454          223360.820798491           294455.82792945           35478.787278946          411258.878280649          422702.346106183          511111.835129998          633910.342810244          780340.236378774           619555.13471543          633756.837052416           788803.53039827          777517.401792431          773693.310713579          844333.304097923          150471.645079803          20517.8541413645          28176.4764184613            29779.92201853          584153.255246797          757318.123914154          984123.525585508
ydata=[473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,473.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,483.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000,493.150000000000];
xdata=[1.08768000000,4.56728200000,4.99472000000,5.392799100000,6.791900000,8.11234546400000,11.3288295200000,15.2542138200000,23.7895322600000,38.9844910700000,70.9384511000000,100.361732800000,139.249301900000,168.187221200000,275.463489000000,323.697079700000,380.152481600000,511.387441300000,718.018141700000,852.617257800000,927.814899300000,1311.46673800000,1594.53133400000,2.170800000,3.2341306200000,4.161204900000,6.9362766600000,7.13110839100000,12.3316076700000,15.1135745800000,26.7895322600000,37.2653961400000,50.5485060700000,80.8837167700000,89.4791805100000,115.361732800000,129.025620400000,135.094336600000,219.59504900000,258.623830500000,333.908380200000,328.757892400000,444.190385600000,630.380891400000,751.682901300000,871.212537500000,1157.09126200000,1521.09905800000,4.75969298600000,7.91421685000000,10.4898360300000,13.7855398600000,23.1085904900000,30.23414500000,44.5823445700000,67.94350200000,92.278451600000,120.278265300000,122.956453700000,281.618343800000,333.757892400000,419.290040000000,555.402383500000,628.192380300000,831.755300000];
size(xdata)
ans = 1×2
     1    65
size(ydata)
ans = 1×2
     1    61
N = min([length(xdata), length(ydata), length(zdata)]);
data = [xdata(1:N); ydata(1:N)].';  %x y must be columns
fun= @(a,data)(1./a(1)).*(((data(:,1).*exp(a(2)./(8.314.*data(:,2)))./a(3)).^(1./a(4)))+sqrt((((data(:,1).*exp(a(2)./(8.314.*data(:,2)))./a(3)).^(2./a(4)))+1)));
a0 = [3; 66666; 2e+10; 2];
lb=[0; 5555; 1e+10 ; 0];
ub=[258; 205588; 1e+15 ; 4];
a = lsqcurvefit(fun, a0, data, reshape(zdata(1:N), [], 1),lb ,ub)
Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.
a = 4×1
      7.36857657624095e-06
          76932.5702990955
          18995165527.3496
          2.83794699233607

MATLAB: Regularized hypergeometric function 1F2 within matlab

Hi Wozciech,

if you add a prod function,

part3 = prod(gamma([(1-n/2),(3/2 - n/2)]))

then you're good to go. Or, you could use the gamma duplication formula to replace your current part1/part3

part1/part3 =  x.^(1-n) * 2.^(n-1).*sqrt(pi) / (gamma(1-n/2)*gamma(3/2 - n/2))

with

part1/part3 = x.^(1-n) / gamma(2-n)

Best Answer

Related Solutions

MATLAB: Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

MATLAB: Regularized hypergeometric function 1F2 within matlab

Related Question