7 Εξίσωση μεταφοράς

7.2 Μέθοδος των Lax–Wendroff

7.2.1 Μια μη ευσταθής μέθοδος

Στη συνέχεια, θα θεωρήσουμε ένα αριθμητικό σχήμα που να προσεγγίζει με μεγαλύτερη ακρίβεια τη ux, από τις δύο προηγούμενες μεθόδους upwind και downwind. Έτσι, αν στην (7.2) θεωρήσουμε για την προσέγγιση της ux(xi,tj) την δhcu(xi,tj) και για την ut(xi,tj) την δk+u(xi,tj), τότε έχουμε για i=1,,N, και j=0,,M-1,

u(xi,tj+1)-u(xi,tj)k+αu(xi+1,tj)-u(xi-1,tj)2h=ηij, (7.17)

όπου, αν utt,uxxxC([a,b]×[0,T]), λόγω των (2.3) και (2.4),

|ηij|k2maxt[0,T]|utt(x,t)|+h26maxx[a,b]|uxxx(x,t)|. (7.18)

Συνεπώς, ισχύει το ακόλουθο λήμμα.

Λήμμα 7.2.

Έστω u η λύση του (7.1) με utt,uxxxC([a,b]×[0,T]). Τότε για την ηij που δίνεται στην (7.17), έχουμε ότι υπάρχει σταθερά C ανεξάρτητη των k και h, τέτοια ώστε

max0jM-1max1iN|ηij|C(k+h2). (7.19)

Θεωρούμε, λοιπόν, την ακόλουθη μέθοδο για να κατασκευάσουμε προσεγγίσεις Uij των τιμών u(xi,tj),

Uij+1-Uijk+αUi+1j-Ui-1j2h=0,i=1,,N,j=0,,M-1Ui0=g(xi),i=0,,N+1U0j=ϕ0(tj),j=0,,M (7.20)
Σχήμα 7.8: Σημεία ενός πλέγματος για τον υπολογισμό της λύσης με τη μέθοδο (7.20). Με άσπρο σημειώνουμε τα σημεία που αντιστοιχούν σε άγνωστες τιμές της προσέγγισης και με μαύρο αυτά που αντιστοιχούν σε γνωστές τιμές.

Ξεκινώντας από το χρονικό επίπεδο t=0 όπου γνωρίζουμε ότι Ui0=g(xi), i=0,,N+1, και χρησιμοποιώντας την (7.20), μπορούμε να βρούμε με άμεσο τρόπο την προσέγγιση Ui1, i=0,,N, στο επόμενο χρονικό επίπεδο t1=k, όμως παρατηρούμε ότι δεν μπορούμε να υπολογίσουμε την UN+11. Επομένως, όπως και στη μέθοδο (7.11), για να βρούμε όλες τις τιμές Uij στα επόμενα χρονικά επίπεδα tj, j=1,,M, χρειάζεται να γνωρίζουμε την u(b,t), για t[0,T]. Έτσι, υποθέτουμε ότι στο πρόβλημα (7.1), η ακριβής λύση u ικανοποιεί επιπλέον και τη συνοριακή συνθήκη u(b,t)=ϕ1(t), t[0,T], οπότε θέτουμε UN+1j=ϕ1(tj), j=0,,M.

Συνεπώς, αν συμβολίσουμε τώρα λ=αkh, τότε η (7.20) μπορεί να γραφεί

Uij+1 =-λ2Ui+1j+Uij+λ2Ui-1j, i=1,,N,j=0,,M-1, (7.21)
Ui0 =g(xi), i=0,,N+1,
U0j =ϕ0(tj) και UN+1j=ϕ1(tj), j=0,,M.

Οπότε τώρα χρησιμοποιώντας την (7.21), βρίσκουμε με άμεσο τρόπο όλες τις τιμές Uij στα επόμενα χρονικά επίπεδα tj, j=1,,M.

Εύκολα παρατηρούμε ότι χρησιμοποιώντας τη μέθοδο (7.21), η τιμή της προσέγγισης στο σημείο (xi,tj), εξαρτάται τελικά από τις τιμές στα σημεία (xi-j,0), , (xi,0), , (xi+j,0), βλ. Σχήμα 7.9. Για να ικανοποείται τώρα η συνθήκη CFL αρκεί το διάστημα [xi-j,xi+j], το οποίο περιέχει τα σημεία του χωρίου υπολογιστικής εξάρτησης της μεθόδου, να περιέχει το x¯0=xi-αtj, το οποίο ισχύει ανεξάρτητα από το πρόσημο της σταθεράς α, αν |λ|=|αkh|1. Επίσης, ακολουθώντας τα βήματα για τη μέλετη της ευστάθιας von Neumann για τις μεθόδους upwind και downwind, έχουμε ότι, αν και πάλι θέσουμε Uij=wjerxiI, με r,

