Pi sayısını hesaplamak için Monte Carlo tümlev alma yöntemi de kullanılabilir.

Birim çemberin alanının pi olduğunu ve çeyreğinin alanının da pi/4 olduğunu bildiğimiz için çeyrek çemberi kullanıyoruz.
Burada rastgele noktalar koyuyoruz. Daha sonra çemberin içine düşen noktalarla toplam noktaların oranını buluyoruz. Örneğin 10 nokta koymuşsak ve bunlardan 9 tanesi çemberin içine düştüyse
çemberin içine düşen nokta / toplam nokta = 0.9
olmuş oluyor.
Peki noktanın çemberin içinde olduğunu nasıl anlıyoruz.
(x,y) noktamız olsun. Çemberin denklemi x^2 + y^2 idi hatırlayacağınız üzere. Yarıçapı 1 olan çember kullandığımız için x^2+y^2 <= 1 ise nokta çemberin içinde değilse dışında oluyor.
Bu şekilde örneğin 1.000.000 kadar nokta oluşturduğumuzda çemberin içine dışına düşen noktaların oranı 0.786054 gibi bir rakam çıkıyor. Bunu dörtle çarptığımızda
pi = 4 * 0.786054 = 3.144216 gibi bir rakam çıkıyor.
Bahsettiğim gibi Monte Carlo bir tümlev alma yöntemi olduğu için tümlevi alınabilen tüm işlevlerle kullanabiliyoruz. Diğer yöntemlerle kolayca tümlevi alınamayan, ya da hesaplanması güç katlı tümlev hesabında Monte Carlo integrali kullanılıyor.
Merak edenler için Python programını da yazayım.
import math
import random
# from time import perf_counter
# burayı alttaki satırla değiştirdim
from timeit import default_timer as perf_counter
def monte_carlo(start, stop, num_points, f):
hits = 0
upp_bound = 0
for i in range(num_points):
x = random.uniform(start, stop)
y = random.uniform(0, 4) # this is cheating since I already know it's 4
if y <= f(x):
hits += 1
if y > upp_bound:
upp_bound = y
ans = (hits / num_points) * ((stop - start) * upp_bound)
return ans
def riemann(start, stop, step, f):
return sum(f(start+step*m)*step for m in range(int((stop - start)/step)))
f = lambda x: math.sqrt(16 - 16*x**2)
start = perf_counter()
a = monte_carlo(0, 1, 10000000, f)
m_time = perf_counter() - start
start = perf_counter()
b = riemann(0, 1, 1/15000, f)
r_time = perf_counter() - start
print('Monte Carlo: %s\nTime elapsed: %s' % (a, m_time))
print('Riemann sum: %s\nTime elapsed: %s' % (b, r_time))