# Numerical convergence depending on summation order

I'm looking for an example of convergent series such that the numerical convergence depends on the order of summation?

Or perhaps a series of positive terms where the partial sums value depend on the order of summation?

Or perhaps a divergent sequence which numerically converges, because the crucial part gets eaten by numerical resolution?

An example which I dlislike might be $$S_n = \sum_{k=1}^{n}(-1)^k2^{64}+\frac{1}{k}$$

It diverges, but straightforward calculation using 64-bit floating point will say it is bounded.

I'm looking for something more sensible, any help will be appreciated!

• "numerical convergence"? – David Mitra Jan 16 '15 at 17:13
• @DavidMitra Should I call it something else? What I mean is that if you were to calculate the sum using a computer in an ordinary, straightforward way, then the errors would stack up and change the result significantly. – dtldarek Jan 16 '15 at 17:26
• It's fine, I thought that's what you meant. Just wanted to make sure. – David Mitra Jan 16 '15 at 17:32
• @dtldarek: Why do you dislike your example? I think it's fine; though just a straight harmonic series might be more pithy. – user14717 Jan 16 '15 at 20:31
• @NickThompson Because it's ugly and obvious. Straight harmonic series would also be bad, because even if you don't know any numerical analysis, surely the computer will have troubles representing very small terms: from some point the numbers will become $0.0$ and the sum becomes finite, so of course it won't diverge (and the order does not matter much). I would like an example that is more sneaky, something that we could expect the computer to be able to calculate, but it fails. Something along your answer is much better. – dtldarek Jan 16 '15 at 22:57

When doing floating point arithmetic, you should regard each addition as inducing a relative error of approximately the float $\epsilon$, defined as the smallest positive number such that $1+\epsilon == 1$ evaluates to true. This number can be queried using std::numeric_limits<float>::epsilon() in C++, or FLT_EPSILON in C.

So to generate a sequence whose value depends on the order of summation, find a sequence you know is $\approx 1$, and whose later terms are smaller than $\epsilon$, but can nonetheless accumulate to something $>\epsilon$. Here's some code which does just that:

#include <limits>
#include <iostream>
#include <cmath>

int main()
{
float zeta_2_ls = 0;
//Make n strictly > 1/sqrt(eps), or else the two sums will be identical:
unsigned long long n = std::sqrt(100000000/std::numeric_limits<float>::epsilon());
for(unsigned long long ii = 1; ii <= n; ++ii)
{
zeta_2_ls += 1.0f/(ii*ii);
}

float zeta_2_sl = 0;
for(unsigned long long ii = n; ii > 0; --ii)
{
zeta_2_sl += 1.0f/(ii*ii);
}

std::cout << std::setprecision(12);
std::cout << std::fixed;
std::cout << "Zeta(2), summed largest to smallest: " << zeta_2_ls << std::endl;
std::cout << "Zeta(2), summed smallest to largest: " << zeta_2_sl << std::endl;
std::cout << "Zeta(2), using pi^2/6 to float prec: " << M_PI*M_PI/6.0f << std::endl;
}


The output:

Zeta(2), summed largest to smallest: 1.644725322723
Zeta(2), summed smallest to largest: 1.644933938980
Zeta(2), using pi^2/6 to float prec: 1.644934066848


You can calculate that the ULP distance from the exact answer is only 1 for the small-to-large summed series, but 1751 for the large-to-small.

Here's another example, where the terms go to zero, the series converges, but the value of the sum depends on the number of terms

#include <stdio.h>
#include <math.h>

float saw_tooth(float x, unsigned int terms)
{
float y = 0;
for(unsigned int ii = 1; ii <= terms; ++ii)
{
if(ii & 1)
{
y += sinf(ii*x)/ii;
}
else
{
y -= sinf(ii*x)/ii;
}
}
y *= 2.0/M_PI;
return y;
}

int main()
{
float x = M_PI;
unsigned int terms = 10;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 1000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 10000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));

x = M_PI/2.0;
terms = 10;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 1000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 10000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));
terms = 100000000;
printf("Saw tooth(%g) with %d terms = %g\n", x, terms, saw_tooth(x,terms));

printf("Why is the so far from zero? sinf(M_PI) = %g\n", sinf(M_PI));
}


Output:

Saw tooth(3.14159) with 10 terms = -5.30531e-07
Saw tooth(3.14159) with 100000 terms = -0.0055615
Saw tooth(3.14159) with 1000000 terms = -0.0555674
Saw tooth(3.14159) with 10000000 terms = -0.491704
Saw tooth(3.14159) with 100000000 terms = -0.6344
Saw tooth(1.5708) with 10 terms = 0.531527
Saw tooth(1.5708) with 100000 terms = 0.499998
Saw tooth(1.5708) with 1000000 terms = 0.500003
Saw tooth(1.5708) with 10000000 terms = 0.500004
Saw tooth(1.5708) with 100000000 terms = 0.616448
Why is the so far from zero? sinf(M_PI) = -8.74228e-08


Of course, the range reduction is causing the evaluation of sinf(ii*M_PI) to become ~FLT_EPSILON, and so the terms just sorta . . . drift. But to be honest, I cannot explain why they drift so far, given the mitigating effect of the division by ii.

• This is a great answer, thank you! – dtldarek Jan 17 '15 at 22:08
• Sure thing; by the way, the reason I suspected the sawtooth would be so bad is due to randomascii.wordpress.com/2014/10/09/…, but in fact more pedestrian mechanisms came into play. Still the link will give you tons of insight into floating point arithmetic. – user14717 Jan 18 '15 at 0:44