Numerische Methoden II, 170 024, und Übungen, 170 025

Vorlesungs-Notizen C. Brand, 7. März 2018

1 Finite-Differenzen-Methode

Unterlagen, Informationsquellen Für die ersten Einheiten zum Thema Finite-Differenzen-

Methode verwenden wir:

Kapitel 1 und Teile von Kapitel 2 aus

Randall J. LeVeque, Finite Difference Methods for Ordinary an Partial Differential
Equations, SIAM 2007.

Zu diesem Buch gibt es eine Webseite mit Aufgaben und Matlab-Code:

[Link]

Und ergänzend Kapitel 9.1 und 9.2 aus

Robert Plato, Numerische Mathematik kompakt, Grundlagenwissen für Studium und
Praxis 3., aktualisierte und verbesserte Auflage, Vieweg, Wiesbaden 2006

Als erstes sieht man heutzutage in Wikipedia nach. Dort steht zu diesem Stichwort:

Finite-Differenzen-Methoden sind eine Klasse numerischer Verfahren zur

Lösung gewöhnlicher und partieller Differentialgleichungen.

Zunächst wird das Gebiet, für das die Gleichung gelten soll, durch eine
endliche Zahl von Gitterpunkten diskretisiert. Eindimensionale Intervalle
werden dazu in gleich lange Teilintervalle zerlegt, mehrdimensionale Gebie-
te in Rechteckgitter. Die Ableitungen der gesuchten Funktion an den Git-
terpunkten werden dann durch Differenzenquotienten approximiert [. . . ].
Die Differentialgleichung wird auf diese Weise durch ein System von Diffe-
renzengleichungen angenähert, die mittels verschiedener Algorithmen zur
numerischen Lösung von Gleichungssystemen gelöst werden können.

Einleitung, Kontext

• Physikalische Gesetze sind zumeist als differentielle Beziehungen (Kraft=Masse

×Beschleunigung, Maxwell-Gleichungen, Navier-Stokes-Gleichungen. . . ) formu-
liert.

Aber: näher an der physikalischen Realität sind oft Integral-Formulierungen: Ge-

samtimpuls für einen Volumsbereich bleibt erhalten,. . . oder Extremal-Prinzipien:
Potentielle Energie wird minimal,. . . Finite Volums- beziehungsweise Finite-
Elemente-Methoden nehmen diese Formulierungen als Ausgangspunkt.

• Analytische Lösungen gibt es nur in den einfachsten Fällen. Ansonsten ist man
auf numerische Lösungen angewiesen.

1
Aufgabenstellung

Gegeben

• Eine Differentialgleichung (gewöhnlich oder partiell, erster oder höherer Ord-

nung, in einer oder mehreren unabhängigen Variablen; oder auch ein System
solcher Differentialgleichungen)

• Ein Gebiet (räumlich ein- oder mehrdimensional, Zeitintervall)

• Rand- und/oder Anfangsbedingungen

Gesucht Eine Funktion, die

• im gesamten Rechengebiet die Differentialgleichung und

• am Rand die Rand- und/oder Anfangsbedingungen erfüllt.

Aber: Diese Forderung ist zu streng, es gibt sinnvolle Aufgabenstellungen, wo eine

solche Lösung (eine Lösung im klassischen Sinn) nicht existiert, obwohl Integral- oder
Extremal-Formulierungen sogenannte schwache Lösungen zulassen.

Grundidee, Lösungsansatz: Diskretisierung des kontinuierlichen Problems

• Überziehe das Gebiet mit einem Netz aus Gitterpunkten

• Ersetze Differentialoperatoren durch Differenzen von Gitterwerten

• Aus (nicht-)linearen Differentialgleichungen werden Systeme (nicht-)linearer Glei-

chungen.

Vorteile

• einfache Herleitung, leichte Implementierung für geometrisch einfache Gebiete.

Nachteile

• auf unregelmäßigen Geometrieen. Nicht so mächtiger Theorie-Unterbau wie bei

Finiten Elementen.

1.1 Differenzenformeln

approximieren Differentialoperatoren auf einem Rechengitter

2
Elementare Formeln für erste und zweite Ableitung

lassen sich an einem Funktionsgraph geometrisch anschaulich interpretieren: Sekanten-

Steigung als Näherung an Tangenten-Steigung. Siehe Figure 1.1 in [Link].

Erste Ableitung

u(x + h) − u(x)
u0 (x) ≈ Vorwärtsdifferenz (1)
h
u(x) − u(x − h)
u0 (x) ≈ Rückwärtsdifferenz (2)
h
0 u(x + h) − u(x − h)
u (x) ≈ Zentrale Differenz (3)
2h

Zweite Ableitung

u(x − h) − 2u(x) + u(x + h)

