Convergence to \(\pi\)

Finding numerical approximations to \(\pi\) has fascinated people for millenia. One famous formula is

$ \displaystyle \frac{\pi^2}{6} = 1 + \frac{1}{2^2} + \frac{1}{3^2} + \cdots. $

Say \(s_k\) is the sum of the first \(k\) terms of the series above, and \(p_k = \sqrt{6s_k}\). Here is a fancy way to compute these sequences in a compact code.

a = [1/k^2 for k=1:100] 
s = cumsum(a)        # cumulative summation
p = @. sqrt(6*s)

using Plots
plot(1:100,p,m=:o,leg=:none,xlabel="k",ylabel="p_k",title="Sequence convergence")

This graph suggests that \(p_k\to \pi\) but doesn’t give much information about the rate of convergence. Let \(\epsilon_k=|\pi-p_k|\) be the sequence of errors. By plotting the error sequence on a log-log scale, we can see a nearly linear relationship.

ep = @. abs(pi-p)    # error sequence
plot(1:100,ep,m=:o,l=nothing,
    leg=:none,xaxis=(:log10,"k"),yaxis=(:log10,"error"),title="Convergence of errors")

This suggests a power-law relationship where \(\epsilon_k\approx a k^b\), or \(\log \epsilon_k \approx b (\log k) + \log a\).

k = 1:100
V = [ k.^0 log.(k) ]     # fitting matrix
c = V \ log.(ep)         # coefficients of linear fit
2-element Array{Float64,1}:
 -0.1823752497283019
 -0.9674103233127926

In terms of the parameters \(a\) and \(b\) used above, we have

@show (a,b) = exp(c[1]),c[2];
(a, b) = (exp(c[1]), c[2]) = (0.8332885904225771, -0.9674103233127926)

It’s tempting to conjecture that \(b\to -1\) asymptotically. Here is how the numerical fit compares to the original convergence curve.

plot!(k,a*k.^b,l=:dash)