5 Εξίσωση θερμότητας: Πεπερασμένες διαφορές

5.2 Πεπλεγμένη μέθοδος του Euler

Σε αυτή την παράγραφο θα μελετήσουμε τη μέθοδο που προκύπτει αν αντί για την προσέγγιση, δk+ της ut που χρησιμοποιήσαμε στην άμεση μέθοδο του Euler, θεωρήσουμε την δk-.

Έτσι, θα θεωρήσουμε και πάλι μια διαμέριση του [0,L]×[0,T] όπως και πριν και θα κατασκευάσουμε προσεγγίσεις Uji των τιμών u(xi,tj), i=0,,N+1, j=0,,M, της ακριβούς λύσης του (5.1), βλέπε Σχήμα 5.1. Λόγω των αρχικών και συνοριακών συνθηκών του προβλήματος (5.1) θέτουμε και πάλι U0j=UN+1j=0, j=0,,M, και Ui0=g(xi), i=1,,N. Στη συνέχεια, προσεγγίζουμε την ut(xi,tj) στις εξισώσεις (5.3) χρησιμοποιώντας τη διαφορά δk-u(xi,tj) που θεωρήσαμε στην (2.1).

Έτσι, αν υποθέσουμε ότι utt(x,t)C([0,L]×[0,T]), λόγω της (2.3), η (5.3) γίνεται για i=1,,N,

u(xi,tj)-u(xi,tj-1)k=uxx(xi,tj)+η^ij,j=1,,M, (5.25)

όπου

|η^ij|k2maxt[0,T]|utt(xi,t)|. (5.26)

Στη συνέχεια, για να προσεγγίσουμε την uxx(xi,tj), στην (5.25) χρησιμοποιούμε, όπως και στην άμεση μέθοδο του Euler, την δh,2cu(xi,tj). Έτσι, αν υποθέσουμε ότι 4x4u(x,t)C([0,L]×[0,T]), λόγω της (2.9), η (5.25) γίνεται για i=1,,N, και j=1,,M,

u(xi,tj)-u(xi,tj-1)k=u(xi+1,tj)-2u(xi,tj)+u(xi-1,tj)h2+η^ij+η~ij, (5.27)

όπου η η~ij φράσσεται και πάλι όπως στην (5.7).

Συνεπώς, για ηij=η^ij+η~ij, η (5.27) δίνει για i=1,,N, και j=1,,M,

δk-u(xi,tj)-δh,2cu(xi,tj)=ηij, (5.28)

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

Λήμμα 5.2.

Έστω u η λύση του (5.1) με 4x4uC([0,L]×[0,T]) και uttC([0,L]×[0,T]). Τότε, για την ηij που δίνεται στην (5.28), έχουμε ότι υπάρχει σταθερά C ανεξάρτητη των k και h, τέτοια ώστε

max1jMmax1iN|ηij|C(k+h2). (5.29)

Για να κατασκευάσουμε, λοιπόν, προσεγγίσεις Uij των u(xi,tj), i=1,,N, j=1,,M, θεωρούμε την ακόλουθη μέθοδο