wj+1erxiI=-λ2wjerxi+1I+wjerxiI+λ2wjerxi-1I,

Οπότε, εφαρμόζοντας την τριγωνομετρική ταυτότητα (5.18), παίρνουμε

wj+1 =wj(-λ2wjerhI+1+λ2wje-rhI)
=wj(1-Iλsin(rh)).

Συνεπώς, wj=κjw0, με κ=1-Iλsin(rh), και

|κ|2=1+λ2sin2(rh).

Επομένως, η μέθοδος (7.21) δεν είναι von Neumann ευσταθής, γιατί |κ|>1.

Σχήμα 7.9: Χωρίο εξάρτησης της προσεγγιστικής λύσης με τις μεθόδους (7.21) και Lax–Wendroff.

7.2.2 Μέθοδος Lax–Wendroff

Στη συνέχεια, θα κατασκευάσουμε μια νέα μέθοδο για την προσέγγιση της (7.1) η οποία θα χρησιμοποιεί τα ίδια σημεία του πλέγματος όπως η (7.21), θα είναι von Neumann ευσταθής και θα έχει και αυτή τοπικό σφάλμα ηij το οποίο θα τείνει στο μηδέν, καθώς k και h τείνουν στο μηδέν. Θέτουμε, λοιπόν, για i=1,,N, j=0,,M-1,

u(xi,tj+1)=Au(xi+1,tj)+Bu(xi,tj)+Γu(xi-1,tj)+kηij, (7.22)

και θέλουμε να προσδιορίσουμε τα A,B και Γ, έτσι ώστε το τοπικό σφάλμα ηij να είναι όσο το δυνατόν μικρότερο. Χρησιμοποιώντας το γεγονός ότι ut=-αux, utt=-αuxt=-α(ut)x=α2uxx και uttt=-α3uxxx και αναπτύσσοντας σύμφωνα με το Θεώρημα Taylor, έχουμε

u(xi,tj+k)=u(xi,tj)+kut(xi,tj)+k22utt(xi,tj)+k36uttt(xi,tj)+η^ij=u(xi,tj)-αkux(xi,tj)+α2k22uxx(xi,tj)-α3k36uxxx(xi,tj)+η^ij, (7.23)

με

|η^ij|k424maxt[0,T]|t4u(xi,t)|, (7.24)

και όμοια

u(xi±h,tj)=u(xi,tj)±hux(xi,tj)+h22uxx(xi,tj)±h36uxxx(xi,tj)+η~i±j, (7.25)

με

|η~i±j|h424maxx[a,b]|x4u(x,tj)|. (7.26)

Αν αντικαταστήσουμε τις (7.23) και (7.25) στην (7.22), έχουμε

kηij=u(xi,tj+1)-Au(xi+1,tj)-Bu(xi,tj)-Γu(xi-1,tj)=u(xi,tj)-αkux(xi,tj)+α2k22uxx(xi,tj)+α3k36uxxx(xi,tj)+η^ij-A(u(xi,tj)+hux(xi,tj)+h22uxx(xi,tj)+h36uxxx(xi,tj)+η~i+j)-Bu(xi,tj)-Γ(u(xi,tj)-hux(xi,tj)+h22uxx(xi,tj)-h36uxxx(xi,tj)+η~i-j).

Επομένως

kηij=(1-A-B-Γ)u(xi,tj)-(αk+Ah-Γh)ux(xi,tj)+(α2k22-(A+Γ)h22)uxx(xi,tj)+(α3k36-(A-Γ)h36)uxxx(xi,tj)+η^ij-Aη~i+j-Γη~i-j=(1-A-B-Γ)u(xi,tj)-h(λ+A-Γ)ux(xi,tj)+h22(λ2-(A+Γ))uxx(xi,tj)+h36(λ3-A+Γ)uxxx(xi,tj)+η^ij-Aη~i+j-Γη~i-j.

Στη συνέχεια, εξισώνοντας τους συντελεστές των u,ux,uxx στην παραπάνω εξίσωση με το μηδέν, προκύπτουν οι ακόλουθες σχέσεις

1-A-B-Γ=0λ+A-Γ=0λ2-(A+Γ)=0.

Επομένως, έχουμε A=-12λ(1-λ), B=1-λ2 και Γ=12λ(1+λ). Συνεπώς, για αυτήν την επιλογή των A,B και Γ έχουμε

k|ηij|h36(|λ3|+|λ|)|uxxx(xi,tj)|+|η^ij|+λ2max(|η~i+j|,|η~i-j|)αk6(α2k2+h2)|uxxx(xi,tj)|+|η^ij|+α2k2h2max(|η~i+j|,|η~i-j|). (7.27)

Οπότε, λόγω των (7.24) και (7.26), έχουμε ότι υπάρχει σταθερά C ανεξάρτητη των k και h, τέτοια ώστε

max1iN,0jM-1|ηij|C(k2+h2).

Η μέθοδος που προκύπτει για την παραπάνω επιλογή των συντελεστών A,B και Γ ονομάζεται μέθοδος των Lax–Wendroff και για να μπορέσουμε να βρούμε την προσέγγιση σε κάθε χρονικό επίπεδο tj, j=1,,M, χρειάζεται, όπως και στην (7.21) να θεωρήσουμε γνωστή την ακριβή λύση για x=a και x=b. Έτσι, η μέθοδος Lax–Wendroff ορίζεται ως εξης

Uij+1=-12λ(1-λ)Ui+1j+(1-λ2)Uij+12λ(1+λ)Ui-1j,i=1,,N,j=0,,M-1Ui0=g(xi),i=0,,N+1U0j=ϕ0(tj), και UN+1j=ϕ1(tj),j=0,,M (7.28)

Εύκολα μπορούμε να δούμε ότι το χωρίο υπολογιστικής εξάρτησης για το (xi,tj) είναι το ίδιο με τη μέθοδο (7.21) και εξαρτάται από τα σημεία {xi-jh,,xi,,xi+jh}, βλέπε Σχήμα 7.9. Επίσης για να ικανοποείται η συνθήκη CFL πρέπει το σημείο x¯0=xi-αtj να περιέχεται στο διάστημα [xi-j,xi+j], το οποίο ισχύει ανεξάρτητα από το πρόσημο της σταθεράς α, αν |λ|=|αkh|1.

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

Θεωρούμε και πάλι το πρόβλημα των Παραδειγμάτων 7.1 και 7.2, στο διάστημα [-1,9], με g όπως στο Παράδειγμα 1.2, ϕ0(t)=ϕ1(t)=g(x-αt) και α=1. Επειδή γνωρίζουμε ότι u(a,t)=u(b,t)=0, για t[0,7], μπορούμε να θέσουμε U0j=UN+1j=0, j=0,,M. Επίσης, διαμερίζουμε το [-1,9] σε N+2 σημεία, με N=99 και το [0,7] σε M+1, με M=63,77,98. Στα Σχήματα 7.10 και 7.11 βλέπουμε την ακριβή λύση u και τις προσεγγίσεις με τη μέθοδο των Lax–Wendroff (7.28) για t=0,1,5. Παρατηρούμε ότι η προσεγγιστική λύση έχει παρόμοια συμπεριφορά όπως και αυτή της μεθόδου upwind στο Παράδειγμα 7.1. Έτσι αν λ=70/63>1, για M=63, τότε το σφάλμα της προσεγγιστικής λύσης είναι μεγάλο, ενώ το μικρότερο σφάλμα προκύπτει αν λ<1 και είναι πιο κοντά στη μονάδα.

Σχήμα 7.10: Η ακριβής λύση u(x,t) του Παραδείγματος 7.4 και οι προσεγγίσεις της για N=99 και M=63,77,98, με τη μέθοδο Lax-Wendroff για t=0 και 1.
Σχήμα 7.11: Η ακριβής λύση u(x,t) του Παραδείγματος 7.4 και οι προσεγγίσεις της για N=99 και M=63,77,98, με τη μέθοδο Lax-Wendroff για t=5.

7.2.3 Ευστάθεια και σύγκλιση

Για τη μελέτη της ευστάθειας von Neumann για τη μέθοδο των Lax–Wendroff, θα θεωρήσουμε και πάλι ότι η προσεγγιστική λύση είναι της μορφής Uij=wjerxiI, r, οπότε σύμφωνα με την (7.28), έχουμε

wj+1erxiI= -12λ(1-λ)wjerxi+1I+(1-λ2)wjerxiI
 +12λ(1+λ)wjerxi-1I.

Χρησιμοποιώντας τώρα τις (5.23) και (5.18), παίρνουμε

wj+1=(-12λ(1-λ)+12λ(1+λ))wjcos(h)
+(1-λ2)wj-I(12λ(1-λ)+12λ(1+λ))wjsin(h)
=wj(1-λ2+λ2cos(h)-Iλsin(rh)).

Συνεπώς, wj=κw0, με κ=1-λ2+λ2cos(h)-Iλsin(rh). Οπότε

