Finite-Differenzen-Methode

Perspektivische Grafik der numerischen Lösung einer partiellen DGL in 2D nach der FDM in Zylinderkoordinaten für ein physikalisches Problem

Finite-Differenzen-Methoden (FDM), auch Finite-Differenzen-Verfahren, Methoden/Verfahren der endlichen Differenzen oder schlicht Differenzenverfahren, sind eine Klasse numerischer Verfahren zur Lösung gewöhnlicher und partieller Differentialgleichungen (DGL).

Die grundlegende Idee des Verfahrens ist es, die Orts- und/oder Zeitableitungen in der Differentialgleichung in einem vorgegebenen Intervall der unabhängigen Variablen an Gitterpunkten in diesem Intervall durch Differenzenquotienten zu approximieren. Diese Approximationen der Differenzialgleichung in den Gitterpunkten stellen dann ein Gleichungssystem dar, welches zu lösen ist.

Die erste Grafik illustriert die numerische Lösung $u(r,z)$ einer partiellen DGL in 2D nach der FDM in einem krummlinigen $(r,z)$-Koordinatensystem (Zylinderkoordinaten) mit 21 Gitterpunkten (20 Intervallen) in der r-Richtung und 11 Gitterpunkten (10 Intervallen) in der z-Richtung. Der diskretisierte Körper ist ein Zylinder mit einer Höhe von 100 cm und mit einem Radius von 200 cm. Es handelt sich eigentlich um einen (realen) 3D-Zylinder, der Besonderheiten besitzt. Die in die Berechnung eingehende Materialverteilung im Zylinder kann sich in Prinzip von Zelle zu Zelle ändern. In dem vorliegenden Modell hängt jedoch die Materialverteilung nicht von der polaren Koordinate, der $\varphi$-Koordinate, ab. Diese Eigenschaft muss prinzipiell auch jede physikalische Lösungsfunktion besitzen, die eine Observable ist.

Somit wird das Problem zu einem mathematischen 2D-, einem $(r,z)$-Problem. Da die Materialverteilung außerdem in z-Richtung eine Spiegelsymmetrieebene besitzt, wurde für $z=100$ als eine der Randbedingungen der Wert $0$ für den Gradienten vorgegeben und damit die Hälfte der Gitterpunkte eingespart. Der Zylinder ist also eigentlich 200 cm hoch. Als weitere Randbedingung wurde vorgegeben, dass die gesuchte Lösung an den Außenrändern verschwinden soll. Dass die Näherungslösung dem gehorcht, ist auf der Grafik zu erkennen. Es handelt sich bei diesem Zylindermodell um ein zwar einfaches, aber reales physikalisches Problem aus der Neutronendiffusionstheorie für ein spezielles Kernreaktormodell eines BN-Reaktors. Dargestellt wird der schnelle Neutronenfluss.

In den Jahren von 1950 bis 1980 dominierte die FDM in den numerischen Programmen der Reaktorphysik zur Berechnung des Neutronenflusses. Verfahren dieser Art finden verbreitete Anwendung unter anderem bei fluiddynamischen Simulationen, zum Beispiel in der Meteorologie und der Astrophysik. 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.

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.

Geschichte des Namens

In keinem deutsch- oder englischsprachigen Lehrbuch bis etwa zum Jahr 1980 wird man die Termini Finite-Differenzen-Methode bzw. finite difference method finden. Zum Beispiel fehlt dieses Stichwort im Lexikon Mathematik von 1977, wo aber der Differenzenquotient enthalten ist. Man sprach in der englischsprachigen Fachliteratur schon früher von finite-difference calculations, finite difference equations oder finite difference schemes, in der deutschsprachigen von Differenzenrechnungen, Differenzengleichungen und Diskretisierungsschemata. In der deutschsprachigen Fachliteratur tauchten Formen des englischen Adjektivs finite (begrenzt, endlich) in diesem Zusammenhang nicht auf. Streng genommen ist der Wortbestandteil finit überflüssig, denn Differenzen sind immer endlich. Sie sind das Gegenstück zum Differential, das einen unendlich kleinen Abschnitt auf der Achse eines Koordinatensystems beschreibt. Auch handelt es sich nicht um eine spezielle Methode, sondern umfasst eine Vielzahl von Verfahren.

