Κυριακή, 22 Μαρτίου 2015

Υπολογισμός του π

Υπάρχουν διάφοροι αλγόριθμοι που υπολογίζουν τον αριθμό π. Συνήθως ο υπολογισμός γίνεται με διάφορες σειρές ή με κλάσματα. Παρακάτω θα δούμε ένα αλγόριθμο που βασίζεται στην τύχη. Όμως με την βοήθεια της θεωρίας των πιθανοτήτων το αποτέλεσμα δεν είναι τυχαίο.


Αρχικά έστω ότι έχουμε το μοτίβο του σχήματος 1.

σχήμα 1, ένα τετράγωνο διαστάσεων 1x1 και ένας κύκλος εφάπτεται στις πλευρές
του. Συνεπώς ο κύκλος έχει ακτίνα 0.5 .
Στο σχήμα έχουμε ένα τετράγωνο και έναν κύκλο μέσα σε αυτό. Το τετράγωνο έχει διαστάσεις 1x1 και ως εκ τούτου ο κύκλος έχει ακτίνα 0.5 . Επίσης, με βάση ένα σύστημα αναφοράς (του οποίου η αρχή ταυτίζεται με την κάτω αριστερή γωνία του τετραγώνου) φαίνονται τα σημεία Α, Β, Γ και Κ καθώς και οι συντεταγμένες τους.

Φανταστείτε τώρα πως παίρνετε ένα πλήθος από βόλους και τους ρίχνετε τυχαία μέσα στο τετράγωνο. Προφανώς ορισμένοι βόλοι θα πέσουν μέσα στον κύκλο ενώ άλλοι έξω από αυτόν. Έστω n το πλήθος των βόλων που ρίχνετε και s το πλήθος των βόλων που πέφτουν μέσα στον κύκλο. Είναι φανερό πως ο λόγος s/n είναι ίσος με τον λόγο εμβαδόν κύκλου/εμβαδόν τετραγώνου. Δηλαδή
Όμως γνωρίζουμε πως
Συνεπώς

Άρα γνωρίζοντας το s και το n και επιλύοντας την παραπάνω σχέση ως προς π μπορούμε να έχουμε μία εκτίμηση του π. Δηλαδή
Αυτό που πρέπει να κάνουμε είναι να προσομοιώσουμε την παραπάνω διαδικασία. Όσο μεγαλύτερο αριθμό βόλων ρίξουμε μέσα στο τετράγωνο τόσο καλύτερη εκτίμηση για το π θα έχουμε. Έτσι πρέπει να γράψουμε έναν αλγόριθμο ώστε ο υπολογιστής να πραγματοποιήσει την εξής διαδικασία.

Πως θα το κάνουμε αυτό; Πολύ απλά θα φτιάξουμε ένα πλήθος n σημείων με τυχαίες συντεταγμένες σε σχέση με το σύστημα αναφοράς του σχήματος 1. Επειδή θέλουμε όλα τα σημεία να περιέχονται στο τετράγωνο, το εύρος των x και των y τους θα κυμαίνεται από το 0 έως το 1. Έπειτα θα ελέγξουμε ποια σημεία "έπεσαν" μέσα στον κύκλο. Η εξίσωση του συγκεκριμένου κύκλου είναι
Οπότε ένα σημείο με συντεταγμένες (a,b) θα είναι εντός του κύκλου αν ισχύει η ανισότητα
Έτσι θα μετρήσουμε τον αριθμό των σημείων s για τα οποία ισχύει αυτή η ανισότητα και γνωρίζοντας πλέον και το n και το s μπορούμε να υπολογίσουμε μία εκτίμηση του π.

Ο κώδικας που υλοποιεί τον παραπάνω αλγόριθμο είναι ο ακόλουθος:

Private Sub Command1_Click()

'δήλωση των μεταβλητών
Dim x1 As Double, x2 As Double, m As Integer, s As Long

'εισαγωγή του αριθμού m που...
m = Text1.Text
'... σχετίζεται με το n με τον τύπο n=10^m
n = 10 ^ m
'θέτουμε τον μετρητή μας ίσο με μηδέν
s = 0

'ανοίγουμε ένα αρχείο για την εξαγωγή αποτελεσμάτων
Open "points.txt" For Output As #1

'ξεκινάμε μία διαδικασία n φορές
For i = 1 To n
'δημιουργούμε ένα σημείο με τυχαίες συντεταγμένες (x1, x2)
x1 = Rnd
x2 = Rnd
'αν αυτό το σημείο είναι εντός του κύκλου αυξάνουμε τον μετρητή μας κατά ένα
If (x1 - 0.5) ^ 2 + (x2 - 0.5) ^ 2 < 0.5 ^ 2 Then
s = s + 1
End If
'τυπώνουμε το σημείο στο αρχείο
Print #1, i, x1, x2
'πηγαίνουμε στο επόμενο σημείο
Next i

'κλείνουμε το αρχείο
Close

'υπολογίζουμε την εκτίμηση του π
pi_value = 4 * s / n
'υπολογίζουμε την πραγματική τιμή του π (δεν υπάρχει έτοιμη στην vb6)
pi_real = 4 * Atn(1)
'υπολογίζουμε το σχετικό σφάλμα
E = (pi_value - pi_real) / pi_real

'εξάγουμε τα αποτελέσματα στην φόρμα
'για την εκτίμηση του π κρατάμε 8 δεκαδικά ψηφία
Text2.Text = Format(pi_value, "0.00000000")
'και για το σχετικό σφάλμα 3 δεκαδικά ψηφία
Text3.Text = Format(E * 100, "0.000")

End Sub

Η φόρμα εισαγωγής και εξαγωγής δεδομένων στον κώδικα φαίνεται στην εικόνα 1.

εικόνα 1, η φόρμα του κώδικα.
Ξεκινάμε να ρίχνουμε τους βόλους στο τετράγωνο. Επιλέγουμε για αρχή m=3 (οπότε το πλήθος των βόλων είναι 1000) και εκτελούμε τον κώδικα πατώντας το κουμπί "υπολογισμός". Το αποτέλεσμα φαίνεται στην εικόνα 2.

εικόνα 2, το αποτέλεσμα μετά την πρώτη εκτέλεση του κώδικα.
Δύο πράγματα πρέπει να σημειωθούν:
  1. Όσο μεγαλύτερο πλήθος σημείων έχουμε τόσο μικρότερο θα είναι το σφάλμα και η απόκλιση από την πραγματική τιμή του π.
  2. Ακόμη και αν επιχειρήσουμε εκτελέσεις με το ίδιο m κάθε φορά θα πάρουμε διαφορετικά αποτελέσματα αφού η όλη διαδικασία βασίζεται στην τύχη.
Στον πίνακα 1, λοιπόν, παρουσιάζονται τα αποτελέσματα για διάφορες εκτελέσεις του κώδικα. 
πίνακας 1, τα αποτελέσματα του κώδικα για διάφορα m.
Παρατηρούμε πως όσο αυξάνεται το m το σφάλμα μειώνεται. Τέλος ας δούμε πως φαίνονται ορισμένες κατανομές των σημείων για διάφορα m (σχήματα 2, 3 και 4).


σχήμα 2, m=2 π=3.28 error=4.406%

σχήμα 3, m=3 π=3.212 error=2.241%

σχήμα 4, m=4 π=3.16 error=0.586%

Δεν υπάρχουν σχόλια:

Δημοσίευση σχολίου