Σημειωματάριο Πέμπτης 2 Απριλίου 2015

Σχεδιασμός Πασχαλινού αυγού

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

Σχεδιασμός ενός κύκλου

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

Θα ζωγραφίσουμε ένα κύκλο χρησιμοποιώντας την παραμετρική του μορφή στην οποία έχουμε ένα ιδεατό σημείο με συντεταγμένες \((x(t), y(t))\) (τη χρονική στιγμή \(t\)) και αυτό που ζωγραφίζουμε είναι ουσιαστικά η τροχιά του κινητού αυτού. Αν επιλέξουμε \[ x(t) = R \cos{t},\ \ \ y(t) = R \sin{t},\ \ \ \ (0 \le t \le 2\pi), \] τότε το σημείο μας διαγράφει μια πλήρη περιστροφή του κύκλου με ακτίνα \(R\) και κέντρο το \((0,0)\) όταν ο χρόνος κινείται από το 0 έως το \(2\pi\). Μάλιστα η αρχική (\(t=0\)) θέση του σημείου είναι η \((R,0)\), η θέση του σημείου στο χρόνο \(\pi/2\) είναι η θέση \((0,R)\), η θέση του σημείου για χρόνο \(t=\pi\) είναι η \((-R,0)\), κλπ.

Παρακάτω σχεδιάζουμε ένα κύκλο με κέντρο το \((0,0)\) και ακτίνα \(R=10\).

Ξεκινούμε με τη συνάρτηση linspace με την οποία δημιουργούμε ένα διάνυσμα t στο οποίο κρατάμε όλους τους (1000 το πλήθος) χρόνους από 0 έως \(2\pi\) για τους οποίους θα σχεδιάσουμε τη θέση του κινητού.

Έχοντας υπολογίσει τις χρονικές στιγμές που μας ενδιαφέρουν υπολογίζουμε τις αντίστοιχες \(x\) και \(y\) θέσεις (σε αυτές τις χρνικές στιγμές) με τις εντολές

x = R*np.cos(t)
y = R*np.sin(t)

Παρατηρείστε ότι πρόκειται για διανυσματικές πράξεις (τα t, x, y είναι numpy arrays) κατά σημείο.

Έπειτα σχεδιάζουμε τις θέσεις x, y με την εντολή plt.plot(x, y).

Πριν ολοκληρώσουμε την προεργασία και δείξουμε το γράφημα καθορίζουμε μέσω της κλήσης

plt.axes().set_aspect('equal')

ότι θέλουμε το aspect ratio των δύο αξόνων x, y να είναι 1. Με άλλα λόγια θέλουμε ίδιες αποστάσεις στους δύο άξονες να αντιπροσωπεύουν ίδια πραγματικά μήκη στην οθόνη (ώστε, π.χ., ένα τετράγωνο με πλευρές παράλληλες προς τους δύο άξονες να μοιάζει με τετράγωνο και όχι απλά με ορθογώνιο).

In [2]:

import matplotlib.pyplot as plt
import numpy as np

R = 10

t = np.linspace(0, 2*np.pi, 1000)

x = R*np.cos(t)
y = R*np.sin(t)

plt.plot(x, y)
plt.axes().set_aspect('equal')
plt.show()

Όμως το σχήμα του αυγού δεν είναι ακριβώς κυκλικό. Μια πρώτη προσέγγιση είναι να πούμε ότι το αυγό είναι όπως ο προηγούμενος κύκλος που σχεδιάσαμε αλλά με το πάνω ημικύκλιο <<τραβηγμένο>> προς τα πάνω. Ένας τρόπος να το πετύχουμε αυτό είναι να τροποποιήσουμε τις παραμετρικές εξισώσεις του κύκλου \[ x(t) = R \cos t,\ \ \ y(t) = R \sin t,\ \ \ \ (0 \le t \le 2\pi), \] ώστε η ακτίνα \(R\) να μην παραμένει σταθερή στο χρόνο αλλά να αυξάνει από \(t=0\) έως \(t=\pi/2\) (υψηλότερο σημείο), να μειώνεται με συμμετρικό τρόπο από \(t=\pi/2\) έως \(t=\pi\) και να παραμένει σταθερή και ίση με 10 όπως και στον κύκλο για \(\pi \le t \le 2\pi\).