Finite difference method tritt erst in den 1980er Jahren auf, und zwar nach dem „Siegeszug“ der finite element method. Um 1960 wurde der Name finite element method für eine verwandte Methode eingeführt und mit Finite-Elemente-Methode eingedeutscht. Die wahrscheinlich erste Verwendung des Begriffs finite element findet sich 1960 in einem Zeitschriftenaufsatz von Ray W. Clough, einem Professor für Baustatik. Diese Begriffe verbreiteten sich in den 1970er Jahren schnell. Erst deutlich später muss den Mathematikern, Physikern und Ingenieuren aufgefallen sein, dass für die klassischen Diskretisierungsverfahren kein „griffiger“ Name existiert. Sie fanden, angelehnt an den Namen finite element method, dass der Name finite difference method recht passend sei. Ein Standardlehrbuch zu dieser Materie zeigt den Entwicklungsgang: Es trug 1965 noch den Titel Numerical solution of partial differential equations: With exercises and worked solutions. Die Monographie wurde 1970 ins Deutsche übertragen. Eine spätere, erweiterte Auflage der englischen Version hatte 1985 das finite difference methods dann schon im Untertitel. Diese Auflage wurde wahrscheinlich nicht ins Deutsche übersetzt, so dass der Nachweis eines ersten Auftretens von Finite-Differenzen-Methode so nicht erbracht werden kann.

Es dauerte noch eine Weile, bis finite difference method in Finite-Differenzen-Methode eingedeutscht wurde. Man beachte einen kleinen Unterschied. Hätte man den englischen Namen direkt eingedeutscht, müsste es eigentlich Finite-Differenz-Methode heißen. In den 1970er, 1980er Jahren scheint die Beschäftigung mit numerischen Methoden, zumindest in den beiden deutschen Staaten, für Mathematiker kein attraktives Betätigungsfeld gewesen zu sein, wie man durch Recherchen in den Katalogen der deutschsprachigen Bibliotheksverbünde nach entsprechenden Fachbüchern nachweisen kann. Die Anwendung und Weiterentwicklung der Differenzenrechnung verlagerte sich auf Physiker und Ingenieure. Wer sich damit beschäftigte, bezog seine Fachliteratur ohnehin aus der englisch- und russischsprachigen Fachliteratur. In der russischsprachigen verwendet man heute Метод конечных разностей, also Methode der endlichen Differenzen. Eine Recherche in den Katalogen der deutschsprachigen Bibliotheksverbünde ergab, dass wahrscheinlich das erste Fachbuch in deutscher Sprache mit Finite-Differenzen-Methode im Titel im Jahr 1987 erschienen ist, das von Ingenieuren und nicht von Mathematikern verfasst wurde. Es wird im Katalog der deutschen Nationalbibliothek nach den Sachgruppen Bergbau, Bautechnik, Umwelttechnik klassifiziert und ist in vier Auflagen erschienen.

Der Name klassisches Differenzenverfahren wäre sicher treffender gewesen. Klassisch deshalb, weil die Methodik bereits seit Anfang des 19. Jahrhunderts bekannt ist, nachdem Augustin-Louis Cauchy der Differentialrechnung die heute übliche logische Strenge gegeben hatte. Er ging von den infinitesimalen Größen ab und definierte die Ableitung als Grenzwert von Sekantensteigungen. Bereits Leibniz hatte das Wort Differentialquotient eingeführt, das den Grenzwert eines Differenzenquotienten benannte. Richtig Fahrt nahm die Anwendung der Differenzenverfahren erst mit dem Einzug von Großrechnern in Forschungseinrichtungen auf. Das war Mitte der 1950er Jahre, als komplizierte Differentialgleichungen numerisch gelöst wurden. Motor der Differenzenrechnung und der Entwicklung iterativer Methoden zur Lösung der dabei entstehenden großen linearen Gleichungssysteme war die Reaktorphysik, insbesondere die Neutronendiffusionstheorie.

