Finite-Differenzen-Methode

Finite-Differenzen-Methoden (kurz: FDM), auch Methoden der endlichen (finiten) Differenzen sind eine Klasse numerischer Verfahren zur Lösung gewöhnlicher und partieller Differentialgleichungen.

Die grundlegende Idee des Verfahrens ist es, die Ortsableitungen in der Differenzialgleichung an endlich vielen (= „finiten“), äquidistanten Gitterpunkten durch Differenzenquotienten zu approximieren. Die approximierten Lösungen der Differenzialgleichung an den Gitterpunkten lassen sich dann durch das entsprechende Gleichungssystem berechnen.

Verfahren dieser Art finden verbreitete Anwendung unter anderem bei fluiddynamischen Simulationen, zum Beispiel in der Meteorologie und der Astrophysik. In den Jahren von 1950 bis 1980 dominierte die FDM in den numerischen Programmen der Reaktorphysik zur Berechnung des Neutronenflusses. Eine gewisse Verbreitung findet das Differenzenverfahren in der Baustatik. Schon 1904 analysierte Friedrich Bleich den Durchlaufträger; 1909 untersuchte Lewis Fry Richardson elastische Scheiben und 1919 Henri Marcus elastische Platten mit dem Differenzenverfahren.

Eine spezielle Finite-Differenzen-Methode zur numerischen Lösung der Wärmeleitungsgleichung ist das Crank-Nicolson-Verfahren.

Zu den Pionieren des Finite-Differenzen-Verfahrens für partielle Differentialgleichungen zählen Lewis Fry Richardson, Richard Southwell, Richard Courant, Kurt Friedrichs, Hans Lewy, Peter Lax und John von Neumann.

Beispiel zur numerischen Lösung einer gewöhnlichen DGL

Gegeben sei das Randwertproblem

$u''(x)=2$  für  $\quad x\in [0;1]$,

$u(0)=u(1)=3$.

Die Lösungsfunktion $u\colon [0;1]\to \mathbb {R}$ lässt sich hier exakt berechnen zu $u(x)=3+x(x-1)$.

Zur Lösung mit der Differenzenmethode wird das Intervall $[0;1]$ diskretisiert durch die Gitterpunkte $x_{i}=i\cdot h$ für $i=0,\ldots ,n+1$ mit der Maschenweite $h={\tfrac {1}{n+1}}$. Die Diskretisierung der zweiten Ableitung erfolgt mit den zentralen Differenzenquotienten der zweiten Ableitung

$u''(x)\approx {\frac {u(x-h)-2u(x)+u(x+h)}{h^{2}}}\ .$

Dies ergibt an den inneren Gitterpunkten die Differenzengleichungen

${\frac {1}{h^{2}}}(u_{i-1}-2u_{i}+u_{i+1})=2$  für  $i=1,\ldots ,n$

für die numerischen Näherungswerte $u_{i}$ der Lösungswerte $u(x_{i})$. Unter Verwendung der gegebenen Randwerte $u_{0}=u(0)=3$ und $u_{n+1}=u(1)=3$ ist dies ein lineares Gleichungssystem mit $n$ Gleichungen für die $n$ Unbekannten $u_{1},\ldots ,u_{n}$.

In Matrixform lautet das zu lösende System hier:

${\frac {1}{h^{2}}}{\begin{pmatrix}-2&1&0&\dots &0\\1&-2&1&\ddots &\vdots \\0&1&\ddots &\ddots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\dots &0&1&-2\end{pmatrix}}{\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\\\vdots \\u_{n}\end{pmatrix}}={\begin{pmatrix}2-{\frac {3}{h^{2}}}\\2\\2\\\vdots \\2-{\frac {3}{h^{2}}}\end{pmatrix}}$

Da in jeder Zeile maximal nur drei Unbekannte vorkommen, handelt es sich um ein System mit dünnbesetzter Koeffizientenmatrix, genauer um ein System mit Tridiagonal-Toeplitz-Matrix.

Beispiel zur numerischen Lösung einer partiellen DGL

Im Folgenden wird die numerische Lösung der Wärmeleitungsgleichung auf einem beschränkten Gebiet $\Omega$ betrachtet:

$\partial _{t}u-\Delta u=f~~{\text{in }}(0,T)\times \Omega ,~f\in C^{0}((0,T)\times \Omega ),$

$u(t,\cdot )=0~{\text{ auf }}\partial \Omega ,$

$u(0,\cdot )=u_{0}{\text{ in }}\Omega .$

Numerische Lösung im 1D

