5.1 Άμεση μέθοδος του Euler

Όπως και στην προηγούμενη παράγραφο θα θεωρήσουμε ορισμένα σημεία στο χωρίο $[0,L]\times [0,T]$ και θα κατασκευάσουμε προσεγγίσεις της ακριβούς λύσης του (5.1) σε αυτά τα σημεία. Θεωρούμε λοιπόν έναν φυσικό αριθμό $N$ και τη διαμέριση του διαστήματος $[0,L]$ από $N+2$ ισαπέχοντα σημεία , όπου $h=x_{i+1}-x_i$, $i=0,\mathrm{\dots },N$. Επίσης, θεωρούμε και έναν φυσικό αριθμό $M$ και τη διαμέριση του $[0,T]$ από $M+1$ ισαπέχοντα σημεία , όπου $k=t^{j+1}-t^j$, $j=0,\mathrm{\dots },M-1$, βλ. Σχήμα 5.1. Τότε σε κάθε σημείο $(x_i,t^j)$ του διαμερισμού του $[0,L]\times [0,T]$, θα ισχύει:

$$u_t(x_i,t^j)=u_{xx}(x_i,t^j),i=1,\mathrm{\dots },N,j=0,\mathrm{\dots },M.$$

(5.3)

Θα θεωρήσουμε λοιπόν προσεγγίσεις των τιμών της $u$ στα παραπάνω σημεία $(x_i,t^j)$ τις οποίες θα συμβολίσουμε με $U_i^j$, $i=0,\mathrm{\dots },N+1$, $j=0,\mathrm{\dots },M$. Λόγω των συνοριακών συνθηκών $u(0,t^j)=u(L,t^j)=0$, $j=0,\mathrm{\dots },M$, θέτουμε $U_0^j=U_{N+1}^j=0$, $j=0,\mathrm{\dots },M$, και επειδή η λύση $u$ της (5.1) είναι γνωστή για $t=0$, θέτουμε $U_i^0=g(x_i)$, $i=1,\mathrm{\dots },N$. Οι υπόλοιπες τιμες $u(x_i,t^j)$, $i=1,\mathrm{\dots },N$, $j=1,\mathrm{\dots },M$, είναι άγνωστες και θα τις προσεγγίσουμε.

Στο Σχήμα 5.1, για ένα παράδειγμα πλέγματος σημείων, σημειώνουμε με άσπρο κύκλο τα σημεία που αναζητούμε στις τιμές $u(x_i,t^j)$ και με έντονο μαύρο τα σημεία που αντιστοιχούν στις γνωστές συνοριακές τιμές της $u$. Τις άγνωστες τιμές $u(x_i,t^j)$ τις προσεγγίζουμε με τις τιμές $U_i^j$, $i=1,\mathrm{\dots },N$, $j=1,\mathrm{\dots },M$, οι οποίες προκύπτουν με τον ακόλουθο τρόπο.

Για να προσεγγίσουμε τις $u_t(x_i,t^j)$, $i=1,\mathrm{\dots },N$, $j=0,\mathrm{\dots },M-1$, στην (5.3) χρησιμοποιούμε τη ${\delta }_k^{+}u(x_i,t^j)$ που θεωρήσαμε στην (2.1). Έτσι, αν υποθέσουμε ότι $u_{tt}\in C([0,L]\times [0,T])$, λόγω της (2.3), η (5.3) γίνεται για $i=1,\mathrm{\dots },N$,

$$\frac{u(x_i,t^{j+1})-u(x_i,t^j)}{k}=u_{xx}(x_i,t^j)+{\widehat{\eta }}_i^j,j=0,\mathrm{\dots },M-1,$$

(5.4)

όπου

$$|{\widehat{\eta }}_i^j|\le \frac{k}{2}\underset{t\in [0,T]}{\max }|u_{tt}(x_i,t)|.$$

(5.5)

Στη συνέχεια, για να προσεγγίσουμε τις $u_{xx}(x_i,t^j)$, $i=1,\mathrm{\dots },N$, $j=0,\mathrm{\dots },M-1$, στην (5.4) χρησιμοποιούμε την προσέγγιση ${\delta }_{h,2}^{c}u(x_i,t^j)$ που θεωρήσαμε στην (2.8). Έτσι, αν υποθέσουμε ότι $\frac{{\partial }^4}{\partial x^4}u\in C([0,L]\times [0,T])$, λόγω της (2.9), η (5.4) γίνεται για $i=1,\mathrm{\dots },N,$ και $j=0,\mathrm{\dots },M-1,$

(5.6)

και

$$|{\tilde{\eta }}_i^j|\le \frac{h^2}{12}\underset{x\in [0,L]}{\max }|\frac{{\partial }^4}{\partial x^4}u(x,t^j)|.$$