Aus heutiger Sicht ist der Name Finite-Differenzen-Methode dennoch nicht so übel, denn es ist immer hilfreich, wenn sich Fachausdrücke in unterschiedlichen Sprachen ähneln. Zusammengefasst: Die Finite-Differenzen-Methode, die Finite-Elemente-Methode wie auch die Finite-Volumen-Verfahren sind Teilgebiete der Differenzenrechnung, der diskreten Entsprechung der Differentialrechnung.

FDM einer gewöhnlichen DGL 2. Ordnung

Wir suchen die Lösung einer Differentialgleichung 2. Ordnung in einer Dimension (1D) in einem geradlinigen Koordinatensystem

$u''(x)=g(x)$.

Zu berechnen sei die Funktion $u$ im Intervall $[a,b]$ der unabhängigen Variablen $x$, $a\leqslant x\leqslant b$. Die sogenannte Inhomogenität $g$ sei gegeben. Ohne weitere Bedingungen ist $u$ nicht eindeutig bestimmt. Erst mit zwei Anfangs- oder Randbedingungen erhält man Eindeutigkeit. Nur im Fall von Randbedingungen favorisiert man zur numerischen Lösung finite Differenzen.

Das Vorgehen der numerischen Lösung nach der Differenzenmethode besteht darin, aus der Differentialgleichung ein Differenzengleichungssystem zu erzeugen und dieses zu lösen. Dazu hat der Anwender im ersten Schritt auf der $x$-Achse Stützpunkte (Gitterpunkte) zu wählen, über deren Lage er frei verfügen kann.

Äquidistante Stützstellen

Die FDM beruht auf der Diskretisierung einer Funktion auf einem Gitter, hier der Funktion einer unabhängigen Variablen mit äquidistanten Stützstellen

Wählt er ein äquidistantes Gitter, so muss er sich entscheiden, mit wie vielen Stützpunkten er starten möchte. Wählt er als Anzahl der „inneren“ Stützpunkte $n$, so sind $n+1$ die Anzahl der Intervalle und $n+2$ die Anzahl der Stützstellen insgesamt. Dann sind die Gitterschrittweite (Maschenweite, Intervallbreite) $h$

$h={\frac {b-a}{n+1}}$ und die Abszissenwerte $x_{i}$ der Stützpunkte $x_{i}=a+i\cdot h$ für $i=0,\ldots ,n+1$.

Drei-Punkte-Stützstellenschema (englisch Three-Points stencil) für zentrale Differenzen in 1D mit äquidistanten Stützstellen

Im zweiten Schritt hat der Anwender eine Näherung für die zweite Ableitung zu wählen. Die entsprechenden Formeln kann er zum Beispiel aus dem Artikel Koeffizienten für Differenzenquotienten entnehmen. Die kleinst mögliche (und meist auch zu empfehlende) Anzahl von Gitterpunkten für eine Näherung der zweiten Ableitung ist drei. Für Stützstellen, die nicht auf dem Enden des Intervalls $a$ oder $b$ liegen, für die inneren Stützstellen, ist der zentrale Differenzenquotient die passende Wahl:

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

In der Grafik ist die rote die vorgegebene Stützstelle, die beiden Nachbarn sind die blauen Punkte.

Für die inneren Stützstellen folgen daraus die Differenzengleichungen

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

für die unbekannten numerischen Näherungswerte $u_{i}$ von $u(x_{i})$ und mit den vorgegebenen Werten $g_{i}$, den Werten der Inhomogenität $g$ an den Stützstellen $i$.

Die Gesamtheit der Differenzengleichungen für innere Stützstellen bilden ein lineares Gleichungssystem der Ordnung $n$ für die $n+2$ unbekannten Funktionswerte von $u$.

Randwertproblem für äquidistante Stützstellen

Sind die Funktionswerte $u(a)$ und $u(b)$ vorgegeben, haben wir es mit einem Randwertproblem zu tun. Dann ist $u$ eindeutig bestimmt.

