Τμήμα Μαθηματικών

Πληροφορική Ι

Εβδομάδα 10 · Αριθμητικός Υπολογισμός, Γραφήματα και Προσομοιώσεις

(Οδηγίες χρήσης του διαδραστικού κώδικα: δείτε την αρχική σελίδα.)

1 Αριθμητικός (Προσεγγιστικός) Υπολογισμός

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

Σε αυτή την εκτενή εβδομάδα θα δούμε πολλά κλασικά παραδείγματα προσεγγιστικού υπολογισμού και προσομοίωσης, θα μάθουμε να σχεδιάζουμε πραγματικά γραφήματα με το matplotlib, και θα κλείσουμε με το σύνολο Mandelbrot και μερικούς αλγορίθμους Monte Carlo.

1.1 Προσέγγιση του π με τη μέθοδο Monte Carlo

Θεωρούμε το τετράγωνο \([-1,1] \times [-1,1]\) και τον εγγεγραμμένο μοναδιαίο κύκλο \(x^2+y^2 \le 1\). Ο λόγος του εμβαδού του κύκλου προς το εμβαδόν του τετραγώνου είναι

\[ \frac{\pi \cdot 1^2}{2 \cdot 2} = \frac{\pi}{4}. \]

Αν “ρίξουμε” \(n\) τυχαία σημεία μέσα στο τετράγωνο και μετρήσουμε πόσα από αυτά πέφτουν μέσα στον κύκλο (hits), το ποσοστό hits/n προσεγγίζει το \(\pi/4\), άρα \(\pi \approx 4 \cdot \text{hits}/n\). Αυτή είναι μια προσομοίωση Monte Carlo: χρησιμοποιούμε τυχαιότητα για να προσεγγίσουμε μια ντετερμινιστική ποσότητα.

Τρέξτε τον κώδικα αρκετές φορές με το ίδιο n: θα πάρετε κάθε φορά μια ελαφρώς διαφορετική εκτίμηση, αφού τα σημεία είναι τυχαία. Αυξήστε το n σε 200000 και παρατηρήστε ότι το σφάλμα τείνει να μικραίνει — αλλά όχι μονότονα, και όχι πολύ γρήγορα (ο ρυθμός σύγκλισης μιας προσομοίωσης Monte Carlo είναι της τάξης \(1/\sqrt{n}\), πολύ πιο αργός από αιτιοκρατικές μεθόδους προσέγγισης του \(\pi\), όπως θα δούμε αμέσως παρακάτω).

1.2 Οπτικοποίηση με το matplotlib

Το matplotlib το γνωρίσαμε ήδη στην Εβδομάδα 1: plt.plot(x, y) για γράφημα γραμμής, plt.scatter για διάγραμμα διασποράς, plt.bar για ραβδόγραμμα, και οι xlabel / ylabel / title / legend για επεξηγήσεις.

Εδώ πάμε λίγο βαθύτερα — πολλαπλές σειρές στο ίδιο γράφημα, λογαριθμικές κλίμακες, και εικόνες (imshow) — τα εργαλεία που θα χρειαστούμε σε όλη την εβδομάδα:

Ας ξαναδούμε τη Monte Carlo προσομοίωση του \(\pi\) παραπάνω, αυτή τη φορά βλέποντας τα σημεία: μέσα στον κύκλο με ένα χρώμα, έξω με άλλο — ένα διάγραμμα διασποράς με δύο σειρές:

Από εδώ και πέρα θα χρησιμοποιούμε το matplotlib κάθε φορά που χρειαζόμαστε πραγματικό γράφημα.

Ας δούμε και τη σύγκλιση: το σφάλμα της εκτίμησης καθώς αυξάνεται το n, σε διπλά λογαριθμικό άξονα. Η μέθοδος Monte Carlo συγκλίνει σαν \(1/\sqrt{n}\) — δηλαδή για \(100\) φορές περισσότερα σημεία, μόλις \(10\) φορές μικρότερο σφάλμα:

