8 Εξίσωση κύματος 8 Εξίσωση κύματος 8.2 Άλλες μέθοδοι και προβλήματα

8.1 Μέθοδος των Courant, Friedrichs και Lewy

Θεωρούμε το ακόλουθο πρόβλημα: Ζητείται μια συνάρτηση $u:[0,L]\times[0,T]\to\mathbb{R}$ , τέτοια ώστε

\begin{split}\displaystyle u_{tt}(x,t)&\displaystyle=\alpha^{2}u_{xx}(x,t),% \qquad x\in[0,L],\ t\in[0,T],\\ \displaystyle u(0,t)&\displaystyle=u(L,t)=0,\quad t\in[0,T],\\ \displaystyle u(x,0)&\displaystyle=g_{0}(x),\qquad x\in[0,L],\\ \displaystyle u_{t}(x,0)&\displaystyle=g_{1}(x),\qquad x\in[0,L],\end{split}

(8.1)

όπου $\alpha,L>0$ , $g_{0},g_{1}\in C[0,L]$ . Είναι γνωστό, βλ. π.χ. (Ακρίβης και Αλικάκος, (2012)), ότι η λύση του (8.1) δίνεται από τη σειρά

u(x,t)=\sum_{n=1}^{\infty}[A_{n}\cos(\alpha\lambda_{n}t)+B_{n}\sin(\alpha% \lambda_{n}t)]\sin(\lambda_{n}x),\text{ με }\lambda_{n}=\dfrac{n\pi}{L},

και $A_{n}=\dfrac{2}{L}\displaystyle\int_{0}^{L}g_{0}(x)\sin(\lambda_{n}x)$ , $B_{n}=\dfrac{2}{\alpha n\pi}\displaystyle\int_{0}^{L}g_{1}(x)\sin(\lambda_{n}x)$ .

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

u_{tt}(x_{i},t^{j})=\alpha^{2}u_{xx}(x_{i},t^{j}),\quad i=1,\dots,N,\ j=0,% \dots,M.

(8.2)

Θα προσεγγίσουμε καταρχήν τις $u_{xx}$ και $u_{tt}$ στα σημεία του διαμερισμού και θα χρησιμοποιούμε τις προσέγγισεις $\delta^{c}_{h,2}u(x_{i},t^{j})$ και $\delta^{c}_{κ,2}u(x_{i},t^{j})$ , αντίστοιχα, που θεωρήσαμε στην (2.8). Έτσι, αν υποθέσουμε ότι $\frac{\partial^{4}u}{\partial x^{4}},\frac{\partial^{4}u}{\partial t^{4}}\in C% ([0,L]\times[0,T])$ , λόγω της (2.11), η (8.2) γίνεται

\begin{split}&\displaystyle\dfrac{u(x_{i},t^{j+1})-2u(x_{i},t^{j})+u(x_{i},t^{% j-1})}{k^{2}}\\ &\displaystyle=\alpha^{2}\dfrac{u(x_{i+1},t^{j})-2u(x_{i},t^{j})+u(x_{i-1},t^{% j})}{h^{2}}+\eta_{i}^{j},\end{split}

(8.3)

για $i=1,\dots,N$ , $j=1,\dots,M-1$ , όπου

|\eta_{i}^{j}|\leq\dfrac{k^{2}}{12}\max_{t\in[0,T]}|\frac{\partial^{4}u}{% \partial t^{4}}(x_{i},t)|+\dfrac{\alpha^{2}h^{2}}{12}\max_{x\in[0,L]}|\frac{% \partial^{4}u}{\partial x^{4}}(x,t^{j})|.