Die beiden „Ränder“ sind die erste und die letzte Stützstelle, die das vorgegebene $x$-Intervall links und rechts begrenzen. Randwertprobleme nennt man in der Mathematik die wichtige Klasse von Problemstellungen, bei denen zu einer vorgegebenen Differentialgleichung Lösungen gesucht werden, die auf dem Rand des Definitionsbereiches vorgegebene Funktionswerte, die Randbedingungen, annehmen sollen.

Die erste Differenzengleichung wird mit den Randwerten dann zu

${\frac {1}{h^{2}}}(u(a)-2u_{1}+u_{2})=g_{1}$ oder ${\frac {1}{h^{2}}}(-2u_{1}+u_{2})=g_{1}-{\frac {1}{h^{2}}}u(a)$,

und die letzte

${\frac {1}{h^{2}}}(u_{n-1}-2u_{n}+u(b))=g_{n}$ oder ${\frac {1}{h^{2}}}(u_{n-1}-2u_{n})=g_{n}-{\frac {1}{h^{2}}}u(b)$.

In Matrixform lautet das zu lösende System, wenn wir beide Seiten aller Differenzengleichungen mit dem Faktor $-h^{2}$ multiplizieren, dann:

${\begin{pmatrix}\quad 2&-1&\quad 0&\dots &\quad 0\\-1&\quad 2&-1&\ddots &\vdots \\\quad 0&-1&\ddots &\ddots &\quad 0\\\vdots &\ddots &\ddots &\ddots &-1\\\quad 0&\dots &\quad 0&-1&\quad 2\end{pmatrix}}{\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\\\vdots \\u_{n}\end{pmatrix}}=-h^{2}{\begin{pmatrix}g_{1}-{\frac {1}{h^{2}}}u(a)\\g_{2}\\g_{3}\\\vdots \\g_{n}-{\frac {1}{h^{2}}}u(b)\end{pmatrix}}$

Die Matrix heißt Koeffizientenmatrix und ist vom Rang $n$. Sie ist regulär und somit ist das erzeugte Gleichungssystem eindeutig lösbar und ihre Inverse ist sogar formelmäßig angebbar. Da in jeder Zeile maximal nur drei Unbekannte vorkommen, gehört sie zu den dünnbesetzten Matrizen. Sie ist eine Bandmatrix, bei der zusätzlich zur Hauptdiagonalen nur eine obere und eine untere Nebendiagonale Elemente ungleich Null aufweisen und ist somit eine Tridiagonalmatrix. Da die Elemente der Hauptdiagonalen und der beiden Nebendiagonalen konstant sind, zählt sie darüber hinaus zu den Tridiagonal-Toeplitz-Matrizen.

Ähnlich spezielle Koeffizientenmatrizen treten auch anderweitig bei der Diskretisierung nach der Finiten-Differenzen-Methode auf. Solche speziellen Matrizentypen verringert die Anzahl der Rechenoperationen bei der Lösung des Gleichungssystems. Insbesondere aber spart man Arbeitsspeicher, um die Elemente der Koeffizientenmatrix zu speichern. Bei einer Anzahl der Unbekannten $n$ und einer voll besetzen Koeffizientenmatrix würde man $n^{2}$ Speicherplätze benötigen, dagegen sind es im Fall einer Tridiagonalmatrix nur $3n$. In der Frühzeit der Anwendung der FDM, in den 1950er Jahren, musste man mit dem Speicherplatz sehr haushalten. Deshalb konnte man vorrangig nur 1D-Probleme behandeln, wobei die Anzahl der Unbekannten in realen Problemstellungen durchaus schon die Millionengrenze erreichte.

Beispiel für äquidistante Stützstellen

Exakte Lösungsfunktion des Beispiels und Werte der Näherungslösungen mit 2, 3 und 4 äquidistanten Intervallen

Gegeben sei das Randwertproblem

$u''(x)=2$  für $0\leqslant x\leqslant 1$

mit den Randbedingungen

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

Die Inhomogenität $g=2$ ist also konstant.