u00 (x) ≈ Zentrale Differenz (4)
h2

Es gibt auch noch andere Differenzenformeln, zum Beispiel die unsymmetrische Vierpunkt-
Formel
1
u0 (x) ≈ [2u(x + h) + 3u(x) − 6u(x − h) + u(x − 2h)] (5)
6h
Beispiel 1: Für u(x) = sin(x) ist u0 (x) an der Stelle x = 1 gleich cos(1). Je nach
gewählter Schrittweite h liefern die Differenzenformeln unterschiedlich genaue Nähe-
rungen. Die Fehler (Näherung minus korrekter Wert) sind hier für einige Schrittweiten
tabellarisch gelistet.

h vorwärts rückwärts zentral vierpunkt

1,0 × 10−1 −4,2939 × 10−2 4,1138 × 10−2 −9,0005 × 10−4 6,8207 × 10−5
5,0 × 10−2 −2,1257 × 10−2 2,0807 × 10−2 −2,2510 × 10−4 8,6491 × 10−6
1,0 × 10−2 −4,2163 × 10−3 4,1983 × 10−3 −9,0050 × 10−6 6,9941 × 10−8
5,0 × 10−3 −2,1059 × 10−3 2,1014 × 10−3 −2,2513 × 10−6 8,7540 × 10−9
1,0 × 10−3 −4,2083 × 10−4 4,2065 × 10−4 −9,0050 × 10−8 6,9959 × 10−11

Interessant ist, wie rasch der Fehler für h → 0 abnimmt. Man erkennt:

• Vorwärts- und Rückwärtsdifferenz haben absolut etwa gleich große Fehler, aber
mit unterschiedlichen Vorzeichen.

• Bei Vorwärts- und Rückwärtsdifferenz nimmt der Fehler pro Zehnerpotenz in h

ebenfalls um eine Zehnerpotenz ab.

• Bei der zentralen Differenz nimmt der Fehler pro Zehnerpotenz in h zwei Zeh-
nerpotenzen ab.

• Bei der Vierpunkt-Formel: pro Zehnerpotenz in h drei Zehnerpotenzen!

3
Ein log-log-Diagramm zeigt besonders gut, dass die Fehler je nach Verfahren propor-
tional zu unterschiedlichen Potenzen von h sind.

Bemerkung: In doppelt logarithmischer Darstellung werden Funktionen der Art y = cxp zu

Geraden. Die Geradensteigung entspricht der Hochzahl.

Das MATLAB-Programm chap1ex1.m1 berechnet die Tabelle und zeichnet das Dia-
gramm.

Fehler verschiedener Differenzenformeln

0
10

−2
10

−4
10
Abschneidefehler

−6
10

−8
10
Vorw.
Rückw.
−10 Zentr.
10
Vierpkt.

−12
10
−3 −2 −1
10 10 10
Schrittweite h

Groß-O Schreibweise verwenden wir, um den Zusammenhang Gitterweite–Fehler

einfach zu beschreiben.

Für eine Funktion g = g(h) bedeutet der Ausdruck g = O(hp ), sprich „g ist groß-Oh
von h hoch p“: Der Betrag von g ist kleiner als eine Konstante mal |h|p , sofern h
genügend nahe bei 0 liegt.

Formaler geschrieben: Es gilt g = O(hp ) für h → 0, wenn es Konstante C und h0 gibt,

so dass |g(h)| < C|h|p ∀h mit |h| < h0 gilt.

Abschneidefehler Der Fehler einer Differenzenformel lässt sich durch Taylorreihen-

entwicklung untersuchen. Wir nehmen an, dass u(x) hinreichend oft differnzierbar ist.
Dann gilt

1 1 3 000 1 4 0000
u(x + h) = u(x) + hu0 (x) + h2 u00 (x) + h u (x) + h u (x) + O(h5 ) (6)
2 6 24
1 1 3 000 1 4 0000
u(x − h) = u(x) − hu (x) + h2 u00 (x) −
0
h u (x) + h u (x) + O(h5 ) (7)
2 6 24
1
Aus dem Begleitmaterial zu R. LeVeque

4
Dann ergeben passende Umformungen

u(x + h) − u(x) 1 00
= u0 (x) + hu (x) + O(h2 ) Vorwärts-Differenz (8)
h 2
u(x) − u(x − h) 1 00
= u0 (x) − hu (x) + O(h2 ) Rückwärts-Differenz (9)
h 2
u(x + h) − u(x − h) 1 2 000
= u0 (x) + h u (x) + O(h4 ) Zentrale Differenz (10)
2h 6

Man spricht von Abschneidefehler (truncation error), weil er durch „Abschneiden“ der
höheren Terme in der Taylorreihe entsteht.