Ας θέσουμε \(R_0 = 10\) για την ακτίνα του κάτω ημικυκλίου και \(R_1 = 13\) για την ακτίνα στο υψηλότερο σημείο του αυγού, και ας γράψουμε \(r(t)\) για την ακτίνα στη χρονική στιγμή \(t\). Η παραμετρική παράσταση του αυγού μας είναι \[ x(t) = r(t) \cos t,\ \ \ y(t) = r(t)\sin t,\ \ \ \ (0 \le t \le 2\pi). \] Η συνάρτηση \(r(t)\) είναι σταθερή και ίση με \(R_0\) για \(\pi \le t \le 2\pi\). Θέλουμε να αυξάνει με συνεχή τρόπο από \(t=0\) έως \(t=\pi/2\) ανάμεσα στις τιμές \(R_0\) και \(R_1\) και να φθίνει συμμετρικά για \(\pi/2 \le t \le \pi\). Ένας απλός τρόπος να το πετύχουμε αυτό είναι να ορίσουμε τη συνάρτηση ως τη λεγόμενη γραμμική παρεμβολή ανάμεσα στις δύο τιμές: \[ r(t) = \frac{2}{\pi}\left(R_1 t + R_0 (\frac{\pi}{2}-t) \right),\ \ \ \ (0 \le t \le \pi/2). \] Εύκολα επιβεβαιώνει κανείς ότι η συνάρτηση αυτή είναι αύξουσα, συνεχής και παίρνει τις δύο τιμές \(R_0\) και \(R_1\) για \(t=0\) και \(t=\pi/2\) αντίστοιχα.

Αντίστοιχα ορίζουμε τη συνάρτηση \(r(t)\) για \(\pi/2 \le t \le 2\pi\): \[ r(t) = \frac{2}{\pi}\left( R_0(t-\frac{\pi}{2}) + R_1 (\pi-t) \right),\ \ \ \ (\pi/2 \le t \le \pi). \] Στο διάστημα αυτό η συνάρτησή μας είναι φθίνουσα και συνεχής ανάμεσα στις δύο τιμές \(R_1\) και \(R_0\).

Για να υλοποιήσουμε τη συνάρτηση αυτή γράφουμε τη συνάρτηση r(t) η οποία είναι απλά η υλοποίηση αυτού του <<κλαδικού>> τύπου που μόλις περιγράψαμε.

Πέρα από το διάνυσμα των χρόνων,t, θα χρειαστούμε τώρα και το διάνυσμα R των ακτίνων (που περιέχει μια ακτίνα για κάθε μια από τις χρονικές στιγμές που εξετάζουμε). Στη θέση R[i] έχουμε την τιμή r(t[i]) αλλά πριν κάνουμε αυτή την ανάθεση θα πρέπει να δημιουργήσουμε ένα numpy array ίδιου σχήματος με το t. Δημιουργούμε με την εντολή

R = np.zeros(t.shape)

ένα διάνυσμα, ίδιου σχήματος (μήκους) με το t, γεμάτο με μηδενικά.

Οι επόμενες εντολές είναι όπως και πριν για το σχεδιασμό του κύκλου (αλλά δεν ξεχνάμε ότι τώρα το R είναι διάνυσμα και όχι σταθερή τιμή).

In [6]:

import matplotlib.pyplot as plt
import numpy as np

R0 = 10; R1 = 13

def r(t):
    if (0 <= t) and (t <= np.pi/2):
        return R1*t*2/np.pi + R0*2*(np.pi/2-t)/np.pi
    if (np.pi/2 <= t) and (t <= np.pi):
        return R1*2*(np.pi-t)/np.pi + R0*2*(t-np.pi/2)/np.pi
    else:
        return R0