(8.4)

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

		$\displaystyle\begin{aligned}\displaystyle\dfrac{U_{i}^{j+1}-2U_{i}^{j}+U_{i}^{% j-1}}{k^{2}}&\displaystyle=\alpha^{2}\dfrac{U_{i+1}^{j}-2U_{i}^{j}+U_{i-1}^{j}% }{h^{2}},\\ &\displaystyle\quad\text{για }i=1,\dots,N,j=1,\dots,M-1,\end{aligned}$		(8.5)
		$\displaystyle\left\{\begin{aligned}\displaystyle U_{0}^{j}&\displaystyle=U_{N+% 1}^{j}=0,\quad j=0,\dots,M,\\ \displaystyle U_{i}^{0}&\displaystyle=g_{0}(x_{i}),\quad i=1,\dots,N.\end{% aligned}\right.$		(8.6)

Συμβολίζουμε τώρα με $\lambda$ τον λόγο $\dfrac{\alpha k}{h}$ , οπότε η (8.5) γράφεται ως

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

(8.7)

για $i=1,\dots,N$ , $j=1,\dots,M-1$ . Παρατηρούμε ότι χρησιμοποιώντας την (8.7) για να βρούμε την προσέγγιση στο χρονικό επίπεδο $t^{j+1}$ , χρειάζεται να γνωρίζουμε τις προσεγγίσεις στα χρονικά επίπεδα $t^{j}$ και $t^{j-1}$ . Έτσι ενώ, λόγω των (8.6), γνωρίζουμε τις $U^{0}_{i}$ , $i=0,\dots,N+1$ , για να βρούμε την προσέγγιση στο χρονικό επίπεδο $t^{2}$ , χρειάζεται να γνωρίζουμε την προσέγγιση και στο χρονικό επίπεδο $t^{1}$ . Η οποία όμως δεν είναι γνωστή και πρέπει να την προσδιορίσουμε με μια άλλη μέθοδο.

Αν υποθέσουμε ότι $\frac{\partial^{3}u}{\partial t^{3}}\in C([0,L]\times[0,T])$ , χρησιμοποιώντας το ανάπτυγμα Taylor της $u$ στα σημεία $(x_{i},t^{1})$ , $i=1,\dots,N$ , έχουμε

\begin{split}\displaystyle u(x_{i},t^{1})&\displaystyle=u(x_{i},t^{0})+ku_{t}(% x_{i},0)+\dfrac{k^{2}}{2}u_{tt}(x_{i},0)+\hat{\eta}_{i}^{0}\\ &\displaystyle=g_{0}(x_{i})+\alpha^{2}kg_{1}(x_{i})+\dfrac{\alpha^{4}k^{2}}{2}% u_{xx}(x_{i},0)+\hat{\eta}_{i}^{0},\end{split}

(8.8)

όπου

|\hat{\eta}_{i}^{0}|\leq\dfrac{k^{3}}{6}\max_{t\in[0,T]}|\frac{\partial^{3}u}{% \partial t^{3}}(x_{i},t)|

(8.9)

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

\begin{split}&\displaystyle u(x_{i},t^{1})=g_{0}(x_{i})+\alpha^{2}kg_{1}(x_{i}% )+\dfrac{\lambda^{2}}{2}(u(x_{i+1},0)-2u(x_{i},0)+u(x_{i-1},0))\\ &\displaystyle\quad+\hat{\eta}_{i}^{0}+\tilde{\eta}_{i}^{0}\\ &\displaystyle\ =g_{0}(x_{i})+\alpha^{2}kg_{1}(x_{i})+\dfrac{\lambda^{2}}{2}(g% (x_{i+1})-2g(x_{i})+g(x_{i-1}))+\eta_{i}^{0},\end{split}

(8.10)

όπου $\eta_{i}^{0}=\hat{\eta}_{i}^{0}+\tilde{\eta}_{i}^{0}$ και

|\tilde{\eta}_{i}^{0}|\leq\dfrac{\alpha^{4}k^{2}h^{2}}{24}\max_{x\in[0,L]}|% \frac{\partial^{4}u}{\partial x^{4}}(x,0)|.

(8.11)

Λόγω τώρα των (8.4), (8.9) και (8.11), έχουμε ότι ισχύει το ακόλουθο λήμμα.

Λήμμα 8.1.

Έστω $u$ η λύση του (8.1) με $\frac{\partial^{4}}{\partial x^{4}}u,\frac{\partial^{4}}{\partial t^{4}}u\in C% ([0,L]\times[0,T])$ . Τότε για την $\eta_{i}^{j}$ που δίνεται στην (8.3) και (8.10), έχουμε ότι υπάρχει σταθερά $C$ ανεξάρτητη των $k$ και $h$ , τέτοια ώστε

\max_{0\leq j\leq M-1}\max_{1\leq i\leq N}|\eta_{i}^{j}|\leq C(k^{2}+h^{2}).

(8.12)

Επομένως, συμπληρώνουμε τις εξισώσεις (8.3) με τις ακόλουθες προσεγγίσεις $U_{i}^{1}$ της λύσης $u$ του προβλήματος (8.1) στα σημεία $(x_{i},t^{1})$ , $i=0,\dots,N+1$ , οι οποίες ικανοποιούν για $i=1,\dots,N$ ,

\begin{split}\displaystyle U_{i}^{1}&\displaystyle=g_{0}(x_{i})+\alpha^{2}kg_{% 1}(x_{i})+\dfrac{\lambda^{2}}{2}(g_{0}(x_{i+1})-2g_{0}(x_{i})+g_{0}(x_{i-1}))% \\ &\displaystyle=\dfrac{\lambda^{2}}{2}g_{0}(x_{i+1})+(1-\lambda^{2})g_{0}(x_{i}% )+\dfrac{\lambda^{2}}{2}g_{0}(x_{i-1})+\alpha^{2}kg_{1}(x_{i}).\end{split}

(8.13)

Στη συνέχεια, αν συμβολίσουμε με $U^{j}\in\mathbb{R}^{N}$ το διάνυσμα με συνιστώσες $U_{1}^{j},\dots,U_{N}^{j}$ , $U^{j}=$ $(U_{1}^{j},$ $\dots$ , $U_{N}^{j})^{T}$ , μπορούμε να γράψουμε το σύστημα των εξισώσεων (8.3) και (8.13) ισοδύναμα ως

\begin{split}\displaystyle U^{j+1}&\displaystyle=AU^{j}-U^{j-1},\quad j=1,% \dots,M-1,\\ \displaystyle U^{1}&\displaystyle=\dfrac{1}{2}AG_{0}+\alpha^{2}kG_{1},\quad% \text{με }U^{0}=G_{0},\end{split}

(8.14)

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

A=\left(\begin{array}[]{ccccc}2(1-\lambda^{2})&\lambda^{2}&0&&0\\ \lambda^{2}&2(1-\lambda^{2})&\lambda^{2}&0&\\ 0&\ddots&\ddots&\ddots&0\\ &&\lambda^{2}&2(1-\lambda^{2})&\lambda^{2}\\ &&&\lambda^{2}&2(1-\lambda^{2})\\ \end{array}\right),

και $G_{0}=(g_{0}(x_{1}),\dots,g_{0}(x_{N}))^{T}$ και $G_{1}=(g_{1}(x_{1}),\dots,g_{1}(x_{N}))^{T}$ .

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

Εύκολα παρατηρούμε ότι, χρησιμοποιώντας τη μέθοδο (8.3), η τιμή της προσέγγισης στο σημείο $(x_{i},t^{j})$ , εξαρτάται τελικά από τις τιμές στα σημεία $(x_{i-j},0)$ , $\dots$ , $(x_{i},0)$ , $\dots$ , $(x_{i+j},0)$ , βλ. Σχήμα 8.2. Άρα, το χωρίο υπολογιστικής εξάρτησης της μεθόδου (8.14) για το $(x_{i},t^{j})$ είναι το διάστημα $[x_{i-j},x_{i+j}]$ του $x-$ άξονα. Επίσης, λόγω της αναπαράστασης d’Alembert της ακριβούς λύσης $u$ , βλ. (1.48), η τιμή της $u$ στο σημείο $(x_{i},t^{j})$ εξαρτάται από τις τιμές στα σημεία $(x_{i}-\alpha t^{j})$ και $(x_{i}+\alpha t^{j})$ . Επομένως για να ικανοποιείται η συνθήκη CFL για τη μέθοδο (8.14) πρέπει τα σημεία $x_{i}\pm\alpha t^{j}\in[x_{i}-jh,x_{i}+jh]$ , δηλαδή $\lambda\leq 1$ .

Σχήμα 8.2: Χωρίο εξάρτησης της προσεγγιστικής λύσης με τη μέθοδο των Courant, Friedrichs και Lewy.

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

Θεώρημα 8.1.

Έστω ότι η λύση $u$ του προβλήματος (8.1) είναι αρκετά ομαλή και $U_{i}^{j}$ , $i=1,\dots,N$ , $j=0,\dots,M$ η λύση της μεθόδου των Courant, Friedrichs και Lewy, (8.14). Τότε, αν $0<\lambda\leq 1$ , υπάρχει μια σταθερά $C$ , ανεξάρτητη των $k, h$ , τέτοια ώστε

\max_{1\leq j\leq M}\max_{0\leq i\leq N+1}|U_{i}^{j}-u(x_{i},t^{j})|\leq C(k^{% 2}+h^{2}).

(8.15)

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

Θεωρούμε το πρόβλημα (8.1), με $\alpha=1$ και $g_{0}$ που δίνεται από την

g_{0}(x)=\begin{cases}\dfrac{1}{2}(1-\cos(\dfrac{2\pi x}{c})),&αν0\leq x\leq c% ,\\ 0,&\text{διαφορετικά},\end{cases}

(8.16)

όπου $0<c<1$ και $g_{1}(x)=-g_{0}^{\prime}(x)$ . Μπορούμε να δούμε ότι η ακριβής λύση στο $[0,1]\times[0,T]$ , με $T<2-c$ δίνεται από τη συνάρτηση

u(x,t)=\begin{cases}g_{0}(x-t),&αν0\leq x\leq 1-c\\ g_{0}(x-t)-g_{0}(2-x-t),&1-c\leq t\leq 2-c.\end{cases}

(8.17)

Υποθέτουμε ότι $c=0.09$ και $T=1.8$ οπότε για $N=199$ , θέτουμε $M=330,360$ και $400$ . Συνεπώς, για $M=330$ ο λόγος $\lambda>1$ , ενώ για $M=360$ είναι $\lambda=1$ και για $M=400$ , $\lambda<1$ . Στο Σχήμα 8.3 βλέπουμε την ακριβή λύση $u$ και τις προσεγγίσεις με τη μέθοδο Courant, Friedrichs και Lewy, για $t=0,0.1,0.5,0.9,1.2$ και $1.8$ . Παρατηρούμε, λοιπόν, ότι για $M=330$ η προσέγγιση δεν δίνει καλά αποτελέσματα, το οποίο οφείλεται στο γεγονός ότι ο λόγος $\lambda>1$ , ενώ στις άλλες δύο περιπτώσεις το σφάλμα είναι μικρό.

Σχήμα 8.3: Η ακριβής λύση $u(x,t)$ του Παραδείγματος 8.1 και οι προσεγγίσεις της για $N=199$ και $M=330,360,400$ , με τη μέθοδο των Courant, Friedrichs και Lewy για $t=0,0.1,0.5,0.9,1.2,1.8$ .

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

Θεωρούμε το πρόβλημα (8.1) στο $[-2,2]\times[0,1]$ , με $\alpha=1$ και

g(x)=\begin{cases}e^{-\dfrac{1}{1-x^{2}}},&\text{αν }-1\leq x\leq 1,\\ 0,&\text{διαφορετικά},\end{cases}

και $g_{1}=0$ . Μπορούμε να δούμε ότι η ακριβής λύση στο $[-2,2]\times[0,1]$ είναι η $u(x,t)=\frac{1}{2}(g_{0}(x+t)+g_{0}(x-t))$ . Χρησιμοποιούμε τη μέθοδο των Courant, Friedrichs και Lewy, όπου θεωρούμε διαμερίσεις του $[0,1]$ με $k=0.05$ , $0.025$ , $0.0125$ , $0.0063$ , $0.0031$ και $0.0016$ και του $[-2,2]$ , με αντίστοιχο βήμα $h$ , τέτοιο ώστε $k/h=0.75$ . Στον Πίνακα 8.1 βλέπουμε το σφάλμα $\bar{E}=\max_{0\leq i\leq N+1}|U_{i}^{M}-u(x_{i},t^{M})|$ στο χρονικό επίπεδο $t^{M}=T=1$ και την προσεγγιστική τάξη ακρίβειας $p$ του σφάλματος.

$k$	$\bar{E}$	$p$
0.05	0.00640
0.025	0.00117	2.450
0.0125	0.00043	1.433
0.0063	0.00011	1.915
0.0031	0.00003	1.933
0.0016	0.00001	2.023

Πίνακας 8.1: Το σφάλμα

\bar{E}

, στο χρονικό επίπεδo

t=T=1

που προκύπτει για το Παράδειγμα 8.2 με τη μεθόδο των Courant, Friedrichs και Lewy με

k/h=0.75

και η κατά προσέγγιση τάξη ακρίβειας

p

Θα μελετήσουμε τώρα την ευστάθεια von Neumann της μεθόδου των Courant, Friedrichs και Lewy. Αν υποθέσουμε, λοιπόν, ότι οι $U_{i}^{j}=\kappa^{j}e^{rx_{i}I}$ ικανοποιούν την (8.7), έχουμε

\displaystyle\kappa^{j+1}e^{rx_{i}I}=

\displaystyle\kappa^{j}(\lambda^{2}e^{rx_{i+1}I}+2(1-\lambda^{2})e^{rx_{i}I})+% \lambda^{2}e^{rx_{i-1}I})-\kappa^{j-1}e^{rx_{i}I}.

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

	$\displaystyle\kappa^{2}=$	$\displaystyle\kappa(\lambda^{2}e^{rhI}+2(1-\lambda^{2})+\lambda^{2}e^{-rhI})-1$
	$\displaystyle=$	$\displaystyle\kappa(2\lambda^{2}\cos(rh)+2(1-\lambda^{2}))-1$
		$\displaystyle=2\kappa(1-2\lambda^{2}\sin^{2}(\frac{rh}{2}))-1.$

Από το οποίο προκύπτει

\kappa^{2}-2\kappa\sigma+1=0,\quad\text{ με }\sigma=(1-2\lambda^{2}\sin^{2}(% \frac{rh}{2})).

(8.18)

Συνεπώς

\kappa=\sigma\pm\sqrt{\sigma^{2}-1}.

Εύκολα βλέπουμε ότι αν $\sigma^{2}>1$ , τότε υπάρχει μια ρίζα της (8.18) με $|\kappa|>1$ , και αν $\sigma^{2}\leq 1$ τότε ότι $|\kappa|^{2}=1$ . Συνεπώς, για να ισχύει ότι $\sigma^{2}\leq 1$ , αρκεί $\lambda\leq 1$ . Άρα, η (8.7) είναι von Neumann ευσταθής, αν $\lambda\leq 1$ .