Bei der zentralen Formel verschwinden übrigens alle ungeraden Potenzen von h, eine
typische Eigenschaft zentraler Differenzen. Für die Vierpunktformel (5) brauchen wir
noch die Reihenentwicklung von u(x − 2h). Tatsächlich verschwinden in Formel (5)
auch noch die h2 -Terme. Der Differenzenausdruck ist gleich

1 3 000
u0 (x) + h u (x) + O(h4 )
12

Verwendet man die Taylorformel mit Restglied, lassen sich die Abschneidefehler auch
ohne O(hp )-Term angeben, zum Beispiel ist bei der Vorwärtsdifferenz (8)

u(x + h) − u(x) 1
= u0 (x) + hu00 (ξ)
h 2

wobei ξ ein Wert x ≤ ξ ≤ x + h ist. Der Fehler ist also proportional der Ableitung u00 ,
ausgewertet an einer Stelle in der Nähe von x. Für eine Abschätzung reicht es oft, den
Maximalwert maxx |u00 (x)| zu verwenden.

In unserem Beispiel 1 ist 21 u00 (1) = − 12 sin(1) = −0, 42073, und das passt (wenn man
mit dem entsprechendem h multipliziert) gut zu den Zahlenwerten in der Tabelle.

Aufgaben Vergleichen Sie die theoretischen Abschneidefehler der anderen Differen-

zenformeln mit den Tabellenwerten.

Verwenden Sie u(x) = exp(x) und bestimmen Sie für u00 (1) die Abschneidefehler für
verschiedene h und Differenzenformeln. Verwenden Sie dazu aus der folgenden Tabel-
le

• symmetrische Drei-, Fünf- und Siebenpunkt-Formel

• unsymmetrische Vier-, Fünf- und Siebenpunkt-Formel

5
Eine Tabelle mit Differenzenformeln

Eine sehr ausführliche Zusammenstellung verschiedener Differenzenformel, (auch so

genannter kompakter Formeln) aus: Lothar Collatz, Numerische Behandlung von Dif-
ferentialgleichungen, Springer 1955.
Weil es mittels Computer-Algebra-Systemen inzwischen einfach ist, solche Formeln je
nach Bedarf herzuleiten, sind solche Tabellen nicht mehr so wichtig. Herleiten von
Differenzenformeln wird im nächsten Abschnitt erklärt.

6
2

7
1.2 Herleiten von Differenzenformeln

über Taylor-Reihe Ansatz mit unbestimmten Koeffizienten.

Beispiel 2: Einseitige Dreipunkt-Formel für u0 , siehe Ausarbeitung Example 1.2 in R.

LeVeque.

über Interpolationspolynom Grundidee: Finde ein Interpolationspolynom p(x), dif-

ferenziere p und verwende p0 (x) als Näherung für u0 (x).

Beispiel 3: Einseitige Dreipunkt-Formel: Wir leiten sie für die drei Wertepaare

(0; u(0)), (−h; u(−h)), (−2h; u(−2h))

her. (Eine Verschiebung des Koordinatenursprungs ändert nicht die Koeffizienten in

den Differenzenformeln. Wir berechnen daher der Einfachheit halber die Differenzen-
formel für u0 (0).) In der Lagrange-Form lässt sich das Polynom direkt hinschreiben:

p(x) = u(−2h)L−2 (x) + u(−h)L−1 (x) + u(0)L0 (x)

mit
(x + h)x (x + 2h)x (x + 2h)(x + h)
L−2 (x) = L−1 (x) = L0 (x) =
(−2h + h)(−2h) (−h + 2h)(−h) (2h)h

Differenzieren der Li (x) und Auswerten für x = 0 liefert

1 2 3
L−2 (0)0 = L−1 (0)0 = − L0 (00 ) =
2h h 2h
Vergleiche mit Formel 1.11 in LeVeque!

über Vandermonde-System in MATLAB Kapitel 1.5 in LeVeque beschreibt, wie sich

die Koeffizienten als Lösung eines Gleichungssystems mit einer Vandermonde-Matrix
finden lassen und wie MATLAB dieses System erstellt (Skript fdcoeffV.m).

1.3 Die 1-D Poisson-Gleichung

Das übliche einführende Standard-Beispiel für Randwertprobleme. Siehe Plato 9.1,

RLeVq 2.3 und 2.4 und das Beispiel beim Stichwort Finite-Differenzen-Methode in
Wikipedia.

Die 1-D Randwertaufgabe lässt sich einfach und direkt durch Integration lösen; dafür
sind Finite Differenzen-Verfahren gar nicht notwendig. Aber die hier erklärten Grundi-
deen und Lösungsansätze lassen sich direkt auf 2-D und 3-D Problemen übertragen.

