Método de arquimedes para Pi

Considere uma circunferência de diâmetro 1. O número $\pi$ é exatamente o diâmetro desta circunferência. O método de Arquimedes consiste em considerar, para cada número natural $n$, os polígonos regulares inscritos e circunscritos, com $2^n$ lados. O perímetro do polígono inscrito chamaremos de $p_n$ e o do circunscrito $P_n$. As relações $p_n< \pi < P_n$ são claras e além disso $P_n,p_n \to \pi$ quando $n\to \infty$.

Usando relações trigonométricas podemos obter as equações: $$ P_{n+1}=\frac{2P_np_n}{(P_n+p_n)} $$ $$ p_{n+1}=\sqrt{p_nP_{n+1}} $$


In [1]:
import math as m
# usaremos a biblioteca math para a raiz quadrada
def pidearquimedes(n):
    ''' Acha a aproximação inferior e superior até o n-esimo termo '''
    x= 4
    y=2*m.sqrt(2)
    k=2
    while k <= n :
        print( k ," - ", y , " < ", x, "\n")
        xnew = 2*x*y/(x+y)
        y = m.sqrt(xnew*y)
        x = xnew
        k += 1
    return (x+y)/2

In [2]:
pidearquimedes(20)


2  -  2.8284271247461903  <  4 

3  -  3.0614674589207187  <  3.3137084989847607 

4  -  3.121445152258053  <  3.1825978780745285 

5  -  3.1365484905459398  <  3.151724907429257 

6  -  3.1403311569547534  <  3.144118385245905 

7  -  3.1412772509327733  <  3.1422236299424577 

8  -  3.141513801144302  <  3.1417503691689674 

9  -  3.141572940367092  <  3.1416320807031823 

10  -  3.14158772527716  <  3.1416025102568095 

11  -  3.1415914215112  <  3.1415951177495893 

12  -  3.141592345570118  <  3.1415932696293076 

13  -  3.141592576584873  <  3.141592807599645 

14  -  3.141592634338563  <  3.1415926920922543 

15  -  3.1415926487769856  <  3.141592663215408 

16  -  3.141592652386591  <  3.1415926559961966 

17  -  3.141592653288993  <  3.141592654191394 

18  -  3.141592653514593  <  3.1415926537401933 

19  -  3.1415926535709926  <  3.141592653627393 

20  -  3.1415926535850924  <  3.1415926535991923 

Out[2]:
3.1415926535903793

O cálculo acima foi feito na aritmética de uma máquina, em ponto flutuante, aqui erros de aproximação são importantes. O programa abaixo usa uma implementação de representação decimal para qualquer precisão.


In [3]:
import decimal

def ArchPi(precision=99):
    # x: circumference of the circumscribed (outside) regular polygon
    # y: circumference of the inscribed (inside) regular polygon

    decimal.getcontext().prec = precision+1
    D=decimal.Decimal
    
    # max error allowed
    eps = D(1)/D(10**precision)
    
    # initialize w/ square
    x = D(4)
    y = D(2)*D(2).sqrt()

    ctr = D(0)
    while x-y > eps:
        xnew = 2*x*y/(x+y)
        y = D(xnew*y).sqrt()
        x = xnew
        ctr += 1
        
    return str((x+y)/D(2))

In [4]:
print(ArchPi(201))


3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381978

In [5]:
print(ArchPi(1002))


3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198902

In [ ]: