Basel problem, numerically

Question

I want to approximate numerically $\frac{\pi^2}{6}$. My initial idea was to use $$\sum_{n=1}^\infty \frac{1}{n^2}$$ but since I'm using a computer, partial sums never seem to get close enough to the result (I only get 10 correct decimal places). I was wondering if there was a faster converging formula, or some sort of way to break through the 10 decimal mark without having to do millions of calculations.

Answer 1

As I mentioned here, if you consider$$u_n=\frac2{2n+1}+\sum_{k=1}^n\frac1{k^2},$$then $\lim_{n\to\infty}u_n=\frac{\pi^2}6$ too (obviously), but this time there is a $K\in(0,1]$ such that$$(\forall n\in\mathbb{N}):\left|\frac{\pi^2}6-u_n\right|\leqslant\frac K{n^3}.$$

Answer 2

In general, if a series is approaching a limit like $x_n=L+C/n$, then one option is to eliminate $C$ and $n$ from three consecutive terms

$$ x = L+\frac C{n-1}\\y=L+\frac Cn\\z=L+\frac C{n+1}\\ x-2y+z=\frac{2C}{n^3-n}\\ 2xz-yx-yz=\frac{2CL}{n^3-n}\\ L\approx \frac{2xz-yx-yz}{x-2y+z}$$ I think the error is $O(1/n^3)$; in this case, the furst ten numbers are

1.650000000000000
1.646825396825399
1.645833333333339
1.645429292929354
1.645235042735142
1.645130385487531
1.645069110977597
1.645030889061199
1.645005826404791
1.644988715715510

Answer 3

One insight can turn this into an alternating series: $$\sum_{n=1}^\infty\frac{(-1)^n}{n^2}=\sum_{n=1}^\infty\frac{-1}{n^2}+\sum_{n=1}^\infty\frac{2}{(2n)^2}=-\frac12\sum_{n=1}^\infty\frac{1}{n^2}$$ So $$\sum_{n=1}^\infty\frac{1}{n^2}=-2\sum_{n=1}^\infty\frac{(-1)^n}{n^2}$$ There is a simple acceleration method for alternating series: $$\begin{align}\sum_{n=1}^{\infty}(-1)^na_n&=\sum_{n=1}^k(-1)^na_n+\sum_{n=k+1}^{\infty}(-1)^na_n\\ &=\sum_{n=1}^k(-1)^na_n+\frac12(-1)^{k+1}a_{k+1}+\sum_{n=k+1}^{\infty}\frac12\left((-1)^na_n+(-1)^{n+1}a_{n+1}\right)\\ &=\sum_{n=1}^k(-1)^na_n+\frac12(-1)^{k+1}a_{k+1}+\sum_{n=k+1}^{\infty}(-1)^n\frac12\left(a_n-a_{n+1}\right)\\ &=\sum_{n=1}^k(-1)^na_n+\frac12(-1)^{k+1}a_{k+1}+\sum_{n=k+1}^{\infty}(-1)^nb_n\\ &=\sum_{n=1}^k(-1)^na_n+\frac12(-1)^{k+1}a_{k+1}+\frac12(-1)^{k+1}b_{k+1}+\sum_{n=k+1}^{\infty}(-1)^n\frac12\left(b_n-b_{n+1}\right)\\ &=\sum_{n=1}^k(-1)^na_n+\frac12(-1)^{k+1}a_{k+1}+\frac12(-1)^{k+1}b_{k+1}+\sum_{n=k+1}^{\infty}(-1)^nc_n\end{align}$$ So we can sum the alternating series directly up to a certain point, then start applying acceleration. I found that $k=45$ terms of the alternating series plus $30$ steps of acceleration, ignoring the final remaining sum, was adequate for quadruple precision.

! euler.f90
program euler
   use ISO_FORTRAN_ENV
   implicit none
   integer n, k, m, i
   real(REAL128), allocatable :: a(:)
   real(REAL128) s
   k = 45
   m = 30
   allocate(a(k+m))
   a = [(1/real(n,REAL128)**2,n=1,k+m)]
   do i = 1, m
      a(k+i:k+m) = [a(k+i)/2,((a(k+n-1)-a(k+n))/2,n=i+1,m)]
   end do
   s = sum([((-1)**n*a(n),n=1,k)])+(-1)**(k+1)*sum(a(k+1:k+m))
   s = -2*s
   write(*,*) s
end program euler

The above yielded $$\frac{\pi^2}6\approx1.64493406684822643647241516664604$$

Answer 4