(5.7)

Αν θέσουμε τώρα ${\eta }_i^j={\widehat{\eta }}_i^j+{\tilde{\eta }}_i^j$, η (5.6) δίνει για $i=1,\mathrm{\dots },N,$ και $j=0,\mathrm{\dots },M-1,$

$${\delta }_h^{+}u(x_i,t^j)-{\delta }_{h,2}^{c}u(x_i,t^j)={\eta }_i^j,$$

(5.8)

και λόγω των (5.5) και (5.7) εύκολα βλέπουμε ότι ισχύει το ακόλουθο λήμμα.

Θα κατασκευάσουμε λοιπόν προσεγγίσεις $U_i^j$ της λύσης $u$ του προβλήματος (5.1) στα σημεία $(x_i,t^j)$, $i=0,\mathrm{\dots },N+1$, $j=0,\mathrm{\dots },M$, θεωρώντας τις ακόλουθες εξισώσεις,

		${\displaystyle \frac{U_i^{j+1}-U_i^j}{k}}={\displaystyle \frac{U_{i+1}^j-2U_i^j+U_{i-1}^j}{h^2}},i=1,\mathrm{\dots },N,j=0,\mathrm{\dots },M-1,$		(5.10)
				(5.11)

Κατ’ αναλογία με το Κεφάλαιο 3, το σφάλμα ${\eta }_i^j$ που δίνεται στις (5.8) καλείται τοπικό σφάλμα διακριτοποίησης για τη μέθοδο (5.10)–(5.11).

Έστω τώρα ότι έχουμε θεωρήσει έναν διαμερισμό σημείων $\{x_i\}$ του $[0,L]$, με βήμα $h$ και $\{t^j\}$ του $[0,T]$, με βήμα $k$, τότε προκύπτει ένας διαμερισμός σημείων $\{(x_i,t^j)\}$ του $[0,L]\times [0,T]$. Αν το τοπικό σφάλμα διακριτοποίησης ${\eta }_i^j$, ενός αριθμητικού σχήματος σε ένα σημείο ενός διαμερισμού $\{(x_i,t^j)\}$, φράσσεται κατά απόλυτη τιμή από το γινόμενο μιας θετικής σταθεράς, που δεν εξαρτάται από τα $h$ και $k$, επί $h^{\kappa }$ και $k^{\mu }$, τότε λέμε ότι έχει τάξη ακρίβειας $\kappa$ ως προς $h$ και $\mu$ ως προς $k$.

Επομένως, λόγω των (5.5) και (5.7), το τοπικό σφάλμα διακριτοποίησης της μεθόδου (5.10)-(5.11) έχει τάξη ακρίβειας δύο ως προς $h$ και ένα ως προς $k$.

Επίσης, κατ’αναλογία του Oρισμού 3.1, αν το τοπικό σφάλμα διακριτοποίησης μιας μεθόδου, ${\eta }_i^j$, τείνει στο μηδέν, καθώς τα βήματα $h$ και $k$ του διαμερισμού τείνουν στο μηδέν, τότε η μέθοδος λέγεται συνεπής.

Συνεπώς σύμφωνα με το Λήμμα 5.1, το τοπικό σφάλμα διακριτοποίησης τείνει στο μηδέν, καθώς τα $h$ και $k$ τείνουν στο μηδέν, οπότε η (5.10)-(5.11) είναι μια συνεπής μέθοδος.

Αν συμβολίσουμε τώρα με $\lambda$ τον λόγο ${\displaystyle \frac{k}{h^2}}$, τότε η (5.10) γράφεται ως

$$U_i^{j+1}=\lambda U_{i+1}^j+(1-2\lambda )U_i^j+\lambda U_{i-1}^j,$$

(5.12)

για $i=1,\mathrm{\dots },N$, $j=0,\mathrm{\dots },M-1$. Ξεκινώντας από το χρονικό επίπεδο $t^0=0$, χρησιμοποιώντας δηλαδή τις αρχικές συνθήκες (5.11), από την (5.12) μπορούμε να υπολογίσουμε άμεσα την προσέγγιση στο χρονικό επίπεδο $t^1$. Συνεχίζοντας με αυτό τον τρόπο, λαμβάνουμε, με άμεσο τρόπο μέσω της (5.12), την προσέγγιση στο χρονικό επίπεδο $t^{j+1}$, από τις ήδη γνωστές τιμές στο χρονικό επίπεδο $t^j$.