1.3 Η προσέγγιση του π από τον Αρχιμήδη

Η ιδέα

Ο Αρχιμήδης (3ος αιώνας π.Χ.) προσέγγισε το \(\pi\) παγιδεύοντάς το ανάμεσα σε δύο τιμές, χωρίς να χρειαστεί καθόλου τριγωνομετρία:

  1. Θεωρούμε έναν κύκλο, και δύο κανονικά πολύγωνα με τον ίδιο αριθμό πλευρών: ένα εγγεγραμμένο (μέσα στον κύκλο, με κορυφές πάνω του) και ένα περιγεγραμμένο (γύρω από τον κύκλο, εφαπτόμενο σε αυτόν).
  2. Η περίμετρος του εγγεγραμμένου πολυγώνου είναι πάντα μικρότερη από την περιφέρεια του κύκλου, και η περίμετρος του περιγεγραμμένου πάντα μεγαλύτερη. Άρα ο λόγος (περίμετρος)/(διάμετρος) των δύο πολυγώνων παγιδεύει το \(\pi\) ανάμεσά του: \(b_n < \pi < a_n\).
  3. Όσο αυξάνουμε το πλήθος των πλευρών \(n\), και τα δύο πολύγωνα πλησιάζουν όλο και περισσότερο τον κύκλο, οπότε τα δύο φράγματα σφίγγουν γύρω από το \(\pi\).
  4. Ο Αρχιμήδης βρήκε έναν τρόπο να υπολογίζει τους λόγους του πολυγώνου με διπλάσιες πλευρές (\(2n\)), από τους λόγους του πολυγώνου με \(n\) πλευρές, χρησιμοποιώντας μόνο δύο είδη μέσων όρων:
    • το νέο άνω φράγμα \(a_{2n}\) είναι ο αρμονικός μέσος των \(a_n\) και \(b_n\)·
    • το νέο κάτω φράγμα \(b_{2n}\) είναι ο γεωμετρικός μέσος του νέου \(a_{2n}\) με το παλιό \(b_n\).

\[ a_{2n} = \frac{2\,a_n\,b_n}{a_n + b_n}, \qquad b_{2n} = \sqrt{a_{2n}\, b_n}. \]

  1. Ξεκινώντας από το κανονικό εξάγωνο (\(n=6\)), όπου οι λόγοι είναι γνωστοί με ακρίβεια (\(a_6 = 2\sqrt{3}\), \(b_6 = 3\)), και διπλασιάζοντας 4 φορές (καταλήγοντας σε πολύγωνο με \(96\) πλευρές), ο Αρχιμήδης απέδειξε ότι

\[ \frac{223}{71} < \pi < \frac{22}{7}. \]

Ο κώδικας

Δοκιμάστε να αυξήσετε το doublings (π.χ. σε 10 ή 20) και παρατηρήστε πόσο γρηγορότερα σφίγγει το εύρος αβεβαιότητας σε σχέση με τη μέθοδο Monte Carlo που είδαμε παραπάνω — η μέθοδος του Αρχιμήδη συγκλίνει εκθετικά γρήγορα, όχι μόνο σαν \(1/\sqrt{n}\). Ας το δούμε και σε γράφημα, με λογαριθμική κλίμακα στον κατακόρυφο άξονα:

Μια σχεδόν ευθεία γραμμή σε λογαριθμική κλίμακα σημαίνει ακριβώς εκθετική σύγκλιση.

1.4 Απώλεια ακρίβειας: float32 έναντι float64

