Finite-Differenzen-Methode

Perspektivische Grafik der numerische 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 approximierten Lösungen der Differenzialgleichung an den Gitterpunkten lassen sich dann durch das entstehende lineare Gleichungssystem berechnen.

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.

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.

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(x)$ im Intervall $[a,b]$ der unabhängigen Variablen $x$, $a\leqslant x\leqslant b$. Die sogenannte inhomogene Funktion $g(x)$ sei gegeben.

Das Vorgehen der numerischen Lösung nach der Differenzenmethode besteht darin, die Differentialgleichung in ein Differenzengleichungssystem zu überführen und dieses zu lösen, in der Hoffnung, dass diese Methode zu einem Ergebnis führt, das die exakte Lösung möglichst gut annähert. 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 Five-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}$ der Lösungsfunktion $u(x)$ und mit den vorgegebenen Werten $g_{i}$, den Werten der inhomogen Funktion $g(x)$ an den Stützstellen $i$.

Die Gesamtheit der Differenzengleichungen für innere Stützstellen bilden ein lineares Gleichungssystem der Ordnung $n$. Es sei angemerkt, dass die Finite-Differenzen-Methode den Zweig der numerischen Mathematik umfasst, der das Handwerkszeug zur Diskretisierung bereitstellt. Die (meist iterativen) Verfahren zur Lösung des Gleichungssystems zählen nicht dazu.

Randwertproblem für äquidistante Stützstellen

Sind die Funktionswerte $u(a)$ und $u(b)$ vorgegeben, haben wir es mit einem Randwertproblem zu tun. 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}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}}=-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 Randwertproblem eindeutig lösbar und ihre Inverse ist sogar formelmäßig anzugeben. 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 berechnete man vorrangig 1D-Probleme, 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 Intervallen

Gegeben sei das Randwertproblem

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

mit den Randbedingungen

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

Die inhomogene Funktion $g(x)=2$ ist also konstant.

Die Lösungsfunktion $u(x)$ lässt sich in diesem Fall exakt angeben, und zwar zu $u(x)=3+x(x-1)$. Das 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}2&-1&0\\-1&2&-1\\0&-1&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 inhomogenen Funktion $g(x)$ 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}u_{i-1}+(h_{i}+h_{i+1})u_{i}-h_{i+1}u_{i+1})=-g_{i}$   für   $i=1,\ldots ,n$

oder

$-h_{i}u_{i-1}+(h_{i}+h_{i+1})u_{i}-h_{i+1}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.$

Untersucht man den Diskretisierungsfehler $\delta u_{i}$ (wie üblich durch Taylorentwicklung), so stellt man folgende Struktur fest:

$\delta u_{i}=u_{i}-u(x_{i})={\frac {1}{3}}u'''(x_{i})(h_{i+1}-h_{i})+{\mathcal {O}}(h^{2}).$

Auf äquidistanten Gittern verschwindet der erste Anteil und der Diskretisierungsfehler ist von der Ordnung 2. Auf beliebigen Gittern jedoch reduziert sich der Diskretisierungsfehler 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.

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_{2}u_{2}=h_{1}u(a)-g_{1}{\frac {(h_{1}+h_{2})h_{1}h_{2}}{2}}$

und die letzte

$-h_{n}u_{n-1}+(h_{n}+h_{n+1})u_{n}=h_{n+1}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 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 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.

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, 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.

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., 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