A Recursive Process

(Note: this post is an extension on the calculating pi with python post from a couple of years back. Also here’s another way to inefficiently calculate pi with Buffon’s Needles.)

We’re currently working with Power series and Taylor series in Calculus.
One particularity pretty derivation is going from the series for to the series for
Even better you can use this formula to calculate pi, since , so .
How quickly does this converge to pi? Let’s find out.
Here’s the first ten partial sums:

n= 0 and the partial sum is 4.0
n= 1 and the partial sum is 2.666666666666667
n= 2 and the partial sum is 3.466666666666667
n= 3 and the partial sum is 2.8952380952380956
n= 4 and the partial sum is 3.3396825396825403
n= 5 and the partial sum is 2.9760461760461765
n= 6 and the partial sum is 3.2837384837384844
n= 7 and the partial sum is 3.017071817071818
n= 8 and the partial sum is 3.2523659347188767
n= 9 and the partial sum is 3.0418396189294032
n= 10 and the partial sum is 3.232315809405594

Ok, that’s slow.
Here’s the first 100 terms:

n= 0 and the partial sum is 4.0
n= 10 and the partial sum is 3.232315809405594
n= 20 and the partial sum is 3.189184782277596
n= 30 and the partial sum is 3.1738423371907505
n= 40 and the partial sum is 3.1659792728432157
n= 50 and the partial sum is 3.1611986129870506
n= 60 and the partial sum is 3.157984995168666
n= 70 and the partial sum is 3.155676462307475
n= 80 and the partial sum is 3.1539378622726155
n= 90 and the partial sum is 3.1525813328751204
n= 100 and the partial sum is 3.1514934010709914

Still not to 3.14 even after adding 100 terms of this sequence?
Here’s the first 10,000 terms:

n= 0 and the partial sum is 4.0
n= 1000 and the partial sum is 3.1425916543395442
n= 2000 and the partial sum is 3.1420924036835256
n= 3000 and the partial sum is 3.1419258758397897
n= 4000 and the partial sum is 3.1418425911015158
n= 5000 and the partial sum is 3.1417926135957908
n= 6000 and the partial sum is 3.1417592924821482
n= 7000 and the partial sum is 3.141735490326666
n= 8000 and the partial sum is 3.1417176379662446
n= 9000 and the partial sum is 3.1417037523562383
n= 10000 and the partial sum is 3.1416926435905346

Two factors of ten more than the last set and still not to 5 digits of precision?
Lastly, here’s the first ten million approximations:

n= 0 and the partial sum is 4.0
n= 1000000 and the partial sum is 3.1415936535887745
n= 2000000 and the partial sum is 3.1415931535894743
n= 3000000 and the partial sum is 3.1415929869229293
n= 4000000 and the partial sum is 3.1415929035895926
n= 5000000 and the partial sum is 3.1415928535897395
n= 6000000 and the partial sum is 3.141592820256488
n= 7000000 and the partial sum is 3.1415927964468655
n= 8000000 and the partial sum is 3.141592778589681
n= 9000000 and the partial sum is 3.141592764700862
n= 10000000 and the partial sum is 3.1415927535897814


I really don’t want to calculate this by hand if after 10,000,000 additions, we only have 7 digits of precision.

(Question: why are all these even terms an over-estimate of pi?)

Graph of the first 1000 partial sums:

Let’s compare that to Chudnovsky’s algorithm:

Chudnovsky
n= 1 3.1415926535897342076684535915782983407622332609156
n= 2 3.1415926535897932384626433835873506884758663459963
n= 3 3.1415926535897932384626433832795028841971676788547
n= 4 3.1415926535897932384626433832795028841971693993750
n= 5 3.1415926535897932384626433832795028841971693993750
n= 6 3.1415926535897932384626433832795028841971693993750
n= 7 3.1415926535897932384626433832795028841971693993750


Wow, this algorithm maxed out the 51 digits of precision after 4 iterations. Crazy.

Here’s the python code for using arctan to approximate pi:

from decimal import *

#Sets decimal to 50 digits of precision
getcontext().prec = 50

#This program uses the power series for arctan to calculate pi
#arctan(x) = sum (n = 0 to infinity) (-1)^n * (x^(2n+1))/(2n+1)
#So to calculate pi, compute (arctan (1)) = pi/4 = 1 - 1/3 + 1/5 - 1/7 +...

quarterpi = Decimal(0)
n = 0

while True:
    quarterpi = quarterpi + Decimal(-1)**n / (Decimal(2*n + 1))
    #this line allows the code to print only every 10 lines, increase this number
    #to speed up the calculations and reduce the printing
    if (n % 10 == 0):
        print((quarterpi*Decimal(4)))
    n += 1
    

And here’s the chudnovsky python code:

from decimal import *

#Sets decimal to 50 digits of precision
getcontext().prec = 50

def factorial(n):
    if n<1:
        return 1
    else:
        return n * factorial(n-1)

def chudnovskyBig(n): #http://en.wikipedia.org/wiki/Chudnovsky_algorithm
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
        k += 1
    pi = pi * Decimal(10005).sqrt()/4270934400
    pi = pi**(-1)
    return pi
print ("Chudnovsky")
for i in range(1,50):
    print ("n=",i," ", chudnovskyBig(i))