8
Randwert-Aufgabe für die Poisson-Gleichung, klassische Formulierung Gesucht ist
eine Funktion u : [0, 1] → R, die für gegebenen Quellterm f : (0, 1) → R erfüllt:

u00 = f in Ω = (0, 1) (11)

u(0) = u(1) = 0

Die Aufgabe ist (in dieser Schreibweise) nur sinnvoll , wenn u in (0, 1) zweimal diffe-
renzierbar ist und bis zum Rand hin, also in [0, 1], stetig ist.

Für eine mathematisch präzise Formulierung muss festgelegt sein, welche Funktionen-
klassen für f und u betrachtet werden. Für eine klassische Lösung der Randwert-
aufgabe 11 fordert man f ∈ C 0 (Ω), u ∈ C 2 (Ω) ∩ u ∈ C 0 (Ω̄). Das ist aber für die
praktische Anwendung viel zu einschränkend: Als Quellterme f können auch Punkt-
quellen oder unstetige Stufenfunktionen sinnvoll sein. Andererseits ist in 2-D oder 3-D
die Forderung f ∈ C 0 (Ω) nicht einmal hinreichend für die Existenz einer Lösung u.

1.4 Diskretisierung der Randwertaufgabe

Äquidistantes Gitter für das Rechengebiet Ω der Randwertaufgabe 11: n innere Punkte
xj im Intervall 0 < x < 1

1
xj = jh, j = 1, . . . , n, Gitterweite h =
n+1
führt auf Gleichungssystem

1
Au = f mit A = − tridiag(−1, 2, −1) und f = [f (xj )] (12)
h2

1.5 Lokaler Diskretisierungsfehler

Es sei uex der Vektor der exakten Lösung an den Gitterpunkten: uex = [u(xj )]. Ein-
setzen von uex erfüllt das Gleichungssystem nicht exakt, es bleibt ein Restvektor δ,
der lokale Diskretisierungsfehler.

Auex − f = δ (13)

Aus Taylorreihenentwicklung folgt für j-te Zeile von 13

1
δj = (u(xj−1 ) − 2u(xj ) + u(xj+1 )) − f (xj ) =
h2
1
u00 (xj ) + h2 u0000 (ξ) − f (xj ) =
12
1 2 0000
h u (ξ) mit xj−1 ≤ ξ ≤ xj+1
12

9
Norm des lokalen Diskretisierungsfehler-Vektors
1 2 0000
||δ||∞ ≤ h ||u ||∞
12
2
Maximal-Betrag in δ ist beschränkt durch Maximalwert der vierten Ableitung mal h12 .
Beachte: links Vektornorm im Rn , rechts Norm im Raum der viermal stetig diff’baren
Funktionen.

1.6 Fehler der Lösung

 = uex − u
heißt globaler Fehler.
Der lokale Diskretisierungsfehler δ lässt sich abschätzen, uns aber interessiert der
globale Fehler . Es gilt

A = Auex − Au = Auex − f = δ

also
A = δ (14)
Der globale Fehler  erfüllt ein Gleichungssystem mit derselben Finite-Differenzen-
Matrix wie im Originalproblem, und mit Quellterm δ auf der rechten Seite.
Man kann auch so formulieren

Au = f Näherung u erfüllt dieses System

Auex = f + δ Exakte Lösung erfüllt System mit zusätzlichem Quellterm

Der lokale Diskretisierungsfehler wirkt wie ein zusätzlicher Quellterm!

Entscheidende Frage: folgt aus kleiner Störung δ im Quellterm auch ein kleiner Fehler
 in der Lösung?
Aus Gleichung 14 folgt
 = A−1 δ
und daraus die Forderung, dass A−1 keine „zu starke“ Vergrößerung des δ-Vektors
bewirken darf.

1.7 Stabilität

Angenommen, eine Finite-Differenzen-Diskretisierung auf zunehmend feineren Gittern

führt zu einer Folge von Matrix-Gleichungen

A(h) u(h) = f (h) mit h → 0 .

10
  −1
Wenn (jedenfalls für genügend kleine h) die Inverse A(h) existiert und deren Norm
durch eine von h unabhängige Konstante C beschränkt ist,
 −1
|| A(h) || ≤ C ∀h ≤ h0 ,

dann heißt die Diskretisierung stabil.

(genauer: stabil in der entsprechenden Norm; wir arbeiten hier mit Stabilität in der
Maximumsnorm)

1.8 Konsistenz

Eine Diskretisierung einer Randwertaufgabe heißt konsistent (in der entsprechenden

Norm), wenn ||δ (h) || → 0 für h → 0.

Speziell: konsistent von p-ter Ordnung, wenn ||δ (h) || = O(hp ).

Beachte: Auch die Diskretisierung der Randbedingungen ist mit zu berücksichtigen.

