Matlab numerical error and how to get the correct answer

162 views Asked by At

I somehow obtain the following expression in Matlab (R2014a on W7, 64b)

1/1034591578977116160000*prod(1:19)*(29576428208904825-17729494921579950*k - 20479697577410832*k^2 + 13867226524449248*k^3 - 836937224095392*k^4 - 869194297188672*k^5 + 163710902234880*k^6 + 2589894827520*k^7 - 2476912838400*k^8 + 144848704000*k^9) 

where k is initially a symbolic variable. Then I set k=10 and get the result 370371188237528 using LONGG output format. But if I put the same expression to Mathematica (substituting prod(1:19) with 19!), I get 370371188237525, which I believe is the correct answer. This seems to be a rounding error described many times on this site (is it correct?). How do I avoid it with or without Matlab symbolic toolbox?

1

There are 1 answers

2
Jommy On BEST ANSWER

There is the High Precision Float class. After adding it to the path you can simply write k = hpf(10); before evaluating your expression. The result will be

370371188237525.0106290979251578332118698122510380699168308638036

With the Symbolic Math Toolbox I would write

syms k
expr = 1/1034591578977116160000*prod(1:19)*(29576428208904825-17729494921579950*k - 20479697577410832*k^2 + 13867226524449248*k^3 - 836937224095392*k^4 - 869194297188672*k^5 + 163710902234880*k^6 + 2589894827520*k^7 - 2476912838400*k^8 + 144848704000*k^9);
subs(expr, k, 10);

which evaluates to 3150006955960150124/8505 = 370371188237525.