Računanje pi
V izhodišče koordinatnega sistema postavimo krog s polmerom 1 in kvadrat s polmerom 1, se pravi s stranico 2.
Ploščina kroga je $\pi r^2$, torej $\pi$, ploščina kvadrata pa $a^2$, torej $4$. Delež kvadrata, ki je pokrit s krogom, znaša $\pi / 4$.
Izberimo si dve naključni koordinati znotraj kvadrata, se pravi dve naključni števili med -1 in 1. Verjetnost, da ta točka leži tudi znotraj kroga je enaka $\pi / 4$.
To ponovimo tisočkrat. Znotraj kvadrata bo ležalo približno $n = 1000\times \pi / 4$ točk. Če ne bi vedeli, koliko je $pi$, bi lahko v resnici naredili tale poskus; z njimi bi dobili gornji $n$, iz tega $n$ pa lahko, če obrnemo formulo naokrog, izračunamo vrednost $\pi$.
Obvezna naloga
Napiši program, ki izžreba 1000 naključnih koordinat znotraj kvadrata, (tiho) šteje, koliko jih je znotraj kroga in na koncu izpiše na ta način izračunani $\pi$.
To bodi vam v pomoč: če na začetku programa napišete from random import *, se pojavi funkcija random(), ki vrne naključno število med 0 in 1. Kako to spremeniti v števila med -1 in 1, razmislite sami. Prav tako se sami spopadite z matematiko.
Rešitev
from random import random
n = 1000
i = 0
v_krogu = 0
while i < n:
x = 2 * random() - 1
y = 2 * random() - 1
if x ** 2 + y ** 2 < 1:
v_krogu += 1
i += 1
print("π = ", v_krogu / n * 4)
Naključno število med -1 in 1 lahko naračunamo tako, da naključno število med 0 in 1, ki ga vrne random() pomnožimo z 2, s čimer ga raztegnemo na naključno število med 0 in 2, nato pa od tega odštejemo 1 in smo na intervalu med -1 in 1. Vsaj en študent je odkril funkcijo uniform(a, b), ki vrne naključno število med -1 in 1. V našem primeru torej potrebujemo `uniform(-1, 1).
Število poskusov smo zaradi lepšega shranili v n. Namesto tega bi lahko na vseh mestih v programu namesto n pisali 1000.
Česa drugega pa o gornjem programu ni povedati.
Dodatna naloga
Spremeni program tako, da doda novih 1000 naključnih izbir. In še 1000. In še 1000. In še 1000. In še... To počne toliko časa, dokler se ne zgodi, da se ob nekem dodajanju novih 1000 poskusov $\pi$ spremeni za manj kot 0.000001.
Rešitev
n = 1000
pi_prej = 0
pi = 1
v_krogu = 0
vseh = 0
while abs(pi - pi_prej) > 0.000001:
pi_prej = pi
i = 0
while i < n:
x = 2 * random() - 1
y = 2 * random() - 1
if x ** 2 + y ** 2 < 1:
v_krogu += 1
i += 1
vseh += 1
pi = v_krogu / vseh * 4
print(pi)
Tule imamo dva pi-ja - pi in pi_prej. Prva je svež približek πja, druga pa prejšnji. Program teče, dokler je razlika med njima večja od dovoljene. V vseh bomo šteli, koliko poskusov smo že naredili, v vseh krogih skupaj.
Čisto dodatna naloga (ne šteje nikamor)
Tisti, ki so se ob gornjih nalogah dolgočasili, lahko izračunajo $\pi$ s simulacijo Buffonove igle. Da bo enostavneje, naj bo igla enako široka kot proge na tleh, tako da velja $p = 2 / \pi$.
Rešitev
Tole je pravzaprav lažje od dodatne naloge - če znamo majčkeno trigonometrije.
Delali se bomo, da je igla dolga 1, prav tako, seveda, širina prog.
Žrebali bomo koordinato x levega konca igle in kot. Ker smo rekli, da gre za levo koordinato, mora biti kot obrnjen v desno, torej med -90 in 90. Desna koordinata je enaka levi + kosinusu kota (krat dolžina igle, ki pa je tako ali tako 1). Ali igla seka črto, bomo preverili tako, da bomo primerjali navzdol zaokroženi koordinati x levega in desnega konca igle. Za zaokrožanje navzdol lahko posrbi kar funkcija int.
from random import *
from math import *
n = 10000
seka = 0
i = 0
while i < n:
x0 = uniform(-10, 10)
kot = uniform(radians(-90, 90))
x1 = x0 + cos(kot)
if int(x0) != int(x1):
seka += 1
i += 1
p = seka / n
print(2 / p)
Namesto uniform(radians(-90, 90)) lahko tisti, ki poznamo radiane, pišemo kar uniform(-pi / 2, pi / 2). Ob tem bi sicer kdo poskočil, da pri računanju pi-ja torej uporabljamo pi?!, vendar ga umirimo, da tako ali tako uporabljamo celo kotne funkcije, torej itak goljufamo na veliko. Sploh pa imamo za računanje pi-ja itak boljše metode od te.
Program lahko poenostavimo tako, da naključno žrebamo le razdaljo od neke črte. Leva koordinata x bo torej med 0 in 1. Igla seka črto, če je desna koordinata večja od 1.
from random import *
from math import *
n = 10000
seka = 0
i = 0
while i < n:
x0 = uniform(0, 1)
kot = pi / 2 * random()
if x0 + cos(kot) > 1:
seka += 1
i += 1
p = seka / n
print(2 / p)