9 Νοέμβρη 2002

1. 9Νοέμβρη 2002

2. Περιεχόμενα • Ένα Απλό Παράδειγμα: Εκθετική Αύξηση Πληθυσμού. • Ένα πιο Δύσκολο Παράδειγμα: Τροχιές Δορυφόρων. • Απλουστεύσεις του Μοντέλου • Αριθμητική Επίλυση του Μοντέλου • Αριθμητική Επίλυση στην Maple • Απογοήτευση! • Ποιο είναι το Πρόβλημα? • Μια Καλλίτερη Αριθμητική Λύση • Καλλίτερα Αποτελέσματα • Μια Ματιά στο Μέλλον

3. Εισαγωγή Οι υπολογισμοί παραγώγων συναρτήσεων είναι απλά ένα παράδειγμα επίλυσης διαφορικών εξισώσεων. Αναγκαζόμαστε να υπολογίσουμε τις παραγώγους προσεγγιστικά μια και πολλές φορές η απ΄ ευθείας χρήση του τύπου της παραγώγου είναι ιδιαίτερα χρονοβόρα. Πολλές φορές όμως αναγκαζόμαστε να κάνουμε το αντίθετο – ξεκινάμε με μια διαφορική εξίσωση, και θέλουμε να βρούμε την λύση της, ή τουλάχιστον να μπορέσουμε να την υπολογίσουμε αριθμητικά. Πολλά φυσικά φαινόμενα μπορούν να μοντελοποιηθούν με διαφορικές εξισώσεις: • Αλλαγές στον πληθυσμό των αρπακτικών και της λείας τους. • Αλλαγές στην οικονομία και οι επιπτώσεις τους στις επενδύσεις, την ανεργία ... • Προβλέψεις πολέμων τοπικής ή ευρείας κλίμακας. • Το μεγαλύτερο μέρος της Φυσικής (πχ τροχιά δορυφόρου).

4. Ένα Απλό Παράδειγμα: Εκθετική Αύξηση Πληθυσμού Η διακύμανση του πληθυσμού των οργανισμών σε ένα περιβάλλον με αρκετούς πόρους μπορεί να περιγραφθεί με την εξής διαφορική εξίσωση: dx(t) / dt = a x(t) Όπου, η x(t) παριστά τον πληθυσμό την χρονική στιγμή t. Η σταθερά a περιγράφει τον ρυθμό αύξησης του πληθυσμού. Η διαφορική αυτή εξίσωση μπορεί να λυθεί συμβολικά. Η Maple ξέρει πως, τουλάχιστον στην περίπτωση που δώσουμε μια συγκεκριμένη τιμή στο a (εδώ a = 3): > dsolve( { x(0) = 5, diff(x(t),t) = 3*x(t) }, { x(t) } ); x(t) = 5 exp(3 t) Όπου, η x(0) = 5 καθορίζει τον αρχικό πληθυσμό. Πόσο αποδεκτή είναι η λύση αυτή σαν μοντέλο αύξησης ρεαλιστικών πληθυσμών?

5. Παράδειγμα: Τροχιές Γύρω από Πλανήτες Οι άνθρωποι έχουν προσπαθήσει να μοντελοποιήσουν την κίνηση των πλανητών και των δορυφόρων του από την αρχαιότητα. Ο J. Kepler βρήκε ένα μοντέλο που περιγράφει καλά την κίνηση των πλανητών. Ο I. Newton εξήγησε πώς το μοντέλο αυτό μπορεί να προκύψει από τους (δικούς του) νόμους κίνησης και τον νόμο της βαρύτητας. Θα προσπαθήσουμε να χρησιμοποιήσουμε τους νόμους του Newton για να μοντελοποιήσουμε την κίνηση ενός αντικειμένου σε τροχιά γύρω από έναν πλανήτη. Υποθέτουμε ότι ξέρουμε : • Την μάζα του πλανήτη. • Την αρχική θέση του αντικειμένου. • Την αρχική ταχύτητα του αντικειμένου. Θέλουμε το μοντέλο μας να μπορεί να μας πει την τροχιά του αντικειμένου για κάποιο χρονικό διάστημα στο μέλλον.

