v2.3 – Finite Differenzen für zweidimensionale Gebiete

 

2.3.1 Poissongleichung im Einheitsquadrat

Wir wollen hier eine einfache Beispielaufgabe besprechen:

-\Delta u = f\quad in\:\:\:\Omega = {\left( {0,1} \right)^2},\quad u = 0\quad auf\:\:\partial \Omega

2.3.1.1 Definition des Gitters

Das Differenzenverfahren definiert man wie im eindimensionalen Fall. Wir wählen als Schrittweite h = \frac{1}{N} und definieren die Gitter:

{\Omega _h} = \left\{ {\left( {{x_i},{y_j}} \right):\left( {\frac{i}{N},\frac{j}{N}} \right),\:\:0 < i,j < N} \right\}

{\Gamma _h} = \left\{ {\left( {{x_i},{y_j}} \right):\left( {\frac{i}{N},\frac{j}{N}} \right),\:\:i \in \left\{ {0,N} \right\},j = 0, \ldots ,N\:\:\: \vee \:\:\:j \in \left\{ {0,N} \right\},i = 0, \ldots ,N} \right\}

num-203-gitter-gebiet-zweidimensional-diskretisierung

2.3.1.2 Aufstellen des Gleichungssystems

Bezeichne {u_{i,j}} = {u_h}\left( {{x_i},{y_j}} \right) den gesuchten Näherungswert für u\left( {{x_i},{y_j}} \right):

num-203-aufstellen-gleichungssystem-gitter

dann wird

\frac{{{\partial ^2}u}}{{\partial {x^2}}} durch \frac{1}{{{h^2}}}\left[ {{u_{i+1,j}}-2{u_{i,j}}+{u_{i-1,j}}} \right]

\frac{{{\partial ^2}u}}{{\partial {y^2}}} durch \frac{1}{{{h^2}}}\left[ {{u_{i,j+1}}-2{u_{i,j}}+{u_{i,j-1}}} \right]

approximiert. Somit lautet das Differenzenverfahren:

-\left( {{\Delta _h}{u_h}} \right): = -\frac{1}{{{h^2}}}\left[ {{u_{i+1,j}}+{u_{i-1,j}}+{u_{i,j+1}}+{u_{i,j-1}}-4{u_{i,j}}} \right] = {f_{i,j}}

für alle \left( {{x_i},{y_j}} \right) \in {\Omega _h}, wobei {f_{i,j}} = {f_h}\left( {{x_i},{y_j}} \right) = \left( {{R_h}f} \right)\left( {{x_i},{y_j}} \right) = f\left( {{x_i},{y_j}} \right) gilt. Für den ersten Punkt {u_{1,1}} ergibt sich zum Beispiel:

{f_{1,1}} = -\frac{1}{{{h^2}}}\left[ {{u_{2,1}}+{u_{0,1}}+{u_{1,2}}+{u_{1,0}}-4{u_{1,1}}} \right]

Es gehen immer die benachbarten Knoten in die Gleichung ein:

num-203-funf-punkte-stern-zweidimensional-gitter

Um diese Gleichungen in ein Gleichungssystem zu schreiben, müssen die Gitterpunkte als Vektor dargestellt werden. Da sie zweidimensional angeordnet sind, müssen sie z.B. zeilenweise nummeriert werden:

num-203-knoten-nummerierung-rand-gebiet

Die Randpunkte müssen in der Nummerierung nur enthalten sein, wenn die Lösung auf dem Rand ungleich 0 ist. Wir haben in unserem Beispiel homogene Dirichletrandbedingungen, daher sind alle Randknoten gleich 0 und fallen im Gleichungssystem weg. Wir nummerieren nur die inneren Knoten. Bei der angeführten Nummerierungsmethode ergibt sich für N = 5 folgendes Gleichungssystem:

