Derivation of formula for $\pi$

Question

There is a well known formula for approximating the value $\pi$, namely

$$\pi = 4 \times arctg 1 = 4 \int^1_0 (arctg x)' dx = 4 \int^1_0 \frac{1}{1+x^2}dx \approx \frac{1}{N} \sum_{i=1}^N \frac{4}{1 + (\frac{i - 0.5}{N})^2}$$

Even though the formula can be encountered in many sources, I have never been able to understand the last step of the derivation. How do we arrive at $\frac{1}{N} \sum_{i=1}^N \frac{4}{1 + (\frac{i - 0.5}{N})^2}$ as an approximation for the integral? It's probably some well-known numerical method that I don't know about?

I think that the sum tries to emulate the calculation of the area as given by the function $f(x) = \frac{1}{1+x^2}$. We divide the integration range $[0; 1]$ into $N$ parts, and calculate the value of integral in those $N$ points. What I don't get is:

• Why do we substract $0.5$ from $i$? I've written two programs to check whether the results are more accurate when the substition is in there (i.e. when we calculate $\sum_{i=1}^N \frac{4}{1 + (\frac{i}{N})^2}$) - and indeed they are, but the disparity is not very big (for $N = 10000$ I get $\pi = 3.1413926444228837$ in the first case vs $\pi = 3.14129265192314$ in the second case).

• Why do we multiply whole sum by $1/N$?

Answer 1

In geometric and graphical terms, we are approximating the area enclosed by a witch of Agnesi and the $x$-axis over a certain interval, by decomposing such region into a union of rectangles with the same base length:

Can you see it now (pretty literally)?

Answer 2

This is the disappointingly–little-known midpoint rule for approximating integrals, namely $$ \int_a^b f(x) \, dx = h \sum_{n=0}^{(b-a)/h-1} f\left( a + \left(n+\frac{1}{2}\right)h \right) + R, $$ where the error $R$ is bounded by $$ \lvert R \rvert \leq \frac{(b-a)}{24}h^2 \max_{[a,b]} \lvert f''(x) \rvert, $$ (see e.g. here for a derivation). This is half that of the trapezium rule. In my opinion, one should teach it instead of the trapezium rule.