Im 1D-Fall ist $\Omega =(a,b)$ ein beschränktes Intervall. Da in diesem Fall nur eine Ortsableitung betrachtet wird, kann die Wärmeleitungsgleichung folgendermaßen geschrieben werden:

$\partial _{t}u-\partial _{xx}u=f$

Diskretisierung

Um die Finite-Differenzen-Methode anwenden zu können, muss das Intervall $\Omega$ zunächst in endlich viele Teilintervalle unterteilt werden. Hierfür werden $N+2$ äquidistante Stützstellen verwendet:

$x_{0}=a,~x_{N+1}=b,~x_{i}=a+i\cdot {\frac {b-a}{N+1}}$, für $i=0,\ldots ,N+1$.

Die Gitterweite dieser Diskretisierung ist also $h={\tfrac {b-a}{N+1}}$. Nach Voraussetzung verschwindet die gesuchte Funktion $u$ an den Randwerten, d. h. $u(x_{0})=u(x_{N+1})=0$, sodass diese Werte nicht weiter betrachtet werden müssen. Damit lassen sich die Funktionsauswertungen von $u$ an Stützstellen als Vektor im $\mathbb {R} ^{N}$ darstellen:

$U_{h}=(u_{1}(t),\dots ,u_{N}(t))^{T}$

Approximation der Ableitung

Die zweite Ableitung von $u$ bzgl. des Orts kann nun an den Stützstellen durch Differenzenquotienten zweiter Ordnung approximiert werden:

$\partial _{xx}u(t,x_{i})\approx {\frac {u(t,x_{i+1})-2u(t,x_{i})+u(t,x_{i-1})}{h^{2}}}$

Wird die Wärmeleitungsgleichung nach $\partial _{t}u$ umgestellt, ergibt sich damit folgendes System gewöhnlicher Differenzialgleichungen erster Ordnung:

${\dot {U}}_{h}(t)=-A_{h}U_{h}+f_{h}$

wobei $f_{h}:=(f(t,x_{1}),\dots ,f(t,x_{N}))^{T}$ und $A_{h}:={\frac {1}{h^{2}}}{\begin{pmatrix}2&-1&0&\dots &0\\-1&2&-1&\dots &0\\0&\ddots &\ddots &\ddots &\vdots \\0&\dots &\dots &-1&2\end{pmatrix}}$.

Dieses System kann nun durch beliebige Verfahren zur Lösung von gewöhnlichen Differenzialgleichungen, wie z. B. das Runge-Kutta-Verfahren oder das Euler-Verfahren, gelöst werden. 

Güte der Approximation

Eine Finite-Differenzen-Methode erzeugt ein lineares Gleichungssystem (analog Gleichung im Kapitel Beispiel)

$A_{h}u_{h}=b_{h},$

wobei $u_{h}$ die numerische Approximation der Lösung ist und $h$ die Abhängigkeit vom Gitter explizit darstellen soll. Sei $u(x)$ die exakte Lösung und $U\in \mathbb {R} ^{N}$ die endliche Darstellung mittels $U_{i}=u(x_{i})$. Sowohl $u_{h}$ als auch $U$ sind sogenannte Gitterfunktionen, sie sind nur definiert in allen Gitterpunkten des verwendeten Gitters.

Eine FDM heißt konsistent von Ordnung $n\in \mathbb {N}$, falls es ein $C_{K}>0$ gibt mit

$\|A_{h}U-b_{h}\|\leq C_{K}h^{n}.$

Dabei ist $\|\cdot \|$ eine Norm für Gitterfunktionen. Oft verwendet wird die Maximumnorm $\|v\|:=\max _{i}|v(x_{i})|$, sie erlaubt punktweise Fehlerabschätzungen.

Eine FDM heißt stabil, falls es ein $C_{S}>0$ gibt, sodass für alle Gitterfunktionen $v\in \mathbb {R} ^{N}$ gilt

$\|A_{h}v\|\geq C_{S}\|v\|.$

Aus Konsistenz der Ordnung $n$ und Stabilität folgt Konvergenz der Ordnung $n$

$\|U-u_{h}\|\leq {\frac {C_{K}}{C_{S}}}h^{n}.$

Stabilität vorausgesetzt, ist Konsistenz der Ordnung $n$ für Konvergenz der Ordnung $n$ hinreichend, aber nicht notwendig.

Eine Diskretisierung auf einem nichtäquidistanten Gitter

Nichtäquidistante Gitter sind angebracht, wenn die Lösung Besonderheiten aufweist, z. B. Singularitäten oder Grenzschichten. Als Beispiel wird betrachtet