Ο τύπος float της Python αναπαριστά πραγματικούς αριθμούς σε διπλή ακρίβεια (double precision, 64-bit, IEEE 754), με περίπου 15–17 σωστά δεκαδικά ψηφία. Οι περισσότερες γλώσσες προσφέρουν και έναν τύπο απλής ακρίβειας (single precision, 32-bit), με πολύ λιγότερα σωστά ψηφία (περίπου 7). Κάθε αριθμητική πράξη σε πεπερασμένη ακρίβεια στρογγυλοποιεί το αποτέλεσμα στην πλησιέστερη αναπαραστάσιμη τιμή — και αυτές οι μικροσκοπικές στρογγυλοποιήσεις μπορούν να συσσωρευτούν.

Σημείωση

Γιατί “περίπου” και όχι έναν ακριβή αριθμό ψηφίων; Η ακρίβεια αποθηκεύεται σε bits (δυαδικό σύστημα), όχι σε δεκαδικά ψηφία, και τα δύο δεν αντιστοιχίζονται 1-προς-1. Η double έχει \(53\) bits significand (mantissa), και κάθε bit “αξίζει” \(\log_{10}(2) \approx 0.301\) δεκαδικά ψηφία, άρα \(53 \times 0.301 \approx 15{,}95\) — γι’ αυτό λέμε “15 έως 17”: \(15\) δεκαδικά ψηφία είναι πάντα σωστά για κάθε double, ενώ \(17\) ψηφία χρειάζονται μόνο στη χειρότερη περίπτωση για να ξαναφτιάξουμε ακριβώς τον αρχικό αριθμό από τη δεκαδική του γραφή (round-trip). Ομοίως, η single έχει \(24\) bits significand, \(24 \times 0.301 \approx 7{,}22\), άρα “περίπου 7”.

Η Python δεν έχει ξεχωριστό τύπο float32 από προεπιλογή, αλλά μπορούμε να προσομοιώσουμε μία πράξη σε single precision στρογγυλοποιώντας το αποτέλεσμα κάθε βήματος στην ακρίβεια των 32-bit, χρησιμοποιώντας το module struct της πρότυπης βιβλιοθήκης (πακετάρουμε τον αριθμό σε 4 bytes IEEE-754 single precision, και τον ξεπακετάρουμε αμέσως πίσω σε κανονικό float):

Ας ξαναγράψουμε τον αλγόριθμο του Αρχιμήδη ώστε να στρογγυλοποιεί κάθε ενδιάμεσο αποτέλεσμα σε μια δεδομένη ακρίβεια, και ας συγκρίνουμε το σφάλμα της προσέγγισης όταν δουλεύουμε με float32 έναντι της κανονικής float64 — αυτή τη φορά με ένα γράφημα αντί για πίνακα:

Παρατηρήστε τη συμπεριφορά: αρχικά και τα δύο σφάλματα μικραίνουν γρήγορα, όπως θα περιμέναμε από τη θεωρία. Μετά όμως από κάποιο σημείο, το σφάλμα της float32 σταματάει να μικραίνει (ή και χειροτερεύει!), ενώ η float64 συνεχίζει να βελτιώνεται για πολλές ακόμα επαναλήψεις.

Ο λόγος: καθώς το πολύγωνο πλησιάζει τον κύκλο, οι τιμές a και b πλησιάζουν πολύ η μία την άλλη. Οι πράξεις a+b και a*b γίνονται τότε ευαίσθητες σε μικροσκοπικά σφάλματα στρογγύλευσης — το φαινόμενο λέγεται απώλεια σημαντικών ψηφίων (catastrophic cancellation).

Πέρα από ένα σημείο, η «θεωρητική» βελτίωση του αλγορίθμου χάνεται μέσα στον θόρυβο της πεπερασμένης ακρίβειας. Η float64, με πολύ περισσότερα ψηφία, φτάνει σε αυτό το όριο πολύ αργότερα — αλλά τελικά και αυτή σταματάει.

1.5 Η μέθοδος Newton-Raphson για εύρεση ριζών

