Υπολογίστε την τιμή του αριθμού \pi δηλαδή τον λόγο της περιφέρειας ως προς την διάμετρο ενός κύκλου.
Ένας πιθανός τρόπος είναι με την μέθοδο Μοντε Κάρλο και τυχαίους αριθμούς.
Στο σχήμα οι κουκκίδες προκύπτουν από Ν τυχαίους αριθμούς στο διάστημα (0,1). Υπολογίζουμε το πλήθος των κόκκινων κουκίδων Ν_c και εύκολα(*) προκύπτει ο τύπος
\pi = \frac{4 Ν_c}{N}
Υπάρχουν και άλλοι τρόποι υπολογισμού, είτε φυσικοί πετώντας οδοντογλυφίδες στο πάτωμα είτε με κάποιες ακολουθίες που εύκολα θα βρεις με λίγο google.
(*) κωδική φράση των μαθηματικών όταν υπονοούν πως μετά από 3 ώρες πράξεων θα βρεις το αποτέλεσμα, οπότε καλύτερα να το βουλώσεις (αν και εδώ είναι πραγματικά εύκολο να βρεις τον τύπο μόνος σου).
Ενδιαφέρον πρόκληση. Σκέφτηκα να ξεκινήσω υλοποιώντας την ιδέα της εκφώνησης για την μέθοδο Μόντε Κάρλο ώστε να προσεγγίσω την τίμη του π σε python
import random
def norm(u):
return (u[0]**2 + u[1]**2)**(1/2)
def pi(n):
points = []
disk = 0
for i in range(n):
x = random.random()
y = random.random()
points.append((x,y))
if norm(points[i]) <= 1 :
disk += 1
return 4*disk/n
Η ιδέα είναι να πάρω τυχαία ζεύγη (x, y) του R^2 που να ανοίκουν στο τετράγωνο όπως φαίνεται στην εικόνα της εκφώνησης. Έπειτα η βοηθητική συνάρτηση norm(x) ελέγχει αν το μέτρο(“νόρμα 2”) είναι μικρότερη ή ίση του 1, καθώς τότε το σημείο βρίσκεται στο εσωτερικό του δίσκου ακτίνας 1 και με κέντρο το Ο(0,0). Άρα μπορούμε να προσδιορίσουμε το πλήθος των σημείων εντώς δίσκου και τελικά να προσεγγίσουμε το π.
Παρατηρούμε οτι αν πάρουμε λίγα σημεία (π.χ. 1000) κάθε φορά που τρέχουμε την pi(1000) το αποτέλεσμα έχει μεγάλη διαφορά. Όμως όσα περισσότερα παίρνουμε τόσο καλύτερα προσεγγίζουμε το π.
Ο πίνακας points μπορεί να παραλειφθεί γιατί μας γεμίζει και την μνήμη. Απο 'κει και πέρα ας τρέξουμε τον κώδικα 10 φορές για n=10,000 και n=1,000,000 να δούμε τί θα βγάλουμε.
#!/usr/bin/python3
from random import random
def norm2(u):
return sum(x**2 for x in u)
def pi(n):
disk = 0
for i in range(n):
u = random(), random()
if norm2(u) < 1:
disk += 1
return 4*disk/n
def main(n=1000):
print(pi(n))
if __name__ == '__main__':
import sys
try:
n = int(sys.argv[1])
main(n)
except IndexError:
main()
except ValueError:
print(':P')
Υλοποίηση με βάση την μέθοδο Monte Carlo σε FORTRAN :
PROGRAM pi
IMPLICIT NONE
INTEGER :: seed,i,fores
REAL :: x,y,circle
CALL SYSTEM_CLOCK(seed)
CALL SRAND(seed)
PRINT*,"Dwse arithmo simeiwn : "
READ*,fores
circle=0
DO i=1,fores
x=RAND()
y=RAND()
IF ( SQRT((x**2)+(y**2)) <= 1)THEN
circle=circle+1
END IF
END DO
PRINT*,4*circle/fores
END PROGRAM pi
Και υλοποίηση με την ακολουθία Gregory-Leibniz, με υπορουτίνες πάλι σε Fortran :
PROGRAM pi2
IMPLICIT NONE
INTEGER :: n
CALL prompt(n)
CALL gregoryleibniz(n)
END PROGRAM pi2
SUBROUTINE prompt(n)
IMPLICIT NONE
INTEGER :: n
PRINT*,"Poses fores na tre3ei i akolouthia : "
READ*,n
RETURN
END SUBROUTINE prompt
SUBROUTINE gregoryleibniz(n)
IMPLICIT NONE
INTEGER :: i,n,x
REAL :: pi
x=0
pi=0
n=n*2
DO i=1,n,2
x=x+1
IF (modulo(x,2)==0)THEN
pi=pi-4/REAL(i)
ELSE
pi=pi+4/REAL(i)
END IF
END DO
PRINT*,pi
RETURN
END SUBROUTINE gregoryleibniz
Αναφέρεται ότι σε 500.000 επαναλήψεις, βρίσκει το π, με ακρίβεια 5 δεκαδικών ψηφίων.
Η ακολουθία είναι η εξής :
Υλοποίηση της ακολουθίας Gregory-Leibniz, αυτή τη φορά με Function και χωρίς ελέγχους σε μια πιο συμπαγή μορφή :
PROGRAM pi2
IMPLICIT NONE
INTEGER :: n
PRINT*,"Poses fores na tre3ei i akolouthia : "
READ*,n
print*,gregoryleibniz(n)
CONTAINS
FUNCTION gregoryleibniz(n) result(pi)
IMPLICIT NONE
INTEGER :: i,n,sign
REAL(8) :: pi
i=1
pi=0
sign=1
DO WHILE (i<=n*2)
pi=pi+sign*4/REAL(i)
sign=-sign
i=i+2
END DO
RETURN
END FUNCTION gregoryleibniz
END PROGRAM pi2
Υλοποίηση με βάση την ακολουθία Nilakantha, η οποία βρίσκει ακόμα πιο γρήγορα τα δεκαδικά ψηφία(500 επαναλήψεις-για 8 δεκαδικά):
PROGRAM pi3
IMPLICIT NONE
INTEGER :: n
PRINT*,"Poses epanalipseis : "
READ*,n
CALL nilakantha(n)
END PROGRAM pi3
SUBROUTINE nilakantha(n)
IMPLICIT NONE
INTEGER :: n,i,x,sign
REAL(8) :: pi
pi=3.0
x=0.0
i=1
sign=1
DO WHILE (i<=n)
x=x+2.0
pi=pi+sign*4/(REAL(x*(x+1)*(x+2)))
sign=-sign
i=i+1
END DO
PRINT*,pi
END SUBROUTINE nilakantha
Υλοποίηση σε FORTRAN, με βάση την Wallis Formula :
PROGRAM wallis
IMPLICIT NONE
INTEGER (8):: n
PRINT*,"Dwse arithmo epanalipsewn : "
READ*,n
CALL walli(n)
END PROGRAM wallis
SUBROUTINE walli(n)
IMPLICIT NONE
INTEGER (8):: n,i
REAL(8) :: pi,semipi
semipi=1
i=1
DO WHILE (i<=n)
semipi=semipi*(2*i*2*i)/((2*i-1)*(2*i+1))
i=i+1
END DO
pi=semipi*2
PRINT*,pi
RETURN
END SUBROUTINE walli
Σε απόδοση κυμαίνεται στα ίδια επίπεδα με την Leibniz formula
Επίσης να προσθέσω υλοποίηση σε FORTRAN, της Viete Formula :
Η οποία βγάζει ακρίβεια 6 δεκαδικών σε μόλις 15 επαναλήψεις!
PROGRAM Vietes
IMPLICIT NONE
INTEGER :: n
PRINT*,"Dwse arithmo epanalipsewn : "
READ*,n
CALL viete(n)
END PROGRAM Vietes
SUBROUTINE viete(n)
IMPLICIT NONE
INTEGER :: n,i
REAL :: x,pi,y
i=1
x=sqrt(REAL(2))
y=x/2
DO WHILE(i<=n)
pi=2/y
PRINT*,pi
x=sqrt(2+x)
y=y*x/2
i=i+1
END DO
END SUBROUTINE viete