Die homogenen Dirichlet-Randbedingungen in unserem Musterbeispiel lassen sich pro-
blemlos einbinden. In allgemeineren Fällen kann die konsistente Diskretisierung der
Randbedingungen nicht trivial sein.

1.9 Konvergenz

Eine Diskretisierung einer Randwertaufgabe heißt konvergent, wenn ||(h) || → 0 für

h → 0.

Es gilt der wichtige Zusammenhang bei linearen Problemen:

Konsistenz + Stabilität ⇒ Konvergenz

Beweis ist einfach:

 −1  −1
||(h) || = || A(h) δ (h) || ≤ || A(h) || · ||δ (h) || ≤ C||δ (h) || → 0

Ebenso folgt: Konvergenz-Ordnung = Konsistenz-Ordnung!

11
1.10 Eigenschaften der (−1, 2, −1)-Tridiagonalmatrix

Um nachzuweisen, dass die Diskretisierung (12) der Randwertaufgabe (11) stabil ist,
und um den globalen Diskretisierungsfehler abzuschätzen, sind Eigenschaften der Ma-
trix A = tridiag(−1, 2, −1) wesentlich: Nichtsingularität und Norm.
Für diese spezielle Matrix gibt es einfache Formeln und exakte Ergebnisse: Sei A die
n × n Tridiagonalmatrix A= tridiag(-1,2,-1), dann gilt

det A = n + 1
 
 1 n(n + 2)  gerade
8
||A−1 ||∞ = für n
 1 (n + 1)2 ungerade
8

Diese Matrix ist aber nur das prototypische Beispiel für eine allgemeinere Klasse von
Matrizen: irreduzibel diagonaldominante Matrizen. Wir wollen die grundsätzlichen
Ideen vorstellen, ohne zu viel technischen Aufwand zu treiben. Die hier am einfachen
Beispiel der tridiag(−1, 2, −1)-Matrix vorgestellten Argumente und Beweismethoden
lassen sich allgemein für irreduzibel diagonaldominante Matrizen formulieren.

Diagonaldominanz X
|aii | ≥ |aij | ∀i = 1, . . . , n
j6=i

und mindestens in eine Zeile i gilt echt > statt ≥.

Irreduzibilität vereinfacht: „das Rechengitter ist zusammenhängend“

Definition: wenn für i 6= j der Matrixeintrag aij 6= 0 ist, dann heißen die Indizes i
und j direkt verbunden. Die Matrix A ist irreduzibel, wenn es ausgehend von jedem
beliebigen i eine Folge direkter Verbindungen i → j → k → · · · → ` zu jedem anderen
` gibt.
Zusammenhang Matrix–Rechengitter: Ein Differenzenoperator greift auf benachbarte
Punkte im Rechengitter zu. Sind die Gitterpunkte i und j betroffen, dann entsteht in
A ein Eintrag aij 6= 0. In der Regel ist dann auch aji 6= 0. Man sagt: Das Besetzungs-
muster der Matrix ist symmetrisch.
In unserem Musterbeispiel Gleichung 12 gilt sogar aij = aji , die Matrix A ist symme-
trisch.
Aufgaben zur Illustration: Gibt es in der Matrix
 
2 0 0 −1 −1
 0
 2 0 −1 0 

B= 0 0 2 0 −1
 
−1 −1 0
 
2 0
−1 0 −1 0 2

12
eine Verbindung von i = 2 nach j = 3? Wie sieht der zugeordnete Graph aus? Ist A
irreduzibel?

Rechengitter und Matrixstruktur Gibt es in den Matrizen

0 −1 −1 0 0 −1 0
   
4 0 4 0 0
 0
 4 0 0 −1 −1

 0
 4 0 0 −1 −1

−1 0 4 0 −1 −1  0 0 4 0 −1 −1
C= und D = 
   
−1 0 0 4 −1 0  −1 0 0 4 0 0
 
   
 0 −1 −1 −1 4 0  0 −1 −1 0 4 0
0 −1 −1 0 0 4 0 −1 −1 0 0 4

jeweils eine Verbindung von i = 1 nach j = 2? Wie sehen hier die zugeordneten
Graphen aus? Sind B und C irreduzibel?
Es ist gar nicht so einfach, ausgehend von der Matrix die Form des zu Grunde liegen-
de Rechengitters zu rekonstruieren. Die Beispielmatrizen A, B und C stammen von
einfachen Gittern, deren Punkte sehr „durcheinander“ nummeriert sind.
4 5 2 4 5 2

3 5 1 4 2
1 3 6

B C D 1 3 6

Bei „vernünftiger“ Nummerierung sieht die Matrixstruktur regelmäßiger aus. Jeden-

