9 Ελλειπτική εξίσωση: Πεπερασμένες διαφορές

9.1 Ελλειπτική εξίσωση

Θεωρούμε το ακόλουθο πρόβλημα: Ζητείται μια συνάρτηση u:Ω=(0,1)×(0,1), τέτοια ώστε

-(uxx(x,y)+uyy(x,y))+q(x,y)u(x,y)=f(x,y),(x,y)Ωu(x,y)=0,(x,y)Ω (9.1)

όπου q,fC(Ω) και q(x,y)0, για κάθε (x,y)Ω. Στη συνέχεια, θα συμβολίζουμε με qmin=min(x,y)Ωq(x,y).

Θεωρούμε και πάλι έναν φυσικό αριθμό N και διαμερίζουμε το διάστημα [0,1] σε N+2 ισαπέχοντα σημεία 0=x0<x1<<xN<xN+1=1, όπου h=xi+1-xi, i=0,,N. Επίσης θέτουμε yi=xi, i=0,,N+1 και τότε τα σημεία (xi,yj), i,j=0,,N+1 δημιουργούν έναν διαμερισμό του Ω. Τότε, σε κάθε σημείο του διαμερισμού (xi,yj), i,j=1,,N, θα ισχύει:

-(uxx(xi,yj)+uyy(xi,yj))+q(xi,yj)u(xi,yj)=f(xi,yj), (9.2)

για i,j=1,,N. Θα κατασκευάσουμε όπως και στα προηγούμενα κεφάλαια προσεγγίσεις Uij των τιμών u(xi,yj), i,j=0,,N+1. Λόγω των ομογενών συνοριακών συνθηκών Dirichlet, θέτουμε U0j=UN+1j=Ui0=UiN+1=0, i,j=0,,N+1. Για την κατασκευή των υπολοίπων Uij θα προσεγγίσουμε τις uxx(xi,yj) και uyy(xi,yj) στην (9.2) με την δh,2cu(xi,yj), θεωρώντας την πεπερασμένη διαφορά τη μία φορά ως προς τη μεταβλητή x και την αλλή ως προς τη μεταβλητή y, βλ. (2.8). Στο Σχήμα 9.1, βλέπουμε ένα παράδειγμα πλέγματος σημείων, όπου σημειώνουμε τα σημεία των άγνωστων ή γνωστών τιμών της λύσης u του (9.1). Έτσι, αν υποθέσουμε ότι 4ux4,4uy4C(Ω), λόγω της (2.11), η (9.2) γίνεται,

-(u(xi+1,yj)-2u(xi,yj)+u(xi-1,yj)h2+u(xi,yj+1)-2u(xi,yj)+u(xi,yj-1)h2)+q(xi,yj)u(xi,yj)=f(xi,yj)+ηij,i,j=1,,N (9.3)

όπου

|ηij|h212(max(x,y)Ω|4ux4u(x,y)|+max(x,y)Ω|4uy4u(x,y)|). (9.4)
Σχήμα 9.1: Τα σημεία (xi,yj) του χωρίου Ω. Σημειώνουμε με άσπρο κύκλο τα σημεία όπου αναζητούμε τις τιμές Uij και με έντονο μαύρο τα σημεία που αντιστοιχούν στις γνωστές συνοριακές τιμές.

Θεωρούμε λοιπόν τις προσεγγίσεις Uij της u(xi,yj), i,j=1,,N, οι οποίες θα ικανοποιούν τις ακόλουθες εξισώσεις

-(Ui+1j-2Uij+Ui-1jh2+Uij+1-2Uij+Uij-1h2)+q(xi,yj)Uij=f(xi,yj),i,j=1,,N (9.5)

Στο Σχήμα 9.2 φαίνεται ένα παράδειγμα πλέγματος όπου σημειώνουμε τις άγνωστες τιμές της προσέγγισης σύμφωνα με τη μέθοδο (9.5).

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

Επομένως, αν συμβολίσουμε με UN2 το διάνυσμα με συνιστώσες U11,, UN1,U12,, UN2,, U1N,, UNN, U=(U11,,UN1,U12,,UN2,,U1N,,UNN)T, μπορούμε να γράψουμε τις εξισώσεις (9.5) ως γραμμικό σύστημα

AU=h2F,

όπου F=(f(x1,y1),,f(xN,y1),f(x1,y2),,f(xN,y2), , f(x1,yN), , f(xN,yN))T.

Για να κατανοήσουμε τη δομή του πίνακα A, θεωρούμε την περίπτωση q=0 και παρατηρούμε ότι για i=1,,N και j=1 οι εξισώσεις (9.5) γίνονται

-(Ui+11-4Ui1+Ui-11+Uij+1)=h2f(xi,yj),i=1,,N. (9.6)

Επομένως, οι N πρώτες γραμμές του πίνακα Α είναι