6. Το Μαθηματικό Μοντέλο Στο μοντέλο μας, η ταχύτητα του αντικειμένου καθορίζει το πως αλλάζει η θέση του ως προς τον χρόνο. Σε καρτεσιανές συντεταγμένες: dpx/ dt = vx, dpυ/ dt = vυ Όπου, t είναι χρόνος, pxκαι pυείναι οι συντεταγμένες της θέσης και όπου vxκαι vυείναι οι συνιστώσες της ταχύτητας. Οι αλλαγές στην ταχύτητα εξαρτώνται από την επιτάχυνση: dvx/ dt = ax, dvυ/ dt = aυ Τέλος, η επιτάχυνση προκύπτει από την δύναμη της βαρύτητας του πλανήτη (ο οποίος ας υποθέσουμε είναι στην αρχή των αξόνων). Χρησιμοποιώντας F = ma και F = GMm / d2 , where d2 = d2χ + d2υ , και λαμβάνοντας υπ’ όψιν ότι η κατεύθυνση της δύναμης είναι προς τον πλανήτη παίρνουμε ότι ax= ( GM (px/d) / d2 ), ay= ( GM (py/d) / d2 ),

7. Απλουστεύσεις του Μοντέλου Είναι αυτό που παράγαμε ένα καλό μοντέλο? Με λίγη σκέψη θα καταλήξουμε ότι δεν είναι ένα άψογο μοντέλο του πώς πραγματικά αντικείμενα περιστρέφονται γύρω από την πραγματική Γη. • Αγνοεί όλα τα άλλα αντικείμενα του σύμπαντος τα οποία έχουν κάποια βαρυτική (ή άλλη) επίδραση. • Υποθέτει ότι ο πλανήτης δεν κινείται – μια ιδιαίτερα άσχημη υπόθεση εάν το αντικείμενο είναι συγκρίσιμης μάζας με την Γη. • Υποθέτει ότι οι νόμοι της βαρύτητας του Νεύτωνα είναι σωστοί. Ο Αϊνστάιν απέδειξε ότι είναι απλές προσεγγίσεις. • Υποθέτει ότι το αντικείμενο είναι άκαμπτο. Αντικείμενα που δυνατόν να παραμορφωθούν απορροφούν μικρές ποσότητες ενέργειας καθώς περιστρέφονται. Είναι μια καλή προσέγγιση? Εξαρτάται ...

8. Αριθμητική Επίλυση του Μοντέλου Παρόλο που ο Νεύτωνας έλυσε συμβολικά το σύστημα διαφορικών εξισώσεων του μοντέλου μας πριν από 300 χρόνια, κάτι τέτοιο είναι πολύ δύσκολο να το κάνει η Maple. Εάν είναι και για εσάς αρκετά δύσκολο, τότε ας προσπαθήσουμε να το λύσουμε αριθμητικά για κάποια συγκεκριμένη περίπτωση. Πρέπει να προκαθορίσουμε συγκεκριμένες τιμές για την αρχική θέση και την ταχύτητα του αντικειμένου, όπως και την μάζα του πλανήτη. Σαν αποτέλεσμα δεν θα πάρουμε έναν κλειστό τύπο – αλλά μια λίστα αριθμών που θα αφορούν τις x και y συντεταγμένεςτου αντικειμένου σε διάφορες χρονικές στιγμές στο μέλλον. Οι αριθμοί αυτοί θα είναι προσεγγίσεις της πραγματικής λύσης. Κάτι τέτοιο φαντάζει πολύ υποδεέστερο της συμβολικής λύσης. Το θετικό είναι ότι μπορούμε πάντοτε να προσπαθήσουμε να λύσουμε κάτι αριθμητικά – αν και μερικές φορές η προσπάθειά μας δεν θα είναι και πολύ επιτυχής.

9. Αριθμητική Επίλυση στην Maple numerical—orbit := proc (GM, start—xpos, start—ypos, start—xvel, start—yvel, step, n—steps) local xpos, ypos, xvel, yvel, xaccel, yaccel, distance, t, i; xpos := array(0..n—steps); ypos := array(0..n—steps); xvel := array(0..n—steps); yvel := array(0..n—steps); xpos[0] := evalf(start—xpos); ypos[0] := evalf(start—ypos); xvel[0] := evalf(start—xvel); yvel[0] := evalf(start—yvel); for t from 1 to n—steps do distance := sqrt(xpos[t­1]“2 + ypos[t­1]“2); xaccel := ­ xpos[t­1] * GM / distance“3; yaccel := ­ ypos[t­1] * GM / distance“3; xvel[t] := xvel[t­1] + step*xaccel; yvel[t] := yvel[t­1] + step*yaccel; xpos[t] := xpos[t­1] + step*xvel[t­1]; ypos[t] := ypos[t­1] + step*yvel[t­1]; od; [ [ xpos[i], ypos[i] ] $ i=0..n—steps ] end:

10. Απογοήτευση! Υποθέτοντας ότι αποθηκεύσαμε την παραπάνω διαδικασία στο αρχείο orbit.mp ας προσπαθήσουμε να την χρησιμοποιήσουμε: > read `orbit.mp`; > n1a := numerical—orbit(1,1,0,0,0.5,0.001,3000): > plot(n1a,scaling=constrained); Χάλια! Ας ξαναπροσπαθήσουμε με μικρότερο βήμα: > n1b:=numerical—orbit(1,1,0,0,0.5,0.0005,6000): > plot(n1b,scaling=constrained);

11. Ποιο είναι το Πρόβλημα? Προφανώς η αριθμητική λύση είναι πολύ ανακριβής. Σε αυτού του τύπου τα προβλήματα, η ανακρίβεια μπορεί να οφείλεται στην: • Ανεπάρκεια του μαθηματικού μοντέλου • Ανακριβή προσέγγιση του μοντέλου – εδώ ως προς το βήμα στον χρόνο. • Προσεγγιστική παράσταση αριθμών στην μηχανή. Το πιθανότερο εδώ φαίνεται να είναι το (2). Μπορούμε να ξεπεράσουμε το πρόβλημα αυτό χρησιμοποιώντας μικρότερα βήματα στον χρόνο – αλλά αυτό ίσως απαιτεί πάρα πολύ χρόνο. Θα μπορούσαμε εναλλακτικά να σκεφτούμε έναν καλλίτερο τρόπο προσέγγισης.

12. Μια Καλλίτερη Αριθμητική Λύση numerical—orbit—2 := proc (GM, start—xpos, start—ypos, start—xvel, start—yvel, step, n—steps) local xpos, ypos, xvel, yvel, xaccel, yaccel, distance, t, i; xpos := array(0..n—steps); ypos := array(0..n—steps); xvel := array(0..n—steps); yvel := array(0..n—steps); xpos[0] := evalf(start—xpos); ypos[0] := evalf(start—ypos); xvel[0] := evalf(start—xvel); yvel[0] := evalf(start—yvel); for t from 1 to n—steps do distance := sqrt(xpos[t­1]“2 + ypos[t­1]“2); xaccel := ­ xpos[t­1] * GM / distance“3; yaccel := ­ ypos[t­1] * GM / distance“3; xvel[t] := xvel[t­1] + step*xaccel; yvel[t] := yvel[t­1] + step*yaccel; xpos[t] := xpos[t­1] + step*xvel[t]; ypos[t] := ypos[t­1] + step*yvel[t]; od; [ [ xpos[i], ypos[i] ] $ i=0..n—steps ] end:

13. Καλλίτερα Αποτελέσματα Τα αποτελέσματα είναι καλλίτερα τώρα. > n2a := numerical—orbit(1,1,0,0,0.5,0.001,3000): > plot(n2a,scaling=constrained); Ενθαρρυντικό, ας προσπαθήσουμε μεγαλύτερα βήματα: > n2b:=numerical—orbit—2(1,1,0,0,0.5,0.03,100): > plot(n2b,scaling=constrained);

14. Μια Ματιά στο Μέλλον Τι θα συμβεί εάν προσπαθήσουμε να υπολογίσουμε τροχιές πιο μακριά στον χρόνο? > n2bb:=numerical—orbit—2(1,1,0,0,0.5,0.03,1000): > plot(n2bb,scaling=constrained); Όχι απόλυτα ακριβές – η σωστή λύση (υποθέτοντας την ορθότητα του μοντέλου), ο αντικείμενο πρέπει να ακολουθεί την ίδια έλλειψη ακριβώς. Από την άλλη μεριά, δεν είναι και τόσο άσχημα. Το αντικείμενο δεν κινείται σπειροειδώς προς τον (η απομακρυνόμενο από τον) πλανήτη.