Die Lösung $u$ lässt sich in diesem Fall exakt angeben, und zwar durch die Festlegung $u(x)=3+x(x-1)$. Dass das tatsächlich eine Lösung ist, kann man durch Einsetzen in die Differentialgleichung leicht nachweisen. Sie erfüllt auch die Randbedingungen.

Wählen wir vier Intervalle, berechnen die Näherungswerte von $u(x)$ an drei inneren Stützstellen, ist die Gitterschrittweite $h={\frac {1}{4}}=0{,}25$. Es ergibt sich das Gleichungssystem in Matrixform

${\begin{pmatrix}\quad 2&-1&\quad 0\\-1&\quad 2&-1\\\quad 0&-1&\quad 2\\\end{pmatrix}}{\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\end{pmatrix}}=-h^{2}{\begin{pmatrix}2-{\frac {3}{h^{2}}}\\2\\2-{\frac {3}{h^{2}}}\end{pmatrix}}$

mit der Lösung

${\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\end{pmatrix}}={\begin{pmatrix}2{,}8125\\2{,}7500\\2{,}8125\end{pmatrix}}$.

Die Grafik zeigt den Funktionsverlauf der exakten Lösung des Beispiels und Werte der Näherungslösungen mit 2, 3 und 4 Intervallen. Die diskreten Näherungslösungen an den Gitterpunkten sind durch farbige Kreise markiert, die durch Geraden in gleicher Farbe verbunden sind. Es ist anzumerken, dass alle Näherungslösungen dieses Beispiels Funktionswerte liefern, die genau gleich den Funktionswerten der exakten Lösung an den Stützstellen sind. Das liegt daran, dass für ein quadratisches Polynom als exakte Lösung der Konsistenzfehler bei dieser Diskretisierung gleich Null ist und damit auch der Fehler in den Gitterpunkten. Dies ist natürlich ein Ausnahmefall.

Nichtäquidistante Stützstellen

Drei-Punkte-Stützstellenschema für zentrale Differenzen in 1D mit nichtäquidistanten Stützstellen

Nichtäquidistante Stützstellen sind angebracht, wenn die Lösung Besonderheiten aufweist, z. B. Singularitäten oder Grenzschichten. Wir betrachten wieder die Differentialgleichung

$u''(x)=g(x)$

mit der vorgegebenen Inhomogenität $g$ und mit gegebenenfalls vorgegebenen Randwerten $u(a)$ und $u(b).$

Die unabhängige Variable $x$ werde mit folgenden beliebigen, nichtäquidistanten Stützstellen diskretisiert:

$a=x_{0}<x_{1}<\cdots x_{i}<\cdots <x_{n+1}=b$

mit $h_{i}=x_{i}-x_{i-1}$ und $h=\max _{i}h_{i}.$

Da die Nummerierungen $i$ der Gitterschrittweiten $h_{i}$ von Autor zu Autor unterschiedlich sind, könnte folgende anschauliche „Merkregel“ nützlich sein: Betrachten wir den Gitterpunkt $x_{i}$, so bedeuten bei der Nummerierung in diesem Abschnitt $h_{i}$ die Gitterschrittweite vor und $h_{i+1}$ die Gitterschrittweite nach dem Gitterpunkt.

Es ist folgende zentrale Drei-Punkte-Diskretisierung zweckmäßig:

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

Ordnen wir die Formel wieder in der Reihenfolge aufsteigender Stützstellen und multiplizieren beide Seiten mit $-1$, so ergibt sich

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

oder

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

Die sich daraus ergebende Koeffizientenmatrix ist wie im Fall äquidistanter Stützstellen eine Tridiagonalmatrix. Da die Elemente auf der Hauptdiagonalen und der beiden Nebendiagonalen nicht konstant sind, ist es keine Matrix vom Tridiagonal-Toeplitz-Typ.

Setzen wir $h_{i}=h$ für alle $i$, dann ergibt sich wieder die Formel für äquidistante Stützstellen

$-u_{i-1}+2u_{i}-u_{i+1}=-g_{i}h^{2}$   für   $i=1,\ldots ,n.$

