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

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].$$

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]:

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]:

You can generate an explicit interval with:

In [5]:

```
x = iv.mpf([2,3])
x
```

Out[5]:

In [6]:

```
x.delta
```

Out[6]:

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]:

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

and the upper endpoint from `iv.mpf(1)/5`

.

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)}')
```

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]:

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]:

The actual value is $\frac{\pi^2}{6}$ which is `1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890061975870...`

according to the The On-Line Encyclopedia of Integer Sequences.

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]:

In [16]:

```
53/(m.log(10)/m.log(2))
```

Out[16]:

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)}.')
```

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.