\begin{equation*}%\label{A-matrix-partI} \begin{aligned} &\left( \begin{array}{cccccc} 4 & -1 & 0 & &0 &|\\ -1& 4 & -1 & 0 &0 &|\\ 0 &\ddots&\ddots&\ddots&0 &|\\ & &-1 & 4 &-1&|\\ & & & -1 & 4&|\\ \end{array}\right.\!\!\!\!\! &&\begin{array}{cccccc} 1 &0&0 & &0 &|\\ 0 &1&0 & &0 &|\\ 0 & &\ddots& &0 &|\\ 0 & & 0 &1&0 &|\\ 0 & & &0&1 &|\\ \end{array}\!\!\!\!\! &&\left.\begin{array}{ccc} 0& &0\\ 0&\dots &0\\ 0&\dots &0\\ 0& &0\\ 0& &0\\ \end{array} \right)\\ &\hspace{1cm}\underbrace{\hspace{3cm}} &&\hspace{0.1cm}\underbrace{\hspace{3cm}} &&\hspace{0.1cm}\underbrace{\hspace{1.5cm}}\\ &\hspace{2cm}\text{$N$ στήλες} &&\hspace{1cm}\text{$N$ στήλες} &&\text{$N^2-2N$ στήλες} \end{aligned}\end{equation*}

Συμβολίζουμε με B τον N×N πίνακα

\[B= \left( \begin{array}{rrrrr} 4 & -1 & 0 & &0 \\ -1& 4 & -1 & 0 &0 \\ 0 &\ddots&\ddots&\ddots&0 \\ & &-1 & 4 &-1\\ & & & -1 & 4\\ \end{array}\right) \]

και με \(\mathcal{I}\) τον μοναδιαίο N×N πίνακα, και λόγω των εξισώσεων (9.5) για i=1,,N και j=2, οι N+1 έως 2N γραμμές του πίνακα Α είναι

\begin{equation*} \begin{aligned} \left( \begin{array}{cccccccccc} & & &|& & & |& & &|\\ & \mathcal{I}& & |& B& &|& \mathcal{I} & &|\\ & & & |& & &|& & &|\\ % \multicolumn{9}{c}{}& &\multicolumn{3}{c}{\myunderbrace{}{\text{$N^2-3N$ στήλες}}} \end{array} \right. &\left. \begin{array}{ccc} 0&0&0\\ 0&0&0\\ 0&0&0\\ % \multicolumn{9}{c}{}& &\multicolumn{3}{c}{\myunderbrace{}{\text{$N^2-3N$ στήλες}}} \end{array} \right)\\ &\hspace{.1cm}\underbrace{\hspace{1.5cm}}\\ &\hspace{.1cm}{\text{$N^2-3N$ στήλες}} \end{aligned}\end{equation*}

Οπότε, ο πίνακας A έχει τη μορφή

\begin{equation}\label{9.9} A= \left( \begin{array}{rrrrr} B & \mathcal{I} & 0 & &0 \\ \mathcal{I}& B & -1 & 0 &0 \\ 0 &\ddots&\ddots&\ddots&0 \\ & &\mathcal{I} & B &\mathcal{I}\\ & & & \mathcal{I} & B\\ \end{array}\right)\tag{9.9} \end{equation}

ο οποίος δεν έχει αυστηρά κυριαρχική διαγώνιο, είναι όμως αντιστρέψιμος. Αν υποθέσουμε τώρα ότι q0 και ότι qmin0, τότε το γραμμικό σύστημα που προκύπτει από τις (9.5) είναι

(A+h2Q)U=h2F, (9.10)

A,U και F όπως παραπάνω και Q ένας διαγώνιος N2×N2 πίνακας με στοιχεία στη διαγώνιο q(x1,y1), , q(xN,y1), q(x1,y2), , q(xN,y2), , q(x1,yN), , q(xN,yN). Είναι απλό να δούμε ότι ο πίνακας A+h2Q έχει όμοια μορφή με αυτή της (\ref{9.9}), διότι αλλάζουν μόνο τα διαγώνια στοιχεία και έχει αυστηρά κυριαρχική διαγώνιο, αν qmin>0.

Επομένως, μπορούμε να επιλέξουμε μια άμεση μέθοδο, όπως είναι αυτή της απαλοιφής Gauss, για να λύσουμε το γραμμικό σύστημα (9.10), βλ. π.χ. (Ακρίβης και Δουγαλής, (2015)). Ο πίνακας του γραμμικού συστήματος (9.10) είναι μεγάλος σε διάσταση και είναι και αραιός, έχει δηλαδή πολλά μηδενικά στοιχεία. Λόγω της (9.6), μπορούμε να δούμε ότι κάθε γραμμή του πίνακα A έχει το πολύ 5 μη-μηδενικά στοιχεία. Στην πράξη, δεν επιλέγουμε τη μέθοδο απαλοιφής Gauss για την αριθμητική επίλυση μεγάλων αραιών γραμμικών συστημάτων, αλλά άλλες μεθόδους όπως είναι οι επαναληπτικές μέθοδοι, Jacobi, Gauss-Seidel, SOR και άλλες, βλ. π.χ. (Ακρίβης και Δουγαλής, (2015), Strikwerda, (2004)).