Going back to Micheal's answer if we assume $$x_{0n}=\sum_{k=1}^n\frac1{k^2}\approx\frac{\pi^2}6+\frac{C_1}n$$ Then we can extract $C_1$ by $$\begin{align}x_{1n}&=n(n+1)\left(x_{0n}-x_{0,n+1}\right)=n(n+1)\left(-\frac1{(n+1)^2}\right)=-\frac n{n+1}\\ &\approx n(n+1)\left(\frac{C_1}n-\frac{C_1}{n+1}\right)=C_1\end{align}$$ And then we could use the value of $C_1$ obtained to get a better approximation of $\frac{\pi^2}6$. If we assumed the error of the new approximation was roughly $\frac{D_2}{n^2}$ then we could attempt, sort of like LutzL suggested in his comment, to extract $D_2$ similarly. But we are going to use a slightly different ansatz: $$x_{0n}=\sum_{k=1}^n\frac1{k^2}=\frac{\pi^2}6+\sum_{k=1}^m\frac{(n-1)!C_k}{k!(n+k-1)!}+R_{0nm}$$ Where we are using factorial polynomials instead of powers of $n$. Now we have $$\begin{align}x_{1n}&=-\frac n{n+1}=C_1+\sum_{k=1}^{m-1}\frac{(n+2-1)!C_{k+1}}{k!(n+2+k-1)!}+R_{1n,m-1}\\ &\cdots\\ x_{in}&=(n+2i-2)(n+2i-1)\left(x_{i-1,n}-x_{i-1,n+1}\right)=C_i+\sum_{k=1}^{m-i}\frac{(n+2i-1)!C_{k+i}}{k!(n+2i+k-1)!}+R_{in,m-i}\\ &\cdots\\ x_{mn}&=(n+2m-2)(n+2m-1)\left(x_{m-1,n}-x_{m-1,n+1}\right)=C_m+R_{mn0}\end{align}$$ So now if we approximate each $R_{in,m-i}=O\left(\frac1{n^{m-i+1}}\right)\approx0$ we can solve the equations from the bottom up to get $C_m,\dots,C_i,\dots,C_1,\frac{\pi^2}6$. There is a complication in that $x_{in}=(n+2i-2)(n+2i-1)\left(x_{i-1,n}-x_{i-1,n+1}\right)$ is an error magnifying formula so we are going to have to take special precautions to prevent roundoff error from completely overwhelming our results. In this case it was possible to observe that $y_{in}=\frac{(n+i-1)!}{n!}x_{in}$ is an integer for $i>1$ so we used $y_{in}$ for intermediate calculations and then divided in the end to recover $x_{in}$ and then used $x_{in}$ to get $C_i$ and eventually $\frac{\pi^2}6$. With $n=60$ terms of partial sum and then computing up to $C_m=C_{37}$ we got to about $\text{IEEE-754}$ quadruple precision accuracy.

program euler3
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL128), allocatable :: terms(:)
   real(REAL128) diff, den
   integer n, k, m, i, j
   k = 60
   m = 37
   allocate(terms(k+m))
   terms(1:k) = [(1/real(n,REAL128)**2,n=1,k)]
   terms(k+1) = -k
   terms(k+2:k+m) = [(n+3,n=k,k+m-2)]
   do i = 3, m
      terms(k+i:k+m) = [((terms(n+i-1)*(n+i-1)-terms(n+i)*(n+1))* &
         (n+2*i-2)*(n+2*i-1),n=k,k+m-i)]
   end do
   terms(k+1) = terms(k+1)/(k+1)
   den = 1
   do n = k+2, k+m
      den = den*(n-1)
      terms(n) = terms(n)/den
   end do
   do i = k+m, k, -1
      diff = 0
      do j = k+m, i+1, -1
         diff = (diff+terms(j))/((j-i)*(i+j-k-1))
      end do
      terms(i) = terms(i)-diff
   end do
   write(*,*) sum(terms(1:k))
end program euler3

This program produced $$\frac{\pi^2}6\approx1.64493406684822643647241516664604$$ EDIT: Whoa! It occurred to me that you can prove readily by mathematical induction that $C_i=(i-1)!(i-2)!$ for $i>1$. This means that most of the hard work of the above program can be replaced by the relatively compact expression $$\frac{\pi^2}6\approx\sum_{k=1}^n\frac1{k^2}+\frac1n-\sum_{i=2}^m\frac{(i-1)!(i-2)!(n-1)!}{i!(n+i-1)!}$$ It seemed to converge to $\text{IEEE-754}$ quadruple precision accuracy at about $n=50$, $m=45$

program euler4
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL128) sum
   integer n, k, m
   n = 50
   m = 45
   sum = 0
   do k = m, 3, -1
      sum = (sum+1)*(k-1)*(k-2)/(k*(n+k-1))
   end do
   sum = (sum+1)/(2*(n+1))
   sum = (1-sum)/n
   do k = n, 1, -1
      sum = sum+real(1,REAL128)/(k**2)
   end do
   write(*,*) sum
end program euler4

Result was $$\frac{\pi^2}6\approx1.64493406684822643647241516664603$$