|κ|2=(1-λ2+λ2cos(h))2+λ2sin2(h)=(1-λ2)2+λ4cos2(h)+2λ2(1-λ2)cos(h)+λ2sin2(h)=(1-λ2)2+λ4(1-sin2(h))+2λ2(1-λ2)cos(h)+λ2sin2(h)=1-2λ2(1-λ2)+λ2(1-λ2)sin2(h)+2λ2(1-λ2)cos(h)=1-2λ2(1-λ2)[1-cos(h)]+λ2(1-λ2)sin2(h)=1-4λ2(1-λ2)sin2(h/2)+λ2(1-λ2)sin2(h).

Χρησιμοποιώντας τώρα την τριγωνομετρική ιδιότητα sin2(h)=4sin2(h/2)-4sin4(h/4), η παραπάνω σχέση γίνεται

|κ|2=1-4λ2(1-λ2)sin4(h/4).

Είναι προφανές ότι, αν λ21, τότε |κ|1, και, άρα, η μέθοδος Lax–Wendroff είναι von Neumann ευσταθής, αν λ[-1,1].

Για να δείξουμε τη σύγκλιση της μεθόδου Lax–Wendroff χρειάζεται να εφαρμόσουμε παρόμοια επιχειρήματα όπως αυτά της μεθόδου ενέργειας στην Παράγραφο 3.4, βλ. παραδείγματος χάριν (Larsson and V. Thomée, (2009), Παράγραφος 12.3), η απόδειξη του οποίου είναι πολύπλοκη και ξεφεύγει από τους σκοπούς αυτών των σημειώσεων. Διατυπώνουμε, επομένως, το επόμενο θεώρημα χωρίς απόδειξη.

Θεώρημα 7.3.

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

max1jMmax0iN+1|Uij-u(xi,tj)|C(k2+h2). (7.29)

Στο ακόλουθο παράδειγμα βλέπουμε και πειραματικά ότι το σφάλμα με τη μέθοδο Lax–Wendroff έχει μεγαλύτερη τάξη προσέγγισης ως προς k και h από το σφάλμα της μεθόδου upwind ή downwind.

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

Βασική υπόθεση στα θεωρήματα σύγκλισης για τις μεθόδους που μελετήσαμε σε αυτό το κεφάλαιο είναι ότι η ακριβής λύση u, και κατά επέκταση η συνάρτηση της αρχικής συνθήκης g, είναι αρκετά ομάλη. Είναι προφανές ότι η g που θεωρήσαμε στα προηγούμενα παραδείγματα δεν είναι ομαλή συνάρτηση, για αυτό τον λόγο τώρα θεωρούμε το πρόβλημα (7.13) στο διάστημα [-11,9]×[0,7], με

g(x)={e-11-x2,αν -1x1,0διαφορετικά.

Επειδή η σταθερά α του (7.13) είναι θετική για να ικανοποιείται η συνθήκη CFL θα χρησιμοποιήσουμε τις μεθόδους upwind και Lax–Wendroff. Θεωρούμε τώρα διαμερίσεις του [0,7] με k=0.05, 0.025, 0.0125, 0.0063, 0.0031, 0.0016, 0.0008 και του [-11,9], με αντίστοιχο βήμα h, τέτοιο ώστε k/h=0.9. Θέτουμε E¯upwind=max0N+1|UiM-u(xi,tM)|, με tM=T=7 το σφάλμα της προσεγγιστικής λύσης που προκύπτει με τη μεθόδο upwind και με ανάλογο τρόπο ορίζουμε το σφάλμα E¯lax με τη μέθοδο των Lax–Wendroff. Στον πίνακα 7.1 βλέπουμε τα σφάλματα αυτών των δύο μεθόδων και την προσεγγιστική τάξη ακριβείας p.

Upwind Lax–Wendroff
k E¯upwind p E¯lax p
0.05 0.0483 0.01965
0.025 0.0325 0.574 0.00960 1.034
0.0125 0.0208 0.642 0.00413 1.216
0.0063 0.0127 0.713 0.00151 1.456
0.0031 0.0073 0.785 0.00045 1.743
0.0016 0.0040 0.852 0.00012 1.848
0.0008 0.0021 0.906 0.00003 1.992
Πίνακας 7.1: Τα σφάλματα των μεθόδων upwind και Lax–Wendroff στο χρονικό επίπεδo t=T=7, E¯upwind και E¯lax, αντίστοιχα, του Παραδείγματος 7.5, καθώς και η κατά προσέγγιση τάξη ακρίβειά τους p, αν θεωρήσουμε διαμερίσεις όπου k/h=0.9.