Interval arithmetic and series

This notebook was finished live during class on Sept. 23, 2020.

The estimate

In notes that were worked through in class, we showed that $$\sum_{k=1}^\infty \frac{1}{k^2}$$ lies in the interval obtained by interval arithmetic as below: $$\sum_{k=1}^N \frac{1}{k^2} + \left[\frac{1}{N+1},\frac{1}{N}\right].$$

Reviewing interval arithmetic

Interval arithmetic was covered in a special notebook in class. We will just briefly review what we need for this notebook.

You can import an interval arithmetic package with:

In [1]:
from mpmath import iv

You set the number of bits of precision with:

In [2]:
iv.prec = 53 # Match the precision of Python floats

You can convert integers to (degenerate) intervals with:

In [3]:
iv.mpf(7)
Out[3]:
mpi('7.0', '7.0')

If you do arithmetic and one of the numbers is an interval, then the operation will be done with interval arithmetic:

In [4]:
iv.mpf(1)/5
Out[4]:
mpi('0.19999999999999998', '0.20000000000000001')

You can generate an explicit interval with:

In [5]:
x = iv.mpf([2,3])
x
Out[5]:
mpi('2.0', '3.0')
In [6]:
x.delta
Out[6]:
mpi('1.0', '1.0')

You can use intervals in this construction. So, this is something like what we want for our remainder:

In [7]:
iv.mpf([iv.mpf(1)/6, iv.mpf(1)/5])
Out[7]:
mpi('0.16666666666666666', '0.20000000000000001')

The new interval formed above uses the lower endpoint from iv.mpf(1)/6 and the upper endpoint from iv.mpf(1)/5.

Constructing the interval

Below we write a function which returns the interval:

$$\sum_{k=1}^N \frac{1}{k^2} + \left[\frac{1}{N+1},\frac{1}{N}\right].$$
In [8]:
def interval_bound(N):
    total = iv.mpf(0)
    for k in range(1, N+1):
        total += iv.mpf(1) / (k**2)
    remainder = iv.mpf([iv.mpf(1)/(N+1), iv.mpf(1)/N])
    return total +  remainder
In [9]:
for N in [1, 2, 3, 10, 100, 1000, 1000000]:
    print(f'With N={N} we get {interval_bound(N)}')

With N=1 we get [1.5, 2.0]
With N=2 we get [1.5833333333333332593, 1.75]
With N=3 we get [1.6111111111111109384, 1.6944444444444446418]
With N=10 we get [1.640676822075630481, 1.6497677311665410738]
With N=100 we get [1.6448848902838917319, 1.6449839001849027031]
With N=1000 we get [1.6449335676804512918, 1.6449345666816703737]
With N=1000000 we get [1.6449340667367555735, 1.6449340669597960485]

Computing the series to a given accuracy

Now we want to confine the value of the series $\sum_{k=1}^\infty \frac{1}{k^2}$ to an interval of length less than $\epsilon>0$. Ignoring round off error, this amounts to making the interval $$\left[\frac{1}{N+1},\frac{1}{N}\right]$$ have length less than $\epsilon$.

Observe that the length of this interval is $$\frac{1}{N}-\frac{1}{N+1}=\frac{1}{N(N+1)} < \frac{1}{N^2}.$$ So, we can guarantee that the interval has length less than $\epsilon$ if $1/N^2 < \epsilon$. This is equivalent to saying that $N > 1/\sqrt{\epsilon}$. Since $N$ should be an integer, we take the smallest integer larger than $1/\sqrt{\epsilon}$. This is the ceiling of $1/\sqrt{\epsilon}$. Our $\epsilon$ will be a float so we can use the ceil function from the math library.

In [10]:
import math as m
In [11]:
m.ceil(m.pi)
Out[11]:
4
In [12]:
def series_sum(epsilon):
    N = m.ceil(1/m.sqrt(epsilon))
    return interval_bound(N)

We can use this to print the sum to a given number of digits.

In [13]:
d = 11 # Number of digits of precision
x = series_sum(10**(1-d))
iv.nstr(x, d)
Out[13]:
'[1.6449340668, 1.6449340669]'

The actual value is $\frac{\pi^2}{6}$ which is 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890061975870... according to the The On-Line Encyclopedia of Integer Sequences.

Looking at round off error

The round off error of a calculation depends on the precision (number of bits used in the mantissa). For this problem we have a very good understanding of the truncation error, and round off error is the remaining contribution to error.

Here we are interpreting error as the width of the interval returned. The following returns the error in our computed interval when $N$ is given.

In [14]:
def actual_error(N):
    x = interval_bound(N)
    return x.delta

We'll look into this for $N=10^5$ which should lead to a truncation error of $\frac{1}{N(N+1)} \approx 10^{-10}$. We should get 10 digits of precision.

In [15]:
m.log(10)/m.log(2)
Out[15]:
3.3219280948873626
In [16]:
53/(m.log(10)/m.log(2))
Out[16]:
15.954589770191003
In [17]:
for dps in range(10, 22, 1):
    iv.dps = dps # Set the number of digits of precision
    error = actual_error(10**5)
    print(f'When at {dps} DPS, error is {iv.nstr(error)}.')

When at 10 DPS, error is [1.4551e-6, 1.4551e-6].
When at 11 DPS, error is [1.8197e-7, 1.8197e-7].
When at 12 DPS, error is [2.2834e-8, 2.2834e-8].
When at 13 DPS, error is [1.5209e-9, 1.5209e-9].
When at 14 DPS, error is [2.7761e-10, 2.7761e-10].
When at 15 DPS, error is [1.222e-10, 1.222e-10].
When at 16 DPS, error is [1.0277e-10, 1.0277e-10].
When at 17 DPS, error is [1.0017e-10, 1.0017e-10].
When at 18 DPS, error is [1.0002e-10, 1.0002e-10].
When at 19 DPS, error is [1.0e-10, 1.0e-10].
When at 20 DPS, error is [9.9999e-11, 9.9999e-11].
When at 21 DPS, error is [9.9999e-11, 9.9999e-11].

The first thing to observe is that for large enough digits of precision, the total error matches our truncation error. For values returned close to $10^{-10}$ there is almost no round off error contribution.

Beginning at about $15$ digits or so, the error is dominated by truncation error (contribuiting about $10^{-10}$), so there is no point to increasing the precision further. For values smaller than $15$, the error is dominated by round off error.

Playing with this some more, it seems that adding $10^k$ numbers (at a similar scale) leads to a loss of about $k$ digits of precision. This makes sense since adding adding one number to another leads to an error on the order of the smallest significant digit.

To compute the quantity to about $D$ digits of precision, we'd use $N=10^{D/2}$. Adding this many numbers would cause us to loose $D/2$ digits of precision, and so we need to use $3D/2$ digits of precision for the calculation of $D$ digits.