falls illustrieren diese Aufgaben auch:
Die Matrixstruktur hängt von der Nummerierung der Gitterpunkte ab. Umnumme-
rieren bedeutet gleichzeitiges Vertauschen entsprechender Zeilen und Spalten in der
Matrix.

A = tridiag(−1, 2, −1) ist nicht singulär Für diese Matrix lässt sich das am ein-
fachsten durch Gauß-Elimination nachrechnen. Aber wir zeigen das in der Form: aus
Ax = 0 folgt x = 0. Das ist nämlich eine Beweismethode, die auch allgemein für
irreduzible diagonaldominante Matrizen gilt.
Wenn (indirekt angenommen) ein Vektor x 6= 0 die Gleichung Ax = 0 erfüllt, dann
gibt es mindestens eine Komponente xi 6= 0. Wir wählen ein xi mit maximalem Betrag.
Das kann nicht x1 sein, weil es gilt (erste Zeile von Ax = 0 umgeformt)

2x1 = x2 .

Das wäre ein Widerspruch zu |x1 | ≥ |x2 |. Und es kann auch nicht (gleiches Argument)
xn sein. Also ist es ein xi mit i ∈ {2, . . . , n − 1}. Dann gilt (Gleichung in Zeile i von
Ax = 0 umgeformt)

1
xi = (xi−1 + xi+1 )
2

13
Diese Gleichung sagt: xi ist der Mittelwert seiner beiden Nachbarn. Es können also
nicht beide Nachbarn kleiner als xi sein. Ebenso wenig können beide Nachbarn größer
sein. Es kann aber schon gar nicht ein Nachbar größer und der andere kleiner sein,
weil wir |xi | ≥ |xj | ∀j angenommen haben.
Daher müssen xi+1 = xi = xi−1 alle gleich sein. Dann gilt aber auch (in gleicher Weise
argumentierend)
xi = xi−1 = xi−2 = · · · = x2 = x1 .
(Hier habern wir die offensichtliche Tatsache verwendet, dass im Rechengitter xi mit
xi−1 und weiter bis schließlich x1 verbunden ist. Im allgemeinen Beweis geht hier die
Irreduzibilität ein.)
Gleichzeitig gilt auch 2x1 = x2 (siehe oben). Daraus folgt schließlich:

xi = 0 ∀i

Bemerkung: Dieser Beweis lässt sich mit ganz ähnlicher Argumentation, aber etwas
mehr technischem Aufwand allgemein für irreduzible, diagonaldominante Matrizen
führen.
Deswegen sind Diagonaldominanz und Irreduzibilität wichtige Eigenschaften: sie ga-
rantieren die Lösbarkeit des Finite-Differenzen-Gleichungssystems.

A ist invers-positiv Eine Matrix X heißt invers-positiv, wenn X −1 existiert und

(elementweise) X −1 ≥ 0 gilt.
Für die (-1;2;-1)-Tridiagonalmatrix A gilt sogar: A−1 > 0.
Dass A−1 existiert, haben wir gerade gezeigt. Der Beweis von A−1 > 0 argumentiert
ähnlich. Für die k-te Spalte der Inversen gilt:

Ax = ek (k-ter Einheitsvektor auf der rechten Seite) (15)

Der k-te Einheitsvektor hat Komponenten

ek = [δki ],

(mit Kronecker-δ-Symbol; nicht zu verwechseln mit lokalem Diskretisierungsfehler

δ)
Zeilenweise ausgeschrieben:
1
x1 = x2 + δk1 (16)
2
1
xj = (xj−1 + xj+1 ) + δkj (j = 2, . . . , n − 1)
2
1
xn = xn−1 + δkn (17)
2

14
Weglassen des δkj -Terms macht daraus die Ungleichungen

1
xj ≥ (xj−1 + xj+1 ) für j = 2, . . . , n − 1
2
Es können nicht alle Komponenten in x gleich sein, denn dann würde aus einer der bei-
den Randgleichungen 16,17 folgen: xj = 0 ∀j = 1, . . . , n. Dann kann aber Gleichung 15
nicht erfüllt sein; Widerspruch.
Sei nun xi die kleinste Komponente in x (oder eine davon, falls es mehrere gibt), und
wenigstens einer der beiden Nachbarn sei echt größer. Angenomen, i ∈ 2, . . . , n − 1,
dann gilt für so ein kleinstes xi
1 1
xi ≥ (xi−1 + xi+1 ) > (xi + xi )
2 | {z } | {z } 2
≥xi ≥xi

