What is the strange pattern in the behaviour of this approximation of pi?

Question

I have been playing around with this approximation of pi recently: $$\lim_{n\to\infty} \sum_{i=0}^{n-1} \frac{n}{n^2+i^2} = \frac{\pi}{4}$$ and although I am perfectly aware that as far as approximations of pi go, this one is pretty terrible, I still thought it might be fun to improve on it a little.

To achieve that, one must come up with an approximation of the difference between $4\cdot\sum_{i=0}^{n-1} \frac{n}{n^2+i^2}$ and $\pi$ (as a function of n). Here are the respective approximations for $n=$ 10, 100, 1 000 and 10 000 ($\pi\approx 3.1415926535$): $$4\cdot\sum_{i=0}^{10-1} \frac{10}{10^2+i^2} \approx 3.2399259889$$ $$4\cdot\sum_{i=0}^{100-1} \frac{100}{100^2+i^2} \approx 3.1515759869$$ $$4\cdot\sum_{i=0}^{1000-1} \frac{1000}{1000^2+i^2} \approx 3.1425924869$$ $$4\cdot\sum_{i=0}^{10000-1} \frac{10000}{10000^2+i^2} \approx 3.1416926519$$ and here are the respective errors of the approximations ($\pi - \pi_{\text{approx}}$) with the same values for n: $$\pi - 4\cdot\sum_{i=0}^{10-1} \frac{10}{10^2+i^2} \approx -0.0983333353 \approx -0.1 = -\frac{1}{10}$$ $$\pi - 4\cdot\sum_{i=0}^{100-1} \frac{100}{100^2+i^2} \approx -0.0099833333 \approx -0.01 = -\frac{1}{100}$$ $$\pi - 4\cdot\sum_{i=0}^{1000-1} \frac{1000}{1000^2+i^2} \approx -0.0009998333 \approx -0.001 = -\frac{1}{1000}$$ $$\pi - 4\cdot\sum_{i=0}^{10000-1} \frac{10000}{10000^2+i^2} \approx -0.0000999983 \approx -0.0001 = -\frac{1}{10000}$$ And as you can see, the error is almost exactly $\frac{1}{n}$ each time.

This means that we can improve on the approximation by correcting and adding $\frac{1}{n}$ to obtain: $$-\frac{1}{n} + 4\cdot\sum_{i=0}^{10-1} \frac{n}{n^2+i^2} \approx \pi$$

For further improvements on this approximation, I will make the (possibly wrong, but at least useful for now) assumption that $\pi - 4\cdot\sum_{i=0}^{n-1} \frac{n}{n^2+i^2}$ can be expressed as a sum of the form: $$\pi - 4\cdot\sum_{i=0}^{n-1} \frac{n}{n^2+i^2} = \sum_{k=1}^{\infty} \frac{a_k}{n^k}$$ I hope it will soon be clear why I have made this assumption.

While the first term ($-\frac{1}{n}$) was found just by looking at the error rate for different values (as it was relatively obvious), I will find the next terms with a piece of code to make my job easier (and more accurate, hopefully). Here is the code (python 3):

from mpmath import mp

mp.dps = 150  # a somwhat arbitrary decision of the precision

pi = mp.pi

def pi_approx(n):
    total = 0
    for i in range(n):
        total += mp.fraction(4*n, n**2+i**2)
    return total

def improve_approx(n_1, n_2):  
    pi_1 = pi_approx(n_1)
    pi_2 = pi_approx(n_2)


    error_1 = pi - pi_1
    error_2 = pi - pi_2

    exp = -mp.log(error_1/error_2, n_1/n_2)  # which corresponds to k in the calculation above

    coeff = 1/(error_1 * (n_1**exp))  # which corresponds to 1/a_k (so I can put it in the denominator)

    return exp, coeff

If we run the code with the numbers n_1 = $10^3$ and n_2 = $10^4$, we obtain the following results: $\text{exp} \approx -0.99993484985549576581$ and $\text{coeff} \approx -1.0006169120237216770$, which is quite reassuring as it means we have $\approx -1\cdot n^{-1}$ as the error, which is exactly what we found previously.

Now, we can update our approximations by correcting them as such: $\pi \approx -\frac{1}{n} + \sum_{i=0}^{n-1} \frac{n}{n^2+i^2}$ and continue finding new terms. By doing this, we obtain: exp $\approx -1.999999999998800345$ and coeff $\approx 6.0000000000490435656$, which I will round to $\frac{1}{6\cdot n^2}$ as the next term, which gives us $\pi \approx -\frac{1}{n} + \frac{1}{6\cdot n^2} + \sum_{i=0}^{n-1} \frac{n}{n^2+i^2}$

And we can continue and find (you can check my results), we find: exp $\approx 5.999999999999792888$ and coeff $\approx -504.00000000096160606$, which we can round as before and get: $$\pi \approx -\frac{1}{n} + \frac{1}{6\cdot n^2} - \frac{1}{504\cdot n^6} + \sum_{i=0}^{n-1} \frac{n}{n^2+i^2}$$

By continuing this process, we see the next exponents are $10$ and $14$, and that the next coefficients are $1056$ and $-384$ (all of them being rounded up or down by very little). We therefore have: $$\pi \approx -\frac{1}{n} + \frac{1}{6\cdot n^2} - \frac{1}{504\cdot n^6} + \frac{1}{1056\cdot n^{10}} - \frac{1}{384\cdot n^{14}} + \sum_{i=0}^{n-1} \frac{n}{n^2+i^2}$$