Η μέθοδος (5.10)-(5.11) καλείται άμεση μέθοδος του Euler. Στο Σχήμα 5.2 φαίνεται ένα παράδειγμα πλέγματος, κοντά στο $(x_i,t^j)$, όπου σημειώνουμε με άσπρο κύκλο το σημείο του χρονικού επιπέδου $t^{j+1}$, όπου σύμφωνα με την (5.12), αντιστοιχεί σε άγνωστη τιμή της προσέγγισης και με μαύρο αυτά του χρονικού επιπέδου $t^j$ που αντιστοιχούν σε γνωστές τιμές της προσέγγισης.

Αν συμβολίσουμε με $U^j\in {ℝ}^N$ το διάνυσμα με συνιστώσες, $U_1^j,\mathrm{\dots },U_N^j$, $U^j=$ $(U_1^j,$ $\mathrm{\dots }$, $U_N^j){}^T$, μπορούμε να γράψουμε το σύστημα των εξισώσεων (5.12) ισοδύναμα ως

(5.13)

όπου $A$ είναι ο $N\times N$ πίνακας

και $G={(g(x_1),\mathrm{\dots },g(x_N))}^T$.

Επομένως, για να υπολογίσουμε την προσέγγιση στο χρονικό επίπεδο $t^{j+1}$, απλώς πολλαπλασιάζουμε τον πίνακα $A$ με τη γνωστή, από το προηγούμενο βήμα, $U^j$. Επειδή για να υπολογίσουμε την προσέγγιση σε κάθε χρονικό βήμα δεν χρειάζεται να αντιστρέψουμε κάποιον πίνακα, δηλαδή να λύσουμε κάποιο γραμμικό σύστημα με αριθμητικές μεθόδους, η μέθοδος καλείται άμεση.

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

Στη συνέχεια, θα δείξουμε το ακόλουθο θεώρημα για την ευστάθεια της άμεσης μεθόδου του Euler.

Για να αποδείξουμε την ευστάθεια της μεθόδου στο Θεώρημα 5.1 υποθέσαμε ότι $\lambda \le 1/2$. Λόγω της ύπαρξης μιας συνθήκης ανάμεσα στις παραμέτρους διακριτοποίησης $k$ και $h$, λέμε ότι η μέθοδος είναι ευσταθής υπό συνθήκες. Το γεγονός ότι αυτή η συνθήκη είναι αναγκαία και δεν οφείλεται σε τυχόν αδυναμία της απόδειξης του Θεώρηματος 5.1 που δώσαμε φαίνεται από το ακόλουθο παράδειγμα.

Είναι απλό να δούμε ότι αν $L=1$ και $g(x)=\sin (n\pi x)$, $n=1,2,\mathrm{\dots },$ τότε το πρόβλημα (5.1) γίνεται

	$u_t(x,t)$	$=u_{xx}(x,t),$		$x\in [0,1],t\in [0,T],$		(5.14)
	$u(0,t)$	$=u(1,t)=0,$		$t\in [0,T],$		(5.15)
	$u(x,0)$	$=\sin (n\pi x),$		$x\in [0,1],$		(5.16)

και η ακριβής λύση $u$ είναι η $u(x,t)=e^{-n^2{\pi }^{2}t}\sin (n\pi x)$, μια δηλαδή από τις θεμελιώδεις συναρτήσεις που εμφανίζονται στην έκφραση (5.2). Παρατηρούμε λοιπόν ότι σε κάθε χρονικό επίπεδο $t\ge 0$ καθεμιά από αυτές τις θεμελιώδεις συναρτήσεις είναι μια ημιτονοειδής συνάρτηση, της οποίας το πλάτος μειώνεται συνεχώς με εκθετικό τρόπο, καθώς το χρονικό επίπεδο αυξάνει. Είναι προφανές ότι, αν θεωρήσουμε για αρχική συνθήκη $g$ μια από τις $\sin (n\pi x)$, $n=1,2,\mathrm{\dots },$ και η προσεγγιστική λύση που λαμβάνουμε με τη μέθοδο (5.10), δεν αποτυπώνει αυτήν την ιδιότητα, ή τουλάχιστον το πλάτος δεν παραμένει φραγμένο, δεν θα προσεγγίζει “καλά” την αντίστοιχη θεμελιώδη συνάρτηση.

Επειδή η ακριβής λύση του (5.14)–(5.16) αποτελείται από το γινόμενο της $e^{-n^2{\pi }^{2}t}$, που εξαρτάται μόνο από τον χρόνο $t$, και της $\sin (n\pi x)$, υποθέτουμε ότι η λύση $U_i^j$ της (5.10) έχει τη μορφή $U_i^j=w_j\sin (n\pi x_i)$. Δηλαδή αποτελείται από δύο παράγοντες, ο ένας από τους οποίους εξαρτάται μόνο από τον χρόνο $t^j$ και ο άλλος είναι ακριβώς η ημιτονοειδής συνάρτηση που βρίσκεται στην έκφραση της ακριβούς λύσης του (5.14). Αντικαθιστούμε λοιπόν την παραπάνω μορφή της $U_i^j$ στην εξίσωση (5.12), οπότε παίρνουμε για $i=1,\mathrm{\dots },N,j=0,\mathrm{\dots },M-1$,

