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)
Out[2]:
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))
In [5]:
print(ArchPi(1002))
In [ ]: