Matlab sum is wrong: Double symsum gives incorrect result
Mia Lopez
I'm trying to calculate the double sum
$$ \frac{1}{10} \sum_{x=1}^{10} \left( \frac{1}{x} \sum_{n=0}^{floor(log_{10}x)} 10^n \right).$$
In MATLAB, my result is
>> syms n x
>> vpa(symsum(symsum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10)
ans =
0.29289682539682539682539682539683This is incorrect. Am I using wrong syntax? Is it a MATLAB bug? Wolfram Alpha gives
sum(sum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10)
ans = 0.392897WxMaxima gives
float((sum(sum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10));
0.3928968253968254The results from Wolfram Alpha and WxMaxima are correct, or at least they match my hand calculations.
NOTE: I'm using 10 here for the upper limit on the first sum since it's the smallest value for which this discrepancy appears. I'm guessing it has something to do with the upper limit on the inner sum, but when I test $floor(\log_{10}x) = floor(\log(x)/\log(10))$ for different values of $x$ in MATLAB, I get the expected results.
$\endgroup$ 43 Answers
$\begingroup$You're using symbolic math, but the log(10) in the middle of your code is being evaluating numerically before any symbolic calculation, which leads to imprecision. If you evaluate just your upper bound, you'll see that
floor(log(x)/log(10))returns floor((1125899906842624*log(x))/2592480341699211) rather than floor(log(x)/log(10)). Order of operations matter when coding for symbolic math and the log function has no way of knowing that the surrounding context is symbolic so it defaults to numeric evaluation.
Instead, use:
syms n x
vpa(symsum(symsum(10^n, n, 0, floor(log(x)/log(sym(10))))/x, x, 1, 10)/10)or just:
syms n x
vpa(symsum(symsum(10^n, n, 0, floor(log10(x)))/x, x, 1, 10)/10)This can also be solved numerically via:
s = 0;
for x = 1:10 s = s+sum(10.^(0:floor(log10(x))))/x;
end
s = s/10Finally, WxMaxima and Wolfram Alpha are both computer algebra systems that work symbolically by default. Matlab is numeric by default.
$\endgroup$ 2 $\begingroup$The funny thing about your sum is that it has an exact closed form that is trivial to compute:
Sum[10^n/x, {x, 1, 10}, {n, 0, Floor[Log[10, x]]}]/10gives the result $9901/25200$ in Mathematica, because for all $1 \le x < 10$, $\lfloor \log_{10} x \rfloor = 0$, thus the inner sum is simply $10^0 = 1$ for all but the last value of $x$, in which case the sum is $10^0 + 10^1 = 11$. Therefore, the given double sum is simply $$\frac{1 + H_{10}}{10},$$ where $H_n = \sum_{k=1}^n 1/k$ is the $n^{\rm th}$ harmonic number.
$\endgroup$ $\begingroup$This is likely a numerical error coming from the fact that $\log_{10}(10) = 1$ and the approximation $$\log_{10}(x) \approx \frac{\log(x)}{\log(10)}$$ Feeding this function $x=10$ would make the slightest numerical floating point rounding error (downwards) would give us a value just slightly lower than 1 ( even if by as little as $10^{-15}$ ) which then the floor function would round down to 0 instead of 1.
$\endgroup$ 2