(5.17)

Στη συνέχεια, εφαρμόζοντας την τριγωνομετρική ταυτότητα

(5.18)

η (5.17) γίνεται

(5.19)

Επίσης, λόγω της τριγωνομετρικής ιδιότητας

$$\begin{array}{c}\hfill 2{\sin }^2(\frac{\theta }{2})=1-\cos (\theta ),\theta \in [0,2\pi ]\end{array}$$

(5.20)

η (5.19) δίνει

$$w_{j+1}=(1-4\lambda {\sin }^2(\frac{n\pi h}{2}))w_j.$$

(5.21)

Συνεπώς, $w_j={\kappa }^j{w}_0$, με $\kappa =1-4\lambda {\sin }^2(\frac{n\pi h}{2}).$ Άρα, η μέγιστη τιμή της $U_i^j$ είναι φραγμένη, αν και μόνο αν $|\kappa |\le 1$. Θα πρέπει, λοιπόν,

$$-1\le 1-4\lambda {\sin }^2(\frac{n\pi h}{2}))\le 1\text{ ή }2\lambda {\sin }^2(\frac{n\pi h}{2}))\le 1.$$

Επομένως, επειδή $|\sin (\frac{n\pi h}{2})|\le 1$, για $n=1,2,\mathrm{\dots }$, αρκεί $\lambda \le {\displaystyle \frac{1}{2}}$. Βλέπουμε ότι η συνθήκη ευστάθειας στο Θεώρημα 5.1 είναι αναγκαία για να έχουμε “καλές” προσεγγίσεις με την άμεση μέθοδο του Euler.

Μπορούμε να γενικεύσουμε την παραπάνω απόδειξη θεωρώντας τώρα ότι

$$U_i^j=w_j{e}^{r{x}_{i}I},r\in ℝ,$$

με $I=\sqrt{-1}$ τη φανταστική μονάδα, ικανοποιούν την εξίσωση (5.12). Αν οι τιμές $|U_i^j|$ παραμένουν φραγμένες για κάθε $r\in ℝ$, θα καλούμε αυτή την ιδιότητα ευστάθεια von Neumann. Είναι φανερό ότι για να είναι φραγμένες οι $|U_i^j|$, αρκεί να φράσσονται οι $|w_j|$ για κάθε $j$. Μια πιο διεξοδική μελέτη της ευστάθειας von Neumann ο ενδιαφερόμενος αναγνώστης μπορεί να βρεί στο (Strikwerda, (2004)).

Μπορεί να αποδειχθεί ότι για την εξίσωση της θερμότητας η ευστάθεια von Neumann συνεπάγεται την ευστάθεια μιας μεθόδου, βλ. π.χ. (Morton and Mayers, (2005)). Στη συνέχεια, θα δείξουμε ότι για να είναι η άμεση μέθοδος του Euler von Neumann ευσταθής πρέπει και πάλι να ισχύει $\lambda \le {\displaystyle \frac{1}{2}}$.

Υποθέτουμε λοιπόν ότι οι

$$U_i^j=w_j{e}^{r{x}_{i}I},r\in ℝ,$$

με $I=\sqrt{-1}$ τη φανταστική μονάδα, ικανοποιούν την εξίσωση (5.12) και, άρα,

(5.22)

Λόγω τώρα της ταυτότητας

$$e^{\theta I}=\cos (\theta )+I\sin (\theta ),$$

έχουμε

$$\cos (\theta )=\frac{1}{2}(e^{\theta I}+e^{-\theta I}),\theta \in [0,2\pi ].$$

(5.23)

Οπότε η (5.22) λόγω της (5.20) γίνεται

	$w_{j+1}$	$=(1-2\lambda +2\cos (rh))w_j{e}^{r{x}_{i}I}$
		$=(1-4\lambda {\sin }^2({\displaystyle \frac{rh}{2}}))w_j.$

Συνεπώς $w_j={\kappa }^j{w}_0$, με $\kappa =1-4\lambda {\sin }^2(\frac{rh}{2})$. Για να παραμένει η λύση $U_i^j$ φραγμένη κατ’απόλυτη τιμή θα πρέπει $|\kappa |\le 1$, για κάθε $r\in ℝ$. Καταλήγουμε δηλαδή στο ίδιο αποτέλεσμα όπως και παραπάνω, δηλαδή ότι $\lambda \le {\displaystyle \frac{1}{2}}$.

Στη συνέχεια, δείχνουμε το ακόλουθο για τη σύγκλιση της μεθόδου.