Στην Εβδομάδα 7 είδαμε τη μέθοδο Newton ειδικά για τον υπολογισμό τετραγωνικών ριζών. Η ίδια ιδέα γενικεύεται άμεσα στην εύρεση ρίζας οποιασδήποτε παραγωγίσιμης συνάρτησης \(f\), δηλαδή μιας λύσης της εξίσωσης \(f(x) = 0\): ξεκινώντας από μια αρχική εκτίμηση \(x_0\), ακολουθούμε την εφαπτομένη της \(f\) στο τρέχον σημείο μέχρι να τμήσει τον άξονα \(x\), και επαναλαμβάνουμε:

\[ x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}. \]

Επειδή στην Python οι συναρτήσεις είναι απλά αντικείμενα, μπορούμε να γράψουμε μία γενική συνάρτηση newton που δέχεται ως ορίσματα την ίδια τη \(f\) και την παράγωγό της \(f'\):

Γεωμετρικά: από το σημείο \((x_k, f(x_k))\) κατεβαίνουμε στην εφαπτομένη μέχρι να τμήσει τον άξονα — αυτό το \(x\) είναι το επόμενο \(x_{k+1}\). Λίγα βήματα φτάνουν:

Η γενικότητα της μεθόδου φαίνεται καθαρά αν τη χρησιμοποιήσουμε για μια εξίσωση που δεν έχει καμία σχέση με ρίζες αριθμών — π.χ. την εξίσωση \(\cos(x) = x\), δηλαδή \(f(x) = \cos(x) - x = 0\):

Η λύση αυτή (\(x \approx 0.739085\)) λέγεται αριθμός Dottie, και είναι το μοναδικό σταθερό σημείο του συνημιτόνου. Δοκιμάστε τη συνάρτηση newton με τη δική σας \(f\) και \(f'\) — αρκεί η εξίσωση να έχει (πραγματική) ρίζα κοντά στην αρχική σας εκτίμηση \(x_0\).

Μια εναλλακτική: η μέθοδος της διχοτόμησης (bisection)

Η μέθοδος Newton συγκλίνει πολύ γρήγορα, αλλά έχει ένα μειονέκτημα: χρειάζεται την παράγωγο \(f'\), η οποία δεν είναι πάντα εύκολο (ή καν δυνατό) να υπολογιστεί. Η μέθοδος της διχοτόμησης δεν χρειάζεται καθόλου παράγωγο — στηρίζεται μόνο στο θεώρημα ενδιάμεσης τιμής: αν η \(f\) είναι συνεχής και τα πρόσημα των \(f(a)\), \(f(b)\) είναι αντίθετα, τότε υπάρχει σίγουρα μια ρίζα ανάμεσα στο \(a\) και το \(b\). Η ιδέα είναι ολόιδια με τη δυαδική αναζήτηση που είδαμε στην Εβδομάδα 1 και στην Εβδομάδα 4: εξετάζουμε το μέσο \(m=(a+b)/2\), και κρατάμε το μισό διάστημα όπου σίγουρα υπάρχει ρίζα, πετώντας το άλλο μισό:

Συγκρίνετε το πλήθος επαναλήψεων με αυτό της newton για το ίδιο πρόβλημα (x**3 - 2). Η διχοτόμηση χρειάζεται πολύ περισσότερα βήματα: σε κάθε επανάληψη απλώς «μισαίνει» το εύρος αβεβαιότητας (γραμμική σύγκλιση), ενώ η Newton διπλασιάζει τα σωστά ψηφία (τετραγωνική σύγκλιση).

Το αντιστάθμισμα: η διχοτόμηση είναι εγγυημένα ευσταθής (αρκεί η αρχική αλλαγή προσήμου), ενώ η Newton μπορεί να αποκλίνει για «κακή» αρχική εκτίμηση ή προβληματική παράγωγο.

Σημείωση

Μην μπερδεύετε αυτή τη μέθοδο διχοτόμησης με το module bisect της πρότυπης βιβλιοθήκης της Python — εκείνο αφορά αναζήτηση/εισαγωγή σε ήδη ταξινομημένες λίστες, όχι εύρεση ριζών συναρτήσεων. Η ιδέα της διχοτόμησης του χώρου αναζήτησης είναι κοινή, αλλά η εφαρμογή διαφορετική.

1.6 Σειρές Taylor

Η ιδέα

Η σειρά Taylor μιας συνάρτησης \(f\) γύρω από ένα σημείο \(a\) την προσεγγίζει με ένα πολυώνυμο, χτισμένο από τις παραγώγους της στο \(a\):

\[ f(x) \approx \sum_{k=0}^{n} \frac{f^{(k)}(a)}{k!}(x-a)^k. \]

Θα δουλέψουμε με σειρές Maclaurin (\(a=0\)). Η εκθετική συνάρτηση έχει την ιδιαίτερα απλή μορφή \(e^x = \sum_{k=0}^{\infty} x^k/k!\), αφού όλες οι παράγωγοί της είναι ίδιες με την ίδια τη συνάρτηση:

Το ημίτονο και το συνημίτονο, και οπτικοποίηση της σύγκλισης

\[ \sin(x) = \sum_{k=0}^{\infty} \frac{(-1)^k x^{2k+1}}{(2k+1)!}, \qquad \cos(x) = \sum_{k=0}^{\infty} \frac{(-1)^k x^{2k}}{(2k)!}. \]

Παρατηρήστε ότι κοντά στο μηδέν, ακόμα και λίγοι όροι δίνουν εξαιρετική προσέγγιση, ενώ μακριά από το μηδέν το πολυώνυμο “ξεφεύγει” γρήγορα — χρειάζονται πολύ περισσότεροι όροι για να ακολουθήσει το \(\sin(x)\) σε μεγαλύτερο εύρος. Οι σειρές Taylor προσεγγίζουν καλύτερα κοντά στο σημείο ανάπτυξης.

1.7 Αλγόριθμοι Monte Carlo, ξανά

Είδαμε ήδη πώς η τυχαιότητα προσεγγίζει το \(\pi\). Ας δούμε τρία ακόμα κλασικά παραδείγματα Monte Carlo.

Το πρόβλημα των γενεθλίων (Birthday Paradox)

Σε μια τάξη \(k\) ατόμων, ποια είναι η πιθανότητα τουλάχιστον δύο να έχουν τα ίδια γενέθλια (365 ισοπίθανες ημέρες); Υπάρχει ακριβής τύπος:

\[ P(\text{σύγκρουση}) = 1 - \prod_{i=0}^{k-1} \frac{365-i}{365}. \]

Το αποτέλεσμα είναι διάσημο ακριβώς επειδή είναι αντι-διαισθητικό: χρειάζονται μόλις \(23\) άτομα για πιθανότητα \(>50\%\).

Ολοκλήρωση Monte Carlo

Η ίδια ιδέα των τυχαίων σημείων γενικεύεται στον υπολογισμό οποιουδήποτε ορισμένου ολοκληρώματος: αν \(\bar{f}\) είναι η μέση τιμή της \(f\) σε τυχαία σημεία του \([a,b]\), τότε \(\int_a^b f(x)\,dx \approx (b-a)\cdot\bar{f}\):

Σε αντίθεση με τις προηγούμενες μεθόδους αυτής της εβδομάδας, η ολοκλήρωση Monte Carlo δεν χρειάζεται καμία ιδιότητα ομαλότητας της \(f\), και γενικεύεται εύκολα σε πολλές διαστάσεις — με το τίμημα της αργής σύγκλισης \(1/\sqrt{n}\).

Το πρόβλημα Monty Hall

Σε ένα τηλεπαιχνίδι υπάρχουν \(3\) πόρτες: πίσω από τη μία ένα αυτοκίνητο, πίσω από τις άλλες δύο τίποτα. Διαλέγετε μια πόρτα· ο παρουσιαστής, που ξέρει πού είναι το αυτοκίνητο, ανοίγει μία άλλη, άδεια πόρτα, και σας δίνει την επιλογή να αλλάξετε ή να μείνετε. Συμφέρει να αλλάξετε; Η διαίσθηση λέει “δεν έχει σημασία” — είναι λάθος:

Η αλλαγή πόρτας κερδίζει περίπου \(2/3\) των φορών, ενώ η παραμονή μόνο \(1/3\) — ο παρουσιαστής, ανοίγοντας επίτηδες μια άδεια πόρτα, “συμπυκνώνει” την πιθανότητα των δύο πορτών που δεν διαλέξατε σε μία.

1.8 Μιγαδικοί Αριθμοί και το Σύνολο Mandelbrot

Μιγαδικοί αριθμοί στην Python

Η Python έχει ενσωματωμένο τύπο complex: γράφουμε 3 + 4j (το j είναι η φανταστική μονάδα). Οι πράξεις +, -, *, / δουλεύουν κατευθείαν, και η abs() δίνει το μέτρο:

Ορισμός και αλγόριθμος χρόνου διαφυγής

Για κάθε μιγαδικό \(c\), ορίζουμε \(z_0=0,\ z_{n+1}=z_n^2+c\). Το σύνολο Mandelbrot είναι το σύνολο όλων των \(c\) για τα οποία η ακολουθία παραμένει φραγμένη. Μόλις \(|z_n|>2\), η ακολουθία είναι καταδικασμένη να αποκλίνει, άρα αρκεί να μετρήσουμε το χρόνο διαφυγής:

Οπτικοποίηση του συνόλου

Εφαρμόζουμε την mandelbrot_escape σε ένα πλέγμα σημείων, και ζωγραφίζουμε το αποτέλεσμα με plt.imshow:

Το χαρακτηριστικό σχήμα “καρδιάς με φυσαλίδες” είναι το σύνολο Mandelbrot· το σύνορό του είναι ένα φράκταλ — όσο κι αν μεγεθύνουμε, αποκαλύπτει νέα λεπτομέρεια. Δοκιμάστε να αλλάξετε τα όρια x_min, x_max, y_min, y_max για να “μεγεθύνετε” σε μια περιοχή (π.χ. γύρω από \(c=-0.75+0.1i\)), αυξάνοντας και το max_iter για περισσότερη λεπτομέρεια.

2 Σύνοψη

Σε αυτή την (εκτενή) εβδομάδα είδαμε:

  • Την τεχνική της προσομοίωσης Monte Carlo: τυχαία δειγματοληψία για αριθμητική εκτίμηση μιας ποσότητας.
  • Γραφήματα με matplotlib: plot, scatter, imshow, ετικέτες, υπόμνημα, λογαριθμική κλίμακα.
  • Επαναληπτικές μεθόδους σύγκλισης με άνω/κάτω φράγματα που “σφίγγουν” γύρω από μια τιμή.
  • Την αναπαράσταση αριθμών κινητής υποδιαστολής (float32 έναντι float64) και τα σφάλματα στρογγύλευσης που συσσωρεύονται.
  • Δύο αλγορίθμους εύρεσης ρίζας συνάρτησης: Newton-Raphson και διχοτόμηση.
  • Την τεχνική προσέγγισης μιας συνάρτησης με ένα πεπερασμένο άθροισμα όρων.
  • Την προσομοίωση Monte Carlo εφαρμοσμένη σε προβλήματα πιθανότητας, και τη σύγκρισή της με έναν ακριβή, κλειστό τύπο.
  • Έναν αλγόριθμο χρόνου διαφυγής (escape-time) πάνω σε πλέγμα σημείων, χρησιμοποιώντας μιγαδικούς αριθμούς ως τύπο δεδομένων.

⤺ Πίσω στην αρχική σελίδα / πρόγραμμα μαθήματος  ·  Διαφάνειες αυτής της εβδομάδας