This approximation is much better than the original one: if we plug in $100$ in the original approximation we get: $$\pi \approx 3.1515759869231285559 \text{ meaning an error of around } -0.0099833333333353174603$$ and with the new approximation: $$\pi \approx 3.1415926535897932385 \text{ meaning an error of around } 2.3859012707969532008 \cdot 10^{-38}$$

But if we continue, a weird thing happens: the computer returns exp $\approx 18$, which is expected, as the exponents were jumping by steps of four each time ($2, 6, 10, 14$), but it returns coeff $\approx 41.912873$, which is much farther to the closest integer than the previous ones.

My question is: where do the numbers $6, -504, 1056, -384, \text{~}42.91$ (the numbers in the denominators) come from, and why is the last one (out of the first few) not an integer? Additionally, why do the exponents jump by steps of four each time? Or have I just missed something and all of this is completely wrong?

Answer 1

$$S_n=\sum_{k=0}^{n-1} \frac{n}{n^2+k^2} =\sum_{k=0}^{n-1} \frac{n}{(n+i k)(n-i k)} =\frac 12 \Bigg[\sum_{k=0}^{n-1}\frac{1}{n+i k} +\sum_{k=0}^{n-1}\frac{1}{n-i k}\Bigg]$$ $$S_n=\frac{i}{2} (\psi (-i n)-\psi (i n)-\psi ((1-i) n)+\psi((1+i) n))$$ Using four times the asymptotics of the digamma function and continuing with Taylor seris $$4\,S_n=\pi +\frac{1}{n}-\frac{1}{6 n^2}+\frac{1}{504 n^6}-\frac{1}{1056 n^{10}}+\frac{1}{384 n^{14}}-\frac{43867}{1838592 n^{18}}+$$ $$\frac{77683}{141312 n^{22}}-\frac{657931}{24576 n^{26}}+\frac{1723168255201}{703954944 n^{30}}-\frac{151628697551}{393216 n^{34}}+O\left(\frac{1}{n^{38}}\right)$$ and then your result if you truncate to $O\left(\frac{1}{n^{18}}\right)$ which aloow to have all numerators equal to $1$ (which is nice looking).

$$\frac 1{\frac{43867}{1838592 }}=\frac{1838592 }{43867}=41.9129$$

Answer 2

There is a very fishy step here, and you might need to work a little to correct it, but the approach has the advantage of only using standard calculus and Taylor series, if it can be fixed.

The power series at $0\leq x<1$ for $\arctan(x+\epsilon)$ gives: $$\arctan\left(x+\frac1n\right)=\arctan x +\frac1{n(x^2+1)}-\frac{x}{n^2(1+x^2)^2}+O\left(\frac1{n^3}\right).$$

This works for $n(1-x)\leq 1.$

The. For $k=0,\dots,n-1,$ you set $x=\frac{k}{n}$ and get:

$$\arctan\frac{k+1}n-\arctan \frac{k}{n}=\frac{n}{k^2+n^2}-\frac{nk}{(n^2+k^2)^2}+O\left(\frac1{n^3}\right)$$

Summing, we get $$\frac{\pi}4=\arctan(1)-\arctan(0)=\sum_{k=0}^{n-1}\left(\frac n{n^2+k^2}-\frac{nk}{(n^2+k^2)^2}\right)+O(1/n^3)$$

The summing of the $n$ cases of $O(1/n^3)$ is fishy, however.

So the real question, other than the fishy step, is how to estimate:

$$\sum_{k=0}^{n-1}\frac{nk}{(n^2+k^2)^2}=\frac1{n^2}\sum\frac{k/n}{(1+(k/n)^2)^2}$$

But that is $\frac1n$ times the Reimann sum for $$\int_{0}^{1}\frac{x\,dx}{(x^2+1)^2}=\left.-\frac12\frac1{x^2+1}\right|_0^1=\frac1{4}$$

So you get $n$ times the difference approaches $1/4.$

So, that is the result you want.

Answer 3

Notice that $$\frac n{n^2+i^2} = \frac1{1+(i/n)^2}\cdot\frac 1n = f(i/n)\cdot\frac1n\tag 1$$ with a function $f(x) = \frac1{1+x^2}$.

The right-hand side of (1) can be understood as a rectangle of width $1/n$ and height $f(i/n)$ so that $$S_n=\sum_{i=0}^{n-1} f(i/n)\cdot \frac1n$$ is the sum of $n$ such rectangles where $i/n$ is ranging in $n$ equal steps from 0 to 1.

The $S_n$ are called "Riemann sums", and if $f$ is well-behaved (Riemann integrable) over $[0,1]$, then the Riemann sums converge to the Riemann integral: $$S = \lim_{n\to\infty}S_n = \int_0^1 f(x)\,dx \tag 2$$ which is the area between $f(x)$ and the $x$-axis from 0 to 1. Solving the integral for the above function is straight forward, just substitute $x = \tan t$ and $dx = (1+\tan^2t) dt$: $$S = \int_0^1\frac1{1+x^2}\,dx = \int_{\arctan 0}^{\arctan 1} \!\!\!\!\!\!\!dt = \arctan t\Big|_{t=0}^{t=1} = \frac\pi 4 $$

---

Thus, what you are doing is basically trying to approximate this definite integral by approximating and digging into the anatomy and behaviour of the integrand, and exploiting the equi-distant decomposition of $[0,1]$.

There are plenty of methods for approximating such integrals, however there are much better / faster approximations of $\pi$ as you certainly know.

If you want to investigate numerical integration further, the Newton-Codes formulas are a good point to start and to understand the degree of the approximation on a more generic ground.