Das hieße aber xi > xi ; Widerspruch. Also ist i = 1 oder i = n. Wir nehmen an i = 1;
der andere Fall lässt sich analog behandeln.
Wäre x1 < 0, dann wäre laut Randgleichung 16 x2 = 2x1 − δk1 < x1 , wiederum ein
Widerspruch. Ebenso ergibt x1 = 0 einen Widerspruch zu den Voraussetzungen.
Es bleibt schließlich als einzige Möglichkeit: die kleinste Komponente xi muss > 0
sein.

Monotonie vereinfacht: „Quellterm ≥ 0 bewirkt Lösung ≥ 0“. Genauer: Eine Matrix

M heißt monoton, wenn aus M u ≥ 0 folgt u ≥ 0 (jeweils zeilenweise).
Für die Monotonie der (-1;2;-1)-Tridiagonalmatrix A gilt sogar: Aus Au ≥ 0 und Au 6=
0 folgt u > 0. In Worten: selbst wenn nur eine einzige Komponente des Quellterms
positiv ist und alle anderen gleich Null, sind alle Komponenten der Lösung positiv.
Diese Aussage folgt direkt aus derEigenschaft A−1 > 0.

Fehlerabschätzung Der Vektor v = [vj ] mit vj = 21 j(n − j + 1) erfüllt

Av = 1 (der lauter-Einsen-Vektor) (18)

Die maximale Komponente von v liegt (bei ungeradem n) „in der Mitte“, für j =
(n + 1)/2, Wert vj = 81 (n + 1)2 .
Die Idee der nun folgenden Argumentation ist: Wenn auf der rechten Seite betrags-
mäßig nicht mehr als 1 steht, ist (wegen Monotonie oder, gleichbedeutend, A−1 ≥ 0)
auch die Lösung betragsmäßig nicht größer als v. Und weiter: Wenn gilt

Ax = b

mit rechter Seite kbk∞ ≤ 1, dann gilt auch für die Lösung kxk∞ ≤ 81 (n + 1)2

15
Gleichbedeutend mit
1
kA−1 k∞ ≤ (n + 1)2
8
für die (-1;2;-1)-Tridiagonalmatrix A. (Bei ungeradem n gilt Gleichheit.) Für die Ma-
trix A = − h12 tridiag(−1, 2, −1) in Gleichung (12) kürzt sich wegen h = n+1 1
der
2
(n + 1) -Term. Für diese skalierte Matrix A gilt dann

1
kA−1 k∞ ≤
8

Entschuldigung, das ist nicht klug, sowohl die mit h2 skalierte als auch die unskalierte
Tridiagonalmatrix mit demselben Symbol A zu bezeichnen. . .
Es gilt dann (siehe Abschnitt Konvergenz) die Abschätzung

1 2 1 0000 1
kk∞ ≤ kA−1 k∞ kδk∞ ≤ ·h ku k∞ = h2 ku0000 k∞
8 12 96

2 Deriving Difference Approximations for Derivatives

Consider a smooth function f (x). (“Smooth” here means that as many derivatives of
f exixt as we may need.)
We will illustrate, by way of example, a method to find approximations to derivatives.
It is sometimes called the method of undetermined coefficients.
Example: Derive an approximation for f 00 at gridpoint x0 using values at gridpoints
x−1 , x0 and x1 of an equidistant grid with spacing h.
Expand f−1 and f1 in Taylor series around x0

h2 00 h3 000 h4 iv
f−1 = f0 − hf00 + f − f + f + O(h5 )
2 0 6 0 24 0
h2 h3 000 h4 iv
f1 = f0 + hf00 + f000 + f + f + O(h5 )
2 6 0 24 0

We now seek a linear combination of the three values f−1 , f0 and f1 approximating f000
in the form
f000 ≈ af−1 + bf0 + cf1 .
To determine the coefficients a, b and c, we insert the Taylor series and collect terms
in f0 , f00 , . . .. We find

h2
f000 ≈ (a + b + c)f0 + h(−a + c)f00 + (a + c)f000 +
2
h3 h4
(−a + c)f0000 + (a + c)f0iv + O(h5 )
6 24

16
Since we have three degrees of freedom (three undetermined coefficients), we can plan
on setting the coefficients of f0 and f00 equal to zero and the coefficient of f000 equal to
one. We are left with the following system of equations.

a+b+c = 0
−a + c = 0
2
a+c =
h2

Solving the above system of equations gives

1 2 1
a= , b=− , c= .
h2 h2 h2
As an additional bonus, this solution makes the coefficient of f0000 equal to zero too.
The coefficient of f0iv turns out to be h2 /12. Hence, we get the approximation

1 2 1 00 h2 iv
f−1 − f0 + f1 = f0 + f + O(h3 ).
h2 h2 h2 12 0
This is the well-known centered difference approximation that we have used already.
As an additional exercise, you may derive a five-point centered approximation to f0000 .
The system of equations for five coefficients a, b, c, d and e in this case will be

