Matlab sum is wrong: Double symsum gives incorrect result

Question

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.29289682539682539682539682539683

This 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.392897

WxMaxima gives

float((sum(sum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10));
0.3928968253968254

The 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.

Answer 1

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/10

Finally, WxMaxima and Wolfram Alpha are both computer algebra systems that work symbolically by default. Matlab is numeric by default.

Answer 2

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.

Answer 3

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]]}]/10

gives 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.