$-u''(x)=f(x),\quad u(0)=u(1)=0$

und eine Diskretisierung auf einem beliebigen Gitter $0=x_{0}<x_{1}<\cdots x_{i}<\cdots <x_{N}=1$ mit $h_{i}=x_{i}-x_{i-1}$ und $h=\max _{i}h_{i}$.

Dann ist folgende Diskretisierung zweckmäßig:

${\frac {2}{h_{k}+h_{k+1}}}({\frac {u_{k}-u_{k+1}}{h_{k+1}}}+{\frac {u_{k}-u_{k-1}}{h_{k}}})=f_{k},\quad u_{0}=u_{N}=0.$

Untersucht man den Konsistenzfehler (wie üblich durch Taylorentwicklung), so stellt man folgende Struktur des Konsistenzfehlers fest:

${\frac {1}{3}}u'''(x_{k})(h_{k+1}-h_{k})+O(h^{2}).$

Auf äquidistanten Gittern verschwindet der erste Anteil und der Konsistenzfehler ist von der Ordnung 2. Auf beliebigen Gittern jedoch reduziert sich der Konsistenzfehler auf die Ordnung Eins! Trotzdem ist der Fehler in der Maximumnorm wieder von der Ordnung Zwei, das kann man mit Hilfe der diskreten Greenschen Funktion oder verbesserten Stabilitätsaussagen verifizieren.

Lösen von Randwertproblemen

Ein Beispiel für ein solches Randwertproblem ist folgende partielle Differentialgleichung zweiter Ordnung:

$A(x,y)\cdot {\frac {\partial ^{2}f(x,y)}{\partial x^{2}}}+B(x,y)\cdot {\frac {\partial ^{2}f(x,y)}{\partial x\partial y}}+C(x,y)\cdot {\frac {\partial ^{2}f(x,y)}{\partial y^{2}}}+D(x,y)\cdot {\frac {\partial ^{2}f(x,y)}{\partial x}}+E(x,y)\cdot {\frac {\partial ^{2}f(x,y)}{\partial y}}+F(x,y)\cdot f(x,y)=0$

Die unbekannte Funktion $f(x,y)$ soll an vorgegebenen Punkten des zweidimensionalen Bereichs bestimmte feste Werte (Randbedingungen) annehmen. Die Funktionen $A(x,y)$, $B(x,y)$, $C(x,y)$, $D(x,y)$, $E(x,y)$ und $F(x,y)$ sind Koeffizienten und sollten daher beschränkt und im gegebenen Integrationsbereich stetig sein. Bei vielen Gleichungen, die spezielle technische Probleme beschreiben, nehmen diese Funktionen konstante Werte an. Diese Gleichungen vom werden oft nach dem Wert der Diskriminante $\Delta (x,y)=B^{2}(x,y)-4\cdot A(x,y)\cdot C(x,y)$ in eine der folgenden Gruppen eingeteilt: hyperbolisch für $\Delta (x,y)>0$, parabolisch für $\Delta (x,y)=0$ und elliptisch für $\Delta (x,y)<0$.

Beispiele für partielle Differentialgleichungen, die für physikalische Probleme formuliert werden, sind:

• Laplace-Gleichung:

$\nabla ^{2}U(x,y,z)=0$

• Poisson-Gleichung:

$\nabla ^{2}U(x,y,z)+f(x,y,z)=0$

• Helmholtz-Gleichung:

$\nabla ^{2}U(x,y,z)+k^{2}U(x,y,z)=0$

• Wellengleichung:

$\nabla ^{2}U(x,y,z)-a^{2}\cdot {\frac {\partial ^{2}U(x,y,z,t)}{\partial t^{2}}}=0$

• Wärmeleitungsgleichung:

$\nabla ^{2}U(x,y,z)-b^{2}\cdot {\frac {\partial U(x,y,z,t)}{\partial t}}=0$

Literatur

• Christian Großmann, Hans-Görg Roos: Numerische Behandlung partieller Differentialgleichungen. 3. Auflage. B. G. Teubner Verlag, Wiesbaden 2005, ISBN 3-519-22089-X.
• Stig Larsson, Vidar Thomée: Partielle Differentialgleichungen und numerische Methoden. Springer-Verlag, Berlin 2005, ISBN 3-540-20823-2.
• Claus-Dieter Munz, Thomas Westermann: Numerische Behandlung gewöhnlicher und partieller Differenzialgleichungen. 3. Auflage. Springer-Verlag, ISBN 978-3-642-24334-9

Weblinks

Commons: Finite-Differenzen-Methode – Sammlung von Bildern, Videos und Audiodateien