a+b+c+d+e = 0
−2a − b + d + 2e = 0
4a + b + d + 4e = 0
−8a − b + d + 8e = 6/h3
16a + b + d + 16e = 0

(Using an algebraic manipulator such as Maple or Mathematica makes solving such

systems much easier!) The solution is
1 1 1 1
a=− , b= , c = 0, d=− , e= .
2h3 h3 h3 2h3

An alternative method: find an interpolating polynomial p(x) for the data points
(xi , fi ). Differentiate the polynomial with respect to x and evaluate at x0 .
If you do this by hand, you should use the Lagrange interpolation formula:
The polynom p(x) that interpolates the n + 1 pairs of values

(x0 , f0 ), (x1 , f1 ), . . . , (xn , fn )

is given by
p(x) = f0 L0 (x) + f1 L1 (x) + · · · + fn Ln (x) ,

17
where
(x − x0 )(x − x1 ) · · · (x − xk−1 )(x − xk+1 ) · · · (x − xn )
Lk (x) =
(xk − x0 )(xk − x1 ) · · · (xk − xk−1 )(xk − xk+1 ) · · · (xk − xn )
for each k = 0, 1, . . . , n.
In Mathematica, this method conveniently finds difference approximations for higher
derivatives.
The commands
data={{-2h,f[-2h]},{-h,f[-h]},{0,f[0]},{h,f[h]},{2h,f[2h]}};
poly=InterpolatingPolynomial[data,x];
appr=D[poly,{x,3}]/.x->0 //Together
find an interpolating polynomial and evaluate its third derivative at x = 0. Mathema-
tica gives
−f (−2 h) + 2 f (−h) − 2 f (h) + f (2 h)
appr = .
2 h3
We may find the truncation error by the command
Series[appr,{h,0,3}]
which produces the output
f (5) (0) h2
f (3) (0) + + O(h)4 .
4

Since many tables provide difference approximations for f 0 , f 00 , . . ., it is normally not

necessary to derive basic formulae. However, for irregular grids, complicated equations
or some specific boundary conditions, these methods are of great value.

3 Non-Equidistant Grid

Write
f − = f (x − h− ), f 0 = f (x), f + = f (x + h+ ),
h− = ri − ri−1 , h+ = ri+1 − ri .

Then − +
0 (f + − f 0 ) hh+ + (f 0 − f − ) hh− h+ h− 000
f (x) = − f
h− + h+ 6
!
00 2 f+ − f0 f0 − f− h+ − h− 000 (h+ )2 − h+ h− + (h− )2 iv
f (x) = − − − f − f
h + h+ h+ h− 3 12

Note that in the second approximation, for h+ 6= h− now there is a low-order trunca-
tion error term involving f 000 . Equidistant grids usually provide approximations with
higher-order truncation errors. So you have to balance the gain by finer meshsize
against possibly higher truncation error.

18
4 Linear solvers

4.1 The Thomas Algorithm for tridiagonal Systems

The so-called Thomas Algorithm is just a form of elimination for solving tridiagonal
systems of linear equations. Llewellyn H. Thomas used it around 1950 to solve elliptic
partial differential equations. The attribution to Thomas seems to be more common
in some engineering disciplines than it is in numerical analysis.

Let A be a tridiagonal matrix and b the right-hand side of a system Ax = b.

   
d1 e1 b1
 c
 2 d2 e2 b2
  
  
c3 d3 e3 b3
   
   
   
   
A= , b= . .
   
.. .. ..  ..

 . . . 
 


   
   
cn−1 dn−1 en−1   bn−1
   
 
cn dn bn

Transform it to upper triangular form, by Gaussian elimination. Use backsubstitution

to solve the resulting equivalent system Ax0 = b0

1 e01 b01
   

 1 e02 


 b02 

1 e03 b03
   
   
   
   
0 0
A = , b = . .
   
.. ..  ..

 . . 
 


   
   
e0n−1 
 0
1  bn−1
  
 
1 b0n

The expressions for the entries in A0 and b0 follow from a step-by-step application of
the usual elimination procedure.

e1 b1
e01 = , b01 = ;
d1 d1

for i = 2, . . . , n:
ei bi − ci b0i−1
e0i = , b0i = ;
di − ci e0i−1 di − ci e0i−1
Backsubstitution: xn = b0n , und für i = n − 1, . . . , 1:

xi = b0i − e0i xi+1 .

19
This algorithm is stable if


 di > 0, i = 1, 2, . . . , n,

 d > |e |,
1 1

 di ≥ |ci | + |ei | und ci 6= 0, ei 6= 0, i = 2, . . . , n − 1,


dn ≥ |cn |.

Essentially these conditions require A to be (weakly) diagonally dominant. They are

sufficient but not necessary.

20