Μέθοδος Όιλερ

Η Μέθοδος Όιλερ (Euler method), αποτελεί μία από τις βασικότερες μεθόδους επίλυσης των συνήθων διαφορικών εξισώσεων, όπου υφίσταται ως δεδομένη η αρχική τιμή. Στη πράξη για παράδειγμα με την μέθοδο αυτή μπορεί να γίνει η περιγραφή ενός πεδίου ροής

Συγκεκριμένα με την μέθοδο Όιλερ ένας παρατηρητής καταγράφει τις ιδιότητες του ρευστού σε μια στοιχειώδη περιοχή του χώρου κατά τον χρόνο που διέρχεται αυτού το ρευστό. Η συνάρτηση του πεδίου ροής είναι της μορφής u = u(x, y, z, t), όπου οι x, y, z, t , οι συντεταγμένες του χώρου και ο χρόνος, που τυγχάνουν όλες ανεξάρτητες μεταβλητές.

Το όνομα της μεθόδου αυτής δόθηκε προς τιμή του Ελβετού μαθηματικού και φυσικού Λεονάρδου Όιλερ (Leonard Euler 1707-1783).

Άτυπη γεωμετρική περιγραφή

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

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

Κάντε ένα μικρό βήμα κατά μήκος της εφαπτομένης γραμμής μέχρι το σημείο A_1. Αν υποθέσουμε ότι το A_1 είναι ακόμα στην καμπύλη, μπορεί να χρησιμοποιηθεί το ίδιο σκεπτικό όπως και για το σημείο A_0 παραπάνω. Μετά από αρκετά βήματα υπολογίζεται μία πολυγωνική καμπύλη A_0A_1A_2A_3\dots. Σε γενικές γραμμές, αυτή η καμπύλη δεν αποκλίνει πολύ από την αρχική άγνωστη καμπύλη, και το σφάλμα μεταξύ των δύο καμπυλών μπορεί να γίνει μικρό, αν το μέγεθος βήματος είναι αρκετά μικρό και το διάστημα του υπολογισμού είναι πεπερασμένο.
Παραγωγή
Απεικόνιση της αριθμητικής ολοκλήρωσης για την εξίσωση y'=y, y(0)=1. Το μπλε είναι η μέθοδος του Όιλερ, το πράσινο η μέθοδος ενδιάμεσης τιμής, το κόκκινο η ακριβής λύση y=e^t. Το μέγεθος του βήματος είναι h = 1.0.
Η ίδια απεικόνιση για h = 0.25. Φαίνεται ότι η μέθοδος ενδιάμεσης τιμής συγκλίνει γρηγορότερα απ’ ότι η μέθοδος του Όιλερ.

Θέλουμε να προσεγγίσουμε την λύση της αρχικής τιμής του προβλήματος

y'(t) = f(t,y(t)), \qquad \qquad y(t_0)=y_0,

χρησιμοποιώντας τους πρώτους δύο όρους από την επέκταση της σειράς Taylor του y, η οποία παρουσιάζει τη γραμμική προσέγγιση γύρω από το σημείο (t0,y(t0)) . Ένα βήμα της μεθόδου Όιλερ από το tn στο tn+1 = tn + h είναι

\( y_{n+1} = y_n + hf(t_n,y_n). \qquad \qquad \)

Η μέθοδος του Όιλερ είναι ρητή, π.χ. η λύση y_{n+1} είναι μία ρητή συνάρτηση του y_i για i \leq n.

Ενώ η μέθοδος του Όιλερ ενσωματώνει πρώτου βαθμού Συνήθεις Διαφορικές Εξισώσεις, κάθε Συνήθης Διαφορική Εξίσωση βαθμού N μπορεί να αναπαρασταθεί ως μία πρώτου βαθμού Συνήθης Διαφορική Εξίσωση:

\( y^{(N)}(t) = f(t, y(t), y'(t), \ldots, y^{(N-1)}(t)) , \)

Παρουσιάζουμε βοηθητικές μεταβλητές \( z_1(t)=y(t), z_2(t)=y'(t),\ldots, z_N(t)=y^{(N-1)}(t) \) και παίρνουμε την αντίστοιχη εξίσωση

\( \mathbf{z}'(t) = \begin{pmatrix} z_1'(t)\\ \vdots\\ z_{N-1}'(t)\\ z_N'(t) \end{pmatrix} = \begin{pmatrix} y'(t)\\ \vdots\\ y^{(N-1)}(t)\\ y^{(N)}(t) \end{pmatrix} = \begin{pmatrix} z_2(t)\\ \vdots\\ z_N(t)\\ f(t,z_1(t),\ldots,z_N(t)) \end{pmatrix} \)

Αυτό είναι ένα σύστημα πρώτου βαθμού στην μεταβλητή \( \mathbf{z}(t) \) και μπορεί να λυθεί εφαρμόζοντας την μέθοδο του Όιλερ ή κάθε άλλο σχήμα για πρώτου βαθμού συστήματα.
Παράδειγμα

Δίνεται η διαφορική εξίσωση y'=y και το αρχικό σημείο y(0)=1, θα θέλαμε να χρησιμοποιήσουμε τη μέθοδο του Όιλερ για να προσεγγίσουμε το y_3 χρησιμοποιώντας μέγεθος βήματος h=1.

Η μέθοδος Όιλερ είναι

\( y_{n+1} = y_n + hf(t_n,y_n). \qquad \qquad \)

οπότε πρώτα πρέπει να υπολογίσουμε f(t_0, y_0). Αυτή η απλή διαφορική εξίσωση εξαρτάται μόνο από το y.

\( f(y_0)=1. \qquad \qquad \)

Κάνοντας το παραπάνω βήμα βρίσκουμε την κλίση της εφαπτομένης της καμπύλης που είναι η λύση στο σημείο (0,1). Θυμηθείτε ότι η κλίση ορίζεται ως η μεταβολή του y διαιρούμενη από την μεταβολή του t, ή \( \operatorname{d}y/\operatorname{d}t. \)

Το επόμενο βήμα είναι να πολλαπλασιάσουμε την παραπάνω τιμή με το μέγεθος του βήματος h.

\( h \cdot f(y_0) = 1 \cdot 1 = 1. \qquad \qquad \)

Δεδομένου ότι το μέγεθος του βήματος είναι η μεταβολή του t, όταν πολλαπλασιάσουμε το μέγεθος του βήματος και την κλίση της εφαπτομένης, θα έχουμε μεταβολή του y value. Η τιμή αυτή, στη συνέχεια, προστίθεται στην αρχική τιμή του y και προκύπτει η επόμενη τιμή που θα χρησιμοποιηθεί για υπολογισμούς.

\( y_0 + hf(y_0) = y_1 = 1 + 1 \cdot 1 = 2. \qquad \qquad \)

Τα παραπάνω βήματα πρέπει να επαναληφθούν για να υπολογιστεί το \( y_2 \) και το \( y_3 \) .

\( y_2 = y_1 + hf(y_1) = 2 + 1 \cdot 2 = 4 \qquad \qquad \)
\( y_3 = y_2 + hf(y_2) = 4 + 1 \cdot 4 = 8 \qquad \qquad \)

Λόγω του επαναληπτικού χαρακτήρα αυτού του αλγορίθμου, μπορεί να είναι χρήσιμο να οργανώσουμε τους υπολογισμούς σε μορφή πίνακα, όπως φαίνεται παρακάτω, ώστε να αποφευχθούν τα λάθη.

\( y_n t_n y'(t) h dy y_{n+1} \)
1 0 1 1 1 2
2 1 2 1 2 4
4 2 4 1 4 8
Σφάλμα

Αν υποθέσουμε ότι οι f(t) και y(t) είναι γνωστές ακριβώς τη χρονική στιγμή t_0, τότε η μέθοδος του Όιλερ δίνει την προσεγγιστική λύση την χρονική στιγμή t_0+h:

\( y(t_0 + h) = y(t_0) + h f(t_0,y(t_0)) \, \qquad \qquad \)

Συγκριτικά, η επέκταση της σειράς Taylor στο h στο t_0 δίνει:

\( y(t_0 + h) = y(t_0) + h y'(t_0) + \frac{1}{2}h^2 y''(t_0) + O(h^3). \)

Εφόσον ξέρουμε ότι y'=f(t,y) έπεται ότι:

\( y''={\partial f\over\partial t}(t_0, y(t_0)) + {\partial f\over\partial y}(t_0, y(t_0)) f(t_0, y(t_0)) \)

Αυτό, μαζί με την f(t_0, y(t_0)) μπορούν να δοθούν ως είσοδος στην επέκταση της σειράς Taylor σε h στο t_0.

\( y(t_0 + h) = y(t_0) + h f(t_0, y(t_0)) + \frac{1}{2}h^2[{\partial f\over\partial t}(t_0, y(t_0)) + {\partial f\over\partial y}(t_0, y(t_0)) f(t_0, y(t_0))] + O(h^3). \)

Το σφάλμα που παρουσιάζεται από την μέθοδο του Όιλερ δίνεται από την διαφορά αυτών των εξισώσεων:

\( \frac{1}{2}h^2 [{\partial f\over\partial t}(t_0, y(t_0)) + {\partial f\over\partial y}(t_0, y(t_0)) f(t_0, y(t_0))] + O(h^3). \)

Για μικρό h, το κυρίαρχο λάθος ανά βήμα, ή τοπικό σφάλμα αποκοπής, είναι ανάλογο του h^2. Για την επίλυση του προβλήματος σε ένα δεδομένο σύνολο των t, ο αριθμός των βημάτων που απαιτούνται είναι ανάλογος με το 1/h οπότε πρέπει να περιμένουμε ότι το συνολικό σφάλμα στο τέλος του ορισμένου χρόνου, ή το συνολικό σφάλμα αποκοπής, θα είναι ανάλογο του h (σφάλμα ανά βήμα επί τον αριθμό των βημάτων). Επειδή το συνολικό σφάλμα αποκοπής είναι ανάλογο του h, η μέθοδος του Όιλερ λέγεται ότι είναι πρώτου βαθμού. Το γεγονός αυτό καθιστά τη μέθοδο του Όιλερ λιγότερο ακριβή (για μικρά h) σε σχέση με άλλες τεχνικές ανώτερου βαθμού όπως μέθοδοι Runge-Kutta και γραμμικές μέθοδοι πολλαπλών βημάτων.

Η μείωση του μεγέθους του βήματος μπορεί να βοηθήσει ώστε να κάνει την προσέγγιση πιο ακριβή και να μειωθεί το σφάλμα μεταξύ των δύο καμπυλών. Το σφάλμα με τρία δεκαδικά ψηφία για το παραπάνω παράδειγμα (με μέγεθος βήματος h=1) είναι το ακόλουθο:

\( \text{error} = |e^3 - 8| = 12.086 \qquad \qquad \)

Όταν το μέγεθος βήματος αλλάξει σε h=0.1 η τιμή του \( y_3 \)γίνεται 17.449. Το σφάλμα, τότε, για μέγεθος βήματος h=0.1 είναι το ακόλουθο:

\( \text{error} = |e^3 - 17.449| = 2.637 \qquad \qquad \)

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

Όπως και με άλλες μεθόδους, υπάρχει τρόπος για εμάς να αποφασίσουμε ένα φράγμα σφάλματος για ένα συγκεκριμένο πρόβλημα. Το φράγμα σφάλματος στο συνολικό σφάλμα δίνεται από [1]:

\( |\epsilon_{n+1}| \le \frac{hM}{2L}(e^{L(t-t_0)}-1) \qquad \qquad \_

όπου h είναι το μέγεθος του βήματος, M είναι το άνω όριο στη δεύτερη παράγωγο του y στο δεδομένο διάστημα (το οποίο πρέπει να υπολογιστεί), και L είναι η σταθερά Lipschitz .

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

Για παράδειγμα, ας υπολογίσουμε το μέγεθος βήματος που χρειάζεται ώστε η ολική αποκοπή σφάλματος να είναι \( \epsilon_{n+1}=0.01, \) υποθέτοντας ότι η μέγιστη τιμή της δεύτερης παραγώγου είναι M=10, η Lipschitz σταθερά είναι L=1, και t από το μηδέν στο τέσσερα. Χρησιμοποιώντας την συγκεκριμένη εξίσωση, παίρνουμε τα εξής:

\( \frac{2L}{M(e^{L(t-t_0)}-1)}|\epsilon_{n+1}|\le h \qquad \qquad \)

\( \frac{2}{10(e^{(4)}-1)}|0.01|\le h \qquad \qquad \)

\( h \le 0.0037314721 \qquad \qquad \)

Το οποίο σημαίνει ότι το h πρέπει να είναι μικρότερο από το παραπάνω ώστε να έχουμε το επιθυμητό σφάλμα ή μικρότερο , και 4/h, ή θα χρειαστούν περίπου 1072 επαναλήψεις για να γίνει αυτό. Ο τεράστιος αριθμός των βημάτων, και συνεπώς το υψηλό υπολογιστικό κόστος, ενισχύουν την χρήση εναλλακτικών, μεγαλύτερου βαθμού μεθόδων όπως η μέθοδος Runge–Kutta ή γραμμικές μέθοδοι πολλαπλών βημάτων.

Πηγές

Ascher, Uri M.; Petzold, Linda Ruth. Computer methods for ordinary differential equations and differential-algebraic equations. 1998. SIAM. ISBN 0898714125

↑ Simple Euler method and its modifcations

Εξωτερικοί Σύνδεσμοι
Sister project Τo Βικιβιβλίο Calculus έχει μια σελίδα σχετικά με
Euler's Method

Euler's Method for O.D.E.'s