Uij-Uij-1k=Ui+1j-2Uij+Ui-1jh2,i=1,,N,j=1,,M, (5.30)
{U0j=UN+1j=0,j=0,,M,Ui0=g(xi),i=1,,N. (5.31)

Είναι φανερό ότι, λόγω της (5.28), το τοπικό σφάλμα διακριτοποίησης της μεθόδου (5.30)–(5.31), δίνεται από την ηij. Συνεπώς, σύμφωνα με το Λήμμα 5.2, το τοπικό σφάλμα διακριτοποίησης τείνει στο μηδέν, καθώς τα h και k τείνουν στο μηδέν, οπότε η (5.30) είναι μια συνεπής μέθοδος.

Συμβολίζουμε και πάλι με λ τον λόγο kh2, οπότε η (5.30) γράφεται ως

-λUi+1j+(1+2λ)Uij-λUi-1j=Uij-1, (5.32)

με i=1,,N και j=1,,M. Ξεκινώντας από το χρονικό επίπεδο t0=0 και χρησιμοποιώντας τις αρχικές συνθήκες (5.31), παρατηρούμε ότι, σε αντίθεση με την άμεση μέθοδο του Euler (5.12), δεν μπορούμε άμεσα να υπολογίσουμε την προσέγγιση στο χρονικό επίπεδο t1 από την (5.32). Στο Σχήμα 5.4 φαίνεται ένα παράδειγμα πλέγματος όπου με άσπρο κύκλο σημειώνουμε το σημείο του χρονικού επιπέδου tj+1, όπου σύμφωνα με (5.32), αντιστοιχεί σε άγνωστες τιμές της προσέγγισης και με μαύρο αυτά του χρονικού επιπέδου tj που αντιστοιχούν στη γνωστή τιμή της προσέγγισης.

Σχήμα 5.4: Σημεία ενός πλέγματος για τον υπολογισμό της λύσης με την πεπλεγμένη μέθοδο του Euler. Με άσπρο σημειώνουμε τα σημεία που αντιστοιχούν σε άγνωστες τιμές της προσέγγισης και με μαύρο αυτά που αντιστοιχούν σε γνωστές τιμές.

Αν θέσουμε UjN το διάνυσμα με συνιστώσες U1j,,UNj, U=(U1j, , UNj)T, μπορούμε να γράψουμε το σύστημα των εξισώσεων (5.32) ισοδύναμα ως το γραμμικό σύστημα

BUj=Uj-1,j=1,,M, με U0=G, (5.33)

όπου B είναι ο N×N πίνακας

B=(1+2λ-λ00-λ1+2λ-λ000-λ1+2λ-λ-λ1+2λ),

και G=(g(x1),,g(xN))T.

Αφού γνωρίζουμε τη U0 για να υπολογίσουμε τώρα την προσέγγιση στο χρονικό επίπεδο t1 χρειάζεται να λύσουμε το γραμμικό σύστημα BU1=U0. Είναι φανερό ότι ο B έχει αυστηρά κυριαρχική διαγώνιο και, άρα, είναι αντιστρέψιμος. Επομένως, το γραμμικό σύστημα BU1=U0 έχει μοναδική λύση. Επειδή ο πίνακας B είναι τριδιαγώνιος για τον υπολογισμό της U1, μπορούμε να χρησιμοποιήσουμε τη μέθοδο της Παραγράφου 3.2. Συνεχίζουμε με αυτό τον τρόπο, οπότε σε κάθε χρονικό βήμα γνωρίζουμε την προσέγγιση στο χρονικό επίπεδο tj-1 και για να υπολογίσουμε την προσέγγιση στο χρονικό επίπεδο tj λύνουμε το γραμμικό σύστημα (5.33). Επειδή, κάθε φορά χρειάζεται να λύσουμε ένα γραμμικό σύστημα, η μέθοδος καλείται πεπλεγμένη ή έμμεση. Συγκεκριμένα, τη μέθοδο (5.30) ή (5.32) την ονομάζουμε πεπλεγμένη μέθοδο του Euler.

Ο πίνακας B στο γραμμικό σύστημα (5.33) είναι τριδιαγώνιος και έχει αυστηρά κυριαρχική διαγώνιο, οπότε είναι αντιστρέψιμος και άρα η ύπαρξη και η μοναδικότητα της προσεγγιστικής λύσης σε κάθε χρονικό βήμα είναι δεδομένη.

Για την ευστάθεια της μεθόδου μπορούμε να δείξουμε το ακόλουθο θεώρημα.

Θεώρημα 5.3.

Έστω UjN, j=0,,M, τα διανύσματα που ικανοποιούν την (5.33), και U0j=UN+1j=0, j=0,,M. Τότε, ισχύει η ακόλουθη ανισότητα

max0iN+1|Uij|max0iN+1|Ui0|,για j=0,1,2,M.
Απόδειξη.

Από τη σχέση (5.32) εύκολα παίρνουμε

(1+2λ)|Uij|λ(|Ui+1j|+|Ui-1j|)+|Uij-1|,1iN,j=1,,M.

Αν θέσουμε λοιπόν U¯j=max0iN+1|Uij|, j=0,,M, παίρνουμε

(1+2λ)|Uij|2λU¯j+U¯j-1,j=1,,M,

και άρα

U¯jU¯j-1U¯0,

από όπου προκύπτει η ζητούμενη εκτίμηση. ∎

Στη συνέχεια, θα μελετήσουμε την ευστάθεια von Neumann για την πεπλεγμένη μέθοδο του Euler. Όπως και στην άμεση μέθοδο του Euler, υποθέτουμε ότι οι

Uij=wjerxiI,r,

ικανοποιούν την (5.32). Οπότε

-wj(λerxi+1I-(1+2λ)erxiI+λerxi-1I)=wj-1erxiI. (5.34)

Χρησιμοποιούμε τώρα και πάλι τις ιδιότητες (5.23) και (5.20) και έχουμε

-wj(2λcos(rh)-(1+2λ))=wj(1+4λsin2(rh2))=wj-1. (5.35)

Συνεπώς, wj=κjw0, με κ=1/(1+4λsin2(rh2)). Άρα, επειδή για κάθε λ και r, ισχύει |κ|1, η απόλυτη τιμή της Uij είναι φραγμένη, για κάθε σημείο του διαμερισμού, και άρα δεν είναι αναγκαία καμία συνθήκη για την ευστάθεια της πεπλεγμένης Euler.

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

Θεώρημα 5.4.

Έστω ότι η λύση u του προβλήματος (5.1) είναι αρκετά ομαλή και Uij, i=0,,N+1, j=0,,M, η λύση της μεθόδου (5.32). Τότε, υπάρχει μια σταθερά C, ανεξάρτητη των k,h, τέτοια ώστε

max0jMmax0iN+1|Uij-u(xi,tj)|C(k+h2). (5.36)
Απόδειξη.

Θέτουμε Eij=Uij-u(xi,tj), i=0,,N+1, j=0,,M, όπου λόγω των σχέσεων U0j=u(0,tj)=0 και UN+1j=u(L,tj)=0, έχουμε E0j=EN+1j=0. Αφαιρούμε τώρα κατά μέλη τις (5.32) και (5.27), οπότε παίρνουμε

-λEij+(1+2λ)Eij-λEi-1j=Eij-1+kηij,i=1,,N,j=1,,M.

Θέτουμε, στη συνέχεια, E¯j=max1iN|Eij|, η¯=max1iN,1jM|ηij| και έχουμε

(1+2λ)Eij2λE¯j+E¯j-1+kη¯,

οπότε

E¯jE¯j-1+kη¯E¯0+jkη¯.

Από εδώ, επειδή jkT και λόγω του Λήμματος 5.2, προκύπτει η ζητούμενη εκτίμηση. ∎


Παράδειγμα 5.2.

Αν στο Παράδειγμα 5.1 χρησιμοποιήσουμε την πεπλεγμένη μέθοδο του Euler, παρατηρούμε ότι οι αντίστοιχες προσεγγίσεις Uij της u(xi,tj) είναι καλύτερες. Στο Σχήμα 5.5 απεικονίζονται οι γραφικές παραστάσεις της u για t=0, καθώς και οι προσεγγίσεις για t=0.025,0.05,0.0625, με M=24,32,128. Επίσης, στον Πίνακα 5.2 βλέπουμε το μέγιστο σφάλμα στα χρονικά επίπεδα t=0.025,0.05,0.0625,0.075. Στο Σχήμα 5.5 παρατηρούμε ότι όλες οι προσεγγιστικές λύσεις φαίνεται να προσεγγίζουν την ακριβή λύση, καθώς το χρονικό επίπεδο t μεγαλώνει, το οποίο φαίνεται και στα αποτελέσματα του Πίνακα 5.1. Αυτή η συμπεριφορά των προσεγγιστικών λύσεων οφείλεται στην ευστάθεια της μεθόδου για όλες τις τιμές της παράμετρου λ, σε αντίθεση με τη συμπεριφορά της άμεσης μεθόδου του Euler, βλ. Παράδειγμα 5.1.

Σχήμα 5.5: Η ακριβής λύση u(x,t)=e-4π2tsin(2πx) του Παραδείγματος 5.1 και οι προσεγγίσεις της για N=23 και M=24,32,128, με την πεπλεγμένη μέθοδο του Euler για t=0,0.025,0.05 και 0.0625.

t=0.025 t=0.05 t=0.0625 t=0.075
M j E¯j j E¯j j E¯j j E¯j
24 6 0.0303 12 0.0235 15 0.0183 18 0.0137
32 8 0.0236 16 0.0181 20 0.0140 24 0.0105
128 32 0.0077 64 0.0058 80 0.0044 96 0.0033
Πίνακας 5.2: Το σφάλμα E¯j=max1iN|Uij-u(xi,tj)| στα χρονικά επίπεδα t=0.025,0.05 και 0.0625 του Παραδείγματος 5.2 για N=23 και M=24,32,128. Στη στήλη αριστερά του κάθε σφάλματος δίνεται το χρονικό βήμα j, τέτοιο ώστε t=jk, με k=1/M.