\frac{1}{{{h^2}}}\left( {\begin{array}{*{20}{c}}  4 & {-1} & {} & {} &\vline & {-1} & {} & {} & {} &\vline & {} & {} & {} & {} &\vline & {} & {} & {} & {} \\{-1} & 4 & {-1} & {} &\vline & {} & {-1} & {} & {} &\vline & {} & {} & {} & {} &\vline & {} & {} & {} & {} \\{} & {-1} & 4 & {-1} &\vline & {} & {} & {-1} & {} &\vline & {} & {} & {} & {} &\vline & {} & {} & {} & {} \\{} & {} & {-1} & 4 &\vline & {} & {} & {} & {-1} &\vline & {} & {} & {} & {} &\vline & {} & {} & {} & {} \\ \hline{-1} & {} & {} & {} &\vline & 4 & {-1} & {} & {} &\vline & {-1} & {} & {} & {} &\vline & {} & {} & {} & {} \\{} & {-1} & {} & {} &\vline & {-1} & 4 & {-1} & {} &\vline & {} & {-1} & {} & {} &\vline & {} & {} & {} & {} \\{} & {} & {-1} & {} &\vline & {} & {-1} & 4 & {-1} &\vline & {} & {} & {-1} & {} &\vline & {} & {} & {} & {} \\{} & {} & {} & {-1} &\vline & {} & {} & {-1} & 4 &\vline & {} & {} & {} & {-1} &\vline & {} & {} & {} & {} \\ \hline{} & {} & {} & {} &\vline & {-1} & {} & {} & {} &\vline & 4 & {-1} & {} & {} &\vline & {-1} & {} & {} & {} \\{} & {} & {} & {} &\vline & {} & {-1} & {} & {} &\vline & {-1} & 4 & {-1} & {} &\vline & {} & {-1} & {} & {} \\{} & {} & {} & {} &\vline & {} & {} & {-1} & {} &\vline & {} & {-1} & 4 & {-1} &\vline & {} & {} & {-1} & {} \\{} & {} & {} & {} &\vline & {} & {} & {} & {-1} &\vline & {} & {} & {-1} & 4 &\vline & {} & {} & {} & {-1} \\ \hline{} & {} & {} & {} &\vline & {} & {} & {} & {} &\vline & {-1} & {} & {} & {} &\vline & 4 & {-1} & {} & {} \\{} & {} & {} & {} &\vline & {} & {} & {} & {} &\vline & {} & {-1} & {} & {} &\vline & {-1} & 4 & {-1} & {} \\{} & {} & {} & {} &\vline & {} & {} & {} & {} &\vline & {} & {} & {-1} & {} &\vline & {} & {-1} & 4 & {-1} \\{} & {} & {} & {} &\vline & {} & {} & {} & {} &\vline & {} & {} & {} & {-1} &\vline & {} & {} & {-1} & 4 \\   \end{array} } \right)\left( {\begin{array}{*{20}{c}}{{u_1}} \\{{u_2}} \\{{u_3}} \\{{u_4}} \\ \hline{{u_5}} \\{{u_6}} \\{{u_7}} \\{{u_8}} \\ \hline{{u_9}} \\{{u_{10}}} \\{{u_{11}}} \\{{u_{12}}} \\ \hline{{u_{13}}} \\{{u_{14}}} \\{{u_{15}}} \\{{u_{16}}} \\   \end{array} } \right) = \left( {\begin{array}{*{20}{c}}{{f_1}} \\{{f_2}} \\{{f_3}} \\{{f_4}} \\ \hline{{f_5}} \\{{f_6}} \\{{f_7}} \\{{f_8}} \\ \hline{{f_9}} \\{{f_{10}}} \\{{f_{11}}} \\{{f_{12}}} \\ \hline{{f_{13}}} \\{{f_{14}}} \\{{f_{15}}} \\{{f_{16}}} \\   \end{array} } \right)

2.3.1.3 Eigenschaften und Lösung des Gleichungssystems

Die Steifigkeitsmatrix

  • ist großdimensioniert (Anzahl der Gleichungen ist O\left( {{h^{-2}}} \right))
  • ist schwach besetzt (maximal 5 Nicht-Nullelemente pro Zeile
  • hat eine Bandweite O\left( N \right) = O\left( {{h^{-1}}} \right)
  • ist symmetrisch und positiv definit (weil die Hauptdiagonalelemente größer als die Summe der restlichen Elemente einer Zeile sind)
  • ist schlecht konditioniert (\kappa \left( {{K_h}} \right) = \left\| {{K_h}} \right\|\left\| {K_h^{-1}} \right\| = O\left( {{h^{-2}}} \right))

Große Dimension, schwache Besetztheit und schlechte Kondition sind typische Eigenschaften von Matrizen, die bei der Diskretisierung von Randwertproblemen entstehen. Die Symmetrie folgt aus den Eigenschaften des Differentialoperators. Wäre ein Term erster Ordnung vorhanden, wäre die Matrix nicht symmetrisch. Die Bandweite folgt aus der Knotennummerierung.

Gleichungssysteme mit den aufgelisteten Eigenschaften löst man am besten mit

  • CG-Verfahren mit guter Vorkonditionierung
  • Mehrgitterverfahren
  • FFT (Fast Fourier Transform)

2.3.1.4 Nicht-rechteckige Gebiete

Nicht-rechteckige Gebiete können auch behandelt werden. Am Rand muss das Gitter entsprechend angepasst werden:

num-203-nicht-rechteckiges-gebiet-rand-diskretisierung

Trotz Konsistenzordnung 1 in randnahen Punkten ergibt sich globale Konvergenzordnung 2.