Προγραμματιστική πρόκληση No 6: Υπολογισμός της τιμής του 𝜋

Υπολογίστε την τιμή του αριθμού \pi δηλαδή τον λόγο της περιφέρειας ως προς την διάμετρο ενός κύκλου.

Ένας πιθανός τρόπος είναι με την μέθοδο Μοντε Κάρλο και τυχαίους αριθμούς.

%CE%B5%CE%B9%CE%BA%CF%8C%CE%BD%CE%B1

Στο σχήμα οι κουκκίδες προκύπτουν από Ν τυχαίους αριθμούς στο διάστημα (0,1). Υπολογίζουμε το πλήθος των κόκκινων κουκίδων Ν_c και εύκολα(*) προκύπτει ο τύπος

\pi = \frac{4 Ν_c}{N}

Υπάρχουν και άλλοι τρόποι υπολογισμού, είτε φυσικοί πετώντας οδοντογλυφίδες στο πάτωμα είτε με κάποιες ακολουθίες που εύκολα θα βρεις με λίγο google.

(*) κωδική φράση των μαθηματικών όταν υπονοούν πως μετά από 3 ώρες πράξεων θα βρεις το αποτέλεσμα, οπότε καλύτερα να το βουλώσεις (αν και εδώ είναι πραγματικά εύκολο να βρεις τον τύπο μόνος σου).

2 Likes

Ενδιαφέρον πρόκληση. Σκέφτηκα να ξεκινήσω υλοποιώντας την ιδέα της εκφώνησης για την μέθοδο Μόντε Κάρλο ώστε να προσεγγίσω την τίμη του π σε 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) το αποτέλεσμα έχει μεγάλη διαφορά. Όμως όσα περισσότερα παίρνουμε τόσο καλύτερα προσεγγίζουμε το π.

1 Like

Εξαιρετικό, και πολύ καθαρός κώδικας! Μπορείς να δώσεις και κάποια αριθμητικά αποτελέσματα; Με κάποιο γράφημα ίσως; Σε τι εξυπηρετεί ο πίνακας points;

1 Like

Ο πίνακας points μπορεί να παραλειφθεί γιατί μας γεμίζει και την μνήμη. Απο 'κει και πέρα ας τρέξουμε τον κώδικα 10 φορές για n=10,000 και n=1,000,000 να δούμε τί θα βγάλουμε.

Α/Α    n=1,000   n=1,000,000
0      3.212    3.14044
1      3.084    3.141628
2      3.16     3.142376
3      3.204    3.142348
4      3.192    3.141492
5      3.228    3.141336
6      3.192    3.144912
7      3.172    3.139024
8      3.192    3.138416
9      3.212    3.140508

Από τις μετρήσεις αυτές μπορούμε να υπολογίσουμε την μέση τιμή και το στατιστικό σφάλμα.
Για το pi(1000) έχουμε :

pi = 3.1848 +- 0.0128

Για το pi(1000000):
pi = 3.141248 +- 0.000580

Προφανώς έχοντας ένα εκατομ. σημεία όχι μόνο είμαστε πολύ κοντα στο π αλλά και το σφάλμα μας είναι πολύ μικρό. (Άφησα αρκετά δεκαδικά σκόπιμα).

Υ.Γ.: Ας δούμε και μία γραφική παράσταση, με μπλε το pi(1000) και με κόκκινο το pi(1000000)

2 Likes

απλά επειδή δεν μού πολυάρεσε η χρήση τής append:

#!/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')
1 Like

Υλοποίηση με βάση την μέθοδο 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 δεκαδικών ψηφίων.
Η ακολουθία είναι η εξής :

π = (4/1) - (4/3) + (4/5) - (4/7) + (4/9) - (4/11) + (4/13) - (4/15)...
2 Likes

Και λίγη numpy

#!/usr/bin/python3

import numpy as np

def count(c):
    x, y = np.random.random(c), np.random.random(c)
    return sum(x**2 + y**2 < 1.)


def pi(n):
    return 4*count(n)/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')
3 Likes

Υλοποίηση της ακολουθίας 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 δεκαδικά):

π = 3 + 4/(2*3*4) - 4/(4*5*6) + 4/(6*7*8) - 4/(8*9*10) + 4/(10*11*12) - 4/(12*13*14) ...
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
4 Likes

Ας βάλω και εγώ την nilakantha σε numpy

#!/usr/bin/python3

import numpy as np

def nilakantha(n):
    a2 = 2.*np.arange(n) + 2.
    a3 = a2 + 1.
    a4 = a2 + 2.
    s = np.ones(n)
    s[1::2] -= 2.
    return 3. + 4.*sum(s/a2/a3/a4)

def main(n=500):
    print(nilakantha(n))
    
if __name__ == '__main__':
    import sys
    try:
        n = int(sys.argv[1])
        main(n)
    except IndexError:
        main()
    except ValueError:
        print(':P')
3 Likes

Υλοποίηση σε FORTRAN, με βάση την Wallis Formula :
wallisformula

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

2 Likes

Άσκηση: φτιάξτε σε numpy την ίδια φόρμουλα. (δείτε παραπάνω παράδειγμα για Νιλακάνθα, το οποίο είναι και λίγο πιο δύσκολο).

Γιατί numpy; Για να παιδευτούμε λίγο με την τεχνική τής ανυσματοποίησης. :stuck_out_tongue:

2 Likes

Φαντάζομαι θα είναι κάπως έτσι :

#!/usr/bin/python3

import numpy as np

def wallis(n):
    a2 = 2.*np.arange(n)+1.
    a3 = a2 + 1.
    a4 = a3 + 1.
    return 2*np.prod((a3**2)/(a2*a4))

def main(n=10000000):
    print(wallis(n))
    
if __name__ == '__main__':
    import sys
    try:
        n = int(sys.argv[1])
        main(n)
    except IndexError:
        main()
    except ValueError:
        print(':P')

Επίσης να προσθέσω υλοποίηση σε FORTRAN, της Viete Formula :
vietesformula
Η οποία βγάζει ακρίβεια 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
1 Like