Zur Untersuchung des Konsistenzfehlers wird die exakte Lösung in die Diskretisierung eingesetzt und man hat folgenden Term abzuschätzen:

${\frac {-2}{h_{i}+h_{i+1}}}({\frac {u(x_{i})-u(x_{i}-h_{i})}{h_{i}}}+{\frac {u(x_{i})-u(x_{i}+h_{i+1})}{h_{i+1}}})-u''(x_{i}).$

Taylorentwicklung liefert

$\delta ={\frac {1}{3}}u'''(x_{i})(h_{i+1}-h_{i})+{\mathcal {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 selbst in der Maximumnorm wieder von der Ordnung Zwei. Das kann man mit Hilfe der diskreten Greenschen Funktion oder verbesserten Stabilitätsaussagen verifizieren. Es gilt also

$|u(x_{i})-u_{i}|={\mathcal {O}}(h^{2}).$

Randwertproblem für nichtäquidistante Stützstellen

Sind wie oben wieder die Randwerte $u(a)$ und

$u(b)$ gegeben, wird die erste Differenzengleichung zu

$(h_{1}+h_{2})u_{1}-h_{1}u_{2}=h_{2}u(a)-g_{1}{\frac {(h_{1}+h_{2})h_{1}h_{2}}{2}}$

und die letzte

$-h_{n+1}u_{n-1}+(h_{n}+h_{n+1})u_{n}=h_{n}u(b)-g_{n}{\frac {(h_{n}+h_{n+1})h_{n}h_{n+1}}{2}}.$

Ansonsten gilt auch im Fall nichtäquidistanter Stützstellen das im Abschnitt Randwertproblem für äquidistante Stützstellen gesagte.

Beispiel für nichtäquidistante Stützstellen

Exakte Lösungsfunktion des Beispiels und Werte der Näherungslösungen mit 2, 3 und 4 nichtäquidistanten Intervallen

Wie im Fall äquidistanter Intervalle soll wieder folgendes Randwertproblem numerisch gelöst werden

$u''(x)=2$  für $0\leqslant x\leqslant 1$

mit den Randbedingungen

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

Die Inhomogenität $g=2$ ist also wieder konstant.

Die Lösungsfunktion $u$ kennen wir schon. Sie ist gegeben durch $u(x)=3+x(x-1)$.

Bevor wir die Intervalle wählen, eine Vorbemerkung zu physikalischen Größen. Das ist wichtig für Beispiele, wo Zahlenwerte für Größen eingesetzt werden. In mathematischen Lehrbüchern findet man zum Beispiel solche Texte wie: Die Gitterschrittweite $h$ beträgt $h=1/4$. Das ist, vom physikalischen Standpunkt (und mit Differentialgleichungen werden bevorzugt physikalische Probleme gelöst), nicht korrekt. Die Gitterschrittweite ist eine physikalische Größe und besteht immer aus dem Produkt einer Zahl und einer Einheit, also zum Beispiel $h=1/4\;\mathrm {cm}$. Darauf sollten auch Mathematiker hinweisen: Die Gitterschrittweite beträgt $h=1/4\;\mathrm {cm}$ (mit welcher Längeneinheit auch immer) oder Der Zahlenwert der Gitterschrittweite beträgt $h=1/4.$ Die nachfolgend verwendeten Symbole $h_{i}$ stehen nicht für physikalische Größen, sondern für Zahlenwerte der physikalischen Größen.

Behandeln wir für dieses Beispiel wieder den Fall mit vier Intervallen, berechnen folglich die Näherungswerte der Funktion $u(x)$ an drei inneren Stützstellen. Als Zahlenwerte der nicht gleich breiten Gitterschrittweiten sollen vorgegeben werden

$h_{1}={\frac {1}{6}},\;h_{2}={\frac {1}{3}},\;h_{3}={\frac {1}{3}},\;h_{4}={\frac {1}{6}}.$

Es ergibt sich das Gleichungssystem in Matrixform

${\begin{pmatrix}\quad {\frac {1}{2}}&-{\frac {1}{6}}&\quad 0\\-{\frac {1}{3}}&\quad {\frac {2}{3}}&-{\frac {1}{3}}\\\quad 0&-{\frac {1}{6}}&\quad {\frac {1}{2}}\\\end{pmatrix}}{\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\end{pmatrix}}={\begin{pmatrix}\quad 0{,}97222\\-0{,}07407\\\quad 0{,}97222\end{pmatrix}}$

mit der Lösung

${\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\end{pmatrix}}={\begin{pmatrix}2{,}86111\\2{,}75000\\2{,}86111\end{pmatrix}}.$

Die Grafik zeigt den Funktionsverlauf der exakten Lösung des Beispiels und Werte der Näherungslösungen mit zwei, drei und vier Intervallen. Bei zwei Intervallen wurden die Schrittweiten

$h_{1}={\frac {1}{3}},\;h_{2}={\frac {2}{3}},$

gewählt, bei drei Intervallen

$h_{1}={\frac {1}{4}},\;h_{2}={\frac {1}{2}},\;h_{3}={\frac {1}{4}}.$

Die diskreten Näherungslösungen an den Gitterpunkten der Grafik sind wieder durch farbige Kreise markiert, die durch Geraden in gleicher Farbe verbunden sind. Ergebnis: Für diese Differentialgleichung liefern alle Näherungslösungen in Fall nichtäquisistanter Stützstellen ebenfalls Werte, die genau gleich den Funktionswerten der exakten Lösung an den Stützstellen sind, wie auch bei äquidistanten Gitterschrittweiten. Das war nach der oben gegebenen Formel für den Konsistenzfehler $\delta$ auch zu erwarten.

Upwind finite Differenzen für ein Konvektion-Diffusionsproblem

Betrachtet wird das Konvektions-Diffusionsproblen

$-\varepsilon u''+u'=0,\quad u(0)=0,\,\,u(1)=1$

mit einem kleinen Parameter $0<\varepsilon <<1$. Durch den kleinen Parameter dominiert in diesem Problem der Konvektionsterm $u'$ und nicht wie oft üblich der Diffusionsterm $-\varepsilon u''$. Das normale Differenzenverfahren mit der äquidistanten Schrittweite $h$ lautet

$-\varepsilon {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}+{\frac {u_{i+1}-u_{i-1}}{2h}}=0.$

Mit $u_{0}=0$ und $u_{N}=1$ kann man das diskrete Problem exakt lösen und erhält mit $r=-{\frac {2\varepsilon +h}{h-2\varepsilon }}$

$u_{i}={\frac {r^{i}-1}{r^{N}-1}}.$

Wenn $\varepsilon$ extrem klein ist, verbietet es sich, mit Schrittweiten zu arbeiten, die noch kleiner als $2\varepsilon$ sind, insbesondere für mehrdimensionale Probleme. Ist aber $h>2\varepsilon$, so ist $r$ negativ, und die numerische Lösung oszilliert und approximiert die exakte Lösung extrem schlecht!

Deshalb ist es üblich, eine upwind Finite-Differenzen-Methode einzusetzen. Dabei wird $u'$ statt durch den obigen zentralen Differenzenquotienten durch einen einseitigen ersetzt, also durch $(u_{i+1}-u_{i})/h$ oder $(u_{i}-u_{i-1})/h$. Dies geschieht in Abhängigkeit von der Stromrichtung (deshalb upwind), hier in Abhängigkeit vom Vorzeichen von $u'$. Im eindimensionalen Fall wird so erreicht, dass die Koeffizientenmatrix des diskreten Problems positive Diagonalelemente und negative Außerdiagonalelemente besitzt (dies garantiert Stabilität des diskreten Problems). Das upwind-Verfahren lautet also

$-\varepsilon {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}+{\frac {u_{i}-u_{i-1}}{h}}=0,$

es liefert eine robuste Approximation der Lösung auch im Fall $h>2\varepsilon$.

Beispiel zur numerischen Lösung der Wärmeleitungsgleichung

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 Differentialgleichungen 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&\ddots &\vdots \\0&-1&\ddots &\ddots &0\\\vdots &\ddots &\ddots &\ddots &-1\\0&\dots &0&-1&2\end{pmatrix}}$.

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

Beispiel zur numerischen Lösung der zweidimensionalen Poissongleichung

Betrachtet wird exemplarisch

$-\triangle u=-\nabla ^{2}u=f\,\,{\rm {in}}\,\,\Omega =(0,1)^{2}\,\,{\rm {mit}}\,\,u|_{\partial \Omega }=0.$

Auf dem Gitter mit den Gitterpunkten $(x_{i},y_{j})$ mit $x_{i}=ih,\,y_{j}=jh$ mit $h=1/(n+1)$ und $i,j=0,1,\cdots ,n+1$ wird eine Näherung $u_{ij}$ für $u(x_{i},y_{j})$ gesucht.

Analog wie im eindimensionalen Fall diskretisiert man die zweiten Ableitungen in jeder Koordinatenrichtung in den inneren Gitterpunkten und erhält für $i,j=1,2,\cdots ,n$ das diskrete Problem

${\frac {1}{h^{2}}}(4u_{ij}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1})=f(x_{i},y_{j}).$

Dabei setzt man natürlich die Werte auf dem Rand $\partial \Omega$ gleich Null. Diese Diskretisierung heißt Fünf-Punkte-Stern für die Poissongleichung. Man kann zeigen, dass das erzeugte Gleichungssystem eine eindeutige Lösung besitzt. Zudem gilt für den Fehler

$|u(x_{ij})-u_{ij}|={\mathcal {O}}(h^{2}),$

falls die exakte Lösung des Problems ausreichend glatt ist. Dies ist im zweidimensionalen Fall insbesondere in Gebieten mit Ecken keinesfalls trivial.

Konsistenz, Stabilität und Konvergenz

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, so dass man den Konsistenzfehler folgendermaßen abschätzen kann:

$\|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.

Klassen von partiellen Differentialgleichungen zweiter Ordnung

Betrachtet werden folgende partielle Differentialgleichung zweiter Ordnung:

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

Gesucht ist die Funktion $u$, von der oft noch Anfangs- oder Randbedingungen vorgeschrieben werden.

Die Funktionen $A(\cdot ,\cdot )$, $B(\cdot ,\cdot )$, $C(\cdot ,\cdot )$, $D(\cdot ,\cdot )$, $E(\cdot ,\cdot )$ und $F(\cdot ,\cdot )$ sind gegebene 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$.

Typische Zusatzbedingungen für elliptische Probleme sind Randbedingungen, für parabolische und hyperbolische Probleme jedoch Anfangsbedingungen und Randbedingungen.

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

• Clarence Hudson Richardson: An introduction to the calculus of finite differences. Van Nostrand, New York 1954 (142 S.). 
• Christian Großmann, Hans-Görg Roos: Numerische Behandlung partieller Differentialgleichungen. 3., völlig überarb. und erw. Auflage. B. G. Teubner, Wiesbaden 2005, ISBN 978-3-519-22089-3 (572 S., eingeschränkte Vorschau in der Google-Buchsuche). }
• Stig Larsson, Vidar Thomée: Partielle Differentialgleichungen und numerische Methoden. Springer, Berlin, Heidelberg, New York 2005, ISBN 978-3-540-20823-5 (272 S., eingeschränkte Vorschau in der Google-Buchsuche). 
• Claus-Dieter Munz, Thomas Westermann: Numerische Behandlung gewöhnlicher und partieller Differenzialgleichungen: Ein anwendungsorientiertes Lehrbuch für Ingenieure. 4., verbesserte und überarbeitete Auflage. Springer Vieweg, Berlin 2019, ISBN 978-3-662-55886-7 (394 S., eingeschränkte Vorschau in der Google-Buchsuche). 

Weblinks

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

Siehe auch

• Differenzenrechnung
• Differenzenquotient
• Koeffizienten für Differenzenquotienten
• Numerische Differentiation