t = np.linspace(0, 2*np.pi, 1000)
R = np.zeros(t.shape)

for i in range(len(R)):
    R[i] = r(t[i])

x = R*np.cos(t)
y = R*np.sin(t)

plt.plot(x, y)
plt.axes().set_aspect('equal')
plt.show()

Ένα σχεδόν οφθαλμοφανές μειονέκτημα του προηγούμενου σχήματος είναι η έλλειψη ομαλότητας, ή, με άλλα λόγια, το γεγονός ότι η καμπύλη πουτ σχεδιάσαμε δε φαίνεται (με το μάτι) να έχει παράγωγο (εφαπτόμενη) στα σημεία που αντιστοιχούν στους χρόνους \(t=0, \pi/2, \pi\). Εύκολα επιβεβαιώνουμε κοιτώντας τον τύπο για τη συνάρτηση \(r(t)\) ότι όντως η συνάρτηση αυτή δεν είναι παραγωγίσιμη για αυτές τις τιμές του \(t\), πράγμα όχι τόσο απρόσμενο μια και στα σημεία αυτά η συνάρτηση δίνεται από διαφορετικό τύπο αριστερά και δεξιά του σημείου.

Μια λύση σε αυτό το πρόβλημα είναι να βρούμε μια συνάρτηση \(r(t)\) που να μην έχει τέτοιο πρόβλημα έλλειψης παραγώγου και ταυτόχρονα να έχει τα χαρακτηριστικά μονοτονίας και σταθερότητας που ζητάμε.

Μια λύση είναι η συνάρτηση \[ r(t) = \begin{cases} R_0 + (R_1-R_0) \sin^4 t & 0 \le t \le \pi\\ R_0 & \pi \le t \le 2\pi \end{cases} \] (ο εκθέτης 4 θα μπορούσε να είναι και 2 εδώ, απλά το 4 παράγει αισθητικά καλύτερα αποτελέσματα). Η συνάρτηση αυτή, παρά το ότι δίνεται κι αυτή από κλαδικό τύπο, έχει παραγώγους και στα σημεία \(t=0\) και \(t=\pi\) (εκεί όπου αλλάζει ο κλαδικός τύπος) και αυτό οφείλεται στην ύψωση σε δύναμη. Εύκολα αποδεικνύονται όλα αυτά καθώς και το ότι η συνάρτηση αυτή είναι αύξουσα στο \([0,\pi/2]\), φθίνουσα στο \([\pi/2, \pi]\) και σταθερή στο \([\pi, 2\pi]\).

Όλα τα υπόλοιπα γίνονται με τον ίδιο τρόπο όπως και πριν με την εξαίρεση του χρώματος σχεδιασμού που τώρα είναι το κόκκινο. Αυξήσαμε επίσης (και πάλι για αισθητικούς λόγους) το \(R_1\) σε 15 από 13.

In [8]:

import matplotlib.pyplot as plt
import numpy as np

R0 = 10; R1 = 15

def r(t):
    if 0 <= t <= np.pi:
        return R0 + (R1-R0)*(np.sin(t)**4)
    else:
        return R0

t = np.linspace(0, 2*np.pi, 1000)
R = np.zeros(t.shape)

for i in range(len(R)):
    R[i] = r(t[i])

x = R*np.cos(t)
y = R*np.sin(t)

plt.plot(x, y, 'r')
plt.axes().set_aspect('equal')
plt.show()

Στο επόμενο θα ζωγραφίσουμε κι ένα πράσινο <<φυλλαράκι>> πάνω στο κόκκινο αυγό (δείτε το σχήμα παρακάτω). Το φυλλαράκι αποτελείται από δύο συμμετρικά κυκλικά τόξα.

Οι δύο παράμετροι που ρυθμίζουν το σχήμα αυτών των κυκλικών τόξων είναι το \(R_2\) (τα κέντρα των αντίστοιχων κύκλων είναι στις θέσεις \((-R_2,0)\) και \((R_2,0)\)) και το \(\delta\) (οι ακτίνες των δύο κύκλων είναι \(R_2+\delta\)). Έτσι τα σημεία τομής αυτών των κυκλικών τόξων με τον άξονα των \(x\) είναι τα \((-\delta,0)\) και \((\delta,0)\). Με λίγη προσεκτική γεωμετρία βρίσκουμε ότι τα δύο αυτά κυκλικά τόξα αντιστοιχούν σε γωνία \(2\phi\) με \[ \cos\phi = \frac{R_2}{R_2+\delta}, \] και άρα μπορούμε να υπολογίσουμε το \(\phi\) χρησιμοποιώντας την αντίστροφη τριγωνομετρική συνάρτηση \(\cos^{-1}\) που στην python γράφεται ως acos. Αυτό και κάνουμε αφού δώσουμε τις τιμές \(R_2=20\) και \(\delta=1\) (μετά από διάφορες δοκιμές για να δούμε ποιο αποτέλεσμα μας αρέσει περισσότερο).

Επίσης ξαναγράφουμε τη συνάρτηση r(t) με τρόπο ώστε να μπορέσουμε να αποφύγουμε το for loop που είχαμε μέχρι τώρα (το οποίο κοστίζει σε χρόνο αρκετά). Γραμμένη με τον τρόπο αυτό η συνάρτηση r(t) μπορεί να εφαρμοστεί κατ' ευθείαν σε ένα διάνυσμα t χωρίς να χρειαστεί να γράψουμε loop. Για να γίνει όμως αυτό θα πρέπει η συνάρτηση να είναι γραμμένη κατάλληλα και η δυσκολία εδώ είναι ότι άλλος τύπος εφαρμόζεται για άλλες τιμές του t λόγω του κλαδικού χαρακτήρα της συνάρτησης.

Στη νέα μορφή της συνάρτησης r(t) ορίζουμε το βοηθητικό διάνυσμα R, ίδιας διάστασης με το t, το οποίο έχει ως αρχική τιμή το R0 σε όλες τις θέσεις. Με την εντολή

R[t<np.pi] += ((R1-R0)*(np.sin(t)**4))[t<np.pi]

λέμε ότι στις θέσεις i όπου t[i]<np.pi θα πρέπει να γίνει η ανάθεση της τιμής (R1-R0)*(np.sin(t[i])**4) στη θέση R[i]. Αυτή είναι ακριβώς η τροποποίηση που χρειάζεται να γίνει στις τιμές του R ώστε να <<επιμηκυνθεί>> το αυγό. Επιστρέφουμε το διάνυσμα R.

Ζωγραφίζουμε έπειτα το αυγό όπως και πριν και μετά ζωγραφίζουμε τα πράσινα κυκλικά τόξα. Γι' αυτό παίρνουμε ένα νέο διάνυσμα χρόνων t, αυτή το φορά από \(-\phi\) έως \(\phi\) (έχουμε δηλ. και εδώ όπως και πριν μια ταύτιση της έννοιας του χρόνου με την έννοια της γωνίας). Τα x, y τώρα ορίζονται με τις παραμετρικές εξισώσεις τωμν κύκλων που περιέχουν αυτά τα κυκλικά τόξα που είναι οι κύκλοι ακτίνας \(R_2+\delta\) με κέντρο το σημείο \((-R_2, 0)\) και το σημείο \((R_2, 0)\). Αυτά τα τόξα ζωγραφίζονται πράσινα. Μάλιστα εκμεταλλευόμαστε τη συμμετρία και αφού ζωγραφίσουμε το ένα από αυτά παίρνουμε τις συντεταγμένες του άλλου αλλάζοντας απλώς τις τιμές του x σε -x.