Numerik IIB

— Numerische Verfahren für Differentialgleichungen —

Bernhard Schmitt

Sommersemester 1997

Inhaltsverzeichnis

1 Gewöhnliche Differentialgleichungen 5 1.1 Theoretische Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Einschrittverfahren für Anfangswertprobleme . . . . . . . . . . . . . . . . . . . . 9

1.2.1 Herleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2 Konsistenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2.3 Stabilität . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.4 Konvergenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2.5 Schrittweitensteuerung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

1.3 Mehrschrittverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.3.1 Adams-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.3.2 Lineare Mehrschrittverfahren und Stabilität . . . . . . . . . . . . . . . . . 24

1.4 Extrapolationsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.5 Steife Anfangswertprobleme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1.6 Schießverfahren für Randwertprobleme . . . . . . . . . . . . . . . . . . . . . . . . 36 1.7 Differenzenverfahren für Randwertprobleme . . . . . . . . . . . . . . . . . . . . . 42

2 Gemeinsame Prinzipien von Diskretisierungsverfahren 46 2.1 Konvergenz, Konsistenz, Stabilität . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.2 Allgemeine Verfahrens-Ansätze . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2 INHALTSVERZEICHNIS

3 Partielle Differentialgleichungen 52 3.1 Allgemeine Eigenschaften . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.2 Differenzenverfahren für elliptische Randwertprobleme . . . . . . . . . . . . . . . 54

3.2.1 Die Poissongleichung auf einfachen Gebieten . . . . . . . . . . . . . . . . 54 3.2.2 Eigenschaften der Differenzen-Matrizen . . . . . . . . . . . . . . . . . . . 59 3.2.3 Diskretisierungsfehler, Konvergenz . . . . . . . . . . . . . . . . . . . . . . 62 3.2.4 Allgemeinere Gebiete und Gleichungen . . . . . . . . . . . . . . . . . . . . 66

3.3 Finite-Elemente-Verfahren für elliptische Probleme . . . . . . . . . . . . . . . . . 69 3.3.1 Variationsformulierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.2 Rayleigh-Ritz-Galerkin-Verfahren . . . . . . . . . . . . . . . . . . . . . . . 72 3.3.3 Spline-Räume, Finite Elemente . . . . . . . . . . . . . . . . . . . . . . . . 73

3.4 Partielle Anfangs-Randwert-Probleme . . . . . . . . . . . . . . . . . . . . . . . . 77

Index 82

Einleitung

In den (Natur-) Wissenschaften lassen sich viele der beobachteten Effekte durch Differentialgleichungen unterschiedlicher Art beschreiben. Ein gutes aktuelles Beispiel ist der Komet Hale-Bopp, da man sich hier tatsächlich für unterschiedliche Aspekte interessiert.

Zur Bestimmung seiner Bewegung durch das Sonnensystem reicht es aus, ihn als punktförmige Masse im Schwerefeld von Sonne, Planeten usw. zu betrachten. Die Bahn (in Abhängigkeit von der Zeit) wird dann durch ein System von gewöhnlichen Differentialgleichungen beschrieben, wie sie im ersten Teil der Vorlesung behandelt werden.

Die sichtbare Gestalt des Kometen am Himmel kommt dagegen durch die Einwirkung der Sonnenstrahlung und des Sonnenwinds zustande. Der Schweif bildet sich aus verdampfter Kometenmaterie, seine Form wird durch die Strömung im Sonnenwind bestimmt. Die Modellierung der Materiedichte im dreidimensionalen Raum erfordert partielle Differentialgleichungen. Einfache Fälle solcher Dgln werden im zweiten Teil der Vorlesung besprochen.

Einige durch Differentialgleichungen beschriebene Problemstellungen Einfache Modelle physik., biolog., chem.,.. Prozesse; z.B.,

a) Massenbewegung unter Gravitation.

Vereinfachende Annahmen: Vakuum, Massen starr, punktförmig, konzentriert,... a1) Wurf auf Erde: Flugbahn u(t) = ( \begin{matrix}x(t) \\ y(t) \end{matrix} ), Newtonsches Gesetz:

mu′′(t) = m ( \begin{matrix}x′′(t) \\ y′′(t) \end{matrix} ) = −m ( \begin{matrix}0 \\ g \end{matrix} ). Explizite Lösung Wurfparabel

’Anfangswertproblem’: Start-Ort &-Geschwindigkeit gegeben ⇒ ∃_{1} Lösung! ’Randwertproblem’: Start-& Zielpunkt gegeben, Lösung?

a2) Astronomie (Raumflug, n-Körper-Problem): nicht explizit lösbar, numerische

Verfahren erforderlich (hoher Genauigkeit)

b) Wachsende Population in Nährmedium (Bakterien, Protozoen,..)

b1) Kleiner Behälter, gute Durchmischung:

Population p(t) zur Zeit t wächst nach Gesetz

\frac{dp}{dt} = αp − βp^{2} = α(1 − \frac{β}{α}p)p.

Dabei: α = Vermehrungsrate, βp^{2} = Wettbewerb. Lösung: logistische Kurve

Gewöhn- liche Dgln

b2) Rohrförmiger Behälter, p abhängig von x und t: ^{a b \overline{}x}

Wanderung (Diffusion) im Rohr prop. zu p_{xx}, für p(t, x) gilt

\frac{∂p}{∂t} = αp − βp^{2} + γ \frac{∂^{2}p}{∂x^{2}}

Zusatzinformation: Startpopulation p(0, x), Rohrende: p_{x}(t, a) = p_{x}(t, b) = 0, d.h., Anfangswerte bzgl. t, Randwerte bzgl. x.

<<<< >>>> Partielle Dgln 2.Ordnung, Einteilung nach Typen (hier in Orts-Gebiet D ⊂ R^{2})

a) Parabolische Dgln: Ausgleichsvorgänge u(t, x, y), z.B., Wärmeleitungs-Gl.

u_{t} = α(u_{xx} + u_{yy}).

Sinnvolle Aufgabe: Anfangswerte in t = 0, Randwerte auf ∂D

b) Elliptische Dgln: Gleichgewichtszustände u(x, y), z.B., Potentialgleichung

u_{xx} + u_{yy} = 0.

Sinnvolle Aufgabe: Randwerte auf ∂D

c) Hyperbolische Dgln: Schwingungen/Wellen u(t, x, y), z.B., Wellen-Gl.

u_{tt} = α(u_{xx} + u_{yy}).

Sinnvolle Aufgabe: Anfangswerte u, u_{t} in t = 0, Randwerte auf ∂D Parti- elle Dgln

1 Gewöhnliche Differentialgleichungen

1.1 Theoretische Grundlagen

Im allgemeinsten Fall sind gewöhnliche Differentialgleichungen Systeme von nichtlinearen Gleichungen, die verschiedene Ableitungen einer Vektor-Funktion u : R → R^{n} miteinander verknüpfen. Meist betrachtet man aber zur Standardisierung explizite Gleichungen der Form

u′(x) = \frac{du}{dx} = f(x, u(x)), kurz u′ = f(x, u),  (1.1.1)

wobei f : R × R^{n} → R^{n} eine gegebene stetige Funktion ist. Als Lösung einer solchen Differentialgleichung (Dgl) ist eine stetig differenzierbare Vektorfunktion u : [a, b] → R^{n}, also eine Raumkurve gesucht, die für jeden x-Wert aus [a, b] im Punkt u(x) den Ableitungswert f(x, u(x)) annimmt, sich also in das durch f(x, y), x ∈ [a, b], y ∈ R^{n} gegebene Richtungsfeld einpaßt.

Die ausführliche Bezeichnung für (1.1.1) ist System von Dgln erster Ordnung, da nur die erste Ableitung u′ der gesuchten Funktion auftaucht. Jede Dgl m-ter Ordnung,

z^{(m)}(x) = g(x, z(x), z′(x), …, z^{(m−1)}(x))

für eine Funktion z : R → R (zur Vereinfachung) kann durch Einführung von Hilfsfunktionen

u_{1}(x) := z(x), u_{2}(x) := z′(x), …, u_{m}(x) := z^{(m−1)}(x),

in die Standardform (1.1.1) mit n = m Komponenten gebracht werden, denn dann ist

u′ = ( u′_{1} ; ⋮ ; u′_{m−1} ; u′_{m} ) = ( u_{2} ; ⋮ ; u_{m} ; g(x, u_{1}, …, u_{m}) ) .

Aus diesem Grund werden im folgenden meist Systeme erster Ordnung behandelt, und nur in wenigen Fällen der häufig auftretende Fall von Dgln zweiter Ordnung.

Die Dgl (1.1.1) definiert durch jeden Punkt (x, y) ∈ R × R^{n} eine eigene Lösungskurve. Zur Festlegung einer eindeutigen Lösung sind weitere Bedingungen erforderlich. Dazu werden zwei Standardfälle behandelt. Beim Anfangswertproblem (AWP) wird der Funktionswert an einer Stelle x_{0} ∈ [a, b] vollständig vorgeschrieben (im folgenden oBdA x_{0} = a)

u(a) = u_{0} ∈ R^{n}.  (1.1.2)

Die Lösung ist dann in einem (evtl. unbeschränkten) Intervall [a, b] gesucht. Für dieses Problem gibt es recht allgemeine Existenz- und Eindeutigkeitsaussagen. Sowohl aus theoretischer als auch praktischer Sicht schwieriger ist das Randwertproblem (RWP), wo nach Lösungen gesucht wird, für die eine bestimmte Beziehung der Funktionswerte am Anfangs- und Endpunkt eines Intervalls [a, b] gefordert wird, etwa in Form eines zusätzlichen (evtl. nichtlinearen) Gleichungssystems

r(u(a), u(b)) = 0, r : R^{n} × R^{n} → R^{n}.  (1.1.3)

Es folgen einige grundlegende theoretische Aussagen zu den genannten Problemen.

Für das Anfangswertproblem (1.1.1),(1.1.2) gibt es eine sehr einfache und allgemeine Existenzaussage im Satz von Picard-Lindelöf. Er beruht auf dem Übergang von der Dgl durch Integration zur äquivalenten Integralgleichung

u(x) = u_{0} + ∫_{a}^{x} f(t, u(t)) dt,  (1.1.4)

der auch für die Konstruktion numerischer Verfahren grundlegend ist.

Satz 1.1.1 Die Funktion f sei auf dem Streifen Ω := {(x, y) : x ∈ [a, b], y ∈ R^{n}} über dem endlichen Intervall [a, b] definiert und stetig. Außerdem gelte dort die Lipschitzbedingung

∥f(x, u) − f(x, v)∥ ≤ L∥u − v∥ ∀x ∈ [a, b], u, v ∈ R^{n},  (1.1.5)

mit einer Lipschitzkonstanten L ≥ 0. Dann existiert zu jedem u_{0} ∈ R^{n} genau eine Lösung u ∈ (C^{1}[a, b])^{n} des Anfangswertproblems (1.1.1),(1.1.2). Diese Lösung hängt stetig vom Anfangswert ab. Für zwei Lösungen u, v ∈ (C^{1}[a, b]))^{n} mit

\begin{matrix}u′ = f(x, u) \\ u(x_{0}) = u_{0} \end{matrix} , \begin{matrix}v′ = f(x, v) \\ v(x_{0}) = v_{0} \end{matrix} ,

gilt dabei die Schranke

∥u(x) − v(x)∥ ≤ e^{L|x−a|} ∥u_{0} − v_{0}∥ ∀x ∈ [a, b].  (1.1.6)

Beweis a) Bei (1.1.4) handelt es sich um ein Fixpunktproblem u = T u mit der Abbildung

T : \{ \begin{matrix}(C[a, b])^{n} → (C[a, b])^{n}, \\ v ↦ T v, \end{matrix} mit (T v)(x) := u_{0} + ∫_{a}^{x} f(t, v(t)) dt.

Für y ∈ (C[a, b])^{n} wird die gewichtete Norm

∥y∥_{L} := max{∥y(x)∥e^{−L|x−a|} : x ∈ [a, b]}

eingeführt, mit der (C[a, b])^{n} ein Banachraum wird. Aus Voraussetzung (1.1.5) läßt sich direkt die Kontraktivität von T ableiten, denn für beliebige u, v ∈ (C[a, b])^{n} gilt

∥T u − T v∥_{L} = max_{x∈[a,b]} e^{−L|x−a|}∥ ∫_{a}^{x} (f(t, u(t)) − f(t, v(t)))dt∥
≤ max_{x∈[a,b]} e^{−L|x−a|} ∫_{a}^{x} ∥f(t, u(t)) − f(t, v(t))∥ dt
≤ L max_{x∈[a,b]} e^{−L|x−a|} ∫_{a}^{x} e^{L|t−a|}e^{−L|t−a|}∥u(t) − v(t)∥ dt
≤ L max_{x∈[a,b]} e^{−L|x−a|} ∫_{a}^{x} e^{L|t−a|} dt ∥u − v∥_{L}
≤ (1 − e^{−L|b−a|})∥u − v∥_{L}.

Da der Vorfaktor kleiner als eins ist, ist T eine Kontraktion und nach dem Banachschen Fixpunktsatz existiert ein eindeutiger Fixpunkt u = T u ∈ (C[a, b])^{n}. Mit u ist auch die Funktion f(x, u(x)) = u′(x) stetig.

b) Aus der Differenz der beiden für u und v geltenden Gleichungen (1.1.4) folgt wie eben

∥u(x) − v(x)∥ = ∥u_{0} − v_{0} + ∫_{a}^{x} (f(t, u(t)) − f(t, v(t))) dt∥
≤ ∥u_{0} − v_{0}∥ + L ∫_{a}^{x} ∥u(t) − v(t)∥ dt.

Für die skalare Funktion δ(x) := ∥u(x) − v(x)∥ gilt also die Ungleichung δ(x) ≤ ∥u_{0} − v_{0}∥ + L ∫_{a}^{x} δ(t) dt. Mit dem folgenden Gronwall-Lemma ergibt sich daraus die Schranke (1.1.6).

Im Beweis wurde die folgende Aussage benutzt.

Lemma 1.1.2 (Gronwall-Lemma) Die Funktion δ ∈ C[a, b] erfülle die Ungleichung

δ(x) ≤ α + L ∫_{a}^{x} δ(t) dt ∀x ∈ [a, b],

mit α, L ≥ 0. Dann folgt

δ(x) ≤ α e^{L(x−a)} ∀x ∈ [a, b].

Beweis Es sei ε > 0. Für die Funktion γ(x) := (α + ε)e^{L(x−a)} gilt γ(x) = α + ε + L ∫_{a}^{x} γ(t) dt. Offensichtlich ist δ(a) < γ(a). Sei nun a < x_{1} ≤ b die erste Stelle mit δ(x_{1}) = γ(x_{1}), insbesondere gelte also δ(x) ≤ γ(x) ∀a ≤ x ≤ x_{1}. Dann folgt in x_{1} mit

δ(x_{1}) ≤ α + L ∫_{a}^{x_{1}} δ(t) dt < α + ε + L ∫_{a}^{x_{1}} γ(t) dt = γ(x_{1})

aber ein Widerspruch. Daher gilt δ(x) < (α + ε)e^{L(x−a)} für jedes ε > 0.

Besonders einfach ist das Studium von linearen Differentialgleichungen,

u′(x) = A(x)u(x) + g(x), d.h., f(x, y) = A(x)y + g(x),  (1.1.7)

mit einer Matrixfunktion A(x) ∈ R^{n×n} und einer Vektorfunktion g(x) ∈ R^{n}, da hier explizite Lösungsdarstellungen gelten. Mit Hilfe eines Fundamentalsystems der Dgl (1.1.7), d.h. einer regulären Matrixfunktion W (x) ∈ R^{n×n}, die die Dgl

W′(x) = A(x)W (x)

erfüllt, ist die Lösung des Anfangswertproblems (1.1.7), (1.1.2) explizit gegeben durch

u(x) = W (x)W (a)^{−1}u_{0} + W (x) ∫_{a}^{x} W (t)^{−1}g(t) dt.

Da mit W (x) auch W (x)M (M regulär) ein Fundamentalsystem ist, vereinfacht sich diese Darstellung durch Übergang zum ausgezeichneten System Y (x) := W (x)W (a)^{−1} mit Y (a) = I.

Satz 1.1.3 Die Lösung des linearen Anfangswertproblems (1.1.7),(1.1.2) mit einer beschränkten, stetigen Matrixfunktion A(x) ist

u(x) = Y (x)u_{0} + Y (x) ∫_{a}^{x} Y (t)^{−1}g(t) dt,  (1.1.8)

wobei Y (x) das Fundamentalsystem mit Y′ = AY, Y (a) = I, ist.

Explizit angebbar ist das Fundamentalsystem bei konstanter Matrix A. Über den Ansatz u(t) = e^{λx}y läßt sich die homogene Dgl auf das Eigenwertproblem (A−λI)y = 0 zurückführen. Bei einer diagonalisierbaren Matrix A mit Eigen-Wert-/-Vektor-Paaren (λ_{j}, y_{j}) ist u(x) = ∑_{j=1}^{n} α_{j}e^{λ_{j}x}y_{j}. Äquivalent dazu ist analog zum skalaren Fall (n = 1) die Darstellung

Y (x) = e^{(x−a)A} := ∑_{k=0}^{∞} \frac{(x − a)^{k}}{k!} A^{k}.

Dies vereinfacht die Lösungsformel (1.1.8) weiter zu

u(x) = e^{(x−a)A}u_{0} + ∫_{a}^{x} e^{(x−t)A}g(t) dt.

Die Beziehung dieser expliziten Formeln zu Satz 1.1.1 ist folgende. Bei (1.1.7) gilt die Lipschitzbedingung (1.1.5) mit L = max_{x} ∥A(x)∥, aus der für die Differenz von Lösungen folgt

∥u(x) − v(x)∥ = ∥Y (x)(u_{0} − v_{0})∥ ≤ e^{L|x−a|}∥u_{0} − v_{0}∥,

also ∥Y (x)∥ ≤ e^{L|x−a|}. Im Fall konstanter Koeffizienten reduziert sich dies auf die Ungleichung ∥e^{tA}∥ ≤ e^{t∥A∥}, t ≥ 0, die auch direkt aus der Reihendarstellung folgt.

Mit Hilfe der expliziten Lösungsdarstellungen soll kurz die Situation bei linearen Randwertproblemen diskutiert werden. Gesucht ist eine Lösung von

u′(x) = A(x)u(x) + g(x), x ∈ (a, b), R_{0}u(a) + R_{1}u(b) = d,  (1.1.9)

mit festen n × n-Matrizen R_{0}, R_{1} und d ∈ R^{n}. Dabei ist rang(R_{0}, R_{1}) = n vorauszusetzen. Im Gegensatz zum Anfangswertproblem können hier alle Lösbarkeitsfälle auftreten. Die Lösung u wird durch die explizite Formel (1.1.8) mit einem unbekannten Startvektor η ∈ R^{n} dargestellt,

u(x) = Y (x)η + γ(x), γ(x) := ∫_{a}^{x} Y (x)Y (t)^{−1}g(t) dt.

Nach Einsetzen in die Randbedingung ergibt sich das lineare Gleichungssystem

(R_{0}Y (a) + R_{1}Y (b))η = d − Rγ(b)  (1.1.10)

für η. Wie aus der Linearen Algebra bekannt, gilt die Fallunterscheidung,

Wenn R_{0}Y (a) + R_{1}Y (b) \{ \begin{matrix}regulär: & es existiert eine eindeutige Lösung \\ singulär: & \{ \begin{matrix}es gibt keine Lösung \\ es gibt unendlich viele Lösungen \end{matrix} \end{matrix} .

Es gibt bei Randwertproblemen also keine zum Satz 1.1.1 vergleichbare allgemeine Existenzoder Eindeutigkeitsaussage.

1.2 Einschrittverfahren für Anfangswertprobleme

1.2.1 Herleitung

Die Lösung der Dgl u′ = f(x, u) ist diejenige Kurve (x, u(x)), die dem Richtungsfeld (x, y, f(x, y)) folgt. Ein elementares Approximationsverfahren beruht darauf, daß man vom Startwert u_{0} = u(x_{0}) ’ein Stück weit’ in Richtung der Ableitung u′(x_{0}) = f(x_{0}, u_{0}) geht und dann dort dieses Verfahren wiederholt (Eulerscher Polygonzug, s.u.). Einen allgemeinen Zugang zu solchen schrittweise vorgehenden numerischen Integrationsverfahren bekommt man über die Formel (1.1.4). Dazu wird im Punkt x̄ eine Schrittweite h > 0 gewählt und die Dgl von x̄ bis x̄ + h integriert. Dies ergibt

\frac{1}{h}[u(x̄ + h) − u(x̄)] = \frac{1}{h} ∫ x̄+h ; x̄ f(t, u(t)) dt = ∫_{0}^{1} f(x̄ + hs, u(x̄ + hs)) ds.  (1.2.1)

In numerischen Verfahren wird das Integral auf der rechten Seite approximiert, z.B., durch eine geeignete Quadraturformel. In diesem Abschnitt werden Näherungen der Form

∫_{0}^{1} f(x̄ + hs, u(x̄ + hs)) ds = \underbrace{f_{h}(x̄, u(x̄))}_{Näherung} + \underbrace{T_{h}(x̄)}_{Fehler} .  (1.2.2)

betrachtet, in denen die Verfahrensfunktion f_{h} nur vom Wert (x̄, u(x̄)) abhängt. Im Verfahren wird der lokale Fehler T_{h} vernachlässigt und die Differentialgleichung (1.1.1) durch eine Differenzengleichung ersetzt. Somit ergeben sich Näherungswerte u_{Δ}(x) an u(x) aus der Formel

\frac{1}{h}[u_{Δ}(x̄ + h) − u_{Δ}(x̄)] := f_{h}(x̄, u_{Δ}(x̄))

nur an endlich vielen Stellen x_{i}, also auf einem Gitter Δ.

Definition 1.2.1 Ein Einschrittverfahren zur Lösung des Anfangswertproblems (1.1.1), (1.1.2) besteht aus der Wahl einer Verfahrensfunktion f_{h}(x, y) und eines Gitters

Δ : a := x_{0} < x_{1} < … < x_{N} = b, h_{i} := x_{i+1} − x_{i}, H := max_{i=0}^{N−1} h_{i}, |Δ| := H.

Damit berechnet sich die Näherungslösung u_{Δ}(x_{i}) = y_{i}, i = 0, …, N auf dem Gitter rekursiv

\begin{aligned}y_{0} & := & u_{0}, \\ y_{i+1} & := & y_{i} + h_{i}f_{h_{i}}(x_{i}, y_{i}), i = 0, …, N − 1. \end{aligned}  (1.2.3)

Die Bezeichnung Einschrittverfahren bezieht sich auf die Tatsache, daß in (1.2.3) nur Näherungswerte y_{i}, y_{i+1} aus einem Schritt x_{i} → x_{i+1} verknüpft werden. Es folgen vier einfache Beispiele.

Verfahren 1, Euler-Cauchy-Polygonzug: Das Integral in (1.2.2) wird durch den Funktionswert am linken Rand ersetzt,

∫_{0}^{1} f(x̄ + hs, u(x̄ + hs)) ds ≅ f(x̄, u(x̄)) =: f_{h}(x̄, u(x̄)),

das Verfahren ist also

y_{i+1} := y_{i} + h_{i}f(x_{i}, y_{i}),

i = 0, 1, . . .. Klartext: ‘In x_{i} geht man einen Schritt der Länge h_{i} in Richtung des Richtungsfeldes’

[figure: y; u_{0}; x_{0}; y_{1}; x_{1}; y_{2}; x_{2}; u(x)]

Verfahren 2, von Runge: Integralapproximation durch die Rechteckregel,

∫_{0}^{1} f(x̄ + hs, u(x̄ + hs)) ds ≅ f(x̄ + \frac{h}{2}, u(x̄ + \frac{h}{2})).

Der unbekannte Wert u(x̄+ \frac{h}{2}) wird durch Verfahren 1 angenähert,

u(x̄ + \frac{h}{2}) ≅ u(x̄) + \frac{h}{2}f(x̄, u(x̄)) =: ỹ,

die Verfahrensfunktion ist hier also

f_{h}(x, y) := f(x + \frac{h}{2}, y + \frac{h}{2}f(x, y)). [figure: x_{i}; x_{i} + \frac{h}{2}; x_{i+1}; y_{i}; y_{i+1}; ỹ]

Verfahren 3, von Heun: Integralapproximation durch die Trapezregel,

∫_{0}^{1} f(x̄ + hs, u(x̄ + hs)) ds ≅ \frac{1}{2}[f(x̄, u(x̄)) + f(x̄ + h, u(x̄ + h))].

Wieder wird u(x̄ + h) durch Verfahren 1 angenähert. Jetzt ist

f_{h}(x, y) := \frac{1}{2}[f(x, y) + f(x + h, y + hf(x, y))].

Verfahren 4, Klassisches Runge-Kutta-Verfahren: Integralapproximation durch eine modifizierte Simpsonregel,

\begin{aligned}f_{h}(x, y) & := & \frac{1}{6}[k_{1} + 2k_{2} + 2k_{3} + k_{4}], & mit \\ k_{1} & := & f(x, y) & ≅ u′(x) \\ k_{2} & := & f(x + \frac{h}{2}, y + \frac{h}{2}k_{1}) & ≅ u′(x + \frac{h}{2}) \\ k_{3} & := & f(x + \frac{h}{2}, y + \frac{h}{2}k_{2}) & ≅ u′(x + \frac{h}{2}) \\ k_{4} & := & f(x + h, y + hk_{3}) & ≅ u′(x + h) \end{aligned}

Beispiel 1.2.2 Dgl u′ = u =: f(x, u), x_{0} = 0, u_{0} = 1 ⇒ u(x) = e^{x}.

Verfahren Näherung  lokal.Fehler
1  y_{1} = y_{0} + hy_{0} = 1 + h  h^{2}/2
2  y_{1/2} = 1 + \frac{h}{2}, y_{1} = 1 + h(1 + \frac{h}{2}) = 1 + h + \frac{1}{2}h^{2}  h^{3}/6
3  ỹ_{1} = 1 + h, y_{1} = 1 + \frac{h}{2}(1 + 1 + h) = 1 + h + \frac{1}{2}h^{2}  h^{3}/6
4  k_{1} = 1, k_{2} = 1 + \frac{h}{2}, k_{3} = 1 + \frac{h}{2} + \frac{h^{2}}{4} , k_{4} = 1 + h + \frac{h^{2}}{2} + \frac{h^{3}}{4} ⇒
y_{1} = 1 + h + \frac{1}{2}h^{2} + \frac{1}{6}h^{3} + \frac{1}{24}h^{4}  h^{5}/120

Also erzeugen alle Verfahren den Anfang der Exponentialreihe unterschiedlich weit. Sie sind Spezialfälle einer allgemeinen Verfahrensklasse, der Runge-Kutta-Verfahren. Diese bestehen aus m ∈ N Stufen und werden durch eine Tabelle von Koeffizienten definiert,

\frac{c A}{b} = \begin{aligned}c_{1} & a_{11} & ⋯ & a_{1m} \\ ⋮ & ⋮ & & ⋮ \\ c_{m} & a_{m1} & ⋯ & a_{mm} \\ & b_{1} & ⋯ & b_{m} \end{aligned}

wobei c_{i} = ∑_{j=1}^{m} a_{ij}, i = 1, …, m gilt. Die Verfahrensfunktion beim allgemeinen Runge-Kutta-Verfahren ist damit

f_{h}(x, y) := ∑_{i=1}^{m} b_{i} k_{i} ,  (1.2.4)

k_{i} := f(x + hc_{i}, y + h ∑_{j=1}^{m} a_{ij}k_{j}), i = 1, …, m.

Wenn die Koeffizientenmatrix A strikt untere Dreieckgestalt besitzt, also für a_{ij} = 0 ∀j ≥ i, können die Hilfsgrößen k_{i} der Reihe nach explizit aus diesen Gleichungen berechnet werden. Dann liegt ein explizites Runge-Kutta-Verfahren vor. Die Verfahren 1–4 gehören zu dieser Klasse, das klassische Runge-Kutta-Verfahren 4 etwa hat die Tabelle

0 0 ⋯
\frac{1}{2} \frac{1}{2} 0 · ·
\frac{1}{2} 0 \frac{1}{2} 0 ·
1 0 0 1 0
\frac{1}{6} \frac{1}{3} \frac{1}{3} \frac{1}{6}

Im allgemeinen Fall (∃a_{ij} ≠ 0, j ≥ i) ist (1.2.4) ein implizites Runge-Kutta-Verfahren, da zur Berechnung der k_{i} i.a. nichtlineare Gleichungssysteme zu lösen sind. Solche Verfahren können sehr gute Eigenschaften besitzen und werden daher später noch einmal betrachtet.

Die entscheidende Frage bei den numerischen Verfahren (1.2.3) ist die nach der Konvergenz der Näherungswerte u_{Δ}(x) → u(x), |Δ| → 0, (x fest) bei feiner werdendem Gitter. Der Nachweis dieser Tatsache erfolgt nach einem sehr allgemeinen Prinzip zweistufig in der Form

Konsistenz + Stabilität ⟹ Konvergenz

1.2.2 Konsistenz

Eine Verfahrensfunktion f_{h} sollte sinnvollerweise so gewählt sein, daß die exakten Lösungen u der Differentialgleichung auch die Differenzengleichung (1.2.3) näherungsweise erfüllen.

Definition 1.2.3 Das Verfahren (1.2.3) heißt konsistent, wenn mit jeder Lösung u der Dgl (1.1.1) bei stetig differenzierbarer rechter Seite f für den lokalen Fehler

T_{h}(x) := \frac{1}{h}[u(x + h) − u(x)] − f_{h}(x, u(x)), x ∈ [a, b − h],  (1.2.5)

gilt T_{h} → 0 für h → 0. Das Verfahren heißt konsistent mit Ordnung p > 0, wenn

T_{h}(x) = O(h^{p}), h → 0, x ∈ [a, b),  (1.2.6)

für genügend oft differenzierbare rechte Seiten f.

Beispiel 1.2.4 Verfahren 2, f_{h} = f(x + \frac{h}{2}, y + \frac{h}{2}f(x, y)). Die Taylorentwicklung einer Lösung u in x ist

\frac{1}{h}[u(x + h) − u(x)] = u′(x) + \frac{h}{2}u′′(x) + O(h^{2}).

Für f ∈ C^{2} erhält man die Ableitungen u^{(j)}(x) der Lösung aus u′ = f(x, u), u′′ = f_{x} + f_{y}u′ = f_{x} + f_{y}f. Daher stimmt die Taylorentwicklung von f_{h} (nach h),

f(x + \frac{h}{2}, u + \frac{h}{2}f) = f(x, u) + \frac{h}{2}(f_{x} + f_{y}f) + O(h^{2})

mit der obigen in den ersten beiden Gliedern überein, das Verfahren hat also Konsistenzordnung p = 2.

Konsistenznachweise sind zwar elementar durchführbar, aber sehr arbeitsaufwendig, Standardmethode ist dabei die im Beispiel angewendete Taylorentwicklung. Umgekehrt können durch Vergleich der Taylorentwicklungen von u und f_{h} in (1.2.4) die erforderlichen Bedingungen an die Koeffizienten hergeleitet und damit geeignete Verfahren konstruiert werden. Für ein dreistufiges Verfahren mit Ordnung drei sind, z.B., folgende Ordnungsbedingungen zu erfüllen:

\begin{aligned}b_{1} + b_{2} + b_{3} & = & 1 \\ b_{2}c_{2} + b_{3}c_{3} & = & \frac{1}{2} \\ b_{2}c_{2}^{2} + b_{3}c_{3}^{2} & = & \frac{1}{3} \\ b_{3}a_{32}c_{2} & = & \frac{1}{6} \end{aligned} .  (1.2.7)

Dabei entsprechen die ersten drei Bedingungen denen an die Quadraturformel mit Stützstellen c_{i} und Gewichten b_{i} (Parabeln werden exakt integriert). Da das Verfahren 6 freie Parameter besitzt, gibt es also eine zweiparametrige Schar von dreistufigen Verfahren mit Ordnung drei. Analoges gilt bei 4-stufigen (8 Bedingungen für 10 Parameter). Keines dieser Verfahren hat aber eine höhere Ordnung als drei bzw. vier. Die Zahl der Ordnungsbedingungen wächst zu höheren Ordnungen stark an, sodaß dann (im Unterschied zu den Quadraturformeln) m > p Stufen für Ordnung p erforderlich werden. Die mit einer bestimmten Stufenzahl m erreichbare Ordnung p ist in der folgenden Tabelle aufgeführt:

\begin{matrix}m = & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & .. & 17 \\ p = & 1 & 2 & 3 & 4 & 4 & 5 & 6 & 6 & 7 & 7 & 8 & .. & 10 \end{matrix}

Die vorher behandelten Verfahren 1–4 besitzen jeweils die höchsten erreichbaren Ordnungen. Die höchste bisher durch explizite Konstruktion erreichte Ordnung ist p = 10 mit m = 17 Stufen. Später werden zwar Einschrittverfahren beliebig hoher Ordnung auftreten, doch werden bei diesen noch sehr viel höhere Stufenzahlen benutzt.

Um die Konvergenz eines Verfahrens bzw. der durch es berechneten Näherungen nachzuweisen, reicht die Konsistenz i.a. nicht aus. Das Verfahren muß eine weitere zentrale Eigenschaft erfüllen.

1.2.3 Stabilität

Stabilität bedeutet allgemein, daß sich die von einem Verfahren berechneten Näherungen bei kleinen Störungen im Verfahren nur wenig (Lipschitz-stetig) änderen. Wegen der einfachen Struktur des diskreten Anfangswertproblems (1.2.3) bietet sich hier aber eine etwas speziellere Formulierung an.

Definition 1.2.5 Ein Einschrittverfahren (1.2.3) heißt stabil, wenn Konstanten S, r, H_{0} existieren so, daß für je zwei gestörte Lösungen {y_{i}}, {z_{i}} mit

\begin{matrix}\frac{1}{h_{i}} [y_{i+1} − y_{i}] − f_{h_{i}}(x_{i}, y_{i}) & = & η_{i} \\ \frac{1}{h_{i}} [z_{i+1} − z_{i}] − f_{h_{i}}(x_{i}, z_{i}) & = & ζ_{i} \end{matrix} \} i = 0, 1, …, N − 1,  (1.2.8)

gilt

max_{i=0}^{N} ∥y_{i} − z_{i}∥ ≤ S(∥y_{0} − z_{0}∥ + max_{i=0}^{N−1} ∥η_{i} − ζ_{i}∥)  (1.2.9)

für alle Gitter Δ mit |Δ| ≤ H_{0} und solange ∥y_{0} − z_{0}∥ ≤ r, ∥η_{i} − ζ_{i}∥ ≤ r, i = 0, …, N − 1.

Im Gegensatz zu den später behandelten Mehrschrittverfahren gilt die Stabilität von Einschrittverfahren unter einfachen Voraussetzungen. In dieser Beziehung und auch in Bezug auf die verwendeten Beweistechniken liegen starke Ähnlichkeiten zu den Aussagen aus §1.1 vor. Umgekehrt kann übrigens auch der Existenzsatz 1.1.1 mit Hilfe des Euler-Cauchy-Verfahrens bewiesen werden. Als Beweishilfsmittel im folgenden Stabilitätssatz wird das diskrete Analogon zum Gronwall-Lemma benötigt.

Lemma 1.2.6 Für Werte δ_{i} ∈ R_{+} gelte mit α, ℓ, h_{i} ≥ 0 rekursiv

δ_{i+1} ≤ (1 + h_{i}ℓ)δ_{i} + h_{i}α, i = 0, …, N − 1.

Dann folgt (mit x_{i} := x_{0} + h_{0} + … + h_{i−1}) die explizite Schranke

δ_{i} ≤ e^{ℓ(x_{i}−x_{0})}(δ_{0} + (x_{i} − x_{0})α), i = 0, …, N.

Beweis Da 1 + ℓh_{i} ≤ e^{ℓh_{i}} ist und δ_{0} die Behauptung erfüllt, folgt induktiv

δ_{i+1} ≤ e^{ℓh_{i}}δ_{i} + h_{i}α ≤ e^{ℓh_{i}}(e^{ℓ(x_{i}−x_{0})}[δ_{0} + (x_{i} − x_{0})α]) + h_{i}α
≤ e^{ℓ(x_{i+1}−x_{0})}[δ_{0} + (x_{i} − x_{0})α + h_{i}α].

Als wesentliche Voraussetzung wird beim Stabilitätsbeweis von Einschrittverfahren wieder eine Lipschitzbedingung benötigt.

Satz 1.2.7 (Stabilitätssatz) Die das Einschrittverfahren (1.2.3) definierende Funktion f_{h} erfülle in dem in Satz 1.1.1 definierten Streifen Ω eine Lipschitzbedingung

∥f_{h}(x, y) − f_{h}(x, z)∥ ≤ ℓ∥y − z∥ ∀x ∈ [a, b], y, z ∈ R^{n}.  (1.2.10)

Dann ist das Verfahren stabil, es gilt (1.2.9) mit der (von Δ unabhängigen) Konstanten S = e^{(b−a)ℓ} max{1, b − a} und H_{0} = b − a, r = ∞.

Beweis Die Subtraktion der beiden Gleichungen in (1.2.8) ergibt

y_{i+1} − z_{i+1} = y_{i} − z_{i} + h_{i}[f_{h_{i}}(x_{i}, y_{i}) − f_{h_{i}}(x_{i}, z_{i})] + h_{i}[η_{i} − ζ_{i}].

Mit der Lipschitzbedingung folgt daraus die Ungleichungskette

∥y_{i+1} − z_{i+1}∥ ≤ ∥y_{i} − z_{i}∥ + h_{i}ℓ∥y_{i} − z_{i}∥ + h_{i}∥η_{i} − ζ_{i}∥ = (1 + h_{i}ℓ)∥y_{i} − z_{i}∥ + h_{i}∥η_{i} − ζ_{i}∥.

F’ur die Differenzen δ_{i} := ∥y_{i} − z_{i}∥ entspricht dies den Voraussetzungen von Lemma 1.2.6 mit α = max_{j} ∥η_{j} − ζ_{j}∥ und liefert mit (x_{i} − x_{0})e^{ℓ(x_{i}−x_{0})} ≤ (b − a)e^{ℓ(b−a)} die Behauptung.

Die im Satz benötigte Lipschitzbedingung der Verfahrensfunkton f_{h} folgt bei expliziten Runge-Kutta-Verfahren aus derjenigen von f. Beim Euler-Cauchy-Verfahren ist dies trivialerweise der Fall, elementar ist auch

Beispiel 1.2.8 Verfahren 2, f_{h}(x, y) = f(y + \frac{h}{2}f(y)) die Abhängigkeit von x ist unwesentlich. Aus der Bedingung ∥f(y) − f(z)∥ ≤ L∥y − z∥ folgt

∥f_{h}(y) − f_{h}(z)∥ = ∥f(y + \frac{h}{2}f(y)) − f(z + \frac{h}{2}f(z))∥
≤ L∥y + \frac{h}{2}f(y) − z − \frac{h}{2}f(z)∥ ≤ L(∥y − z∥ + \frac{h}{2}L∥y − z∥),

d.h. ℓ = L + \frac{h}{2}L^{2}.

1.2.4 Konvergenz

Die doppelte Bezeichnung der Näherungen durch y_{i} und u_{Δ}(x_{i}) wurde eingeführt, um einerseits die Schreibweise im Verfahren zu vereinfachen, andererseits aber auch die Konvergenzaussage klarer formulieren zu können. Konvergenz der Näherungslösungen kann natürlich nur bei feiner werdendem Gitter Δ, |Δ| → 0, also auch N → ∞ erwartet werden. Bezogen auf einen festen Punkt x̄ im Integrationsintervall, z.B. x̄ = b, betrifft die Aussage u_{Δ}(x̄) → u(x̄) (|Δ| → 0) Näherungswerte y mit wachsendem Index auf verschiedenen Gittern.

Satz 1.2.9 Ein konsistentes und stabiles Verfahren ist konvergent: Das Einschrittverfahren (1.2.3) sei konsistent mit der Ordnung p > 0, d.h. es gelte (1.2.6). Das Verfahren sei außerdem stabil, etwa aufgrund der Lipschitzbedingung (1.2.10) für f_{h}. Dann konvergieren die Näherungswerte {y_{i} = u_{Δ}(x_{i})} bei feiner werdendem Gitter Δ (|Δ| → 0) gegen die Lösung u des AWPs (1.1.1),(1.1.2). Die Konvergenzordnung ist mindestens p, d.h. es gibt eine Konstante K so, daß für alle |Δ| ≤ H_{0} gilt

max_{x∈Δ} ∥u_{Δ}(x) − u(x)∥ = max_{i=0}^{|Δ|} ∥y_{i} − u(x_{i})∥ ≤ K|Δ|^{p}.  (1.2.11)

Beweis Nach Definition erfüllen die Näherungen y_{i} die Beziehung

\frac{1}{h_{i}} (y_{i+1} − y_{i}) − f_{h_{i}}(x_{i}, y_{i}) = 0.

Für die Gitterwerte z_{i} := u(x_{i}) der Lösung der Dgl gilt andererseits die Konsistenzaussage (1.2.6)

\frac{1}{h_{i}} (z_{i+1} − z_{i}) − f_{h_{i}}(x_{i}, z_{i}) = T_{h_{i}}(x_{i}) = O(h_{i}^{p}).

Diese Gleichungen entsprechen denen aus der Stabilitätsdefinition, (1.2.8), mit η_{i} = 0 und ζ_{i} = T_{h_{i}}(x_{i}). Wegen y_{0} = u_{0} und h_{i}^{p} ≤ |Δ|^{p} folgt daher die Schranke (1.2.11).

Der Beweis des Satzes ist natürlich trivial, nach erfolgter Begriffbildung, da einfach zwei Definitionen zusammengefügt werden. Die eigentliche Arbeit liegt im Einzelnachweis von Konsistenz und Stabilität.

Bemerkung: Für stabile Verfahren gilt Konvergenzordnung ≥ Konsistenzordnung.

Von theoretischer Bedeutung ist die Konvergenzaussage des Satzes, da sie zeigt, daß die Näherungslösung u_{Δ} genügend schnell konvergiert. Bei einer Halbierung der Schrittweiten im Gitter kann, z.B., eine Verkleinerung des Fehlers um den Faktor 2^{−p} erwartet werden. Dies legt die Verwendung von Verfahren möglichst hoher Ordnung nahe. Allerdings hängt die Fehlerkonstante K in (1.2.11) von immer höheren Ableitungen u^{(p)} ab und begrenzt dadurch die verwendbare Ordnung. Diese Konstante K ist nicht explizit bekannt und in der Praxis kaum berechenbar, daher kann aus der Fehlerschranke keine Angabe über die tatsächliche Größe des Fehlers u_{Δ} bei gegebenem Gitter abgeleitet werden. Genausowenig kann man umgekehrt ein Gitter konstruieren, auf dem ein geforderter Fehler erreicht wird. Diese Frage wird jetzt behandelt.

1.2.5 Schrittweitensteuerung

In vielen realen Problemen kann die Gestalt der Lösungen in unterschiedlichen Bereichen des Integrationsintervalls [a, b] sehr unterschiedlich ausfallen (Beispiel Satellitenbahn: stark/schwach gekrümmt nahe/entfernt von Himmelskörpern). Ein Verfahren kann hier nur dann effizient (und auch mit akzeptablen Rundungsfehlern, s.u.) arbeiten, wenn sich die Schrittweiten an diesen Verlauf anpassen. Durch eine Inspektion von Lemma 1.2.6 sieht man, daß die Fehlerschranke (1.2.11) aus

max_{x∈Δ} ∥u_{Δ}(x) − u(x)∥ ≤ K̃ ∑_{j=0}^{N−1} h_{j}∥T_{h_{j}}∥, K̃ = e^{ℓ(b−a)}  (1.2.12)

hervorgeht. Ein Gesamtfehler der Größenordnung (b − a)K̃ε kann daher durch die Forderung

∥T_{h_{i}}(x_{i})∥ ≤^{!} ε ∀i = 0, …, N − 1  (1.2.13)

erreicht werden, vgl. Numerik I, §5.4 (adaptive Integration). Da der exakte lokale Fehler T_{h} nicht bekannt ist, arbeitet man mit einer Schätzung für diesen und nutzt das bekannte Verhalten T_{h} ≐ γh^{p} mit der Forderung (1.2.13) als Richtlinie zur lokalen, fortlaufenden Schrittweitensteuerung. Hat man nämlich nach einem mit Schrittweite h ausgeführten Integrationsschritt (eine Schätzung für) den aktuellen lokalen Fehler T_{h}, dann läßt sich die eigentlich für (1.2.13) erforderliche Schrittweite ĥ schätzen nach

Tĥ ≐ γĥ^{p} ≐ (\frac{ĥ}{h})^{p}T_{h} =^{!} ε ⇒ ĥ ≐ h ^{p} \sqrt{\frac{ε}{T_{h}}} .

Zur Schätzung des lokalen Fehlers können zwei Verfahren verschiedener Ordnung, etwa p und p + 1, verwendet werden, ihre Verfahrensfunktionen seien f_{h} und f̄_{h}. Dazu berechnet man im i-ten Schritt

\begin{matrix}y_{i+1} & = & y_{i} + hf_{h}(x_{i}, y_{i}) \\ ȳ_{i+1} & = & y_{i} + hf̄_{h}(x_{i}, y_{i}) \end{matrix} \} mit Fehler \{ \begin{matrix}O(h^{p}) \\ O(h^{p+1}) \end{matrix} .  (1.2.14)

Ist v_{i}(x) diejenige Lösung der Dgl, die die Anfangsbedingung v_{i}(x_{i}) = y_{i} erfüllt, dann gilt bezüglich dieser Lösung die Beziehung

T_{h}(x_{i}) = \frac{1}{h}[v_{i}(x_{i} + h) − y_{i}] − f_{h}(x_{i}, y_{i}) = \frac{1}{h}[v_{i}(x_{i} + h) − y_{i+1}]
= \frac{1}{h}[ȳ_{i+1} − y_{i+1}] + O(h^{p+1}) = f̄_{h}(x_{i}, y_{i}) − f_{h}(x_{i}, y_{i}) + O(h^{p+1}). (1.2.15)

Daher ist (für genügend kleine Schrittweite) die Differenz f̄_{h}(x_{i}, y_{i}) − f_{h}(x_{i}, y_{i}) eine Schätzung des lokalen Fehlers.

Die Unterschiede zwischen der Beweismethode zur Fehlerabschätzung in Satz 1.2.9 und der hier verwendeten zur Fehlerschätzung lassen sich graphisch mit den Lösungsfächern von Differentialgleichung (1.1.1) und Differenzengleichung (1.2.3) veranschaulichen. Die gestrichelten Fehlerbeiträge der beiden Ansätze sind dabei vergleichbar, aber nicht identisch.

[figure: y; x; x_{0}; x_{1}; x_{2}; x_{3}; u(x); u_{0}; y_{1}; y_{2}; y_{3}; h_{0}ζ_{0}; h_{1}ζ_{1}]

In Satz 1.2.9: Lösungsfächer Differenzengl.

mit ζ_{j} = T_{h_{j}}(x_{j}), j = 0, 1, ..

[figure: y; x; x_{0}; x_{1}; x_{2}; x_{3}; u(x); u_{0}; y_{1}; y_{2}; y_{3}; v_{1}(x); v_{2}(x)]

In (1.2.15): Lösungsfächer Dgl

Mit einer solchen Schätzung des lokalen Fehlers kann auf sehr einfache Weise eine Schrittweitensteuerung formuliert werden:

Algorithmus 1.2.10 Schrittweitensteuerung beim AWP:

x := a; y := u_{0}; h =?, ε =?

wiederhole:

y0 := y + hf_{h}(x, y); y1 := y + hf̄_{h}(x, y); q := εh/∥y1 − y0∥; falls q > 1 dann {Schritt akzeptieren}

[ x := x + h; y := y0; ] h := h ∗ ^{p} \sqrt{q}; falls x + h > b dann h := b − x;

bis x ≥ b;

In der Praxis wird meist eine geringfügig kleinere Schrittweite h := 0.9 ∗ h ∗ ^{p} \sqrt{q} genommen, da eine Schrittwiederholung wesentlich teurer wird als ein etwas zu kleiner Schritt. Außerdem ist es ratsam, einen zu großen Wechsel bei den Schrittweiten durch Zusatzabfragen zu verhindern. Zusätzlich wird meist auch mit der genaueren Lösung y1 weitergerechnet (y := y1 statt y := y0 nach x := x + h), da die lokale Fehlerschätzung nur zur Steuerung verwendet wird. Zur Wahl der Toleranz ε für den lokalen Fehler T_{h} in Relation zum gewünschten globalen Fehler ∥u_{Δ} − u∥ kann keine allgemeine Aussage gemacht werden, da die Konstante K̃ in (1.2.12) i.a. unbekannt ist.

Bei der Konstruktion von Verfahrenspaaren (1.2.14) gibt es zwei übliche Zugänge unterschiedlicher Effizienz und Allgemeinheit:

a) Eingebettete Runge-Kutta-Verfahren: Analog zum Konstruktionsprinzip der Runge-Kutta-Verfahren wird versucht, die gewünschte Näherung ȳ_{i+1} mit möglichst wenig Funktionsauswertungen zu erreichen. Dazu sucht man nach zwei Verfahren (1.2.4) mit m bzw. m + 1 Stufen, bei denen das zweite die Funktionswerte k_{i} des ersten mitverwendet,

\begin{aligned}y_{i+1} & = & y_{i} + h_{i} ∑_{j=1}^{m} b_{j}k_{j}, \\ ȳ_{i+1} & = & y_{i} + h_{i} ∑_{j=1}^{m+1} b̄_{j}k_{j}, \\ k_{j} = f(x + h_{i}c_{j}, y_{i} + h_{i} ∑_{l=1}^{j−1} a_{jl}k_{l}), j = 1, .., m + 1, \end{aligned}  \begin{aligned}c_{1} & 0 & ⋯ & 0 & 0 \\ ⋮ & & ... & & ⋮ \\ c_{m+1} & a_{jl} & & & 0 \\ & b_{1} & ⋯ & b_{m} & 0 \\ & b̄_{1} & ⋯ & · & b̄_{m+1} \end{aligned}

Die lokale Fehlerschätzung zum i-ten Schritt ist damit (und b_{m+1} := 0)

T_{h_{i}}(x_{i}) ≅ \frac{1}{h_{i}} [ȳ_{i+1} − y_{i+1}] = ∑_{j=1}^{m+1}(b̄_{j} − b_{j})k_{j}.

Mit Hilfe der Ordnungsbedingungen (1.2.7) wird kurz die Konstruktion eines eingebetteten Verfahrens zu Verfahren 3 betrachtet. Hier sind c_{2} = a_{21} = 1 vorgegeben. Für die restlichen Parameter bleiben die Bedingungen

b̄_{1} + b̄_{2} + b̄_{3} = 1, b̄_{2} + b̄_{3}c_{3} = \frac{1}{2}, b̄_{2} + b̄_{3}c_{3}^{2} = \frac{1}{3}, b̄_{3}a_{32} = \frac{1}{6}.

Für c_{3} = \frac{1}{2} ergeben sich daraus die Koeffizienten in der folgenden rechten Tabelle. Die linke enthält ein (2,3)-Paar zu Verfahren 2.

(2,3)-Paar zu Verf. 2

0 0
\frac{1}{2}  \frac{1}{2} 0
1 −1 2  0
0 1  0
\frac{1}{6}  \frac{2}{3}  \frac{1}{6}

(2,3)-Paar zu Verf. 3

0 0
1 1 0
\frac{1}{2} \frac{1}{4} \frac{1}{4}  0
\frac{1}{2} \frac{1}{2}  0
\frac{1}{6} \frac{1}{6}  \frac{2}{3}

Bekannt sind Runge-Kutta-Paare mit

Ordnungen (1,2) (2,3) (3,4) (4,5) (4,5) .. (7,8)
Stufen  (1,2) (2,3) (3,5) (4,6) (5,6) .. (11,13)

Bei den 3-stufigen Verfahren gibt es keine eingebetteten 4-stufigen Verfahren. Allerdings fand Fehlberg ein 5-stufiges Paar der Ordnung (3,4), bei dem der letzte Funktionswert f(x + c_{5}h, ...) als erster Funktionswert im nächsten Intervall, f(x + h + c_{1}h, ..), verwendet wird, das also beinahe diese Anforderungen erfüllt. Auch zum klassischen Runge-Kutta-Verfahren gibt es kein eingebettetes mit 5 Stufen, da überhaupt keine 5-stufigen Verfahren mit Ordnung 5 existieren (s.o.). Somit werden sowieso 6 Stufen nötig und es ist günstiger, das Stufenpaar (5,6) zu verwenden. Sehr effiziente Verfahren (und Implementierungen) sind das (4,5)-Paar DOPRI5 bzw. (3,5,8)-Tripel DOP853 von Dormand und Prince.

In praktischen Anwendungen benötigt man oft (z.B. Graphik) die Näherung y bzw. u_{Δ} tatsächlich als Funktion von t. Zu diesem Zweck betrachtet man stetige Erweiterungen. Dabei sind die Koeffizienten in (1.2.4) Polynome b̄_{i}(t) mit b̄_{i}(0) = 0, b̄_{i}(1) = b_{i}, und durch

u_{Δ}(x_{i} + th_{i}) = y_{i} + h_{i} ∑_{j=1}^{m} b_{j}(t)k_{j},

wird eine stetige Funktion u_{Δ} definiert. Für das Verfahren 3 etwa ergibt sich mit

b_{1}(t) = t − \frac{1}{2}t^{2}, b_{2}(t) = \frac{1}{2}t^{2}

eine stetige Erweiterung u_{Δ}(x) der gleichmäßigen Konvergenzordnung 2. Im allgemeinen kann man aber für stetige Erweiterungen nicht die volle normale Ordnung der Endformel mit b̄_{i}(1) erwarten. Für 3-stufige Verfahren ist die glm Ordnung 2 erreichbar, zu DOPRI5 existiert eine Erweiterung mit glm Ordnung 4.

b) Richardson-Extrapolation: Bei vielen Verfahren, z.B. auch den Runge-Kutta-Verfahren, existiert für äquidistante Gitter Δ = {x_{i} = a + ih} eine asymptotische Entwicklung des Fehlers

u_{Δ}(x) − u(x) = h^{p}g_{p}(x) + O(h^{p+1}), h → 0,  (1.2.16)

vgl. Numerik I, §5.5. Durch Linearkombination von zwei Näherungen zu Schrittweiten h und h/2 kann der dominierende Fehleranteil h^{p}g_{p} eliminiert (→ Richardson-Extrapolation) oder aber geschätzt werden, da für eine zweite Lösung u_{Δ′}, zu Δ′ = \{x_{i} = a + i\frac{h}{2} : i = 0, …, 2N\}, gilt

u_{Δ′}(x) − u(x) = 2^{−p}h^{p}g_{p}(x) + O(h^{p+1}) ⇒
(1 − 2^{−p})h^{p}g_{p}(x) = u_{Δ}(x) − u_{Δ′}(x) + O(h^{p+1}).

Daraus kann eine bessere Approximation ū_{Δ}(x) := u_{Δ′}(x) + \frac{1}{2^{p}−1}[u_{Δ′}(x) − u_{Δ}(x)], x ∈ Δ, oder die Fehlerschätzungen

u_{Δ} − u ≐ \frac{2^{p}}{2^{p} − 1}[u_{Δ} − u_{Δ′}], u_{Δ′} − u ≐ \frac{1}{2^{p} − 1}[u_{Δ} − u_{Δ′}],

berechnet werden. Wie in (1.2.15) ergibt sich

T_{h}(x_{i}) = \frac{1}{h}[v_{i}(x_{i} + h) − u_{Δ}(x_{i} + h)] ≐ \frac{1}{h}[ū_{Δ}(x_{i} + h) − u_{Δ}(x_{i} + h)]
= \frac{1}{1 − 2^{−p}} \frac{1}{h}[u_{Δ′}(x_{i} + h) − u_{Δ}(x_{i} + h)] ≐ −h^{p−1}g_{p}(x_{i} + h).  (1.2.17)

Diese Größe ist im speziellen Fall von der Ordnung O(h^{p}). Da man sich hier in (1.2.16) auf (Näherungs-)Lösungen bezieht mit u_{Δ}(x_{i}) = v_{i}(x_{i}) = y_{i}, gilt g_{p}(x_{i}) = 0 und somit nach dem Mittelwertsatz h^{p−1}g_{p}(x_{i} + h) = O(h^{p}). Daher ist (1.2.17) eine Fehlerschätzung von der in Algorithmus 1.2.5 verwendeten Gestalt. Im Vergleich zu eingebetteten Verfahren werden hier aber wesentlich mehr Zusatzfunktionsauswertungen verwendet (bezogen auf die genauere Lösung u_{Δ′} nach 2 Schritten mit h/2 also m/2 pro Schritt).

Beispiel 1.2.11 Periodische Satellitenbahn um Erde und Mond: In der Ebene des rotierenden Erde-Mond-System bewegt sich ein Satellit mit den Koordinaten u(t) ∈ R^{2} nach der Dgl

u′′ = u + ( \begin{matrix}0 & 2 \\ −2 & 0 \end{matrix} ) u′ − ν \frac{u − E}{∥u − E∥^{3}} − μ\frac{u − M}{∥u − M∥^{3}} ,

wobei μ = 1/82.45, ν = 1 − μ die Massenverhältnisse und E = (−μ, 0)T bzw. M = (ν, 0)T die Positionen von Erde und Mond sind. Es gibt eine periodische Bahn, die nahe an der Erde vorbeiführt und dort sehr kleine Schrittweiten erfordert. Ergebnis bei ε =_{10} −4:

Anz.f-Werte h_{min} h_{max} Endfehler
RuKuFe-2/3  35805 2_{10}-6 8_{10}-2  2_{10}-5
RuKuFe-4/5  972 4_{10}-4 6_{10}-1  9_{10}-3

Schrittweitensteuerung hatte zunächst nur das Ziel, die geforderte Genauigkeit der Näherung mit möglichst wenig Aufwand zu erreichen. Zusätzlich reduziert sie aber auch den Einfluß von Rundungsfehlern. Bei Maschinenrechnung werden in (1.2.3) gestörte Näherungen ỹ_{i} berechnet, die pro Schritt anfallenden Rundungsfehler ε_{i} führen also zu einer verfälschten Rekursion

ỹ_{i+1} = ỹ_{i} + h_{i}f_{h_{i}}(x_{i}, ỹ_{i}) + ε_{i+1}, i = 0, …, N − 1,

ỹ_{0} = u_{0}. Daher kommen bei praktischer Rechnung zum Diskretisierungsfehler y_{i} − u(x_{i}) noch Rundungsfehler hinzu. Unter Voraussetzung (1.2.10) gilt dafür

∥ỹ_{i+1} − y_{i+1}∥ ≤ ∥ỹ_{i} − y_{i}∥ + h_{i}∥f_{h_{i}}(x_{i}, ỹ_{i}) − f_{h_{i}}(x_{i}, y_{i})∥ + ∥ε_{i+1}∥
≤ (1 + h_{i}ℓ)∥ỹ_{i} − y_{i}∥ + ε̂_{i+1}, ε̂_{i+1} := ∥ε_{i+1}∥.

Analog zu Lemma 1.2.6 ergibt sich mit ỹ_{0} = y_{0} die Abschätzung

∥ỹ_{i} − y_{i}∥ ≤ ∑_{j=1}^{i} e^{ℓ(x_{i}−x_{j})}ε̂_{j} ≤ S ∑_{j=1}^{N} ε̂_{j}, S := e^{ℓ(b−a)}.  (1.2.18)

Bei einem äquidistanten Gitter führt dies für den Gesamtfehler auf die Schranke

∥ỹ_{i} − u(x_{i})∥ ≤ S(b − a)\frac{ε}{h} + Kh^{p}, ε := max_{i} ε̂_{i}.

Bei äquidistantem Gitter muß die (zunächst unbekannte) Schrittweite h nach der ungünstigsten Stelle im Intervall gewählt werden und erzwingt daher evtl. eine sehr große Zahl N von Schritten. Da der Rundungsfehler nach (1.2.18) im wesentlichen von dieser Schrittezahl N abhängt, können

[figure: 1/h; Gesamtfehler]

diese bei variablem Gitter, das in leichten Passagen große und daher weniger Schritte verwendet, reduziert werden. Zusätzlich sei daran erinnert (vgl. Numerik I, §7.2), daß der kleinstmögliche Wert in der Schranke für ∥ỹ_{i}−u(x_{i})∥ von der Größe O(ε^{p/(p+1)}) ist und mit wachsender Ordnung kleiner wird.

1.3 Mehrschrittverfahren

1.3.1 Adams-Verfahren

Bei Einschrittverfahren der Ordnung p sind m ≥ p Stufen erforderlich und daher zur Approximation des Integrals in (1.2.1),

∫_{0}^{1} f(x_{i−1} + ht, u(x_{i−1} + ht)) dt = \frac{1}{h}[u(x_{i−1} + h) − u(x_{i−1})]

beim Schritt x_{i−1} → x_{i} jeweils mindestens p Funktionswerte k_{j} ≅ f(x_{i−1}+hc_{j}, u(x_{i−1}+hc_{j})) neu zu berechnen. Eine Alternative ist, stattdessen schon bekannte, vor x_{i−1} gelegene Funktionswerte f(x_{i−j}, y_{i−j}) =: f_{i−j}, j = 1, 2, .., m, heranzuziehen.

Durch Integration des Interpolationspolynoms p_{m−1}(x) vom Grad m − 1 zu den Funktionswerten f_{i−1}, ..., f_{i−m}, das mit den Lagrangepolynomen L_{j}(x) zu den Stützstellen x_{i−1}, .., x_{i−m} dargestellt werden kann (vgl. Numerik I, §3.2), ergibt sich für dieses Integral die Approximation

[figure: f; x_{i−m}; x_{i−1}; f_{i−1}; x_{i}]
f_{h}(x_{i−1}, y_{i−m}, …, y_{i−1}) := ∑_{j=−m}^{−1} f_{i+j} \underbrace{∫_{0}^{1} L_{j}(x_{i−1} + ht) dt}_{=:β_{−j}} .  (1.3.1)

Diese Verfahrensfunktion hängt also von m früheren Lösungswerten y ab. Speziell für äquidistante Punkte (d.h. h_{i} ≡ h) erhält man damit die Methode von Adams-Bashforth (A-B):

\begin{matrix}\frac{1}{h}(y_{i} − y_{i−1}) & = & f_{h}(x_{i−1}, y_{i−m}, .., y_{i−1}) \\ f_{h}(x_{i−1}, y_{i−m}, .., y_{i−1}) & := & ∑_{j=1}^{m} β_{j}f(x_{i−j}, y_{i−j}) \end{matrix} \} i = m, m + 1, .., N.  (1.3.2)

Dieses Verfahren ist ein explizites, lineares, m-Schrittverfahren. Pro Schritt fällt hier tatsächlich nur eine neue Funktionsauswertung, f_{i−1}, an. Explizit heißt das Verfahren wieder, da f_{h} nicht vom unbekannten Wert y_{i} abhängt.

Der Konsistenzfehler entspricht dem Quadraturfehler bei der Funktion u′(x), das Verfahren hat daher Ordnung p = m. Die Koeffizienten der ersten A-B-Verfahren enthält die Tabelle.

m  j = 1 2 3 4 Fehler
1 β_{j} = 1 \frac{1}{2}hu′′
2 2β_{j} = 3 −1  \frac{5}{12}h^{2}u′′′
3 12β_{j} = 23 −16 5  \frac{3}{8}h^{3}u^{(4)}
4 24β_{j} = 55 −59 37 −9 \frac{251}{720}h^{4}u^{(5)}

Bei der Approximation (1.3.1) werden nur Daten auf einer Seite des Integrationsintervalls [x_{i−1}, x_{i}] benutzt. Ohne Erhöhung der Zahl m der Teilintervalle kann das A-B-Verfahren durch Berücksichtigung des Wertes im Endpunkt x_{i} verbessert werden, man interpoliert also in den Stellen x_{i−m}, .., x_{i} mit einem Polynom vom Grad m. Dies führt auf die Methode von Adams-Moulton (A-M)

\begin{matrix}\frac{1}{h}(y_{i} − y_{i−1}) & = & f_{h}(x_{i−1}, y_{i−m}, .., y_{i}) \\ f_{h}(x_{i−1}, y_{i−m}, .., y_{i}) & := & ∑_{j=0}^{m} μ_{j}f(x_{i−j}, y_{i−j}) \end{matrix} \} i = m, m + 1, .., N.  (1.3.3)

Hier liegt ein implizites, lineares, m-Schritt-Verfahren vor, da (1.3.3) für y_{i} ein (nichtlineares) Gleichungssystem darstellt. Die Koeffizienten der ersten A-M-Verfahren lauten

m  j = 0 1 2 3 Fehler
0 μ_{j} = 1 −\frac{1}{2}hu′′
1 2μ_{j} = 1 1  −\frac{1}{12}h^{2}u′′′
2 12μ_{j} = 5 8 −1  −\frac{1}{24}h^{3}u^{(4)}
3 24μ_{j} = 9 19 −5 1 −\frac{19}{720}h^{4}u^{(5)}

Die Konsistenz von beiden Verfahrensklassen folgt aus den Approximationseigenschaften der Interpolationspolynome. Beide Verfahrenstypen werden in einem Satz behandelt. Dazu wird das A-B-Verfahren als A-M-Verfahren mit Gewicht β_{0} = 0 interpretiert.

Satz 1.3.1 Der lokale Fehler bei den Adams-Bashforth-Verfahren (1.3.2), γ_{0} := 0, γ_{j} := β_{j}, 1 ≤ j ≤ m, bzw. Adams-Moulton-Verfahren (1.3.3), γ_{j} := μ_{j}, 0 ≤ j ≤ m, ist gegeben durch

T_{h}(x_{i−1}) := \frac{1}{h}[u(x_{i}) − u(x_{i−1})] − ∑_{j=0}^{m} γ_{j}u′(x_{i−j}).

Er hat die Form

T_{h}(x_{i−1}) = \{ \begin{matrix}c_{m}^{B} h^{m} u^{(m+1)}(x_{i−1}) + O(h^{m+1}) & bei A-B \\ c_{m}^{M}h^{m+1}u^{(m+2)}(x_{i−1}) + O(h^{m+2}) & bei A-M \end{matrix} .  (1.3.4)

Beweis Beim A-B-Verfahren sei p das Interpolationspolynom vom Grad m − 1 zu den Werten (x_{i−j}, u′(x_{i−j})), j = 1, .., m. Dann folgt aus der Fehlerformel der Interpolation (Numerik I, 3.3.11) und durch Taylorentwicklung die Aussage

T_{h}(x_{i−1}) = ∫_{0}^{1}[u′(x_{i−1} + ht) − p(x_{i−1} + ht)] dt
= \frac{1}{m!} ∫_{0}^{1}(x_{i−1} + ht − x_{i−1})..(x_{i−1} + ht − x_{i−m})u^{(m+1)}(ξ(t)) dt
= u^{(m+1)}(x_{i−1})h^{m} \underbrace{\frac{1}{m!} ∫_{0}^{1} t(t + 1)..(t + m − 1) dt}_{=: c_{m}^{B}} +O(h^{m+1}).

Bei A-M betrachtet man das Polynom vom Grad m zu (x_{i−j}, u′(x_{i−j})), j = 0, .., m.

Die Fehlerkoeffizienten c_{m} sind in den obigen Tabellen angegeben. Das A-M-Verfahren mit m Stufen besitzt die gleiche Konvergenzordnung h^{m+1} wie das A-B-Verfahren mit m + 1 Stufen, hat aber kleinere Fehlerkoeffizienten |c_{m}^{M}| < |c_{m+1}^{B}|, m > 0. Daher ist das Adams-Moulton-Verfahren vorzuziehen. Die Auflösung der hier auftretenden nichtlinearen Gleichung für y_{i} in (1.3.3) kann durch Kombination mit dem A-B-Verfahren umgangen werden, indem mit A-B zunächst eine Approximation ỹ_{i} für den unbekannten Wert berechnet wird. So ergibt sich das Prädiktor-Korrektor-Verfahren (P-K)

\begin{matrix}ỹ_{i} & := & y_{i−1} + h ∑_{j=1}^{m} β_{j}f(x_{i−j}, y_{i−j}) \\ y_{i} & := & y_{i−1} + hμ_{0}f(x_{i}, ỹ_{i}) + h ∑_{j=1}^{m} μ_{j}f(x_{i−j}, y_{i−j}) \end{matrix} \} i = m, m + 1, …  (1.3.5)

Hier werden pro Schritt zwei Funktionsaufrufe benutzt. Der lokale Fehler dieses Verfahren ist tatsächlich von der Ordnung h^{m+1}. Denn nach Satz 1.3.1 gilt für Lipschitz-stetiges f mit ũ_{i} := u(x_{i−1}) + h ∑_{j=1}^{m} β_{j}f(x_{i−j}, u(x_{i−j})) hier

∥T_{h}^{P K}∥ = ∥T_{h}^{AM} + μ_{0}(u′(x_{i}) − f(x_{i}, ũ_{i}))∥
≤ ∥T_{h}^{AM}∥ + |μ_{0}|L∥u(x_{i}) − ũ_{i}∥ = ∥T_{h}^{AM}∥ + |μ_{0}|Lh∥T_{h}^{AB}∥
= h^{m+1}(|c_{m}^{M}|∥u^{(m+2)}(x_{i−1})∥ + |μ_{0}Lc_{m}^{B}|∥u^{(m+1)}(x_{i−1})∥) + O(h^{m+2}).

Der Fehler bei Runge-Kutta-Verfahren besitzt keine so einfache Gestalt. Im allgemeinen ist er aber für ein Verfahren gleicher Ordnung kleiner als bei Mehrschrittverfahren. Für das Problem u′ = u, u(0) = 1, etwa ergeben sich bei den Verfahren vierter Ordnung nach einem Schritt

Verfahren: Runge-Kutta Adams-Bashf. Adams-Moulton
lok. Fehler: \frac{1}{120}h^{5} \frac{251}{720}h^{5} \frac{19}{720}h^{5}

Bei den Mehrschrittverfahren sind für gleiche Fehlerbeträge also kleinere Schrittweiten erforderlich. Außerdem ist die Programmierung aufwendiger.

Fragen zum praktischen Einsatz:

• Anlaufrechnung: Die m-Schritt-Verfahren können, wie aus (1.3.2)-(1.3.5) ersichtlich, erstmals nach einer Anlaufphase von m Integrationsschritten eingesetzt werden. Daher sind Startwerte y_{0}, …, y_{m−1} ausreichender Genauigkeit auf andere Weise bereitzustellen, entweder mit Runge-Kutta-Verfahren oder Mehrschrittverfahren niederer (wachsender) Ordnung und sehr kleinen Schrittweiten.

• Schrittweitenwechsel: Die Adams-Formeln (1.3.2)-(1.3.5) beziehen sich auf äquidistante Punkte x_{i−m}, .., x_{i} im konstanten Abstand h. Die Koeffizienten β_{j}, μ_{j} ergaben sich dabei durch Integration der Lagrange-Darstellung des Interpolationspolynoms, bei Adams-Bashforth also von p(x) = ∑_{j=1}^{m} L_{j}(x)f_{i−j} in der Form

f_{h}^{AB} = ∑_{j=1}^{m} β_{j}f_{i−j} = ∫_{0}^{1} ∑_{j=1}^{m} L_{j}(x_{i−1} + ht)f_{i−j} dt.

Wird p dagegen als Newton-Polynom (Numerik I, §3.3) dargestellt, gilt

f_{h}^{AB} = ∑_{j=0}^{m−1} h^{j}δ_{j}f[x_{i−j−1}, .., x_{i−1}].

Ein Wechsel der Schrittweite h → ĥ ist in dieser Form jetzt einfach möglich durch Ersetzung der Faktoren h^{j} → ĥ^{j}. Die schon berechneten Differenzen f[.., x_{i−1}] werden übernommen. In der Praxis sind dabei zusätzliche Vorsichtsmaßnahmen zu treffen.

• Schrittweiten-, Ordnungs-Steuerung: Da die A-B-Verfahren Ordnung m und die A-M-Verfahren Ordnung m+1 besitzen, sind die beiden Teile des Prädiktor-Korrektor-Verfahrens (1.3.5) als Verfahrenspaar zur Schätzung des lokalen Fehlers einsetzbar. Algorithmus 1.2.10 könnte damit also auch hier direkt eingesetzt werden (vgl. letzte Anmerkung). Darüberhinaus läßt sich aber auch noch die Ordnung des Verfahrens steuern, da es Adams-Methoden beliebiger Ordnung gibt. Eine einfache Strategie hierfür erhält man aus der Regel:

Ist ein Schritt in [x_{i−1}, x_{i}] mit einem Verfahren der Ordnung m akzeptiert worden, so wähle für [x_{i}, x_{i+1}] dasjenige der drei Verfahren der Ordnungen m − 1, m, m + 1, das die größte Schrittweite ĥ gestattet.

Die erwähnten Maßnahmen ermöglichen auch bei Mehrschrittverfahren eine selbststeuernden Einsatz. Allerdings ist der damit verbundene Verwaltungsaufwand erheblich höher als bei Einschrittverfahren und muß bei einem Vergleich berücksichtigt werden (vgl. §1.4).

1.3.2 Lineare Mehrschrittverfahren und Stabilität

In den Verfahren (1.3.2,1.3.3) werden zwar die Funktionswerte f_{i−j} aus m zurückliegenden Punkten benutzt (’rechte Seite’), Lösungswerte dagegen nur von zwei Stellen, y_{i}, y_{i−1}. Ohne Erhöhung der Schrittzahl m kann man daher eine Verbesserung durch Einführung zusätzlicher Parameter erwarten. Dazu werden jetzt allgemeinere Verfahren der Form

∑_{j=0}^{m} α_{j}y_{i−j} = h ∑_{j=0}^{m} β_{j}f_{i−j} , i = m, m + 1, …,  (1.3.6)

betrachtet mit x_{j} = a + jh, f_{j} := f(x_{j}, y_{j}) und α_{0} = 1. Für β_{0} = 0 ist dies ein explizites, andernfalls ein implizites, lineares m-Schrittverfahren.

Eines der Ziele bei der Konstruktion von Verfahren ist eine möglichst hohe Ordnung des lokalen bzw. Konsistenzfehlers

T_{h}(x_{i−1}) = \frac{1}{h} ∑_{j=0}^{m} α_{j}u(x_{i−j}) − ∑_{j=0}^{m} β_{j}u′(x_{i−j}).  (1.3.7)

Zur Überprüfung werden die Taylorentwicklungen in beiden Summen verglichen. Mit

u(x_{i} − jh) = ∑_{k=0}^{p} \frac{(−jh)^{k}}{k!} u^{(k)}(x_{i}) + O(h^{p+1})

u′(x_{i} − jh) = ∑_{k=0}^{p−1} \frac{(−jh)^{k}}{k!} u^{(k+1)}(x_{i}) + O(h^{p})

in (1.3.7) ergibt sich

T_{h}(x_{i−1}) = ∑_{j=0}^{m} [ ∑_{k=0}^{p} α_{j} \frac{(−j)^{k}}{k!} h^{k−1}u^{(k)} − ∑_{k=0}^{p−1} β_{j} \frac{(−j)^{k}}{k!} h^{k}u^{(k+1)}] + O(h^{p})
= \frac{1}{h}u(x_{i})( ∑_{j=0}^{m} α_{j}) − ∑_{k=0}^{p−1} \frac{(−h)^{k}}{k!} u^{(k+1)}(x_{i}) ∑_{j=0}^{m} (α_{j} \frac{j^{k+1}}{k + 1} + β_{j}j^{k}) + O(h^{p}).

In dieser Entwicklung nach h-Potenzen liest man die Ordnungsbedingungen direkt ab:

Satz 1.3.2 Das Verfahren (1.3.6) besitzt Konsistenzordnung p genau dann, wenn seine Koeffizienten das folgende lineare Gleichungssystem erfüllen, α_{0} = 1,

\begin{matrix}∑_{j=0}^{m} α_{j} = 0, \\ ∑_{j=0}^{m} (\frac{j^{k+1}}{k+1} α_{j} + j^{k}β_{j}) = 0, & k = 0, …, p − 1. \end{matrix}  (1.3.8)

Ordnungsbedingungen sind hier also wesentlich leichter aufzustellen und Lösungen zu berechnen als bei den Einschrittverfahren. Man könnte also versuchen, durch geeignete Wahl der 2m(+1) freien Parameter α_{j}, β_{j} eine möglichst hohe Ordnung zu erreichen.

Beispiel 1.3.3 Konstruktion eines expliziten (β_{0} = 0) 2-Schrittverfahrens mit maximaler Ordnung p = 3.

\begin{matrix}1 + α_{1} + α_{2} & & = & 0 \\ α_{1} + 2α_{2} & +β_{1} + β_{2} & = & 0 \\ \frac{1}{2}(α_{1} + 4α_{2}) & +β_{1} + 2β_{2} & = & 0 \\ \frac{1}{3}(α_{1} + 8α_{2}) & +β_{1} + 4β_{2} & = & 0 \end{matrix} \} ⟹ \begin{matrix}(α_{0}, α_{1}, α_{2}) = (1, 4, −5) \\ (β_{0}, β_{1}, β_{2}) = (0, 4, 2) \end{matrix}

Das Verfahren lautet also

y_{i} + 4y_{i−1} − 5y_{i−2} = h(4f_{i−1} + 2f_{i−2}), i ≥ 2.

Einfaches Zahlenbeispiel: Bei der Dgl u′ = −u, u(0) = 1, mit exakten Startwerten y_{0} := 1, y_{1} = e^{−h}, Berechnung der Näherungen y_{i} := −4(1 + h)y_{i−1} + (5 − 2h)y_{i−2} mit h = 0.2. Numerisches Ergebnis: Die Fehler in den Approximationen oszillieren und wachsen ungefähr um einen Faktor 5 pro Schritt.

\begin{matrix}i & y_{i} \\ 1 & 0.8187 \\ 2 & 0.6701 \\ 3 & 0.5497 \\ 4 & 0.4438 \\ 5 & 0.3986 \\ 6 & 0.1283 \\ 7 & 1.2177 \\ 8 & -5.2551 \\ 9 & 30.8262 \end{matrix}

Das Verfahren ist instabil, Störungen wachsen extrem an. Im Gegensatz zu den Einschrittverfahren ist die Stabilitätsforderung eine echte Einschränkung für Mehrschrittverfahren!

Es gibt sogar eine sogenannte Ordnungsbarriere (von Dahlquist): Ein stabiles m-Schritt-Verfahren besitzt eine (Konsistenz-)

Ordnung p ≤ \{ \begin{matrix}m + 1, wenn & m ungerade \\ m + 2, wenn & m gerade \end{matrix} .

Gegenüber den Adams-Verfahren kann die Ordnung also höchstens um 1 erhöht werden, wenn m gerade ist. Die zugehörigen Verfahren besitzen eine den Adams-Verfahren ähnliche Struktur,

y_{i} − y_{i−k} = h ∑_{j=0}^{m} β_{j}f_{i−j}, i ≥ m, 1 ≤ k ≤ m.  (1.3.9)

Der folgende Konvergenzsatz enthält für diese eine zu den Einschrittverfahren ähnliche Stabilitätsaussage. Auf die Stabilitätsforderungen an allgemeine lineare Mehrschrittverfahren wird kurz noch einmal in §1.4 eingegangen.

Satz 1.3.4 Die rechte Seite f der Dgl erfülle eine Lipschitzbedingung (1.1.5) mit der Konstanten L. Die Schrittweite im Verfahren sei h ≤ h_{0} := 1/(2L|β_{0}|) (:= ∞ für β_{0} = 0). Für die Startwerte gelte ∥y_{i} − u(x_{i})∥ ≤ φ(h), i = 0, .., m − 1, und für die Konsistenzfehler ∥T_{h}(x_{i−1})∥ ≤ τ(h), i = m, …, N. Dann gilt die Fehlerschranke

∥y_{i} − u(x_{i})∥ ≤ e^{λ(x_{i}−x_{m−1})}[φ(h) + 2(x_{i} − x_{m−1})τ(h)], i = m, …, N,  (1.3.10)

mit λ = L ∑_{j=0}^{m} |β_{j}|/(1 − hL|β_{0}|) ≤ 2L ∑_{j=0}^{m} |β_{j}|. Für genügend oft differenzierbares f und τ(h) = O(h^{p}), φ(h) = O(h^{p}) ist also insbesondere auch max_{Δ} ∥u_{Δ} − u∥ = O(h^{p}), h → 0.

Beweis Für hL|β_{0}| < 1 ist die Gleichung (1.3.9) nach y_{i} eindeutig auflösbar (Banachscher Fixpunktsatz). Wieder durch Subtraktion der Gleichungen (1.3.9) und (1.3.7),

y_{i} − y_{i−k} = h ∑_{j=0}^{m} β_{j}f(x_{i−j}, y_{i−j})
u(x_{i}) − u(x_{i−k}) = h ∑_{j=0}^{m} β_{j}f(x_{i−j}, u(x_{i−j})) + hT_{h}(x_{i−1})

folgt für die Fehler ε_{j} := y_{j} − u(x_{j}) die Beziehung

∥ε_{i} − ε_{i−k}∥ = h∥ ∑_{j=0}^{m} β_{j}[f(x_{i−j}, y_{i−j}) − f(x_{i−j}, u(x_{i−j}))] − hT_{h}(x_{i−1})∥
≤ hL ∑_{j=0}^{m} |β_{j}|∥ε_{i−j}∥ + hτ(h).

Daraus erhält man die Schranke

(1 − hL|β_{0}|)∥ε_{i}∥ ≤ ∥ε_{i−k}∥ + hL ∑_{j=1}^{m} |β_{j}| max_{j=1}^{m} ∥ε_{i−j}∥ + hτ(h).

Da auf der rechten Seite viele ältere Fehlerwerte auftreten, liegt es nahe deren Maximum

δ_{j} := max_{ν≤j} ∥ε_{ν}∥, j = m − 1, .., N,

zu betrachten. Nach Division durch 1 − hL|β_{0}| und Berücksichtigung der Identität 1/(1 − hL|β_{0}|) = 1 + hL|β_{0}|/(1 − hL|β_{0}|) ergibt sich so

∥ε_{i}∥ ≤ (1 + hλ)δ_{i−1} + \frac{h}{1 − hL|β_{0}|}τ(h), i = m, .., N.  (1.3.11)

Da nun δ_{i} = max{∥ε_{i}∥, δ_{i−1}} und die rechte Seite größer als δ_{i−1} ist, kann die linke Seite in (1.3.11) durch δ_{i} ersetzt werden. Dann liegt aber die in Lemma 1.2.6 behandelte Situation vor, die auf die Schranke (1.3.10) führt. Als Startpunkt ist dabei allerdings x_{m−1} zu nehmen.

Die stabilen Verfahren höchster Ordnung m + 2, m gerade, ergeben sich bei der Wahl k = m in (1.3.9) aus den abgeschlossenen Newton-Cotes-Quadraturformeln zum Intervall [x_{i−m}, x_{i}],

\frac{1}{h}[y_{i} − y_{i−m}] = ∑_{j=0}^{m} β_{j}f_{i−j} ≅ \frac{1}{h} ∫_{x_{i−m}}^{x_{i}} f(x, u(x)) dx.  (1.3.12)

Bei diesen Verfahren ist allerdings, im Gegensatz zu den Adams-Verfahren, die mit x wachsende Stabilitätsschranke (1.3.10) scharf, selbst dann, wenn im Problem nur exponentiell fallende Lösungen auftreten können (vgl. §1.4). Daher werden in der Praxis doch meist die Adams-Verfahren verwendet.

Aufgrund ihrer Symmetrie bzgl. \frac{1}{2}(x_{i−m} + x_{i}) besitzen Verfahren der Form (1.3.12) aber noch eine weitere Eigenschaft, die von großer praktischer Bedeutung ist. Das einfachste, explizite Verfahren dieser Form ergibt sich aus der Rechteckregel zum Intervall [−2h, 0], seine Symmetrie ist erkennbar in der Form

y_{i+1} − y_{i−1} = 2hf(x_{i}, y_{i}), i = 1, 2, ..

Diese explizite Mittelpunktregel hat nicht nur Ordnung 2, sondern besitzt sogar eine asymptotische Entwicklung nach Potenzen von h^{2}. Daher können mit Richardson-Extrapolation (vgl. Numerik I, §5.5) Näherungen beliebig hoher Ordnung berechnet werden. Das resultierende Gesamt-Verfahren ist vom Typus her aber wieder ein Einschrittverfahren und wird daher in einem eigenen Abschnitt behandelt.

1.4 Extrapolationsverfahren

Die explizite Mittelpunktregel benötigt zur Durchführung noch einen Näherungswert y_{1}. Dieser kann ohne Verlust der Gesamtordnung 2 mit dem Eulerverfahren (Verfahren 1 in §1.2.1) bestimmt werden. Das betrachtete Verfahren ist daher

y_{1} − y_{0} = hf(x_{0}, y_{0}), y_{0} = u_{0},
y_{i+1} − y_{i−1} = 2hf(x_{i}, y_{i}), i = 1, 2, ….  (1.4.1)

Der erste Schritt führt eine Asymmetrie ein, die von der Mittelpunktregel als oszillierende Störung weitergegeben wird. Daher besitzt die Näherungslösung u_{Δ} eine ungewöhnliche asymptotische Entwicklung der Gestalt

u_{Δ}(x_{j}) − u(x_{j}) = ∑_{k=1}^{q} h^{2k}[g_{2k}(x_{j}) + (−1)^{j}g̃_{2k}(x_{j})] + O(h^{2q+2}), h → 0,  (1.4.2)

mit vom Gitter unabhängigen Koeffizientenfunktionen g_{2k}, g̃_{2k}. Insbesondere haben also Punkte x_{j} mit geradem und ungeradem Index verschiedene Entwicklungen und dürfen daher nicht gleichzeitig zur Extrapolation herangezogen werden.

Beispiel 1.4.1 Beim Problem u′ = λu, u(0) = 1 lautet das Verfahren y_{0} = 1, y_{1} = 1 + hλ, y_{i+1} = y_{i−1} + 2hλy_{i}. Wie bei den linearen Dgln führt der Ansatz y_{j} = cξ^{j} in der Differenzengleichung y_{i+1} − 2hλy_{i} − y_{i−1} = 0 auf die quadratische Gleichung

ξ^{2} − 2hλξ − 1 = 0 ⇒ ξ_{1/2} = hλ ± \sqrt{1 + h^{2}λ^{2}} = \{ \begin{matrix}\sqrt{1 + h^{2}λ^{2}} + hλ & > 0 \\ −(\sqrt{1 + h^{2}λ^{2}} − hλ) & < 0 \end{matrix} .

Die y_{j} besitzen daher die Darstellung y_{j} = C_{1}ξ_{1}^{j} + C_{2}ξ_{2}^{j}, wobei die Koeffizienten C_{1}, C_{2} aus den Anfangsbedingungen y_{0} = 1, y_{1} = 1 + hλ zu bestimmen sind. Um den Bezug auf das Gitter herzustellen, wird der Gitterindex bei y_{j} = u_{Δ}(jh) in der Form j = x/h geschrieben. So ergibt sich die Näherungslösung explizit in der Form

2u_{Δ}(x) = (1 + \frac{1}{\sqrt{1 + h^{2}λ^{2}}} )(\sqrt{1 + h^{2}λ^{2}} + hλ)^{x/h}

+(−1)^{x/h} \frac{h^{2}λ^{2}}{1 + h^{2}λ^{2} + \sqrt{1 + h^{2}λ^{2}}} (\sqrt{1 + h^{2}λ^{2}} − hλ)^{x/h}.

Man sieht leicht, daß mit Ausnahme von (−1)^{x/h} alle Ausdrücke eine Talyorentwicklung nach h um h = 0 besitzen. Daraus folgt hier für den Anfang der Entwicklung (1.4.2)

u_{Δ}(x) = e^{λx} + h^{2}[ \underbrace{−λ^{2}(\frac{1}{6}λx + \frac{1}{4})e^{λx}}_{g_{2}} +(−1)^{x/h} \underbrace{\frac{1}{4}λ^{2}e^{−λx}}_{g̃_{2}} ] + O(h^{4}), h → 0.

Wenn die Lösung der Dgl fällt, λ < 0, kann der oszillierende und wachsende Störterm g̃_{2} die numerische Lösung erheblich beeinträchtigen. Dieses Problem läßt sich etwas durch einen zusätzlichen Glättungsschritt abschwächen.

Zur Extrapolation bei der Mittelpunktregel wählt man eine wachsende Folge gerader Zahlen 2 ≤ N_{0} < … < N_{q} und bildet Schrittweiten H_{j} := H/N_{j} sowie Gitter x_{i}^{[j]} := x_{0} + iH_{j} (i = 0, …, N_{j}). Dies ergibt das Verfahren von Gragg-Bulirsch-Stoer:

für j = 0, …, q :
1 y_{0}^{[j]} := u_{0}, y_{1}^{[j]} := y_{0} + H_{j}f(x_{0}, y_{0}), Startwerte
für i = 1, .., N_{j}:  Mittelpunkt-
2  y_{i+1}^{[j]} := y_{i−1}^{[j]} + 2H_{j}f(x_{i}^{[j]}, y_{i}^{[j]})  -Regel
3  P_{j,0} := \frac{1}{4}(y_{N_{j}−1}^{[j]} + 2y_{N_{j}}^{[j]} + y_{N_{j}+1}^{[j]}),  Glättung
für k = 1, .., j: i := j − k;  Extrapo-
4  P_{i,k} := P_{i+1,k−1} + \frac{1}{(N_{j}/N_{i})^{2}−1}(P_{i+1,k−1} − P_{i,k−1})  -lation

(1.4.3)

Im Verfahren (1.4.3) kann einfach die am langsamsten wachsende Folge {N_{j}} = {2, 4, 6, 8, 10, ..} verwendet werden. Schritt 3 ist der erwähnte Glättungsschritt. In Schritt 4 wird die Richardson-Extrapolation ausgeführt, die dem Neville-Algorithmus (I.3.3.5) zur Polynominterpolation entspricht und der Reihe nach jeweils den führenden Term der asymptotischen Entwicklung (1.4.2) eliminiert (Numerik I, §5.5). Im j-ten Lauf durch Schritt 3-4 wird dabei mit der Näherung P_{j,0} vom j-ten Gitter eine zusätzliche Schrägzeile unten an das Extrapolations-Tableau angehängt.

\begin{matrix}P_{0,0} & P_{0,1} & P_{0,2} & ⋯ & P_{0,j} \\ P_{1,0} & P_{1,1} \\ P_{2,0} & ⋮ \\ ⋮ & P_{j−1,1} \\ P_{j,0} \end{matrix}

Danach gilt die (lokale) Fehleraussage

\frac{1}{H} [P_{0,j} − u(x_{0} + H)] = O(H^{2j+2}), j ≤ q.  (1.4.4)

Dieses Endergebnis P_{0,j} gehört zu einem Einschrittverfahren der Ordnung 2j + 2 für den Schritt x → x + H, da es nicht auf Werte vor der Stelle x zugreift. Damit wurde insbesondere folgende Existenzaussage zur Theorie der Einschrittverfahren geliefert, die sich durch Abzählung der insgesamt in (1.4.3) benötigten Funktionsauswertungen ergibt.

Satz 1.4.2 Für jedes p ∈ 2N gibt es ein Einschrittverfahren mit Ordnung p und \frac{1}{4}p^{2} +1 Stufen.

Das Extrapolationsverfahren benutzt also mehr Funktionsauswertungen als die in §1.2 besprochenen Runge-Kutta-Verfahren. Aufgrund der einfachen und einheitlichen Konstruktion von Verfahren verschiedener Ordnung kann hier aber eine Schrittweiten- und Ordnungssteuerung wie bei den Mehrschrittverfahren durchgeführt werden. Insbesondere wird die aktuell verwendete Ordnung 2j + 2(≤ 2q + 2) in (1.4.4) in der Praxis nicht vorgegeben, sondern durch Inspektion des Tableaus P_{i,k} bestimmt.

Demo-Beispiel: Extrapolationsverfahren mit Ordnungs- und Schrittweitensteuerung (vgl. Ende §1.3.1) beim AWP u′ = f(x, u), u(−\frac{1}{2}) = u_{0} auf [−\frac{1}{2}, 1] mit (a := 800, b = ln(1 + a) − 3)

f(x, y) = \{ \begin{matrix}−2axy^{2}, & x ≤ 0 \\ 2[−\frac{ax}{1 + ax^{2}} + bx]y^{2}, & x > 0 \end{matrix} mit u(x) = \{ \begin{matrix}\frac{1}{1 + ax^{2}} , & x ≤ 0 \\ \frac{1}{1 + ln(1 + ax^{2}) − bx^{2}} , & x > 0 \end{matrix} .

Vergleich der Verfahren:

Einschritt  Extrapolation  Mehrschritt
Konsistenz-Ordnung aufwendig  einfach  ⇐  einfach
Stabilität  immer ⇒  immer  Einschränkung
Aufwand  merkbar  hoch  gering
Fehlerschätzung  möglich  einfach  einfach
Schrittweitenänd.  immer ⇒  immer  aufwendig
Ordnungswechsel  nein  einfach  möglich

Die praktische Erfahrung mit den drei behandelten Verfahrensklassen ergibt ein differenziertes Bild, keines kann als Universalverfahren bezeichnet werden. Zwar schneiden in Bezug auf die Gesamtanzahl der zur Lösung benötigten Funktionsaufrufe von f die Prädiktor-Korrektor-Verfahren am günstigsten ab. Da sie im adaptiven Einsatz beim Schrittweitenwechsel aber den höchsten Verwaltungsaufwand haben, sind sie tatsächlich nur selten schneller (in Bezug auf die Rechenzeiten) als Einschrittverfahren, nämlich in Fällen, wo die numerische Auswertung f(x, y) sehr aufwendig ist. Im Vergleich der Einschrittverfahren, Runge-Kutta und Extrapolation, zeigt sich ein Vorteil der möglichen Ordnungssteuerung bei letzteren erst für (extrem) hohe Genauigkeitsanforderungen. Das Runge-Kutta-Tripel DOP853 wird erst für extreme Toleranzen (< 10^{−12}) von Extrapolationsverfahren übertroffen.

1.5 Steife Anfangswertprobleme

In der Praxis trifft man oft auf Differentialgleichungen, deren Lösungen sehr unterschiedliche Zeitskalen besitzen, d.h., es können sowohl sehr schnell, als auch sehr langsam veränderliche Lösungen auftreten. Dies ist, z.B., bei chemischen Reaktionen der Fall. Von größerer Bedeutung im Rahmen dieser Vorlesung ist aber die Tatsache, daß bestimmte Approximationsverfahren für parabolische, partielle Differentialgleichungen (Linienmethode, vgl. §3.4) auf gewöhnliche Anfangswertprobleme dieser Art mit großen Dimensionen n führen. Die Behandlung solcher Probleme ist ein sehr aktuelles Forschungsgebiet, daher können nur die wichtigsten Themen angerissen werden. Eine ausführlichere Darstellung findet sich in den Büchern von Hairer/Wanner und Strehmel/Weiner. Zur Problematik:

Beispiel 1.5.1 u(x) ∈ R^{2},

u′(x) = Lu(x) := ( \begin{matrix}−500 & 499 \\ 499 & −500 \end{matrix} ) u(x), u(0) = u_{0} := ( \begin{matrix}1 \\ 0 \end{matrix} ) .  (1.5.1)

Da die Eigenwerte der Koeffizientenmatrix L die Werte λ_{1} = −1 und λ_{2} = −999 sind, lautet die exakte Lösung dieses Problems

u(x) = ( \begin{matrix}C_{1}e^{−x} + C_{2}e^{λ_{2}x} \\ C_{1}e^{−x} − C_{2}e^{λ_{2}x} \end{matrix} )  (1.5.2)

mit C_{1} = C_{2} = \frac{1}{2}. Bis auf eine sehr kurze Startphase 0 ≤ x ≪ 1/1000 verhält sich die Lösung wie die Funktion \frac{1}{2} \binom{1}{1}e^{−x}, trivialerweise gilt auch u(x) → 0, x → ∞. Wird das Problem, z.B., mit dem (expliziten) Euler-Verfahren und konstantem h behandelt, ergibt sich die numerische Lösung aus

y_{i+1} = (I + hL)y_{i} ⇒ y_{i} = ( \begin{matrix}C_{1}(1 − h)^{i} + C_{2}(1 + hλ_{2})^{i} \\ C_{1}(1 − h)^{i} − C_{2}(1 + hλ_{2})^{i} \end{matrix} ) .

Hier gilt offensichtlich y_{i} → 0 (i → ∞) nur für |1 + hλ_{2}| < 1, d.h., wenn

h < \frac{2}{−λ_{2}} = \frac{1}{500}.

Im anderen Fall wächst y_{i} oszillierend stark an. Obwohl also (z.B.) e^{λ_{2}x} ≤ 10^{−8} für x ≥ 0.02 ist und die Lösung für solche Argumente vollkommen glatt aussieht, auch durch einen Polygonzug mit größeren Schrittweiten gut approximiert werden kann, erzwingt das Vorhandensein des Eigenwerts λ_{2} im Problem die Verwendung kleiner Schrittweiten bei diesem Verfahren.

Beim impliziten Euler-Verfahren wird das Integral ∫_{0}^{1} f(x_{i} +th, u(x_{i} +th)) dt in (1.2.1) durch den Funktionswert am rechten Rand ersetzt. Bei diesem Verfahren,

y_{i+1} = y_{i} + h_{i}f(x_{i+1}, y_{i+1}) ⟺ \{ \begin{matrix}k_{1} & = & f(x_{i} + h_{i}, y_{i} + h_{i}k_{1}) \\ y_{i+1} & = & y_{i} + h_{i}k_{1} \end{matrix}  (1.5.3)

sind k_{1} bzw. y_{i+1} nur implizit beschrieben, müssen also durch Auflösung eines (nicht-)linearen Gleichungssystems bestimmt werden. Im linearen Beispiel (1.5.1) führt dies auf die Beziehungen

y_{i+1} = y_{i} + h_{i}Ly_{i+1} ⟺ (I − h_{i}L)y_{i+1} = y_{i} ⟺ y_{i+1} = (I − h_{i}L)^{−1}y_{i}.

Mit Hilfe der Eigenwertzerlegung von L kann auch hier wieder das Ergebnis für konstante Schrittweite h explizit bestimmt werden,

y_{i} = ( \begin{matrix}C_{1}(1 + h)^{−i} + C_{2}(1 − hλ_{2})^{−i} \\ C_{1}(1 + h)^{−i} − C_{2}(1 − hλ_{2})^{−i} \end{matrix} ) .

Da jetzt 1/(1+h) < 1 und 1/(1−hλ_{2}) < 1 gilt, für alle h > 0, ist zumindestens das asymptotische Verhalten der Lösung u reproduziert, y_{i} → 0, i → ∞ für alle Schrittweiten.

Die beiden Beispielverfahren unterscheiden sich also ganz wesentlich in ihren Stabilitätseigenschaften. Zur Überprüfung solcher Eigenschaften von numerischen Lösungen verschiedener Verfahren betrachtet man deren Verhalten bei bestimmten Testproblemen im Vergleich zur exakten Lösung. Die einfachste Testgleichung ist die skalare, linear-homogene

u′(x) = λu(x), u(0) = 1, λ ∈ C, Re λ ≤ 0,  (1.5.4)

die aber über die Eigenvektor-Zerlegung auch Aussagen über Systeme der Form (1.5.1) gestattet. Für deren Lösungen gilt |u(x + s)| = |e^{λs}u(x)| ≤ |u(x)| ∀s ≥ 0. Runge-Kutta-Verfahren (1.2.4) führen hier auf die skalaren Gleichungen

y_{i+1} = y_{i} + h ∑_{j=1}^{m} b_{j}k_{j},
k_{j} = λy_{i} + hλ ∑_{l=1}^{m} a_{jl}k_{l} ⟺ (I − zA)k = λy_{i}1

wobei z := hλ, k := (k_{1}, .., k_{m})^{T} , A = (a_{jl}), 1 := (1, .., 1)^{T} gesetzt wurde. Analog sei b = (b_{1}, .., b_{m})^{T}. Aus dieser Beziehung ergibt sich die Rekursion y_{i} → y_{i+1} explizit in der Form

y_{i+1} = y_{i} + hb^{T} k = [1 + zb^{T} (I − zA)^{−1}1]y_{i} = φ(z)y_{i},  (1.5.5)
φ(z) := 1 + zb^{T} (I − zA)^{−1}1.  (1.5.6)

Einschrittverfahren führen bei (1.5.4) immer auf eine solche Beziehung y_{i+1} = φ(hλ)y_{i}. In (1.5.5) ist φ(z) eine rationale Funktion der komplexen Variabeln z = hλ, die sogenannte Stabilitätsfunktion. Für explizite Verfahren wurden diese schon in Beispiel 1.2.1 bzw. bei (1.5.3) berechnet:

Verfahren: expl. Euler Runge u.Heun (Nr.2/3) klass. Runge-Kutta (Nr.4) impl. Euler ; φ(z) =  1 + z  1 + z + \frac{1}{2}z^{2}  1 + z + \frac{1}{2}z^{2} + \frac{1}{6}z^{3} + \frac{1}{24}z^{4}  \frac{1}{1−z}

Offensichtlich sind die Stabilitätsfunktionen der expliziten Verfahren ganze rationale Funktionen, d.h., Polynome. Für diese ist ja A eine strikt untere Dreiecksmatrix, also nilpotent mit A^{m} = 0. Daher ist die Neumannreihe ein Polynom, φ(z) = 1 + zb^{T} (I − zA)^{−1}1 = 1 + z ∑_{j=0}^{m−1} z^{j}b^{T} A^{j}1. Wenn ein Verfahren mit Ordnung p konvergiert, bildet die Stabilitätsfunktionen natürlich eine rationale Approximationen an die e-Funktion mit φ(z) = e^{z} + O(z^{p+1}), z → 0.

Bei der Testgleichung (1.5.4), die nur beschränkte Lösungen besitzt, interessieren diejenigen Argumentwerte z = hλ, für die die numerische Lösung y ebenfalls nicht wächst.

Definition 1.5.2 Ein Einschrittverfahren besitze die Stabilitätsfunktion φ, d.h. bei Anwendung auf (1.5.4) führe es auf die Beziehung y_{i+1} = φ(hλ)y_{i}. Das Verfahren heißt absolut stabil bei z ∈ C, falls |φ(z)| ≤ 1 ist. Die Menge

S := {z ∈ C : |φ(z)| ≤ 1}  (1.5.7)

heißt Bereich der absoluten Stabilität des Verfahrens.

Zur stabilen Integration von (1.5.4) muß bei einem gegebenen Verfahren die Schrittweite h so gewählt werden, daß hλ ∈ S gilt, da sonst unsinnige, wachsende Lösungen auftreten können. Für die in R liegenden Stabilitätsintervalle S ∩ R der obigen Beispielverfahren gilt

Verfahren: expl. Euler Runge u.Heun (Nr.2/3) klass.Runge-Kutta (Nr.4) impl. Euler ; Intervall:  [−2, 0]  [−2, 0]  [−2.78, 0]  (−∞, 0]

Die in Beispiel 1.5 für das explizite Euler-Verfahren gefundene Schrittweiten-Einschränkung h ≤ 2/(−λ) entspricht also der Bedingung hλ ∈ [−2, 0] = S ∩ R. Von allen erwähnten Verfahren kommt bei Anwendung auf Probleme der Form (1.5.1) nur das implizite Eulerverfahren ohne Schrittweitenbeschränkung aus, da dessen Stabilitätsbereich die negative reelle Achse umfaßt. Da bei expliziten Verfahren die Stabilitätsfunktion φ ein Polynom ist, mit |φ(z)| → ∞, z → ∞, besitzen diese nur beschränkte Stabilitätsbereiche.

Daher läßt sich allgemein feststellen, daß zur Lösung steifer Probleme, mit möglicherweise schnell ausklingenden Lösungen, nur implizite Verfahren sinnvoll sind, da nur diese sich ausschließlich nach Genauigkeitsforderungen richten müssen. Bei expliziten Verfahren dagegen wird die Schrittweite nicht durch die gewünschte Genauigkeit bestimmt, sondern (sehr restriktiv) durch Stabilitätsschwächen. Eine automatische Schrittweitensteuerung erzeugt hier eine um den Grenzfall schwankende Schrittweitenfolge (vgl. Demo-Beispiel). In Anlehnung an das Lösungsverhalten der Dgl (1.5.4) zeichnet man Verfahren mit analogen Dämpfungseigenschaften aus.

Definition 1.5.3 Ein Einschrittverfahren heißt A-stabil, wenn gilt

S ⊇ C_{−} := {z ∈ C : Re z ≤ 0}, d.h., |φ(z)| ≤ 1 ∀z ∈ C_{−}.

Offensichtlich ist das implizite Eulerverfahren A-stabil, da hier die Bedingung |φ(z)| = |\frac{1}{1−z} | ≤ 1 ⟺ 1 ≤ |1 − z|^{2} = 1 − 2Re z + |z|^{2} für Re z ≤ 0 sicher erfüllt ist.

Implizite Runge-Kutta-Verfahren: Bei den IRK-Verfahren gibt es eine Klasse, bei der gleichzeitig die höchstmöglichen Konvergenzordnungen und gute Stabilitätseigenschaften aufeinandertreffen. Diese Verfahren ergeben sich aus Anwendung der Gaußschen Quadraturformeln auf (1.2.1). Dazu werden die Stützstellen c_{i} in (1.2.4) als Gaußknoten, d.h. die Nullstellen des m-ten Legendre-Polynoms gewählt (vgl. Numerik I, §5.3). Die restlichen Koeffizienten bestimmen sich als Integrale der zugehörigen Lagrangepolynome

L_{i}(x) := ∏_{j≠i^{j=1}}^{m} \frac{x − c_{j}}{c_{i} − c_{j}} , a_{ij} := ∫_{0}^{c_{i}} L_{j}(x) dx, b_{i} := ∫_{0}^{1} L_{i}(x) dx, i, j = 1, .., m.  (1.5.8)

Satz 1.5.4 Für beliebiges m ∈ N sind die Gauß-Runge-Kutta-Verfahren A-stabil. Die Konsistenzordnung des m-stufigen Verfahrens ist 2m.

Die Gauß-Runge-Kutta-Verfahren besitzen also exzellente theoretische Eigenschaften. Es können sogar Stabilitätsaussagen für nichtlineare Probleme gezeigt werden. Allerdings sind bei diesen Verfahren in jedem Verfahrensschritt nichtlineare Gleichungssysteme der Größe mn, z.B. mit dem Newtonverfahren, zu lösen. Daher sind sie für den praktischen Einsatz meist zu teuer. Es gibt verschiedene Ansätze zur effektiveren Implementierung von IRK-Verfahren. Diese schränken aber entweder die Struktur ein, und damit leider auch die erreichbare Ordnung auf m + 1, oder besitzen andere Nachteile.

Eine effiziente Alternative sind Linear-implizite Runge-Kutta-Verfahren. Zur Motivation werde die Lösung der Gleichung (1.5.3) beim impliziten Euler-Verfahren betrachtet mit sog. autonomer rechter Seite f = f(u). Der erste Schritt der Newton-Iteration mit Startwert y_{i+1}^{[0]} := y_{i} und L_{i} := \frac{∂f}{∂y} (y_{i}) ergibt

y_{i+1}^{[1]} = y_{i} + h(I − hL_{i})^{−1}f(y_{i}).

Man kann nun diese Vorschriftt als eigenständiges Verfahren studieren in Bezug auf Stabilität und Konsistenz. Bei mehrstufigen Verfahren sind dabei zusätzliche Parameter γ, γ_{ij} sinnvoll. Die Verfahrensvorschrift lautet dann mit L := \frac{∂f}{∂y} (y_{i}),

(I − hγL)k_{l} = f(y_{i} + h ∑_{j=1}^{l−1} α_{lj}k_{j}) + hL ∑_{j=1}^{l−1} γ_{lj}k_{j}, l = 1, …, m  (1.5.9)
y_{i+1} = y_{i} + h ∑_{j=1}^{m} b_{j}k_{j}.

Die Stabilitätsfunktion dieser Verfahren entspricht der eines IRK-Verfahrens mit Koeffizienten a_{lj} = α_{lj} +γ_{lj}, daher gibt es A-stabile Verfahren. Da in allen Stufen die gleiche Matrix (I −hγL) auftritt, ist pro Schritt y_{i} → y_{i+1} nur eine LR- oder QR-Zerlegung erforderlich. Die Verfahren gehen auf Rosenbrock und Wanner zurück. Man kann auch hier eingebettete Verfahren studieren zur Schrittweitensteuerung, z.B., existieren verschiedene (3,4)-Paare. Diese sind für mittlere Genauigkeitsforderungen konkurrenzfähig zur folgenden Verfahrensklasse.

Implizite Mehrschrittverfahren Im praktischen Einsatz hat sich eine spezielle Verfahrensklasse bewährt, die auf den rückwärtigen Differenzenformeln beruht. Im Gegensatz zu den Adams-Verfahren wird in der Form des allgemeinen Mehrschrittverfahrens (1.3.6) nur ein nichtverschwindender Koeffizient β_{0} ≠ 0 gewählt, während die Koeffizienten α_{j} aus dem linksseitigen Differenzenquotienten abgeleitet werden mit

\frac{1}{h} ∑_{j=0}^{m} α_{j}v(−jh) = β_{0}v′(0) + O(h^{m}), h → 0,  (1.5.10)

für v ∈ C^{m}[−m, 0]. Aufgrund von Satz 1.3.2 besitzt daher die m-stufige Version dieser BDF-Verfahren (backward difference formula) Ordnung m. Die Koeffizienten sind

m  j = 1 2 3 4 β_{0} δ_{m}
1  α_{j} = −1  1 90^{o}
2 3α_{j} = −4 1  \frac{2}{3} 90^{o}
3 11α_{j} = −18 9 −2  \frac{6}{11} 86^{o}
4 25α_{j} = −48 36 −16 3 \frac{12}{25} 73^{o}

Ihre Stabilitätseigenschaften sind allerdings weniger gut als bei Einschrittverfahren. Bei der Testgleichung (1.5.4) haben die BDF-Verfahren die Gestalt (z := hλ)

(1 − β_{0}z)y_{i} + ∑_{j=1}^{m} α_{j}y_{i−j} = 0, i ≥ m.  (1.5.11)

Über den Lösungsansatz y_{i} = Cξ^{i} kommt man wieder auf die allgemeine Lösung y_{i} = ∑_{j=1}^{m} C_{j}ξ_{j}^{i}, wenn das charakteristische Polynom

p_{m}(ξ; z) = (1 − β_{0}z)ξ^{m} + α_{1}ξ^{m−1} + … + α_{m}

der Differenzengleichung (1.5.11) m verschiedene Nullstellen ξ_{j} besitzt. Bei mehrfachen Nullstellen ist die Darstellung durch Polynomanteile i^{l}ξ^{i} zu modifizieren. In Analogie zu den Einschrittverfahren kann man die Menge dieser Nullstellen {ξ_{j}} = {ξ_{j}(z)}_{j=1}^{m} als (mehrwertige) Stabilitätsfunktion des Mehrschrittverfahrens interpretieren. Insbesondere nennt man ein Mehrschrittverfahren absolut stabil an der Stelle z ∈ C, wenn

p_{m}(ξ; z) = 0 ⇒ |ξ| < 1

gilt und A-stabil, wenn es für jedes z ∈ C_{−} absolut stabil ist. Allerdings gibt es eine zweite Ordnungsbarriere von Dahlquist:

Satz 1.5.5 Die Ordnung eines A-stabilen Mehrschrittverfahrens ist höchstens 2.

Daher können die Stabilitätsbereiche S_{m} der BDF-Verfahren höherer Ordnung m ≥ 2 nicht mehr die ganze Halbebene C_{−} umfassen, in der obigen Tabelle ist der (halbe) Öffnungswinkel δ_{m} des größten in S_{m} gelegenen Winkelraums angegeben. Für m > 6 sind die Verfahren unbrauchbar, da dann nicht einmal mehr die negative reelle Achse vollständig in S enthalten ist. Auch für m = 6 ist der Winkel δ_{6} = 17.8^{o} schon sehr klein (vgl. Diagramme).

[figure: Re z; Im z; S_{m}; δ_{m}]

Demo-Beispiel Das nichtlineare System von Dgln

\begin{matrix}u′_{j} = & −jdu_{j}, & j ∈ \{1, n − 1, n\} \\ u′_{j} = & −jdu_{j} + (u_{j−1}u_{j+2} − u_{j+1}u_{j})^{2}, & j = 2, …, n − 2 \end{matrix}

besitzt die Lösung u_{j}(t) = e^{−jdt}. Für große n ist die Lipschitzkonstante der rechten Seite groß (≥ dn). Die Stabilitätsprobleme expliziter Verfahren mit Schrittweitensteuerung sind daran erkennbar, daß Schrittweiten und Rechenzeit i.w. unabhängig von der geforderten Genauigkeit sind.

1.6 Schießverfahren für Randwertprobleme

Hier wird das standardisierte, nichtlineare Randwertproblem (RWP, vgl. §1.1)

u′(x) = f(x, u(x)), x ∈ (a, b), r(u(a), u(b)) = 0,  (1.6.1)

betrachtet mit differenzierbaren Funktionen f(x, y), r(y_{0}, y_{1}). Bei der numerischen Behandlung solcher gewöhnlicher RWPe kann man sowohl Lösungsverfahren von den Anfangswertproblemen übertragen, als auch neue Verfahren untersuchen, die auf partielle Randwertprobleme verallgemeinert werden können. Aus der Reihe der letzteren wird eine spezielle Klasse schon hier besprochen, da deren wesentliche Eigenschaften im gewöhnlichen Fall einfacher zu diskutieren sind. Zunächst wird mit Methoden begonnen, die auf AWP-Verfahren aufbauen. Diese haben den praktischen Vorteil, daß Standardprogramme aus diesem Bereich eingesetzt werden können.

Schießverfahren Aufbauend auf dem letzten Paragraphen wird vorausgesetzt, daß Anfangswertprobleme (mit genügender Genauigkeit) numerisch gelöst werden können. In §1.1 war beim linearen Randwertproblem

u′(x) = A(x)u(x) + g(x), x ∈ (a, b), R_{0}u(a) + R_{1}u(b) = d,  (1.6.2)

mit einer Fundamentallösung Y (x) des Matrix-Anfangswertptoblems Y′(x) = A(x)Y (x), Y (a) = I, aus der Lösungsdarstellung u(x) = Y (x)η + .. der unbekannte Anfangswert η = u(a) ∈ R^{n} durch Einsetzen in die Randbedingung bestimmt worden: R_{0}η + R_{1}[Y (b)η + ..] = d. Damit wurde also die Lösungskurve u(x) durch Variation des Startwertes η = u(a) so angepaßt, daß ihr Wert u(b) am rechten Rand (zusammen mit u(a) = η) die Randbedingung erfüllt (Bildlich Schießen/Werfen: Anpassung der Wurfrichtung so, daß Treffer erzielt wird). Zur Verallgemeinerung dieses Prinzips muß die Abhängigkeit der Lösung vom Anfangswert betont werden. Daher bezeichne jetzt u(x; a, η) die Lösung des AWPs

\frac{d}{dx}u(x; a, η) = f(x, u(x; a, η)), x > a, u(a; a, η) = η.  (1.6.3)

Die spezielle Lösung des RWPs (1.6.1) ist dann diejenige, u(x) = u(x; a, η̂), deren Anfangswert η̂ das nichtlineare Gleichungssystem

F (η) = 0, mit F : \{ \begin{matrix}R^{n} & → & R^{n} \\ η & ↦ & r(η, u(b; a, η)) \end{matrix}  (1.6.4)

erfüllt. Zur Lösung solcher Gleichungssysteme kommt vor allem das Newton-Verfahren in Betracht, vgl. Numerik I, §6.2. Sein Einsatz setzt aber die genügende Differenzierbarkeit der Funktion F voraus. Die Lipschitz-stetige Abhängigkeit der Funktion u(x; a, η) vom Anfangswert η wurde schon in Satz 1.1.1 gezeigt. Darüber hinaus gilt aber auch Satz 1.6.1 Die rechte Seite f sei zweimal stetig differenzierbar, ihre Ableitung f_{y} = \frac{∂f}{∂y} beschränkt auf [a, b] × R^{n}. Dann hängt die Lösung u(x; a, η) von (1.6.3) differenzierbar vom Anfangswert η ab. Ihre Ableitung an der Stelle (x, u(x; a, η)) ist

\frac{∂}{∂η} u(x; a, η) = Y (x; a, η),  (1.6.5)

wobei Y das Fundamentalsystem zur linearisierten Differentialgleichung ist,

\frac{∂}{∂x}Y (x; a, η) = f_{y}(x, u(x; a, η)) · Y (x; a, η), Y (a; a, η) = I.  (1.6.6)

Beweis Die Ableitung wird direkt als in s linearer Anteil der Differenz v(x) := u(x; a, η + s) − u(x; a, η), s ∈ R^{n}, bestimmt. Dies entspricht der Definition der Differenzierbarkeit. Offensichtlich ist v(a) = η + s − η = s. Nach Satz 1.1.1 gilt daher die Schranke ∥v(x)∥ ≤ exp(L(x − a))∥s∥. Da beide Funktionen u Lösungen der Dgl sind, gilt auch für x ∈ (a, b) (mit ′ = ∂/∂x)

v′(x) = u′(x; a, η + s) − u′(x; a, η) = f(x, u(x; a, η) + v(x)) − f(x, u(x; a, η))
= f_{y}(x, u(x; a, η))v(x) + O(∥s∥^{2}).

Die Funktion v̂(x) := Y (x; a, η)s mit Y aus (1.6.6) ist Lösung des AWPs v̂′(x) = f_{y}(x, u(x; a, η))v̂(x), v̂(a) = s. Daher gilt insbesondere v(x) − v̂(x) = O(∥s∥^{2}) und somit

u(x; a, η + s) = u(x; a, η) + Y (x; a, η)s + O(∥s∥^{2}), s → 0.

Damit erfüllt der in s lineare Anteil Y s die Definition der Ableitung von u in Richtung s, also gilt \frac{∂}{∂η} u = Y, wie in (1.6.5) behauptet.

Die im Newtonverfahren zur Lösung des Systems (1.6.4) benötigte Ableitung ergibt sich somit nach der Kettenregel

F′(η) = \frac{d}{dη} r(η, u(b; a, η)) = \frac{∂r}{∂y_{0}} (η, u(b; a, η)) + \frac{∂r}{∂y_{1}} (η, u(b; a, η))\frac{∂u}{∂η}(b; a, η)
= \frac{∂r}{∂y_{0}} (η, u(b; a, η)) + \frac{∂r}{∂y_{1}} (η, u(b; a, η))Y (b; a, η).

Die Bestimmung von F′ ist somit durch Lösung des Matrix-AWPs (1.6.6), also durch (numerische) Berechnung von n Lösungen y_{i}(x) der linearen AWPe y′_{i} = f_{y}(x, u(x; a, η))y_{i}, y_{i}(a) = e_{i}, i = 1, .., n, möglich (e_{i} : i-ter Einheitsvektor). Zur praktischen Durchführung können alle im §1.2-1.4 behandelten Verfahren herangezogen werden, es ist aber einfacher von einer exakten Lösung der AWPe auszugehen. Damit sind die beim Newtonverfahren

η^{[k+1]} := η^{[k]} − (F′(η^{[k]}))^{−1}F (η^{[k]}), k = 0, 1, ..

auftretenden Größen bekannt. Das Verfahren wird noch einmal ausführlich formuliert, wobei bei den Hilfsgrößen u(x; a, η^{[k]}), Y (x; a, η^{[k]}) der Index k und die Abhängigkeit vom Startwert unterdrückt wird.

Algorithmus 1.6.2 Schießverfahren: Gegeben sei η^{[0]} ∈ R^{n}. Für k = 0, 1, … löse man

die AWPe \{ \begin{matrix}u′(x) & = & f(x, u(x)), & u(a) = η^{[k]}, \\ Y′(x) & = & f_{y}(x, u(x))Y (x), & Y (a) = I, \end{matrix}  (1.6.7)
das LGS (r_{y_{0}} + r_{y_{1}}Y (b))(η^{[k+1]} − η^{[k]}) = −r(η^{[k]}, u(b)).  (1.6.8)

Die partiellen Ableitungen von r sind ebenfalls an der Stelle (η^{[k]}, u(b)) auszuwerten. Zunächst sind also in (1.6.7) jeweils n + 1 gekoppelte AWPe zu lösen und danach in (1.6.8) eine lineares n × n-System.

Im linearen Spezialfall (1.6.2) ist r(y_{0}, y_{1}) = R_{0}y_{0} +R_{1}y_{1} −d, d.h. \frac{∂r}{∂y_{0}} = R_{0}, \frac{∂r}{∂y_{1}} = R_{1}. Auch f_{y}(x, y) = A(x) ist unabhängig von u(x; a, η) also auch von η. Somit hängt auch Y nicht von η ab. Ein Schritt des Schießverfahrens mit Startwert η^{[0]} = 0 liefert daher

η^{[1]} = −(R_{0} + R_{1}Y (b))^{−1}(0 + R_{1}u(b) − d).

Dies ist genau der Lösungswert, der bei der expliziten Lösung in (1.1.10) berechnet wurde, (1.6.8) konvergiert also bei linearen Problemen erwartungsgemäß in einem Schritt.

Wie auch sonst beim Newton-Verfahren kann die Konvergenz dieser Iteration für ‘genügend genaue’ Startwerte η^{[0]} gezeigt werden, wenn die zweite Ableitung f_{yy} beschränkt ist, und in der exakten Lösung u des RWPs die Randmatrix F′(u(a)) = r_{y_{0}} + r_{y_{1}}Y (b) regulär ist, wobei Y (b) = Y (b; a, u(a)). In der Praxis ist diese Aussage aber oft wertlos, da die Lösungen u(x; a, η) sich möglicherweise sehr schnell von der gesuchten Lösung u(x; a, u(a)) entfernen. Dann ist der Einzugsbereich der Iteration (1.6.8), und eventuell sogar der Definitionsbereich der Funktion F überhaupt, sehr klein.

Beispiel 1.6.3 Diese Probleme können anhand der allgemeinen Lösung der skalaren Dgl

u′ = \frac{1 + u^{2}}{1 + x^{2}} =: f(x, u) ⇒ u(x; a, η) = \frac{η − a + (1 + aη)x}{1 + aη − (η − a)x}

erläutert werden. Soll das RWP mit der Randbedingung u(0) + u(b) = b, b > 0, (Lösung u(x) ≡ x) mit dem Schießverfahren gelöst werden, ist ein Startwert η < 1/b zu wählen, damit die Lösung u(x; 0, η) = (η + x)/(1 − ηx) überhaupt den rechten Intervallrand erreicht. Mit der explizit bekannten allgemeinen Lösung erhält die Randgleichung (1.6.4) die Form

F (η) = η + u(b; 0, η) − b = η − b + \frac{b + η}{1 − bη} =^{!} 0.

Die Ableitung F′(η) kann hieraus natürlich direkt bestimmt werden. Andererseits gilt auch

Y (x) = \frac{∂u}{∂η}(x; 0, η) = \frac{1 + x^{2}}{(1 − ηx)^{2}} und Y′(x) = 2\frac{η + x}{(1 − ηx)(1 + x^{2})}Y (x) = \underbrace{\frac{2u(x; 0, η)}{1 + x^{2}}}_{f_{y}(x,u)} Y (x).

Bei positiven Startwerten konvergiert das Newtonverfahren ebenfalls für η < 1/b (Ü-Aufgabe).

Die oben erwähnten Schwierigkeiten bei rasch divergierenden Lösungen des Anfangswertproblems können entschärft werden dadurch, daß man die Lösungen jeweils nur über kürzere Intervalle integriert, das Gesamtintervall also in mehrere Teile zerlegt. Dies führt auf die im folgenden beschriebene Mehrzielmethode: Das Gesamtintervall [a, b] wird jetzt unterteilt durch ein Gitter

Δ : a = x_{0} < x_{1} < … < x_{N} = b.

[figure: a = x_{0}; x_{1}; x_{2}; x_{3}; x_{N−1}; x_{N} = b; η_{0}; η_{1}; η_{N−1}; u(x; x_{2}, η_{2})]

In jedem Teilintervall [x_{i}, x_{i+1}] wird das Schießverfahren getrennt angewendet. Dabei werden N Teil-Lösungen u(x; x_{i}, η_{i}) zu

u′(x; x_{i}, η_{i}) = f(x, u(x; x_{i}, η_{i})), u(x_{i}; x_{i}, η_{i}) = η_{i}, i = 0, …, N − 1,  (1.6.9)

berechnet. Eine Lösung u des Randwertproblems (1.6.1) kann nun aus diesen Stücken genau dann zusammengesetzt werden, wenn diese an den Gitterpunkten ‘zusammenpassen’ und außerdem die Randbedingung erfüllen, wenn also gilt

u(x_{i}; x_{i−1}, η_{i−1}) = η_{i}  ≡ u(x_{i}; x_{i}, η_{i}), i = 1, …, N − 1,
r(η_{0}, u(x_{N}; x_{N−1}, η_{N−1})) = 0.

Wegen der Stetigkeit von f ist die so zusammengesetzte Lösung auch stetig differenzierbar. Im Unterschied zum Einfach-Schießverfahren sind nun N unbekannte Startwerte η_{i} zu bestimmen als eine Lösung des nichtlinearen Gleichungssystems für

η⃗ := ( η_{0} ; η_{1} ; ⋮ ; η_{N−1} ) ∈ R^{nN} mit F (η⃗) := ( \begin{matrix}u(x_{i}; x_{i−1}, η_{i−1}) − η_{i}, i = 1, .., N − 1, \\ r(η_{0}, u(x_{N}; x_{N−1}, η_{N−1})),  i = N \end{matrix} ) = 0 (1.6.10)

Hierin ist das einfache Schießverfahren mit N = 1 enthalten. Beim linearen RWP (1.6.2) können die Stück-Lösungen (1.6.9) wieder explizit dargestellt werden (vgl. Satz 1.1.3) durch

u(x; x_{i}, η_{i}) = Y (x)Y (x_{i})^{−1}η_{i} + ∫_{x_{i}}^{x} Y (x)Y (t)^{−1}g(t)dt.

Dabei sind Y (x)Y (x_{i})^{−1} =: Y (x; x_{i}) diejenigen (hier von η_{i} unabhängigen) Fundamentallösungen der linearen Dgl mit Y (x_{i}; x_{i}) = I. Die Integrale entsprechen den speziellen Lösungen u(x; x_{i}, 0). Das System (1.6.10) ist hier dann auch linear und lautet

Y (x_{i}; x_{i−1})η_{i−1} − η_{i} = −u(x_{i}; x_{i−1}, 0), i = 1, .., N − 1,
R_{0}η_{0} + R_{1}Y (b; x_{N−1})η_{N−1} = d − Ru(b; x_{N−1}, 0).

Die (nN) × (nN)-Matrix dieses Systems hat eine spezielle Struktur. Sie besitzt eine zyklische Block-Bidiagonalgestalt (Treppenform) aus n × n-Blöcken,

( \begin{matrix}Y (x_{1}; x_{0}) & −I \\ & Y (x_{2}; x_{1}) & −I \\ & & ... & ... \\ & & & Y (x_{N−1}; x_{N−2}) & −I \\ R_{0} & & & & R_{1}Y (b; x_{N−1}) \end{matrix} )  (1.6.11)

Bei der Lösung des allgemeinen, nichtlinearen Problems (1.6.10) mit dem Newtonverfahren treten lineare Gleichungssysteme der gleichen Struktur auf. Denn die partiellen Ableitungen der Teillösungen nach ihrem Anfangswert sind nach Satz 1.6.1 wieder,

\frac{∂}{∂η_{i}} u(x; x_{i}, η_{i}) = Y (x; x_{i}, η_{i})
Y′(x; x_{i}, η_{i}) = f_{y}(x, u(x; x_{i}, η_{i}))Y (x; x_{i}, η_{i}), Y (x_{i}; x_{i}, η_{i}) = I.  (1.6.12)

Daher besitzt die Ableitungsmatrix F′ = ∂F/∂η⃗ von (1.6.10) genau die Gestalt (1.6.11). Somit ist beim Newtonverfahren zur Lösung des Systems (1.6.10),

\frac{∂F}{∂η⃗} (η^{[k]}) [η^{[k+1]} − η^{[k]}] = −F (η^{[k]}), k = 0, 1, ..

im k-ten Schritt für die Korrekturen

δ_{i} := η_{i}^{[k+1]} − η_{i}^{[k]}, i = 0, .., N − 1,

tatsächlich ein lineares Gleichungssystem der folgenden Form zu lösen,

\begin{aligned}Y (x_{i}; x_{i−1}, η_{i−1})δ_{i−1} − δ_{i} & = & η_{i} − u(x_{i}; x_{i−1}, η_{i−1}), i = 1, .., N − 1, \\ R_{0}δ_{0} + R_{1}Y (b; x_{N−1}, η_{N−1})δ_{N−1} & = & −r(η_{0}, u(b; x_{N−1}, η_{N−1})), \end{aligned}  (1.6.13)

mit R_{0} = \frac{∂r}{∂y_{0}} , R_{1} = \frac{∂r}{∂y_{1}} , jeweils an der Stelle (η_{0}, u(b; x_{N−1}, η_{N−1})). Die Aufstellung dieses Systems erfordert wieder die (numerische) Lösung der Teil-AWPe (1.6.9), (1.6.12).

Die Auflösung der linearen Systeme in (1.6.13) kann mit dem Gauß-Algorithmus durchgeführt werden. Bei Pivotisierung mit Spaltenvertauschung werden dabei Elemente nur in der untersten (Block-) Zeile neu eingeführt und es wird nur mit dem skizzierten Bereich gearbeitet:

(  \begin{matrix}. & . & . \\ ⋯ \end{matrix}  )

Dieses Verfahren ist auch in kritischen Fällen (z.B., wenn die Matrizen Y große Normen besitzen) numerisch recht stabil. In vielen Fällen kann dieses Gleichungssystem aber noch einfacher durch Kondensation gelöst werden. Dazu behandelt man nur δ_{0} als Unbekannte und eliminiert die restlichen δ_{i} mit Hilfe der ersten Gleichungen in (1.6.13):

δ_{1} = Y_{1}δ_{0} + r_{1}, δ_{2} = Y_{2}δ_{1} + r_{2} = Y_{2}Y_{1}δ_{0} + Y_{2}r_{1} + r_{2}, …
δ_{N−1} = Y_{N−1}Y_{N−2}..Y_{1}δ_{0} + r̂,

Y_{i} := Y (x_{i}; x_{i−1}, η_{i−1}). Mit der letzten Gleichung aus (1.6.13) ergibt dies also das System

(R_{0} + R_{1}Y_{N}Y_{N−1}..Y_{2}Y_{1})δ_{0} = −r(η_{0}, ..) − R_{1}Y_{N} r̂,  (1.6.14)

das nur noch von der Größe n × n ist. Diese Kondensation kann übrigens zusammen mit der Berechnung der Y_{i}, etc., erfolgen. Dann sind für einen Newtonschritt jeweils zwei Läufe durch das Intervall nötig, nämlich:

AWPe mit Kondensation, Lösung δ_{0} aus (1.6.14), Berechnung von δ_{1}, …, δ_{N−1}.

Im linearen Fall und in der Nähe der Lösung u unterscheidet sich die Mehrzielmethode mit Kondensation nicht vom einfachen Schießverfahren. Denn, für die Matrix in (1.6.14) gilt mit exakten Startwerten η_{i} ≡ u(x_{i}) die Aussage

R_{0} + R_{1}Y_{n} … Y_{1} = R_{0} + R_{1}Y (b; x_{N−1}, u(x_{N−1})) … Y (x_{1}; a, u(a)) = R_{0} + R_{1}Y (b; a, u(a)).

Letzteres ist aber genau die Matrix aus (1.6.8), da in der Lösung die (Halb-)Gruppenidentität Y (x; s, u(s))Y (s; t, u(t)) = Y (x; t, u(t)) ∀t < s < x gilt. Die lokale Konvergenz (bei u) des Newtonverfahrens hängt also insbesondere nur von der Regularität dieser Matrix ab. Der entscheidende Vorteil der Mehrzielmethode wirkt sich erst im nichtlinearen Fall aus durch (wesentliche) Vergrößerung des Einzugs-, d.h., Konvergenzbereichs beim Newtonverfahren. Die Wahl der Zwischenzielpunkte x_{1}, .., x_{N−1} kann dabei leicht während der Durchführung des Verfahrens modifiziert werden. Z. B. ist es sinnvoll, bei starkem Anwachsen einer Zwischenlösung u(.; x_{i}, η_{i}) einen neuen Zielpunkt einzufügen. Die Auswirkungen der Mehrzielmethode kann anhand des einführenden Beispiels diskutiert werden.

Beispiel 1.6.4 Eine Diskussion der Konvergenz ist hier naturgemäß viel schwieriger als im Fall N = 1. Als Mindestvoraussetzung kann aber die Existenz der Teillösungen in [x_{i}, x_{i+1}] überprüft werden. Beim Problem aus Beispiel 1.6 lautet der Lösungswert in x_{i+1},

u(x_{i+1}; x_{i}, η_{i}) ≡ \frac{η_{i}(1 + x_{i}x_{i+1}) + (x_{i+1} − x_{i})}{1 + x_{i}^{2} − (x_{i+1} − x_{i})(η_{i} − x_{i})}

und existiert insbesondere nur dann, wenn der Nenner positiv ist, also für

η_{i} − x_{i} = η_{i} − u(x_{i}) ≤ \frac{1 + x_{i}^{2}}{x_{i+1} − x_{i}} , i = 0, .., N − 1.

Diese Genauigkeitsforderungen an die lokalen Startwerte η_{i} sind schwächer als beim Einfach-Schießverfahren. Bei äquidistanter Unterteilung, x_{i} = ib/N, etwa ist für η_{0} nur zu fordern η_{0} − u(x_{0}) ≤ N/b, während in Beispiel 1.6 η − u(x_{0}) < 1/b erforderlich war.

Demo-Beispiel: Mehrzielmethode beim Randwertproblem 2. Ordnung u′′ = \frac{3}{2}u^{2}, u(0) = 4, u(1) = 1. Dieses besitzt 2 Lösungen.

1.7 Differenzenverfahren für Randwertprobleme

Die Konstruktion von numerischen Verfahren hatte bisher an der zur Dgl erster Ordnung äquivalenten Integralgleichung angesetzt. Bei allgemeineren Problemstellungen ist aber eine solche Umformung nicht einfach möglich bzw. unbequem, z.B. wenn ursprünglich Dgln höherer Ordnung vorliegen. Hier ist ein Zugang, der auf einer direkten Approximation der Ableitungen aufbaut, naheliegender. Die Standardmethode ist dabei die Approximation von Ableitungen durch Differenzenausdrücke. Diese können explizit konstruiert oder systematisch durch Ableitung von Interpolationspolynomen bestimmt werden (vgl. Numerik I, §5.6). Eine besonders einfache Form besitzen dabei die symmetrischen Differenzenquotienten der Ordnung 2. Hier gilt

Lemma 1.7.1 Für u ∈ C^{m+2}[−ϵ, ϵ], 0 < \frac{m}{2}h ≤ ϵ, ist

\frac{1}{h^{m}} ∑_{j=0}^{m} \binom{m}{j}(−1)^{j}u((j − \frac{m}{2})h) = u^{(m)}(0) + C_{m}h^{2}u^{(m+2)}(ξ), ξ ∈ (−ϵ, ϵ).  (1.7.1)

‘Beweis’: Eine einfache Herleitung der Formel ist durch formale Rechnung mit Operatorreihen bei C^{∞}(R)-Funktionen u (z.B. Polynomen) möglich. Mit h > 0 werden Verschiebungs- und Ableitungsoperatoren punktweise durch

(S_{h}u)(x) := u(x + h), (Du)(x) := u′(x)

definiert. Für r ∈ R ist dann (S_{h}^{r}u)(x) = (S_{rh}u)(x) = u(x + rh). Der Satz von Taylor besagt

S_{h}u(x) = u(x + h) = ∑_{j=0}^{∞} \frac{h^{j}}{j!} u^{(j)}(x) = ∑_{j=0}^{∞} \frac{h^{j}}{j!} D^{j}u(x) = e^{hD}u(x), d.h., S_{h} = e^{hD}.

Durch Differenzbildung erhält man Darstellungen für Differenzenoperatoren, etwa S_{h} − I = e^{hD} − I. Der symmetrische Differenzenausdruck in (1.7.1) folgt aus der binomischen Formel bei

S_{h/2} − S_{−h/2} = e^{\frac{h}{2} D} − e^{− \frac{h}{2} D} = 2 sinh \frac{h}{2}D ⇒  (1.7.2)
(S_{h} − I)^{m}S_{−h/2}^{m} = (S_{h/2} − S_{−h/2})^{m} = (2 sinh \frac{h}{2}D)^{m}
= (hD + \frac{h^{3}}{24}D^{3} + ..)^{m} = h^{m}(D^{m} + \frac{m}{24}h^{2}D^{m+2} + ..).

Dies entspricht (1.7.1), wobei durch einige Zusatzüberlegungen die Darstellung des (Taylor-) Restglieds vereinfacht wurde.

Von besonderem Interesse sind Differenzenverfahren bei (Systemen von) Dgln höherer Ordnung, da sich mit ihnen der Übergang zum äquivalenten System erster Ordnung, mit einer vielfachen Anzahl von Variablen, umgehen läßt. Als Vorbereitung auf die partiellen Dgln wird das häufig auftretende lineare Randwertproblem zweiter Ordnung betrachtet

− u′′(x) + p(x)u′(x) + q(x)u(x) = g(x), u(a) = α, u(b) = β.  (1.7.3)

In den beiden Randpunkten sind also jeweils die Funktionswerte vorgeschrieben. Unter der Voraussetzung p, p′, q, g ∈ C[a, b] und q^{∗} + \frac{1}{4}p^{2} − \frac{1}{2}p′ + (b − a)^{2}/π^{2} > 0, wobei q(x) ≥ q^{∗}, besitzt (1.7.3) eine eindeutige Lösung. Zur Approximation auf einem äquidistanten Gitter

Δ = \{x_{i} := a + ih, i = 0, …, N + 1\}, h := \frac{b − a}{N + 1},

werden die Ableitungen u′′, u′ durch symmetrische Differenzenquotienten ersetzt,

u′′(x_{i}) = \frac{u(x_{i+1}) − 2u(x_{i}) + u(x_{i−1})}{h^{2}} − C_{2}h^{2}u^{(4)}(ξ_{i}),
u′(x_{i}) = \frac{u(x_{i+1}) − u(x_{i−1})}{2h} − 4C_{1}h^{2}u^{(3)}(ζ_{i}).

Die Gitterwerte u_{i} := u(x_{i}), i = 0, .., N, der Lösung (analog p_{i} := p(x_{i}), q_{i} := .., g_{i} := ..) sind am Rand gegeben u_{0} = α, u_{N+1} = β, und erfüllen für i = 1, …, N das Gleichungssystem

\frac{1}{h^{2}} (−u_{i−1} + 2u_{i} − u_{i+1}) + \frac{p_{i}}{2h}(u_{i+1} − u_{i−1}) + q_{i}u_{i} = g_{i} + T_{h}(x_{i}), ⟺
\frac{1}{h^{2}} ( − (1 + \frac{h}{2}p_{i})u_{i−1} + (2 + h^{2}q_{i})u_{i} − (1 − \frac{h}{2}p_{i})u_{i+1}) = g_{i} + T_{h}(x_{i}).  (1.7.4)

Dabei enthält T_{h} die Fehler der Differenzenquotienten und besitzt die Größenordnung O(h^{2}), h → 0 (s.u.). Die Matrix dieses Gleichungsystems ist tridiagonal und hat die Form

A := \frac{1}{h^{2}} ( \begin{matrix}2 + h^{2}q_{1} & −1 + \frac{h}{2}p_{1} \\ −1 − \frac{h}{2}p_{2} & 2 + h^{2}q_{2} & −1 + \frac{h}{2}p_{2} \\ & ... & ... & ... \\ & & −1 − \frac{h}{2}p_{N−1} & 2 + h^{2}q_{N−1} & −1 + \frac{h}{2}p_{N−1} \\ & & & −1 + \frac{h}{2}p_{N} & 2 + h^{2}q_{N} \end{matrix} ) .  (1.7.5)

Durch Vernachlässigung der unbekannten lokalen Fehler T_{h} in (1.7.4) ergibt sich ein Gleichungssystem mit bekannten Koeffizienten für den Vektor y := (y_{1}, …, y_{N})T einer Näherungslösung u_{Δ}, y_{i} = u_{Δ}(x_{i}). Das System besitzt die Gestalt

Ay = g̃,  g̃ = ( g_{1} ; g_{2} ; ⋮ ; g_{N−1} ; g_{N} ) + \frac{1}{h^{2}} ( α(1 + \frac{h}{2}p_{1}) ; 0 ; ⋮ ; 0 ; β(1 − \frac{h}{2}p_{N}) ) ,  (1.7.6)

mit der Matrix A aus (1.7.5). Es ist leicht zu zeigen, daß dieses Gleichungssystem (1.7.6) eindeutig lösbar ist, da die Matrix A viele Eigenschaften besitzt, die auch das ursprüngliche RWP (1.7.3) hat (bei geeigneter Übertragung der Begriffe). Zum Beispiel ist A für p ≡ 0, q ≥ 0, symmetrisch und positiv definit, genauso wie die zum RWP gehörige lineare Abbildung (zwischen Funktionenräumen). Auf dieser Eigenschaft beruht ein Verfahren, das später behandelt wird. Zunächst wird aber nach einem anderen Beweisprinzip vorgegangen.

Satz 1.7.2 Für die Koeffizientenwerte in (1.7.5) gelte q_{i} ≥ q^{∗} > 0, h|p_{i}| ≤ 2, i = 1, …, N. Dann ist die Matrix A streng diagonaldominant, also regulär, und es gilt (elementweise)

A^{−1} ≥ 0 sowie ∥A^{−1}∥_{∞} ≤ \frac{1}{q^{∗}} .

Beweis Nach Voraussetzung sind die Nebendiagonalelemente von A = (a_{ij}) nicht-positiv. Daher gilt A = D − B mit der Diagonalmatrix D = diag(a_{ii}) = diag(\frac{2}{h^{2}} + q_{1}, .., \frac{2}{h^{2}} + q_{N}) und einer nichtnegativen Matrix B. Die Zeilensummen von D^{−1}B für i = 2, .., N − 1 sind

0 ≤ \frac{1}{a_{ii}} (a_{i,i−1} + a_{i,i+1}) = \frac{1 + \frac{h}{2}p_{i} + 1 − \frac{h}{2}p_{i}}{2 + h^{2}q_{i}} = \frac{2}{2 + h^{2}q_{i}} ≤ \frac{2}{2 + h^{2}q^{∗}} < 1.

Für i = 1, N gilt eine analoge Aussage. Damit folgt ∥D^{−1}B∥_{∞} < 1 und die Neumann-Reihe

A^{−1} = (D − B)^{−1} = ∑_{j=0}^{∞} \underbrace{(D^{−1}B)^{j}}_{≥0} D^{−1} ≥ 0

konvergiert. Durch Übergang zu Normen in dieser Reihe ergibt sich mit

∥A^{−1}∥_{∞} ≤ ∑_{j=0}^{∞} ∥D^{−1}B∥_{∞}^{j}∥D^{−1}∥_{∞} = \frac{∥D^{−1}∥_{∞}}{1 − ∥D^{−1}B∥_{∞}} ≤ \frac{h^{2}}{(2 + h^{2}q^{∗})(1 − \frac{2}{2+h^{2}q^{∗}})} = \frac{1}{q^{∗}} .

die zweite Aussage.

Der Satz erfordert ein positives q, hängt aber nicht von der Länge des Intervalls ab. Mit etwas mehr Aufwand läßt sich die Regularität von A unter Voraussetzungen zeigen, die nahe bei der Eindeutigkeitsbedingung an das RWP sind.

Der Satz ist auch von praktischem Interesse, da die Diagonaldominanz die Durchführbarkeit des Gaußalgorithmus ohne Zeilenvertauschung beim Tridiagonalsystem (1.7.6) garantiert mit O(N) Operationen. Aus A^{−1} ≥ 0 folgt, daß die Lösung y monoton von der Inhomogenität g abhängt, und daß daher, z.B., y ≥ 0 gilt, falls g ≥ 0, α, β ≥ 0. Die Normschranke ∥A^{−1}∥ ≤ \frac{1}{q^{∗}} schließlich ist eine Stabilitätsaussage, aus der direkt eine Fehlerschranke für das Verfahren folgt.

Satz 1.7.3 Für die Koeffizientenfunktionen in (1.7.3) gelte p, q, g ∈ C^{2}[a, b] und q(x) ≥ q^{∗} > 0. Die Schrittweite im Verfahren sei so, daß h|p(x)| ≤ 2 ∀x ∈ Δ. Dann gilt für den globalen Fehler der Näherung u_{Δ} aus (1.7.6) mit einer von Δ unabhängigen Konstanten C die Aussage

max_{x∈Δ} |u(x) − u_{Δ}(x)| ≤ Ch^{2}(∥u^{(4)}∥_{∞} + ∥u^{(3)}∥_{∞}).

Beweis Nach den einleitenden Bemerkungen erfüllt der Gittervektor u⃗ der Lösung u das System Au⃗ = g + τ, τ = (T_{h}(x_{1}), .., T_{h}(x_{N}))^{T}. Durch Subtraktion von (1.7.6) folgt

A(u⃗ − y) = τ ⇒ ∥u − y∥_{∞} ≤ ∥A^{−1}∥_{∞}∥τ∥_{∞} ≤ \frac{1}{q^{∗}} ∥τ∥_{∞}.

Unter den Voraussetzungen des Satzes ist u ∈ C^{4}[a, b]. Daher gilt für die Konsistenzfehler T_{h} nach (1.7.4) und Lemma 1.7.1

|T_{h}(x_{i})| = |\frac{1}{h^{2}} (−u_{i−1} + 2u_{i} − u_{i+1}) + \frac{p_{i}}{2h}(u_{i+1} − u_{i−1}) + q_{i}u_{i} − g_{i}|
≤ | \underbrace{−u′′(x_{i}) + p_{i}u′(x_{i}) − q_{i}u(x_{i}) − g_{i}}_{=0} | + h^{2}(C_{2}|u^{(4)}(ξ_{i})| + 4C_{1}|p_{i}u^{(3)}(ζ_{i})|)
≤ C′h^{2}(|p_{i}|∥u^{(3)}∥_{∞} + ∥u^{(4)}∥_{∞}).

Dieses einfache Differenzenverfahren erzeugt also Näherungen an u von der (niedrigen) Konvergenzordnung zwei. Bei gewöhnlichen Randwertproblemen ist dies im Vergleich zu den Möglichkeiten des Schießverfahrens zusammen mit Anfangswertverfahren höherer Ordnung uninteressant. Mit Hilfe der im Beweis von Lemma 1.7.1 gefundenen formalen Identität, \frac{1}{2}(S_{h/2} − S_{−h/2}) = sinh \frac{h}{2}D, aus der, z.B., folgt

D^{2} = (\frac{2}{h}arsinh\frac{1}{2}(S_{h/2} − S_{−h/2}))^{2} = \frac{1}{h^{2}} ((S_{h/2} − S_{−h/2})^{2} − \frac{1}{12}(S_{h/2} − S_{−h/2})^{4} + …),

können aber auch Differenzenapproximationen höherer Ordnung hergeleitet werden. So führen die in der letzten Gleichung angegebenen Terme auf die (wieder symmetrische) Differenzenformel

\frac{1}{h^{2}} (S_{−h} − 2I + S_{h} − \frac{1}{12}(S_{−2h} − 4S_{−h} + 6I − 4S_{h} + S_{2h}))u(x)

= \frac{1}{12h^{2}} ( − u(x − 2h) + 16u(x − h) − 30u(x) + 16u(x + h) − u(x + 2h))
= u′′(x) − \frac{h^{4}}{90}u^{(6)}(ξ).  (1.7.7)

In den randnahen Punkten x_{1}, x_{N} sind allerdings andere, unsymmetrische Formeln zu verwenden. Analoge Approximationen für u′ sind ebenfalls leicht herzuleiten. Da die Formel (1.7.7) Funktionswerte aus 5 Gitterpunkten benutzt, ergibt sich nach Approximation der Dgl anstelle der Tridiagonalmatrix (1.7.5) eine mit fünf besetzten Diagonalen. Auch gehen weitere Eigenschaften, wie Diagonaldominanz, verloren, sodaß auch die einfache Stabilitätsaussage von Satz 1.7.2 nicht mehr möglich ist.

Bei speziellen Gleichungen ohne Beteiligung der ersten Ableitung u′ (d.h. p ≡ 0) kann allerdings ein O(h^{4})-Verfahren ohne diese Nachteile konstruiert werden. Dabei wird ausgenutzt, daß hier bei Kenntnis von u(x) auch u′′(x) = q(x)u(x) − g(x) =: f(x, u(x)) bekannt ist. Daher ersetzt man den h^{2}-Korrekturterm (S_{h/2} − S_{−h/2})^{4}u(x) ≅ h^{4}u^{(4)}(x_{i}) in (1.7.7) durch h^{2}(S_{h/2} − S_{−h/2})^{2}u′′(x_{i}). Dieses Mehrstellenverfahren führt nun auf das Gleichungssystem

\frac{1}{h^{2}} (−y_{i−1} + 2y_{i} − y_{i+1}) + \frac{1}{12}(f(x_{i−1}, y_{i−1}) + 10f(x_{i}, y_{i}) + f(x_{i+1}, y_{i+1})) = 0,

i = 1, .., N, dessen Koeffizientenmatrix wieder tridiagonal und diagonaldominant ist für q^{∗} > 0.

2 Gemeinsame Prinzipien von Diskretisierungsverfahren

2.1 Konvergenz, Konsistenz, Stabilität

Den Fehlerbetrachtungen bei allen bisher behandelten Verfahren liegen gemeinsame Prinzipien zugrunde, die hier noch einmal herausgearbeitet werden sollen. Alle behandelten Probleme können als Gleichungssysteme geschrieben werden, wobei die jeweilige Lösung u implizit

F (u) = 0, u ∈ D, mit F : D ⊂ V → W  (2.1.1)

definiert ist mit Hilfe einer Abbildung F zwischen normierten Räumen V, W. Der Definitionsbereich D von F ist dabei i. a. eine echte Untermenge von V. Bei der Lösung von Differentialgleichungen ist V ein Raum von stetigen, differenzierbaren Funktionen, z.B., (C^{1}[a, b])^{n}. Die Abbildung F besteht hauptsächlich aus einem Ausdruck, der Ableitungen der eingesetzten Funktion enthält, also einem Differentialoperator. Die zusätzlichen Anfangs- bzw. Randbedingungen können entweder in die Definition von F oder in die von V bzw. D eingebaut werden. Bei den Systemen von Dgln erster Ordnung läßt sich

V = (C^{1}[a, b])^{n} mit Norm ∥v∥_{V} := max{∥v(x)∥ : x ∈ [a, b]},
W = R^{n} × (C[a, b])^{n} mit Norm ∥ ( \begin{matrix}d \\ g \end{matrix} ) ∥_{W} := ∥d∥ + max\{∥g(x)∥ : x ∈ [a, b]\}

D = V, wählen. Dabei sei die nicht-indizierte Norm eine beliebige des R^{n}. Man beachte, daß der so gewählte Raum V nicht vollständig ist, da die Norm keine Ableitungen einschließt. Als Abbildung ergibt sich beim

\underline{AWP} F : v ↦ ( \begin{matrix}v(a) − u_{0} \\ v′ − f(., v) \end{matrix} ) ,  \underline{RWP} F : v ↦ ( \begin{matrix}r(v(a), v(b)) \\ v′ − f(., v) \end{matrix} ) .  (2.1.2)

Beim Randwertproblem (1.7.3) zweiter Ordnung können homogene Randbedingungen, α = β = 0, besser in die Definition des Raums V eingebaut werden. Hier bieten sich an

V := {v ∈ C^{2}[a, b] : v(a) = v(b) = 0}, W := C[a, b],  (2.1.3)
F : v ↦ −v′′ + pv′ + qv − g.  (2.1.4)

In V und W kann jeweils die Maximumnorm von v verwendet werden.

Zur numerischen Approximation solcher Probleme wird das Originalproblem durch eine Folge von effektiv lösbaren endlichdimensionalen (diskreten) Problemen ersetzt. Das Ziel bei der Konstruktion solcher Verfahren ist, bei wachsender Dimension im diskreten Problem den Fehler zwischen diskreter und exakter Lösung immer kleiner zu machen. In den bisher behandelten Beispielen entsprachen die endlichdimensionalen Räume jeweils einer Sammlung von Funktionswerten auf einem Gitter Δ ⊂ [a, b], wobei die Raumdimension bei feiner werdendem Gitter anwächst. Stattdessen kommen aber auch etwa Räume von Polynomen oder Splinefunktionen wachsenden Grads in Frage.

Ein Diskretisierungsverfahren besteht also aus einer Schar von endlichdimensionalen Räumen V_{Δ}, W_{Δ}, dim V_{Δ} = dim W_{Δ}, und Abbildungen F_{Δ} : V_{Δ} → W_{Δ}, mit denen diskrete Lösungen durch

F_{Δ}(u_{Δ}) = 0  (2.1.5)

definiert werden. Um diese Näherungen mit der Originallösung u vergleichen zu können, wird mindestens ein Transfer von u nach V_{Δ} (Restriktion) oder von u_{Δ} nach V (Prolongation) benötigt. Bisher trat nur der erste Fall auf. Es sei also eine Restriktion bekannt, d.h., eine (lineare) Abbildung

R_{Δ} : V → V_{Δ}, womit e_{Δ} := u_{Δ} − R_{Δ}u  (2.1.6)

als Fehler der diskreten Lösung definiert wird. Normen in V_{Δ} und W_{Δ} werden beide mit ∥.∥_{Δ} bezeichnet. Diese Normen sollten in Beziehung zu den Normen der Räume V, W stehen. Dazu wird angenommen, daß ein Maß |Δ| für die Feinheit der Gitter Δ existiert mit den Eigenschaften, daß für jedes v ∈ V gilt

|Δ| → 0 ⇒ dim V_{Δ} → ∞, ∥R_{Δ}v∥_{Δ} → ∥v∥_{V} .  (2.1.7)

Für die einfachsten Gitter, die äquidistanten mit Schrittweite h = (b − a)/N, ist |Δ| := h ein solches Maß. Hier liegt sogar eine abzählbare Folge von Gittern vor, die durch die Zahl der Gitterpunkte indiziert werden kann. Die Beziehungen können durch folgendes Diagramm zusammengefaßt werden.

V ⟶^{F} W
R_{Δ} ↓
V_{Δ} ⟶^{F_{Δ}} W_{Δ}

Beispiel 2.1.1 Einschrittverfahren bei Anfangs- und Randwertproblemen. Mit Δ = {x_{i} := a + ih, i = 0, .., N}, h = (b − a)/N, N ∈ N, |Δ| := h, sei

V_{Δ} = {y_{Δ} = ( y_{0} ; ⋮ ; y_{N} ) : y_{i} ∈ R^{n}, i = 0, .., N}, ∥y_{Δ}∥_{Δ} := max_{i=0}^{N} ∥y_{i}∥, R_{Δ}v := ( v(x_{0}) ; ⋮ ; v(x_{N}) ) ,

v ∈ V, wobei dim V_{Δ} = n(N + 1). Die Restriktion R_{Δ}v besteht also aus den Gitterwerten von v. Der Bildraum W_{Δ} kann mit V_{Δ} identifiziert werden. Mit einem durch seine Verfahrensfunktion f_{h} definierten Einschrittverfahren erhält man die diskrete Abbildung

F_{Δ}(y_{Δ}) = ( y_{0} − u_{0} ; \frac{1}{h}(y_{1} − y_{0}) − f_{h}(x_{0}, y_{0}) ; ⋮ ; \frac{1}{h}(y_{N} − y_{N−1}) − f_{h}(x_{N−1}, y_{N−1}) )  (2.1.8)

beim Anfangswertproblem. Beim Randwertproblem ist nur die erste Komponente, y_{0} − u_{0}, von F_{Δ} durch r(y_{0}, y_{N}) zu ersetzen. Beim AWP kann man das Gleichungssystem F_{Δ}(u_{Δ}) = 0 natürlich schrittweise nach y_{0}, y_{1}, … auflösen. Beim RWP entspricht dieses Gleichungssystem dem der Mehrzielmethode (1.6.10) und kann im Rahmen der Newtoniteration auf ein n × n-System kondensiert werden.

Beispiel 2.1.2 Differenzenverfahren beim linearen RWP zweiter Ordnung mit homogenen Randbedingungen. Da die Randwerte durch die Randbedingungen bekannt sind, müssen sie nicht mehr als Unbekannte geführt werden. Mit dem Gitter aus Beispiel 2.1 wählt man hier

V_{Δ} = {y_{Δ} = ( y_{1} ; ⋮ ; y_{N−1} ) : y_{i} ∈ R, i = 1, .., N − 1}, ∥y_{Δ}∥_{Δ} := max_{i=1}^{N−1} |y_{i}|, R_{Δ}v := ( v(x_{1}) ; ⋮ ; v(x_{N−1}) ) ,

v ∈ V, mit dim V_{Δ} = N − 1. Das diskrete Problem hat damit die in §1.7 behandelte Form F_{Δ}(u_{Δ}) = Au_{Δ}−g̃, mit der in (1.7.5) definierten Matrix A (hier mit h = (b−a)/N, Dimensionen (N − 1) × (N − 1)).

Diese Bezeichnungen wurden vor allem eingeführt, um noch einmal die gemeinsame Struktur der alten Beweise hervorzuheben, insbesondere in Bezug auf das Zusammenspiel der Begriffe Konsistenz, Stabilität und Konvergenz, und deren Abhängigkeit von der geeigneten Wahl der Definitionen von F_{Δ} und ∥.∥_{Δ}. Im Beispiel 2.1 hätte man alle Komponenten von F_{Δ} mit h multiplizieren können, ohne das Verfahren zu ändern. In diesem Fall wäre das Verfahren zwar noch konsistent, aber nicht mehr stabil im Sinne der folgenden Definition.

Definition 2.1.3 Ein Diskretisierungsverfahren heißt konsistent (von der Ordnung p > 0), wenn mit der Lösung u ∈ V für den Konsistenzfehler

T_{Δ} := F_{Δ}(R_{Δ}u) gilt ∥T_{Δ}∥_{Δ} → 0, bzw. ∥T_{Δ}∥_{Δ} = O(|Δ|^{p}), |Δ| → 0.

Das Verfahren heißt stabil, wenn Konstanten r, S, H existieren so, daß für alle Gitter mit |Δ| ≤ H für je zwei Elemente v_{Δ}^{[j]} ∈ V_{Δ} mit F_{Δ}(v_{Δ}^{[j]}) = g_{Δ}^{[j]} und ∥g_{Δ}^{[j]}∥_{Δ} ≤ r, j = 1, 2, folgt

∥v_{Δ}^{[1]} − v_{Δ}^{[2]}∥_{Δ} ≤ S∥g_{Δ}^{[1]} − g_{Δ}^{[2]}∥_{Δ} = S∥F_{Δ}(v_{Δ}^{[1]}) − F_{Δ}(v_{Δ}^{[2]})∥_{Δ}.  (2.1.9)

Das wesentliche Merkmal der Stabilitätsforderung ist, daß die Schranke gleichmäßig gilt für die gesamte (unendliche) Folge von Gittern, die Stabilitätskonstante S also nicht von Δ abhängt. Daher bedeutet Stabilität insbesondere auch die gleichmäßig beschränkte Invertierbarkeit der Ableitungen F′_{Δ} für die ganze Folge der Gitter! Die Stabilitätsschwelle r beschränkt die Forderung auf eine nicht zu kleine Umgebung der diskreten Lösung, die der rechten Seite g_{Δ} = 0 entspricht.

Diese Definition stimmt bis auf kleinere Modifikationen mit denen von früher überein, beim AWP war das Def. 1.2.3 und Def. 1.2.5, sowie Satz 1.2.7. Beim RWP zweiter Ordnung wurden die entsprechenden Aussagen direkt nachgewiesen, Konsistenz mit p = 2 in Satz 1.7.3), Stabilität mit S = 1/q^{∗} in Satz 1.7.2. Wie in der obigen Bemerkung gilt die Normschranke dieses Satzes gleichmäßig für die Inversen einer unendliche Folge von Matrizen wachsender Dimension n → ∞.

Die Konvergenz eines Verfahrens folgt nun wieder durch Kombination der beiden Eigenschaften Konsistenz und Stabilität.

Satz 2.1.4 Das diskrete Problem (2.1.5) sei stabil nach (2.1.9) und konsistent von der Ordnung p mit ∥T_{Δ}∥_{Δ} ≤ C_{T} |Δ|^{p}. Dann konvergieren die Näherungslösungen mit Ordnung p, d.h., für alle Gitter mit |Δ| ≤ H_{0} := min{H, (r/C_{T})^{1/p}}, gilt

∥u_{Δ} − R_{Δ}u∥_{Δ} ≤ SC_{T} |Δ|^{p}.

Beweis In (2.1.9) wird gewählt v_{Δ}^{[1]} := u_{Δ} mit g_{Δ}^{[1]} = 0 und v_{Δ}^{[2]} := R_{Δ}u mit g_{Δ}^{[2]} = T_{Δ}. Nach Voraussetzung ist |Δ| ≤ H, ∥g_{Δ}^{[j]}∥_{Δ} ≤ r, j = 1, 2, und die Behauptung folgt aus (2.1.9).

2.2 Allgemeine Verfahrens-Ansätze

In diesem allgemeinen Rahmen sollen noch kurz zwei Prinzipien skizziert werden, die jeweils einen allgemeinen Ansatz zur Konstruktion von Verfahren für unterschiedliche Problembereiche darstellen.

Projektions-Verfahren Näherungen für eine isolierte Lösung von F (u) = 0, u ∈ V, werden hier in einem endlichdimensionalen Unterraum V_{Δ} ⊆ V bestimmt, z.B., einem Raum von Polynomen oder Spline-Funktionen (vgl. Numerik I, §4). Für v_{Δ} ∈ V_{Δ} kann dann F (v_{Δ}) ∈ W gebildet werden, allerdings ist das Gleichungssystem F (u_{Δ}) = 0 in der Regel unlösbar. Hier wird nun ein Projektor P_{Δ} = P_{Δ}^{2} eingeführt, und die Näherung u_{Δ} ∈ V_{Δ} definiert durch

P_{Δ}F (u_{Δ}) = 0, P_{Δ} : W → W_{Δ}.  (2.2.1)

Bei geeigneter Wahl von P_{Δ} und dim W_{Δ} = dim V_{Δ} ist dieses System (lokal) eindeutig lösbar, Stabilitätseigenschaften der Abbildung P_{Δ}F : V_{Δ} → W_{Δ} müssen für das konkrete Verfahren untersucht werden. Aufgrund der speziellen Struktur läßt sich aber die Konvergenz nach einfachen Prinzipien analysieren. Im linearen Fall, F (u) = Lu − r, kann dazu das Hilfsproblem Lv = LR_{Δ}u betrachtet werden. Es hat die Eigenschaft, daß exakte und numerische Lösung übereinstimmen, denn v = v_{Δ} = R_{Δ}u. Durch Subtraktion der Beziehungen P_{Δ}Lu_{Δ} = P_{Δ}r = P_{Δ}Lu und P_{Δ}(Lv_{Δ} − LR_{Δ}u) = 0 folgt die Gleichung P_{Δ}L(u_{Δ} − v_{Δ}) = P_{Δ}L(I − R_{Δ})u. Damit gilt

∥u − u_{Δ}∥ = ∥u − v − (u_{Δ} − v_{Δ})∥ ≤ (1 + K)∥(I − R_{Δ})u∥.

Der globale Fehler kann also bis auf die Konstante 1 + K durch den Approximationsfehler (I − R_{Δ})u bei der Lösung u abgeschätzt werden. Allerdings muß dazu wieder im konkreten Fall die Konstante K als Norm der Abbildung (P_{Δ}L|_{V_{Δ}})^{−1}P_{Δ}L : V → V_{Δ} (≠ Id!) bestimmt werden.

Spezielle Projektionsverfahren sind die in §3.3 behandelten Finite-Element-Verfahren und die Kollokationsverfahren. Ihr Name kommt von der Interpretation, daß man (bei Funktionen u = u(x), x ∈ Ω) für u_{Δ}(x) die Gleichung F (u_{Δ}(x)) = 0 nicht mehr überall fordert, ∀x ∈ Ω, sondern nur noch in endlich vielen Punkten x ∈ Δ̄ ⊆ Ω. Äquivalent dazu ist, daß man P_{Δ} als Interpolationsabbildung mit den x ∈ Δ̄ als Stützstellen definiert und (2.2.1) fordert. Als Beispiel werde eine Dgl k-ter Ordnung betrachtet, etwa (1.7.3) bei k = 2. Dazu wählt man ein Grundgitter Δ : a = x_{0} < x_{1} < … < x_{N} = b und als Ansatzraum einen Raum von Spline-Funktionen

V_{Δ} := S_{Δ,p+k}^{k−1} = {g ∈ C^{k−1}[a, b] : g|_{[x_{i},x_{i+1})} ∈ Π_{p+k}},

die stückweise Polynome vom Grad p+k sind. Dabei ist nur g ∈ C^{m−1} wirklich erforderlich, vgl. Satz von Picard-Lindelöf (T : C^{0} → C^{0}). Als Kollokationspunkte werden in jedem Teilintervall [x_{i}, x_{i+1}) jeweils p affin ähnliche Punkte verwendet Δ̄ := {x_{i}+h_{i}ξ_{j} : j = 1, …, p, i = 0, …, N − 1}. P_{Δ} entspricht der Polynom-Interpolation unabhängig in allen Teilintervallen, Funktionen P_{Δ}v sind also nur stückweise stetig. In der Skizze sind für p = 2 die Kollokationspunkte Δ̄ durch Kreise markiert, der Defekt F (u_{Δ}) verschwindet nur dort.

[figure: x_{i−1}; x_{i+1}; x_{i}; F (u_{Δ})]

Das Verfahren soll nicht weiter diskutiert werden, es sei jedoch erwähnt, daß bei Wahl der Kollokationsstellen ξ_{1}, …, ξ_{p} ∈ (0, 1) wie bei der Gauß-Quadratur (Numerik I, §5.3) sich eine Superkonvergenz auf den Gitterpunkten x_{i} einstellt, d.h., eine Konvergenzordnung, die über der Konsistenzordnung liegt. Dies war auch der Hintergrund für die hohe Ordnung in Satz 1.5.4, die Gauß-Runge-Kutta-Verfahren sind nämlich äquivalent zu Kollokationsverfahren in S_{m}^{0}.

Defektkorrektur-Verfahren Bei Verfahren höherer Ordnung reicht die Kopplung zwischen Unbekannten (Lösungs-Werten, -Koeffizienten) in der Regel weiter als bei solchen niederer Ordnung. Daher ist die Struktur der (bei RWPen) auftretenden Gleichungssysteme komplexer und die Auflösung aufwendiger. Bei den Defektkorrekturverfahren wird dies umgangen dadurch, daß ausgehend von einer Lösung niederer Ordnung schrittweise Näherungen wachsender Ordnung berechnet werden. Wegen der damit auch verfügbaren Fehlerschätzung bieten diese Verfahren als Bonus die Möglichkeit zu Gitter- und Ordnungssteuerung. Als Ansatzpunkt kann eine Newtonartige Iteration

L_{Δ}(u_{Δ}^{[k+1]} − u_{Δ}^{[k]}) = −G_{Δ}(u_{Δ}^{[k]})

für eine nichtlineare Gleichung G_{Δ}(u_{Δ}) = 0 dienen. L_{Δ} ist dabei eine Approximation an G′_{Δ}, da für die (lokale) Konvergenz die Bedingung ∥I − L_{Δ}^{−1}G′_{Δ}(u_{Δ})∥ = ∥L_{Δ}^{−1}(L_{Δ} − G′_{Δ}(u_{Δ}))∥ < 1 erforderlich ist. Diese bietet aber auch den Spielraum merkliche Abweichungen L_{Δ} − G′_{Δ} zuzulassen, wenn L_{Δ} dafür eine einfachere Struktur hat als G′_{Δ}. Z.B., kann L_{Δ} zu einem Verfahren niederer Ordnung gehören als G_{Δ}. Daher werde jetzt konkret folgende Iteration betrachtet,

F′_{Δ}(u_{Δ}^{[0]})(u_{Δ}^{[k+1]} − u_{Δ}^{[k]}) = −F̂_{Δ}(u_{Δ}^{[k]}), k = 0, 1, ….  (2.2.2)

Dabei soll F̂_{Δ} ein Verfahren höherer Ordnung realisieren, während F′_{Δ} die Ableitung zu einem Verfahren niedriger Ordnung darstellt. Als Startwert kann die Lösung u^{[0]} dieses Basisverfahrens dienen, definiert durch F_{Δ}(u^{[0]}) = 0. Dann kann in allen Schritten (2.2.2) die LR-Zerlegung von F′_{Δ}(u_{Δ}^{[0]}) wiederverwendet werden.

Zur Illustration werden Differenzenverfahren beim linearen RWP (1.7.3) mit p ≡ 0 betrachtet. Zum Basisverfahren F_{Δ} der Ordnung 2 gehört die Tridiagonalmatrix (1.7.5). Die Abbildung F̂_{Δ} zum Verfahren 4-ter Ordnung mit (1.7.7) ist für i = 2, …, N − 2

(F̂_{Δ}u_{Δ})(x_{i}) = q(x_{i})y_{i} + \frac{1}{h^{2}} (−y_{i−1} + 2y_{i} − y_{i+1}) + \frac{1}{12h^{2}} (y_{i−2} − 4y_{i−1} + 6y_{i} − 4y_{i+1} + y_{i+2}) − g(x_{i}),

(y_{i} := u_{Δ}(x_{i})), mit Modifikationen in den randnahen Punkten x_{1}, x_{N−1}. Die Ableitungsmatrix F̂′ hat ≥ 5 Diagonalen. Es sei û_{Δ} die genauere Lösung mit F̂(û_{Δ}) = 0. Betrachtet man nun den ersten Schritt der Iteration (2.2.2), ergibt sich für die Fehler dazu die Beziehung

F′_{Δ} · (u_{Δ}^{[1]} − û_{Δ}) = (F′_{Δ} − F̂′_{Δ}) · (u_{Δ}^{[0]} − û_{Δ}),  (2.2.3)

da F̂_{Δ}(u_{Δ}^{[0]}) = F_{Δ}(u_{Δ}^{[0]}) − F̂_{Δ}(û_{Δ}) = F̂′_{Δ} · (u_{Δ}^{[0]} − û_{Δ}). Wegen der Linearität aller Abbildungen ist bei diesen kein Argument angegeben. Die Differenz (F′_{Δ} − F̂′_{Δ})v besteht jetzt nur noch aus dem Differenzenausdruck 4-ter Ordnung

\frac{1}{12h^{2}} (v_{i−2} − 4v_{i−1} + 6v_{i} − 4v_{i+1} + v_{i+2}) = 2h^{2}v[x_{i−2}, …, x_{i+2}].

In (2.2.3) ist v = u_{Δ}^{[0]} − û_{Δ} = O(h^{2}). Für genügend oft differenzierbare Koeffizienten der Dgl kann man (mit einigem Aufwand) nachweisen, daß dann auch für die dividierte Differenz v[x_{i−2}, …, x_{i+2}] = O(h^{2}) gilt, sodaß insgesamt die rechte Seite von (2.2.3) die Größenordnung O(h^{4}) besitzt, d.h., schon ein Defektkorrektur-Schritt erzeugt eine h^{4}-Näherung u_{Δ}^{[1]}.

Wesentlich an dem Verfahren ist, daß der Defekt F̂_{Δ}(u_{Δ}^{[k]}) auf der rechten Seite der Iteration (2.2.2) mit einem Verfahren höherer Ordnung gebildet wird. Wenn allgemein F̂_{Δ} Konsistenz-Ordnung p hat und F_{Δ} die Ordnung 2, ergibt sich bei passender Konstruktion der Fehler u_{Δ}^{[k]} − R_{Δ}u = O(h^{min{2k+2,p}}).

3 Partielle Differentialgleichungen

3.1 Allgemeine Eigenschaften

Im Unterschied zu den gewöhnlichen Differentialgleichungen werden jetzt Funktionen u(x_{1}, .., x_{n}) von mehr als einer Veränderlichen betrachtet. Da partielle Dgln ein sehr umfangreiches Problemgebiet darstellen, das hier nur ansatzweise diskutiert werden kann, werden hier nur lineare Probleme in zwei (Raum-) Dimensionen diskutiert für Funktionen u(x, y). Bei partiellen Differentialgleichungen treten Ableitungen der gesuchten Funktion nach mehr als einer Veränderlichen auf. Daher können sie auch i.a. nicht auf eine Standardform wie (1.1.1) gebracht werden. Geht man von Systemen erster Ordnung für Funktionen u(x, y) := (u_{1}(x, y), .., u_{n}(x, y))T aus, dann ist ein recht allgemeine Form

u_{x} = A(x, y)u_{y} + B(x, y)u + g(x, y).  (3.1.1)

Charakteristische Eigenschaften dieser Dgl, nach denen sich auch die Wahl der zusätzlichen Randbedingungen zu richten hat, werden vor allem durch die n × n-Koeffizientenmatrix A bestimmt. Zwei Sonderfälle erhalten eine Bezeichnung.

Definition 3.1.1 a) Das System (3.1.1) heißt elliptisch im Punkt (x, y), wenn A(x, y) keine reellen Eigenwerte besitzt.

b) Das System (3.1.1) heißt hyperbolisch im Punkt (x, y), wenn A(x, y) reell diagonalisierbar ist.

Für n = 1 ist jede reelle Gleichung hyperbolisch. Für n > 2 sind diese beiden Fälle nicht erschöpfend. Die wichtigsten Beispiele in Dimension n = 2 für u = (v, w)T sind folgende

Beispiel 3.1.2 Die Cauchy-Riemann-Dgln

v_{x} + w_{y} = 0, w_{x} − v_{y} = 0 ⟺ u_{x} = ( \begin{matrix}0 & −1 \\ 1 & 0 \end{matrix} ) u_{y}

sind überall elliptisch (EWe ±i). Wenn die Lösung sogar zweimal differenzierbar ist, führt einmalige Differentiation der Gleichungen auf v_{xx} + w_{yx} = 0, w_{xy} − v_{yy} = 0, woraus die Potentialgleichung v_{xx} + v_{yy} = 0 für v (analog auch für w) folgt. Dies ist eine einzelne partielle Dgl zweiter Ordnung.

Beispiel 3.1.3 Aus dem hyperbolischen System (EWe ±1)

v_{x} + w_{y} = 0, w_{x} + v_{y} = 0 ⟺ u_{x} = ( \begin{matrix}0 & −1 \\ −1 & 0 \end{matrix} ) u_{y}

wird analog die Wellengleichung v_{xx} − v_{yy} = 0 für v und w.

Beispiel 3.1.4 Das System

v_{x} + w_{y} = 0, w_{x} + v = 0 ⟺ u_{x} = − ( \begin{matrix}0 & 1 \\ 0 & 0 \end{matrix} ) u_{y} + ( \begin{matrix}1 & 0 \\ 0 & 0 \end{matrix} ) u

führt auf die Wärmeleitungsgleichungen v_{xx} − v_{y} = 0 und w_{xx} − w_{y} = 0.

Das letzte Beispiel ist durch die Fälle aus Definition 3.1.1 nicht abgedeckt, sein Typ ist parabolisch. Da viele Probleme direkt als Gleichung zweiter Ordnung auftreten, wird eine vollständige Typeneinteilung nur dafür formuliert.

Definition 3.1.5 Die allgemeine lineare Dgl in zwei Variabeln x, y

a(x, y)u_{xx} + 2b(x, y)u_{xy} + d(x, y)u_{yy} + e(x, y)u_{x} + f(x, y)u_{y} = g(x, y)

heißt im Punkt (x, y)

a) elliptisch, wenn a(x, y)d(x, y) − b^{2}(x, y) > 0,
b) hyperbolisch, wenn a(x, y)d(x, y) − b^{2}(x, y) < 0,
c) parabolisch, wenn a(x, y)d(x, y) − b^{2}(x, y) = 0 und Rang ( \begin{matrix}a & b & e \\ b & d & f \end{matrix} ) = 2.

Abhängig vom Typ sind nur bestimmte Randbedingungen sinnvoll. Bei einer (überall) elliptischen Dgln, die auf einem (evtl. unbeschränkten) Gebiet Ω ⊆ R^{2} gilt, sind auf dem ganzen Rand von Ω Bedingungen an u und/oder u_{x}, u_{y} vorzuschreiben. Bei parabolischen und hyperbolischen dagegen nur auf einem Teil davon.

Beispiel 3.1.6 Auf dem L-Gebiet Ω := (−1, 1)^{2} \ [0, 1)^{2} werde in Polarkoordinaten x = r cos ϕ, y = r sin ϕ die Funktion

u(r, ϕ) = r^{2/3} sin \frac{2ϕ − π}{3} , \frac{π}{2} ≤ ϕ ≤ 2π,

definiert. Dann erfüllt u die Potentialgleichung u_{xx}+u_{yy} = 0 auf Ω und besitzt glatte Randwerte, z.B. u(x, y) = 0 für x = 0, y ≥ 0 sowie x ≥ 0, y = 0. Die ersten Ableitungen von u sind jedoch nicht beschränkt in der einspringenden Ecke im Nullpunkt r → 0. In den übrigen Ecken des L treten etwas schwächere Singularitäten auf.

u ≡ 0

Diese Eigenschaft weicht von der Situation bei den gewöhnlichen Dgln ab und führt auf eine neue Rahmenbedingung bei der Diskussion numerischer Verfahren. Wenn eine gewöhnliche Dgl genügend oft differenzierbare Koeffizienten besitzt, kann durch Ableitung der Dgl auf eine entsprechende Differenzierbarkeit der Lösung geschlossen werden. Dies machte den Einsatz von Verfahren hoher Ordnung erst sinnvoll. Bei den partiellen Dgln dagegen kann die Lösung schon bei ganz einfachen Problemen Singularitäten (von Ableitungen) besitzen. Aus diesem und auch anderen, praktischen Gründen macht bei partiellen Dgln der Einsatz numerischer Verfahren hoher Ordnung Schwierigkeiten. Daher werden v.a. Verfahren der Ordnung zwei oder vier behandelt.

3.2 Differenzenverfahren für elliptische Randwertprobleme

3.2.1 Die Poissongleichung auf einfachen Gebieten

Es sei Ω ⊆ R^{2} ein beschränktes Gebiet, d.h. eine offene und zusammenhängende Menge und Γ = ∂Ω sein Rand, der aus stückweise stetig differenzierbaren Kurven bestehen soll. Aus praktischen Gründen werden in diesem Paragraphen zunächst nur Polygone mit achsenparallelen oder diagonalen Kanten behandelt. Ω̄ = Ω ∪ Γ ist der Abschluß von Ω. Mit einer stetigen Funktion g(x, y) lautet die Poissongleichung

− u_{xx} − u_{yy} = g für (x, y) ∈ Ω.  (3.2.1)

Diese Gleichung ist elliptisch, sie stellt sogar die Normalform einer elliptischen Dgl dar. Physikalisch beschreibt sie, u.a., die Wärmeverteilung in einem homogenen Körper mit Wärmequellen g. Die homogene Gleichung (g ≡ 0) heißt Potentialgleichung, ihre Lösungen harmonische Funktionen (→ Funktionentheorie). Zur Definition einer sinnvolle Problemstellung sind bei dieser Dgl auf ganz Γ Randwerte vorzugeben. Historisch unterscheidet man drei Arten von Randbedingungen, die den Funktionswert oder die Ableitung ∂u/∂n in Richtung der äußeren Normalen n enthalten. Dazu wird der Rand disjunkt zerlegt in Γ = Γ_{d} ∪ Γ_{n} ∪ Γ_{c} und auf jedem Teil eine der folgenden Bedingungen gefordert mit stetigen Funktionen ϕ, ψ.

u(x, y) = ϕ_{d}(x, y), (x, y) ∈ Γ_{d} (Dirichlet − RB)  (3.2.2)
\frac{∂u}{∂n}(x, y) = ϕ_{n}(x, y), (x, y) ∈ Γ_{n} (Neumann − RB)  (3.2.3)
\frac{∂u}{∂n}(x, y) + ψ(x, y)u(x, y) = ϕ_{c}(x, y), (x, y) ∈ Γ_{c} (Cauchy − RB).  (3.2.4)

In der obigen physikalischen Interpretation entspricht die Dirichlet-Randbedingung einer Temperaturvorgabe auf Γ_{d}, die Neumann-RB (mit ϕ_{n} ≡ 0) einer isolierten Wand Γ_{n}, während die Cauchy-Bedingung sich als Wärmeabstrahlung durch Γ_{c} interpretieren läßt. In der Praxis kann natürlich auch jeweils nur eine der Randbedingungen auftreten. In diesen Fällen redet man von dem Dirichletschen (Γ = Γ_{d}) bzw. Neumannschen (Γ = Γ_{n}) Randwertproblem.

Definition 3.2.1 Eine Funktion u(x, y) heißt (klassische) Lösung des Randwertproblems (3.2.1), (3.2.2-3.2.4), wenn u ∈ C^{2}(Ω)∩C(Ω̄)∩C^{1}(Ω∪Γ_{n} ∪Γ_{c}) gilt und die Gleichungen (3.2.1), (3.2.2-3.2.4) erfüllt sind.

Während beim Dirichletproblem eine eindeutige Lösung existiert, besitzt das Neumannproblem Lösungen nur unter der Voraussetzung ∫_{Γ} ϕ_{n} ds = ∫ ∫_{Ω} g(x, y) dxdy. Die Lösung ist dann auch nur bis auf eine Konstante eindeutig.

Wie im eindimensionalen Fall (§1.7) lassen sich Näherungsverfahren gewinnen, indem man die Ableitungen in der Dgl durch Differenzenquotienten ersetzt. Mit Schrittweiten h_{x}, h_{y} > 0 gilt für u ∈ C^{4} nach Lemma 1.7.1

u_{xx}(x, y) = \frac{1}{h_{x}^{2}} [u(x − h_{x}, y) − 2u(x, y) + u(x + h_{x}, y)] + O(h_{x}^{2})

u_{yy}(x, y) = \frac{1}{h_{y}^{2}} [u(x, y − h_{y}) − 2u(x, y) + u(x, y + h_{y})] + O(h_{y}^{2})

Bei Addition dieser beiden Näherungen werden zur Approximation der Poissongleichung in (x, y) fünf sternförmig verteilte Funktionswerte miteinander verknüpft. Für den Standardfall h_{x} = h_{y} = h haben die vier Strahlen des Sterns alle die Länge h.

[figure: y − h; y + h; x − h; x + h]

Als erstes Modellproblem wird das Einheitsquadrat Ω = (0, 1) × (0, 1) betrachtet. Wie im eindimensionalen lassen sich mit einer Differenzenapproximation am einfachsten Funktionswerte eines äquidistanten Gitters verknüpfen. Dazu sei mit h := 1/m, m ∈ N, das Gitter

Δ_{h} := {(x_{i}, y_{j}), i, j = 0, .., m}, mit x_{i} = y_{i} := ih, i = 0, .., m.  (3.2.5)

[figure: x_{0}; x_{1}; x_{i}; x_{m}; y_{0}; y_{j}; y_{m−1}; y_{m}; P; W; O; S; N]

Da die Handhabung im zweidimensionalen schwieriger ist als früher, wird diese Indizierung der Gitterpunkte mit zwei Indizes auch für die Näherungswerte v_{ij} = u_{Δ}(x_{i}, y_{j}) ≅ u(x_{i}, y_{j}) benutzt. Zur Beschreibung der Nachbarwerte eines speziellen Punktes P = (x_{i}, y_{j}) ist es jedoch einprägsamer, Himmelsrichtungen zu zu verwenden. Also gelte, zu

v_{P} = v_{i,j} sei v_{N} := v_{i,j+1}, v_{O} := v_{i+1,j}, v_{S} := v_{i,j−1}, v_{W} := v_{i−1,j},  (3.2.6)

wie in der obigen Skizze angedeutet. Die Poissongleichung im Punkt P, −u_{xx}(P)−u_{yy}(P) = g(P) wird nun approximiert durch die Differenzengleichung −\frac{1}{h^{2}} [v_{W} −2v_{P} +v_{O}]− \frac{1}{h^{2}} [v_{S} −2v_{P} +v_{N}] = g_{P} mit g_{P} := g(x_{i}, x_{j}), d.h.,

4 -1 -1 -1  -1·\frac{1}{h^{2}}  \frac{1}{h^{2}} [4v_{P} − v_{N} − v_{O} − v_{S} − v_{W}] = g_{P} .  (3.2.7)

Die Linearkombination in eckigen Klammern wird gerne kompakt durch den links dargestellten Fünf-Punkte-Stern beschrieben. Bei der (homogenen) Potentialgleichung kann diese Beziehung auch in der Form

v_{P} = \frac{1}{4}[v_{N} + v_{O} + v_{S} + v_{W}]

ausgedrückt werden, die aussagt, daß v_{P} der Mittelwert der Funktionswerte in den vier Nachbarpunkten ist (Interpretation: ‘aufgespanntes Gumminetz’).

Dirichletproblem: Hier sind die Funktionswerte der Lösung auf dem ganzen Rand bekannt. Daher können die auf den Rand fallenden Variablen in (3.2.7) durch die entsprechenden Werte von ϕ_{d} ersetzt werden. Aufgrund der einfachen Gestalt des Einheitsquadrats Ω ergeben sich so für die Werte auf den im Bild bei (3.2.5) durch einen Punkt markierten inneren (m − 1)^{2} Gitterpunkte gerade (m−1)^{2} Gleichungen (3.2.7). In Anlehnung an die Formulierung des RWPs (3.2.1),(3.2.2), kann man durch Einführung der Bezeichnungen

Ω_{h} := {(x_{i}, y_{j}) : i, j = 1, .., m − 1}, Γ_{h} := {(x_{i}, y_{j}) : i = 0 ∨ i = m ∨ j = 0 ∨ j = m}, (3.2.8)

d.h., Ω_{h} ∪ Γ_{h} = Δ_{h}, dieses Gleichungssystem ebenfalls in kompakter Form schreiben:

\frac{1}{h^{2}} [4v_{ij} − v_{i,j+1} − v_{i+1,j} − v_{i,j−1} − v_{i−1,j}] = g_{ij},  für (x_{i}, y_{j}) ∈ Ω_{h},
v_{ij} = ϕ_{d}(x_{i}, y_{j})  für (x_{i}, y_{j}) ∈ Γ_{h}.

(3.2.9)

Wie erwähnt, werden die trivialen Randgleichungen dazu benutzt, die entsprechenden Unbekannten aus den Gleichungen auf Ω_{h} zu eliminieren (vgl. auch (1.7.6)). Dann liegt ein lineares Gleichungssystem aus (m − 1)^{2} Gleichungen für ebenso viele Unbekannte vor. Bei den partiellen Dgln ist es angebracht, in der Formulierung wieder stärker zwischen der mit dem Gleichungssystem zusammenhängenden linearen Abbildung L_{Δ} von Gitterfunktionen und den zugehörigen Matrizen zu unterscheiden. Wird zu Gitterfunktionen u_{Δ} : Ω_{h} → R die lineare Abbildung L_{Δ} punktweise definiert durch

(L_{Δ}u_{Δ})(x, y) := \frac{1}{h^{2}} [4u_{Δ}(x, y) − u_{Δ}(x, y + h) − u_{Δ}(x, y − h) − …], (x, y) ∈ Ω_{h},

mit Modifikationen in den randnahen Punkten, entspricht das Gleichungssystem (3.2.9) einem Problem der Form

L_{Δ}u_{Δ} = g_{Δ}, auf Ω_{h},

mit der schon verwendeten Abkürzung v_{ij} = u_{Δ}(x_{i}, y_{j}). ‘Die’ Matrix dieses linearen Systems hängt zumindest von einer Numerierung der bisher doppelt indizierten Unbekannten ab.

Bei zeilenweiser Numerierung (lexikographischer) des Gitters erhalten die drei Unbekannten v_{W} , v_{P} , v_{O} aufeinanderfolgende Nummern bzw. Indizes, während die Nummern von v_{S}, v_{N} genau um m − 1 kleiner bzw. größer sind als die von v_{P}. Die Nummern der Gleichungen (3.2.7)/(3.2.9) seien die von v_{P}, sodaß also das Hauptdiagonalelement immer 4/h^{2} ist. Dann ist die zu dieser Numerierung gehörende Matrix A sehr regelmäßig aufgebaut, sie hat eine Block- und Bandstruktur mit Bandbreite 2m − 1 in der Form

A = \frac{1}{h^{2}} ( \begin{matrix}T & −I \\ −I & T & −I \\ & −I & T & −I \\ . & . & . \\ & & & −I & T \end{matrix} ) , T := ( \begin{matrix}4 & −1 \\ −1 & 4 & −1 \\ & −1 & 4 & −1 \\ . & . & . \\ & & & −1 & 4 \end{matrix} ) .  (3.2.10)

Dabei sind die Matrizen I und T jeweils von der Größe (m − 1) × (m − 1), die Gesamtmatrix A von der Größe (m − 1)^{2} × (m − 1)^{2}.

Stellt man sich andererseits vor, daß die Gitterpunkte wie ein Schachbrett schwarz-weiß eingefärbt sind, erkennt man, daß in (3.2.9) die Nachbarpunkte N, O, S, W immer die zum Punkt P entgegengesetzte Farbe besitzen. Wenn daher zunächst weiße und dann schwarze Punkte numeriert werden, ergibt sich im Großen eine Blockstruktur der Form

\frac{1}{h^{2}} ( \begin{matrix}4I & −M \\ −MT & 4I \end{matrix} ) ,

wobei M in jeder Zeile höchstens vier Einsen und sonst Nullen enthält. Wenn m−1 ungerade ist, haben die Blöcke 4I geringfügig unterschiedliche Größen. Weitere Eigenschaften dieser Matrizen werden später diskutiert. Lösungsverfahren für solche Gleichungssysteme, die die geringe Anzahl der nichttrivialen Matrixelemente nutzen, sind Thema der Vorlesung Numerik IIA.

Neumann-Randbedingung: Hier ist der Wert der Normalenableitung vorgeschrieben. Daher muß auch bei der Randbedingung eine Approximation durch Differenzenquotienten vorgenommen werden. Beim Einheitsquadrat sind die Normalen parallel zu den Achsen. Daher gilt entweder ∂u/∂n = ±u_{x} (auf der Ost- und Westseite), oder ∂u/∂n = ±u_{y} (auf der Nord- und Südseite). Als Beispiel werde die Approximation der Randbedingung u_{x}(1, y) = ϕ_{n}(1, y), y ∈ (0, 1) betrachtet. Dabei bieten sich folgende Möglichkeiten an:

1. Verwendung des Gitters (3.2.5), Approximation der Randbedingung durch den einseitigen Differenzenquotienten mit Schrittweite h: Im Beispiel sei P = (1 − h, y_{j}) ein dem Rand benachbarter Gitterpunkt. Die Randbedingung wird also ersetzt durch \frac{1}{h}(v_{O} − v_{P}) = ϕ_{n}(O). Da die Variable v_{O} sonst nur noch im Differenzenstern von P auftaucht, kann die Bedingung dazu benutzt werden, diese Unbekannte sofort aus dem Stern zu eliminieren. Dann hat das Gleichungssystem wieder die gleiche Struktur wie beim Dirichletproblem, z.B. (3.2.10), nur die Sterne in den Randnachbarpunkten sind abgeändert, im Beispiel ergibt sich statt (3.2.7) jetzt

\frac{1}{h^{2}} (3v_{P} − v_{N} − v_{S} − v_{W}) = g_{p} + \frac{1}{h}ϕ_{n}(O).

Nachteilig ist hierbei die schlechtere Approximation der Randbedingung mit einem O(h)- Fehler wegen der Verwendung unsymmetrischer Differenzenquotienten.

2. Der symmetrische Differenzenquotient \frac{1}{2h}(u(x+h, y)−u(x−h, y)) ist eine O(h^{2})-Approximation für den Wert u_{x}(x, y). Um ihn einsetzen zu können, muß der entsprechende Randpunkt mit ins Berechnungsgitter Ω_{h} aufgenommen werden, dies sei jetzt P = (1, y_{j}) ∈ Γ_{n}. Die approximierte Randbedingung

\frac{1}{2h}(v_{O} − v_{W}) = ϕ_{n}(P)  (3.2.11)

verwendet dann den Punkt O außerhalb von Ω. Die Unbekannte v_{O} kann aber, wie vorher, aus dem Fünfpunktestern im Punkt P eliminiert werden. Dies führt für v_{P} auf den einseitigen Stern

4 -2 -1 -1 ·\frac{1}{h^{2}}  \frac{1}{h^{2}} (4v_{P} − 2v_{W} − v_{N} − v_{S}) = g_{p} + \frac{2}{h}ϕ_{n}(P).

Diese Methode ist besonders naheliegend, da (homogene) Neumann-Randbedingungen oft aufgrund von Symmetrieüberlegungen beim Dirchletproblem auftreten. Zum Beispiel gilt für die Lösung des Problems

−u_{xx} − u_{yy} = 1 auf (−1, 1) × (−1, 1), u = 0 auf dem Rand,

die Aussage u(x, y) = u(±x, ±y), also u_{x}(0, y) ≡ 0, u_{y}(x, 0) ≡ 0. Daher kann dieses Problem auf ein Viertel des usprünglichen Quadrats eingeschränkt werden, etwa das Einheitsquadrat, wobei auf dessen linken und unteren Rand die Neumann-Bedingung (3.2.3) mit ϕ_{n} ≡ 0 gilt. Wird diese mit Hilfe von (3.2.11) diskretisiert, ergibt sich ebenfalls eine Symmetriebedingung für v = u_{Δ}, nämlich in P = (0, y_{j}) die Gleichung \frac{1}{h}(v_{O} − v_{W}) = 0 ⟺ v_{O} = v_{W} und in P = (x_{i}, 0) die Bedingung v_{N} = v_{S}. Diese Approximation der Randbedingung führt daher auf das gleiche Ergebnis, das man durch Berücksichtigung der Symmetrie beim diskreten Problem (3.2.9) erhalten würde. Die Symmetriebetrachtung reduziert dabei die Größe des Gleichungssystems auf ein Viertel. Die Struktur der Matrix zu diesem reduzierten Problem

− u_{xx} − u_{yy} = g in Ω = (0, 1) × (0, 1),  (3.2.12)
u = 0 in Γ_{d} = {(x, y) ∈ Ω̄ : x = 1 ∨ y = 1},  (3.2.13)
\frac{∂}{∂n}u = 0 in Γ_{n} = \{(x, y) : x = 0, y ∈ [0, 1) ∨ y = 0, x ∈ [0, 1)\} (3.2.14)

ist ähnlich zu (3.2.10), wobei allerdings in den Nebendiagonalen auch der Wert −2/h^{2} auftritt. Dies zerstört die Symmetrie der Matrix. Durch Skalierung der Gleichungen (Multiplikation mit \frac{1}{2}, \frac{1}{4}) kann die Symmetrie wiederhergestellt werden. Dies entspricht der Verwendung folgender Differenzensterne auf Γ_{n} (Vorfaktor 1/h^{2}):

2 -\frac{1}{2} -\frac{1}{2} -1  ^{2 -\frac{1}{2}  -\frac{1}{2} -1}  ^{1 -\frac{1}{2} -\frac{1}{2}}

Die Matrix ist auch etwas größer, da jetzt unbekannte Werte auf dem größeren Berechnungsgitter Ω_{h} := {(x_{i}, y_{j}) : i, j = 0, .., m − 1} auftreten. Bei zeilenweiser Numerierung dieser Unbekannten v_{ij} ≅ u(x_{i}, y_{j}), in der Reihenfolge v_{00}, .., v_{m−1,0}, v_{0,1}, .., v_{m−1,1}, …, hat die Systemmatrix zu (3.2.12) die Gestalt

A = \frac{1}{h^{2}} ( \begin{matrix}\frac{1}{2}T & −E \\ −E & T & −E \\ & −E & T & −E \\ . & . & . \\ & & & −E & T \end{matrix} ) ,  (3.2.15)

T := ( \begin{matrix}2 & −1 \\ −1 & 4 & −1 \\ & −1 & 4 & −1 \\ . & . & . \\ & & & −1 & 4 \end{matrix} ) , E := ( \frac{1}{2} 1 1 ... 1 ) .

Die Größe von T und E ist m × m, die der Gesamtmatrix A ist m^{2} × m^{2}.

Die in (3.2.10) und (3.2.15) beobachtete Symmetrie der Matrizen macht den Einsatz effizienterer Verfahren möglich und ist auch bei der Analyse (Stabilität) hilfreich.

3.2.2 Eigenschaften der Differenzen-Matrizen

Zwei Matrix-Eigenschaften und damit verbundene Beweistechniken spielen bei der Behandlung elliptischer Probleme eine wichtige Rolle. Die erste Eigenschaft ist Symmetrie und Definitheit.

Definition 3.2.2 Eine Matrix A = (a_{ij}) ∈ R^{n×n} heißt symmetrisch, wenn A = AT gilt, dh. a_{ij} = a_{ji}∀i, j = 1, .., n. Sie heißt positiv definit bzw. semi-definit, wenn sie symmetrisch ist und ein α > 0 bzw. α ≥ 0 existiert so, daß

vTAv ≥ α vTv ∀v ∈ R^{n}.  (3.2.16)

Bei einer symmetrisch, positiv definiten (’SPD-’) Matrix sind alle Eigenwerte reell positiv, insbesondere ist (3.2.16) mit dem kleinsten Eigenwert, α = min λ(A), erfüllt. Hier ist die Verwendung der Euklidnorm und der dazu gehörigen Matrixnorm,

∥v∥_{2} = \sqrt{vTv}, v ∈ R^{n}, ∥A∥_{2} := \sqrt{ϱ(ATA)},

der Spektralnorm angebracht. Dabei ist ϱ(M) der Spektralradius der Matrix M, d. h. der Maximalbetrag ihrer Eigenwerte. Da jede symmetrische Matrix unitär ähnlich zu einer reellen Diagonalmatrix ist, gilt ∥A∥_{2} = ϱ(A). Zur Abschätzung der Lösungen von linearen Gleichungssystemen gilt bei Definitheit die einfache Aussage,

A SPD ⇒ A^{−1} SPD , ∥A^{−1}∥_{2} = ϱ(A^{−1}) ≤ \frac{1}{α}, α aus (3.2.16).

Diese läßt sich auch auf allgemeinere Matrizen übertragen.

Satz 3.2.3 Für die (i.a. nichtsymmetrischen) Matrizen A_{j} ∈ R^{n×n}, j = 1, …, k, gelte vTA_{j}v ≥ α_{j}vTv ∀v ∈ R^{n}, und α_{1} + … + α_{k} > 0. Dann ist A = A_{1} + … + A_{k} invertierbar und

∥A^{−1}∥_{2} = ∥(A_{1} + … + A_{k})^{−1}∥_{2} ≤ \frac{1}{α_{1} + … + α_{k}} .

Beweis Es sei 0 ≠ v ∈ R^{n} beliebig. Aus der Voraussetzung folgt durch Addition vTAv ≥ αvTv wo α = α_{1} + … + α_{k} > 0. Mit der Cauchy-Schwarz-Ungleichung führt dies auf die Aussage α∥v∥^{2} ≤ vTAv ≤ ∥v∥_{2}∥Av∥_{2} ⇒ ∥Av∥_{2} ≥ α∥v∥_{2}. Für v ≠ 0 ist insbesondere immer Av ≠ 0. Daher ist A regulär und mit v = A^{−1}w folgt ∥A^{−1}w∥_{2} ≤ \frac{1}{α}∥w∥_{2}. Bei der Voraussetzung des Satzes spielt insbesondere nur der symmetrische Anteil

S(A) := \frac{1}{2}(A + AT)

der Matrix A eine Rolle, da vTAv = (vTAv)T = vTATv = vTS(A)v. Die entscheidende Voraussetzung des Satzes ist also die positive Definitheit des symmetrischen Anteils von A, S(A) = S(A_{1}) + … + S(A_{k}). Diese folgt für die Matrizen des letzten Paragraphen einfach aus der Tatsache, daß die Matrix A des Systems (3.2.9) nach Herleitung eine Summe A = A_{1} + A_{2} ist, wobei A_{1}, A_{2} die Differenzenquotienten in x- bzw. y-Richtung enthalten. Beim diskreten Dirichlet-Problem (3.2.9) etwa wird die quadratische Form vTAv betrachtet mit A aus (3.2.10),

\frac{1}{h^{2}} ∑_{i,j=1}^{m−1} v_{ij}(4v_{ij} − v_{i−1,j} − v_{i+1,j} − v_{i,j−1} − v_{i,j+1})

= \frac{1}{h^{2}} ∑_{j=1}^{m−1} ∑_{i=1}^{m−1} v_{ij}(−v_{i−1,j} + 2v_{ij} − v_{i+1,j}) + \frac{1}{h^{2}} ∑_{i=1}^{m−1} ∑_{j=1}^{m−1} v_{ij}(−v_{i,j−1} + 2v_{ij} − v_{i,j+1}).

Randwerte bei v sind hier null zu setzen. Zu den einzelnen Summen gilt

Lemma 3.2.4 Gegeben seien v_{1}, .., v_{m} ∈ R mit v_{m} = 0. Dann gilt für a ∈ {1, 2}

v_{1}(av_{1} − v_{2}) + ∑_{j=2}^{m−1} v_{j}(−v_{j−1} + 2v_{j} − v_{j+1}) = (a − 1)v_{1}^{2} + ∑_{j=2}^{m−1}(v_{j} − v_{j−1})^{2} + v_{m−1}^{2} ≥ C_{a} ∑_{j=1}^{m−1} v_{j}^{2},

mit C_{2} = 4 sin^{2} \frac{π}{2m}, C_{1} > 0.

Beweis Links steht die quadratische Form vTT v der symmetrischen Tridiagonalmatrix

T = ( \begin{matrix}a & −1 \\ −1 & 2 & −1 \\ . & . & . \\ & & −1 & 2 \end{matrix} ) .

Die erste Identität ist elementar, die Ungleichung gilt, da Quadrate summiert werden. Für a = 2 ist die Konstante bekannt, da die Eigenwerte von T gleich 4 sin^{2} \frac{jπ}{2m} sind.

Satz 3.2.5 Die Matrix A (3.2.10) des Gleichungssystems (3.2.9) ist positiv definit mit ∥A^{−1}∥_{2} ≤ C_{h}, C_{h} = \frac{h^{2}}{8} / sin^{2}(hπ/2) ≤ \frac{1}{16} für h ≤ \frac{1}{2}. Insbesondere gilt für die Lösung des diskreten Dirichletproblems (3.2.9) mit ϕ_{d} = 0 die Aussage

∑_{i,j=1}^{m−1} v_{ij}^{2} ≤ C_{h}^{2} ∑_{i,j=1}^{m−1} g_{ij}^{2}.

Die Symmetrie und Definitheit hat auch praktische Bedeutung, da deswegen zur Lösung der linearen Gleichungssysteme effizientere Verfahren als im unsymmetrischen Fall eingesetzt werden können (Cholesky-Zerlegung, spezielle Iterationsverfahren). Ihre Bedeutung im Theoretischen ist nicht so entscheidend, da die Euklidnorm zwar hier die schärfsten Schranken liefert, diese Norm aber bei (diskreten) Funktionen weniger aussagekräftig ist als in der Maximumnorm. Für Abschätzungen in dieser Norm ist eine andere Eigenschaft entscheidend, die mit den Vorzeichen der Matrixelemente zu tun hat. Im folgenden sind Ungleichungszeichen (≤, ≥, >, <) zwischen Vektoren und Matrizen sowie Beträge elementweise zu verstehen, |v| = (|v_{i}|), |A| = (|a_{ij}|).

Definition 3.2.6 Eine Matrix A = (a_{ij}) ∈ R^{n×n} heißt M-Matrix, wenn

a_{ii} > 0 ∀i = 1, .., n, a_{ij} ≤ 0 ∀i ≠ j,  (3.2.17)
A regulär und A^{−1} ≥ 0.  (3.2.18)

Nach Zerlegung der Matrix A = D − N in Diagonale D und Nebendiagonalen −N lautet die Vorzeichenbedingung (3.2.17) kürzer N ≥ 0, a_{ii} > 0∀i. Bei Normbetrachtungen in der Euklidnorm sind Eigenvektoren extremal. Diese Rolle wird hier von den positiven Vektoren übernommen. Wichtige Schlußweisen für nichtnegative Vektoren v ≥ 0 bzw. Matrizen B ≥ 0 sind ∥v∥_{∞} ≤ a ⟺ v ≤ ae mit e = (1, .., 1)T und ∥B∥_{∞} ≤ b ⟺ Be ≤ be. Außerdem übertragen sich Ungleichungen wie v ≤ w ⇒ Bv ≤ Bw. Kriterien für M-Matrizen enthält

Lemma 3.2.7 a) Eine Matrix A = D−N mit der Vorzeichenverteilung (3.2.17) ist genau dann M-Matrix, wenn ϱ(D^{−1}N) < 1 gilt.

b) Eine Matrix A mit der Vorzeichenverteilung (3.2.17) ist (genau dann) M-Matrix, wenn ein Vektor z > 0 existiert mit w := Az > 0. Damit gilt

∥A^{−1}∥_{∞} ≤ \frac{max_{i} z_{i}}{min_{i} w_{i}} .  (3.2.19)

Beweis a) Für ϱ(D^{−1}N) < 1 konvergiert die Neumannreihe

A^{−1} = (I − D^{−1}N)^{−1}D^{−1} = ∑_{j=0}^{∞} \underbrace{(D^{−1}N)^{j}}_{≥0} D^{−1},

die Inverse existiert und ist nichtnegativ. Sei umgekehrt A M-Matrix und v ≠ 0 ein Eigenvektor von D^{−1}N zum Eigenwert λ. Wegen N ≥ 0 gilt |Nv| ≤ N|v|. Aus λDv = Nv folgt daher

|λDv| = |λ|D|v| ≤ N|v| ⇒ (D − N)|v| ≤ (1 − |λ|)D|v|.

Multiplikation mit A^{−1} ≥ 0 führt dann auf die Ungleichung |v| ≤ (1 − |λ|)A^{−1}D|v|, die für |λ| ≥ 1 zum Widerspruch 0 ≤ |v| ≤ 0 mit v ≠ 0 führt.

b) Zu z > 0 wird die reguläre Matrix Z := diag(z_{i}) definiert. Damit ist

Az = Dz − Nz > 0 ⟺ NZe < DZe ⟺ Z^{−1}D^{−1}NZe < e.

Die letzte Ungleichung für diese nichtnegative Matrix bedeutet aber gerade ∥Z^{−1}D^{−1}NZ∥_{∞} < 1, woraus insbesondere ϱ(Z^{−1}D^{−1}NZ) = ϱ(D^{−1}N) < 1 folgt. Nach Teil a) ist A daher M-Matrix. Mit α := min_{i} w_{i} gilt laut Voraussetzung auch Az = w ≥ αe. Aus A^{−1} ≥ 0 folgt damit αA^{−1}e ≤ z ≤ ∥z∥_{∞}e. Dies ist äquivalent zur Behauptung (3.2.19).

Für Matrizen mit (3.2.17) ist die strenge Diagonaldominanz

a_{ii} > ∑_{j≠i} |a_{ij}| ∀i,

ein spezielles M-Matrix-Kriterium und mit z = e in Lemma 3.2.7 in enthalten. Abschwächungen der Diagonaldominanz sind aufwendiger in der Formulierung und liefern nur die Regularitätsaussage, aber keine Schranke wie (3.2.19). Der in Lemma 3.2.7 benötigte Testvektor z kann bei L_{Δ} einfach angegeben werden. Dies führt direkt auf eine Stabilitätsaussage für die betrachteten Modellprobleme auf dem Einheitsquadrat.

Satz 3.2.8 Für das diskrete Dirchletproblem (3.2.9) auf dem Einheitsquadrat mit homogenen Randbedingungen gilt ∥L_{Δ}^{−1}∥_{∞} ≤ \frac{1}{8}, d.h., die Lösung v erfüllt

max_{i,j=1}^{m−1} |v_{ij}| ≤ \frac{1}{8} max_{i,j=1}^{m−1} |g_{ij}|,

Beweis Für z_{Δ}(x, y) := r − (x − \frac{1}{2})^{2} − (y − \frac{1}{2})^{2} ergibt sich sofort L_{Δ}z_{Δ}(x, y) ≥ 4 mit Gleichheit in allen randfernen Punkten. Außerdem ist 0 < z_{Δ} ≤ r in Ω_{h} für r = \frac{1}{2}. Aus Lemma 3.2.7 folgt dann die Aussage.

Für das Problem (3.2.12) mit Neumann-Bedingung auf zwei Seiten des Quadrats gilt eine analoge Aussage mit der größeren Konstanten 1, beim vollständigen Neumann-Problem ist die Matrix natürlich singulär.

3.2.3 Diskretisierungsfehler, Konvergenz

Da sich die Konsistenz des Differenzenverfahrens (für C^{4}-Lösungen) aus Lemma 1.7.1 ergibt und die Stabilität im letzten Abschnitt geklärt wurde, kann sofort auf die Konvergenz der Näherungen geschlossen werden, wenn das Randwertproblem hinreichend glatte Lösungen besitzt. Allerdings ist dies nicht immer der Fall, wie in Beispiel 3.1 gezeigt wurde. Daher wird zunächst eine abgeschwächte Konsistenzaussage behandelt, danach Möglichkeiten zur Beeinflussung des Fehlers durch Gitteranpassung oder Verbesserung des Verfahrens. Bei der folgenden Aussage wird der eigentlich nur auf Gitterfunktionen (Ω_{h} ⊆ Ω) definierte Differenzenoperator L_{Δ} auch auf normale Funktionen angewendet (formal wäre eine Gitter-Restriktion R_{Δ} erforderlich, L_{Δ}R_{Δ}u, vgl. §2).

Satz 3.2.9 Es sei u ∈ C^{k}(Ω), 2 ≤ k ≤ 4. Dann gilt für (x, y) ∈ Ω_{h} die Konsistenz-Aussage

|(L_{Δ}u)(x, y) + u_{xx}(x, y) + u_{yy}(x, y)| ≤  (3.2.20)

Ch^{k−2}( max\{|\frac{∂^{k}}{∂x^{k}} u(ξ, y)| : |ξ − x| < h\} + max\{|\frac{∂^{k}}{∂y^{k}} u(x, η)| : |η − y| < h\}).

Beweis Taylorentwicklung von u um (x, y).

Wird im Satz die Lösung u der Poissongleichung eingesetzt, ergibt sich für den lokalen Fehler T_{Δ} := L_{Δ}u − g_{Δ} die Größenordnung T_{Δ} = O(h^{2}), wenn u ∈ C^{4} gilt. Wegen der Linearität von L_{Δ} führt die Subtraktion von (3.2.9), auf eine Gleichung für den Fehler u_{Δ} − u

L_{Δ}u_{Δ} = g_{Δ} ⇒ −T_{Δ} = g_{Δ} − L_{Δ}u = L_{Δ}u_{Δ} − L_{Δ}u = L_{Δ}(u_{Δ} − u).  (3.2.21)

Auf die tatsächliche Größe dieses globalen Fehlers kann man daraus mit den Stabilitäts-Sätzen 3.2.5 (Euklidnorm) und 3.2.8 (Maximumnorm) schließen. Bei der Quadratsummennorm muß allerdings auf die richtige Skalierung geachtet werden. Denn bei Gittervektoren g_{ij} ist die reine Summe ∑_{i,j} g_{ij}^{2} kein gutes Maß für die Größe von g, da sich wegen der (m − 1)^{2} Summanden bei einer konstanten Funktion der Wert (m − 1)^{2}g^{2} ergibt. Daher werden nun für Gitterfunktionen g_{Δ} folgende skalierten Normen definiert,

∥g_{Δ}∥_{Δ,2} = (h^{2} ∑_{(x,y)∈Ω_{h}} g_{Δ}^{2}(x, y))^{1/2}, ∥g_{Δ}∥_{Δ,∞} = max_{(x,y)∈Ω_{h}} |g_{Δ}(x, y)|.  (3.2.22)

Beide Normen sind Approximationen für die jeweiligen Funktionen-Normen,

∥g∥_{2} = ( ∫_{0}^{1} ∫_{0}^{1} g(x, y)^{2} dxdy)^{1/2}, ∥g∥_{∞} = sup_{(x,y)∈Ω} |g(x, y)|,

insbesondere gilt ∥g_{Δ}∥_{Δ,2} ≤ ∥g_{Δ}∥_{Δ,∞}. Durch Kombination der Ergebnisse aus Formel (3.2.21) und den Sätzen 3.2.5, 3.2.8, 3.2.9, folgt die globale Fehleraussage

Satz 3.2.10 Die Lösung u des Dirichletproblems auf dem Einheitsquadrat erfülle u ∈ C^{4}(Ω̄). Dann gilt für den Fehler der Näherungslösung u_{Δ} aus dem Differenzenverfahren (3.2.9) mit Schrittweite h in jeder der Normen aus (3.2.22) die Abschätzung

∥u_{Δ} − u∥_{Δ,p} ≤ C_{p}h^{2}(∥u_{xxxx}∥_{∞} + ∥u_{yyyy}∥_{∞}),

p ∈ \{2, ∞\}, mit C_{2} = 1/192 (h ≤ \frac{1}{2}), C_{∞} = 1/96.

Beweis Da die Werte von u und u_{Δ} auf dem Rand Γ_{h} übereinstimmen, treten in der Fehlergleichung (3.2.21) nur Beiträge aus dem Konsistenzfehler T_{Δ} auf, der mit Satz 3.2.9 (Konstante C = 1/12) abgeschätzt werden kann. Daher sind die Sätze 3.2.5 und 3.2.8 anwendbar (homogene Randbedingungen) und ergeben die Behauptung. Für das gemischte Randwertproblem (3.2.12) erhält man nach dem selben Prinzip eine ähnliche Aussage mit etwas größeren Konstanten.

Nur auf den ersten Blick ist diese Fehleraussage zufriedenstellend. Denn in wichtigen Fällen erfüllt die Lösung u nicht die Voraussetzungen des Satzes, da dort von g ∈ C^{2}(Ω̄) nur auf u ∈ C^{4}(Ω), (Ω offen!) geschlossen werden kann (vgl. Beispiel 3.1). In diesem Beispiel wächst der lokale Konsistenzfehler (3.2.20) zum Rand hin stark an. Bei den gewöhnlichen Differentialgleichungen wurden solche Effekte durch Verkleinerung der lokalen Schrittweite kompensiert. Überträgt man dieses Prinzip auf ein Produkt-Gitter (x_{i}, y_{j}), indem man in bestimmten Teilen des x- und y-Intervalls die Punkte x_{i}, y_{j} dichter wählt, kann zwar an einer beliebigen Stelle im Gebiet ein feineres zweidimensionales Gitter erzeugt werden. Allerdings ergeben sich dann auch außerhalb dieses Bereichs Regionen mit kleinen Schrittweiten, die in xund y-Richtung außerdem sehr unterschiedlich sind.

Auf einfache Weise läßt sich eine lokale Verfeinerung in einem kleinen Bereich dadurch erreichen, daß man ein Gitter mit halber Schrittweite h/2 in das grobe Gitter ’einhängt’. An einer Übergangsstelle kann das folgendermaßen skizziert werden:

Im groben Gitter (ausgefüllte Punkte) wird der normale Fünfpunktestern mit Schrittweite h verwendet, der an den Übergangsstellen den ersten Punkt des feinen Gitters überspringt. In den Punkten des feinen Gitters (offene Kreise) wird der Stern mit Schrittweite h/2 eingesetzt. In den beiden Punkten des feinen Gitters, die keinen westlichen Nachbarn haben, kann ein diagonaler

[figure]

Fünf-Punkte-Stern mit Schrittweite h/\sqrt{2} benutzt werden (vgl. Lemma 3.2.12). Alle diese Sterne haben die Approximationsordnung 2. Das Verfahren führt nicht mehr auf eine symmetrische Matrix, allerdings läßt sich die Inverse in der ∞-Norm wieder mit Lemma 3.2.7 und dem Testvektor aus Satz 3.2.8 abschätzen. Die lokale Gitterverfeinerung kann in der Nähe einer Problemstelle natürlich mehrfach angewendet werden.

Beispiel 3.2.11 Auf dem L-Gebiet Ω = (0, 1) × (0, 1) \ [\frac{1}{2}, 1) × [\frac{1}{2}, 1) gelte −u_{xx} − u_{yy} = 25, auf dem Rand u = 0. Einsatz ineinandergeschachtelte Gitter mit h = 1/8, 1/16, 1/32. Der steile Anstieg in der einspringenden Ecke wird erst durch die lokal feinere Unterteilung sichtbar.

Wenn die (höheren) Ableitungen der Lösung zu groß werden, kann dies also durch Gitterverfeinerung kompensiert werden. Umgekehrt ist es natürlich sinnvoll, den Aufwand durch Verfahren höherer Konvergenzordnung zu senken, wenn die höheren Ableitungen der Funktion gutartig sind. Die Kenntnis der genauen Fehlerstruktur des oben verwendeten diagonalen Sterns kann dazu verwendet werden.

Lemma 3.2.12 Für eine Funktion u ∈ C^{6}(Ω̄) gilt

\frac{1}{2h^{2}} (4u(x, y) − u(x + h, y + h) − u(x + h, y − h) − u(x − h, y + h) − u(x − h, y − h))
= −(u_{xx} + u_{yy} + \frac{h^{2}}{12} (u_{xxxx} + 6u_{xxyy} + u_{yyyy}))|_{(x,y)} + O(h^{4}),

Beweis Übungsaufgabe. Zum Vergleich lautet der Fehler des normalen Fünf-Punkte-Sterns (Satz 3.2.9)

\frac{h^{2}}{12}(u_{xxxx}(x, y) + u_{yyyy}(x, y)) + O(h^{4}).

Bei der Poissongleichung kann nun durch eine geeignete Linearkombination dieser beiden Sterne ein h^{2}-Fehlerterm erzeugt werden, der sich mit Hilfe der Dgl wieder approximieren läßt. Wendet man nämlich den Laplace-Operator noch einmal auf die Poissongleichung an, folgt

−g_{xx} − g_{yy} = u_{xxxx} + 2u_{xxyy} + u_{yyyy}.

Diejenige Linearkombination der beiden Fünfpunktesterne, deren h^{2}-Fehleranteil den Ausdruck auf der rechten Seite enthält, entspricht dem folgenden Neun-Punkte-Stern

[figure: \frac{1}{6h^{2}} ·; -1  -1; ^{-1  -1}; -4; -4; -4; -4]

Das Gewicht im Mittelpunkt ist 20/(6h^{2}). Der Ausdruck \frac{h^{2}}{12} (g_{xx}+g_{yy}) im führenden Fehlerbeitrag wird natürlich wieder durch den Fünfpunktestern approximiert. Dies führt auf das folgende Verfahren. Eine einzelne Gleichung desselben im Punkt (x, y) wird der Übersicht halber durch Angabe der Gewichte in den Sternen abgekürzt:

\frac{1}{6h^{2}} [ \begin{matrix}−1 & −4 & −1 \\ −4 & 20 & −4 \\ −1 & −4 & −1 \end{matrix} ] u_{Δ}(x, y) = \frac{1}{12} [ \begin{matrix}& 1 \\ 1 & 8 & 1 \\ & 1 \end{matrix} ] g(x, y), (x, y) ∈ Ω_{h}.  (3.2.23)

Dieses Verfahren wird Mehrstellenverfahren genannt, sein lokaler Fehler hat die Form O(h^{4}) für u ∈ C^{6}, bei der Potentialgleichung mit g ≡ 0 besitzt das Verfahren sogar die Ordnung 6. Die Stabilität beim Dirichletproblem kann wieder mit Lemma 3.2.7 wie in Satz 3.2.8 gezeigt werden mit dem gleichen Ergebnis, ∥L_{Δ}^{−1}∥_{∞} ≤ \frac{1}{8}. Der wesentliche Vorteil des Verfahrens ist die Kompaktheit seines Differenzensterns, der bei Anwendung nur geringe Modifikationen gegenüber dem Fünfpunktestern erfordert (Randbedingungen, Matrixstruktur). Da er aber auf der speziellen Form der Poissongleichung beruht, kann er nur in diesem Spezialfall verwendet werden.

Eine auch bei allgemeineren Dgln verwendbare Methode zur Erhöhung der Konsistenzordnung ist der Einsatz von Differenzenausdrücken mit mehr als drei Punkten in jeder Ortsrichtung.

In (1.7.7) wurde eine symmetrische Approximation für u_{xx} mit fünf Funktionswerten angegeben. Die Kombination mit der entsprechenden Approximation für u_{yy} führt bei der Poissongleichung auf den Neun-Punkte-Stern

[figure: \frac{1}{12h^{2}} ·; 60; 1; 1; 1; 1; -16; -16; -16; -16]

mit einem O(h^{4})-Fehler. Da hier jetzt Werte verknüpft werden mit einem maximalen Abstand von 4 Gitterpunkten, ergeben sich erhebliche Probleme am Rand. Insbesondere bei nicht achsenparallelen Ecken können diese Schwierigkeiten auch nicht immer durch Verwendung unsymmetrischer Approximationen behoben werden. Dieses Verfahren wird nicht weiter betrachtet.

3.2.4 Allgemeinere Gebiete und Gleichungen

Der Fünfpunktestern kann in allen Gebieten verwendet werden, dessen Randpunkte auf ein geeignetes quadratisches oder rechteckiges (h_{x} ≠ h_{y}) Gitter Δ_{h} fallen. Für nichtpolygonale Ränder (’krummlinige’) läßt es sich aber meist nicht mehr erreichen, daß alle Nachbarn eines inneren Gitterpunkts entweder wieder Gitter- oder Randpunkte sind. Daher ist nun der Fall zu betrachten, daß ein oder mehrere Nachbarn eines (randnahen) Gitterpunkts P außerhalb von Ω liegen wie in der gezeigten Skizze.

[figure: O’; N’; P; W; S; N; O; P; O; O’; x; x + h; x + d_{O}h]

Der Rand Γ verlaufe also zwischen P und O durch den Punkt O′ und (oder) zwischen P und N durch N′. Als allgemeinen Zugang kann man hier durch Taylor-Abgleich die Gewichte eines Fünfpunktesterns zu den Punkten P, N′, O′, S, W als Approximation an −(u_{xx}+u_{yy}) bestimmen. Auf das gleiche Ergebnis kommt man, wenn Gitterwerte außerhalb von Ω, also v_{N}, v_{O} durch lineare Extrapolation des inneren Wertes v_{P} und der bekannten Randwerte ϕ_{d}(N′), ϕ_{d}(O′) (beim Dirichletproblem) approximiert werden. Wenn d_{N}h, d_{O}h, 0 ≤ d_{N}, d_{O} ≤ 1 die Abstände der Randpunkte von P sind, dann ist

ϕ_{d}(O′) = d_{O}v_{O} + (1 − d_{O})v_{P} ⇒ v_{O} = [ϕ_{d}(O′) − (1 − d_{O})v_{P}]/d_{O}.

Durch Einsetzen in die zweite Differenz in x-Richtung, Anwendung auf u und Taylorentwicklung führt dies auf

−u(x − h, y) + (1 + \frac{1}{d_{O}})u(x, y) − \frac{1}{d_{O}} u(x + d_{O}h, y) =

− u(W) + (1 + \frac{1}{d_{O}})u(P) − \frac{1}{d_{O}} u(O′) = −\frac{h^{2}}{2} (1 + d_{O})u_{xx} + \frac{h^{3}}{6} (1 − d_{O}^{2})u_{xxx} + … (3.2.24)

Division durch h^{2}(1 + d_{O})/2 liefert eine Approximation der zweiten Ableitung. Analoges gilt in y-Richtung. Im oben skizzierten Beispiel ergibt sich so im Punkt P die Gleichung

\frac{1}{h^{2}} [(\frac{2}{d_{O}} + \frac{2}{d_{N}})v_{P} − \frac{2}{1 + d_{O}} v_{W} − \frac{2}{d_{O}(1 + d_{O})}v_{O′} − \frac{2}{1 + d_{N}} v_{S} − \frac{2}{d_{N}(1 + d_{N})}v_{N′}] = g_{P} . (3.2.25)

Da die Gewichte bei den Gitternachbarn v_{S}, v_{W} von den d-Werten abhängen, führt diese Approximation nicht mehr auf eine symmetrische, aber immer noch auf eine M-Matrix, da die Gleichung (3.2.25) diagonaldominant ist. Die Normabschätzung in der ∞-Norm aus §3.2.2 kann hier wiederholt werden, wobei die Konstanten vom Gebiet abängen.

Bei der Diskussion des Gesamtfehlers ist allerdings eine ausführlichere Argumentation notwendig, da die (Konsistenz-) Ordnung des jetzt asymmetrischen Sterns (3.2.25) nur noch eins ist, wie man an der Taylorentwicklung (3.2.24) erkennt. Dennoch kann weiter O(h^{2})-Konvergenz des Gesamtverfahrens gezeigt werden, da der größere lokale Fehler nur am Rand auftritt. Hierbei spielt die M-Eigenschaft der Matrix wieder die entscheidende Rolle, da mit ihr Vergleiche zwischen Vektoren übertragen werden können. Es bezeichne L_{Δ} wieder die lineare Abbildung, die sich beim Dirichletproblem bei Anwendung des normalen Fünf-Punkte-Sterns in inneren Gitterpunkten und einer Gleichung der Form (3.2.25) in randnahen Punkten ergibt. Bei Anwendung dieser Abbildung auf die konstante Funktion e_{Δ}(x, y) ≡ 1 liefert das Bild L_{Δ}e_{Δ} (die Zeilensumme der Matrix) nur in den randnahen Punkten einen nichtverschwindenden Beitrag, da die Dirichletrandwerte im Gleichungssystem L_{Δ}u_{Δ} = g_{Δ} auf der rechten Seite berücksichtigt werden (vgl.(3.2.10)). Speziell in (3.2.25) ergibt dies

(L_{Δ}e_{Δ})(P) = \frac{1}{h^{2}} [\frac{2}{d_{O}} + \frac{2}{d_{N}} − \frac{2}{1 + d_{O}} − \frac{2}{1 + d_{N}}] = \frac{2}{h^{2}} [\frac{1}{d_{O}(1 + d_{O})} + \frac{1}{d_{N}(1 + d_{N})}] > 0. (3.2.26)

Bei einer anderen Anzahl von Randpunkten im Stern von P ist dieser Ausdruck zu modifizieren. Der Konsistenzfehler des Verfahrens hat hier die Form T_{Δ} = L_{Δ}u − g_{Δ} = hr_{Δ} + h^{2}t_{Δ}, wobei der h^{2}-Anteil dem des normalen Sterns aus Satz 3.2.9 entspricht, während hr_{Δ} nur Randfehler enthält. Im Randpunkt P liefert Gleichung (3.2.25) nach (3.2.24) hierzu den Beitrag

hr_{Δ}(P) = \frac{h}{3}[(1 − d_{O})u_{xxx}(x + ξh, y) + (1 − d_{N})u_{yyy}(x, y + ηh)], ξ, η ∈ (0, 1).

Mit K_{j} := ∥∂^{j}u/∂x^{j}∥_{∞} + ∥∂^{j}u/∂y^{j}∥_{∞}, j = 3, 4, kann dieser Beitrag mit (3.2.26) verglichen werden, dabei gilt

h|r_{Δ}(P)| ≤ K_{3} \frac{h}{3}[1 − d_{O} + 1 − d_{N}] ≤ K_{3} \frac{h}{3}[\frac{1}{1 + d_{O}} + \frac{1}{1 + d_{N}}] ≤ K_{3} \frac{h^{3}}{6} (L_{Δ}e_{Δ})(P).

Aus der üblichen Fehlergleichung L_{Δ}(u_{Δ} − u) = −T_{Δ} = −(hr_{Δ} + h^{2}t_{Δ}) kann man damit wegen L_{Δ}^{−1} ≥ 0 komponentenweise Abschätzungen herleiten. In jedem Gitterpunkt gilt

|u_{Δ} − u| = |L_{Δ}^{−1}(hr_{Δ} + h^{2}t_{Δ})| ≤ L_{Δ}^{−1}(h|r_{Δ}| + h^{2}|t_{Δ}|) ≤ L_{Δ}^{−1}(K_{3} \frac{h^{3}}{6} L_{Δ}e_{Δ} + h^{2}|t_{Δ}|)
≤ (K_{3} \frac{h^{3}}{6} + K_{4} \frac{h^{2}}{12}∥L_{Δ}^{−1}∥_{∞})e_{Δ},

also ∥u_{Δ} − u∥_{∞} ≤ Ch^{2} für den globalen Fehler.

Beispiel 3.2.13 Auf dem Viertelkreis Ω = {(x, y) : x^{2}+y^{2} < 1, x > 0, y > 0} sei −u_{xx}−u_{yy} = 4 sowie u = 0 auf dem Kreis und u_{x}(0, y) = 0, u_{y}(x, 0) = 0. Da die Lösung u(x, y) = 1 − x^{2} − y^{2} ist, arbeitet das Verfahren exakt.

Bei einer allgemeinen linearen, elliptischen Dgl 2.Ordnung

au_{xx} + 2bu_{xy} + du_{yy} + eu_{x} + fu_{y} + qu = g  (3.2.27)

(ad − b^{2} > 0) ist gegenüber den schon aus §1.7 und §3.2.1 bekannten Approximationen nur noch die der gemischten Ableitung u_{xy} nachzutragen. Diese kann aus dem Produkt der ersten Differenzen in x- und y-Richtung zusammengesetzt werden. Mit Schrittweite t hat dieser einfache Stern für u_{xy} die Form

\frac{1}{t^{2}} ∗ [figure: 1; 1; −1; −1] u(x, y):= \frac{1}{t^{2}} [u(x + t, y + t) + u(x − t, y − t) − u(x + t, y − t) − u(x − t, y + t)]

Er ist mit Zentrum (x, y) und Schrittweite t = 2h eine O(h^{2})-Approximation an u_{xy}(x, y). Beim Einsatz dieses Sterns besitzen dann aber die Nebendiagonalelemente der entstehenden Matrix (die Gewichte der Sternstrahlen) nicht mehr alle das gleiche Vorzeichen, was für die M-Matrix-Eigenschaft entscheidend war. Dies kann jedoch durch die Verwendung von zwei anderen Sternen in Abhängigkeit vom Vorzeichen von b erreicht werden. Dazu setzt man jeweils zwei der vier Sterne mit Achsenlänge h, die um (x ± h/2, y ± h/2) zentriert sind, durch Mittelung zu einem um (x, y) zentrierten zusammen, der einen Fehler O(h^{2}) besitzt:

[figure: \frac{1}{2h^{2}} ∗; 1; 1; −1; −1; 1; 1; −1; −1; oder; \frac{1}{2h^{2}} ∗; 1; 1; −1; −1; 1; 1; −1; −1]

Von diesen Sternen wird der linke verwendet in Punkten, wo b > 0, der rechte, wenn b < 0 ist:

[figure: u_{xy} ≅ \frac{1}{2h^{2}} ·; 1; -1; -1; 1; 2; -1; -1; · u,; u_{xy} ≅ \frac{1}{2h^{2}} ·; 1; -1; -1; 1; -2; 1; 1; · u.]

Gemeinsam führen diese Approximationen für die allgemeine Dgl (3.2.27) in jedem Gitterpunkt (x, y) auf folgende Differenzengleichung, die wieder in Stern/Matrixform geschrieben wird. Alle Koeffizienten sind ebenfalls in (x, y) auszuwerten. Dabei bedeutet b^{+} := max{b, 0} ≥ 0, b^{−} := min{b, 0} ≤ 0.

\frac{1}{h^{2}} [ \begin{matrix}−b^{−} & d − |b| & b^{+} \\ a − |b| & −2(a + d − |b|) & a − |b| \\ b^{+} & d − |b| & −b^{−} \end{matrix} ] u + \frac{1}{2h} [ \begin{matrix}& f \\ −e & & e \\ & −f \end{matrix} ] u + qu = g.  (3.2.28)

Insbesondere reicht also bei allgemeineren Gleichungen ein Fünfpunktestern nicht mehr aus. Unter einigen Einschränkungen an die Koeffizienten, die einerseits eine Verschärfung der Elliptizität der Dgl darstellen (vgl. Def. 3.1.5) bzw. der Schrittweitenbeschränkung aus Satz 1.7.2 entsprechen, läßt sich auch hier wieder zeigen, daß der Stern (3.2.28) eine M-Matrix erzeugt.

Satz 3.2.14 Für die Koeffizienten der Dgl (3.2.27) gelte in jedem Gitterpunkt die Einschränkung

|b| + \frac{h}{2}|e| < a, |b| + \frac{h}{2}|f| < d, q ≥ 0.

Für Gebiete, bei denen alle relevanten Randpunkte auf dem Gitter liegen, ist die negative Koeffizientenmatrix zu (3.2.28) beim Dirichletproblem eine M-Matrix.

Bei allgemeineren Gebieten oder Randbedingungen sind außerdem in jedem randnahen Punkt spezielle Randsterne analog zu (3.2.25) zu erstellen. Das im nächsten Paragraphen besprochene Verfahren ist im Vergleich dazu wesentlich flexibler in Bezug auf die Geometrie des Problems, allerdings nicht in Bezug auf die Gestalt der Dgl (Ableitungen, Nichtlinearität).

3.3 Finite-Elemente-Verfahren für elliptische Probleme

3.3.1 Variationsformulierung

Ein fundamentales physikalisches Prinzip besagt, daß jedes physikalische System den Zustand einnimmt, in dem seine Gesamtenergie minimal ist. Dieser Zustand kann dann auch durch eine Differentialgleichung beschrieben werden. Im eindimensionalen sei u ∈ C^{2}[0, 1] die Lösung des folgenden Randwertproblems

− (p(x)u′(x))′ + q(x)u(x) = g(x), x ∈ (0, 1), u(0) = u(1) = 0,  (3.3.1)

mit Koeffizienten-Funktionen p ∈ C^{1}[0, 1], q ∈ C[0, 1], p(x) ≥ p^{∗} > 0. Die Randbedingungen wurden zur Vereinfachung homogen gewählt, da dann die Räume

C_{0}^{k}[0, 1] := {v ∈ C^{k}[0, 1] : v(0) = v(1) = 0}, k = 0, 1, 2,

linear sind. In (3.3.1) ist also u ∈ C_{0}^{2}[0, 1] gesucht. Wenn auf beiden Seiten der Dgl in (3.3.1) das Innenprodukt mit einem beliebigen anderen Element v ∈ C_{0}^{2}[0, 1] gebildet wird, folgt durch partielle Integration

(g, v)_{2} := ∫_{0}^{1} g(x)v(x)dx = ∫_{0}^{1}[−(pu′)′v + quv]dx = \underbrace{−[pu′v]_{0}^{1}}_{=0} + ∫_{0}^{1}[pu′v′ + quv]dx.

Der Ausdruck auf der rechten Seite ist schon für einmal differenzierbare Funktionen erklärt, damit läßt sich folgende symmetrische Bilinearform a : C_{0}^{1}[0, 1] × C_{0}^{1}[0, 1] → R definieren

a(u, v) := ∫_{0}^{1}[p(x)u′(x)v′(x) + q(x)u(x)v(x)] dx, u, v ∈ C_{0}^{1}[0, 1].  (3.3.2)

Das Innenprodukt (g, v)_{2} ist ein beschränktes lineares Funktional, das mit f bezeichnet wird, f : C[0, 1] → R, v ↦ (g, v)_{2}. Aus der obigen Umformung folgt also, daß die Lösung u des RWPs (3.3.1) auch folgende Gleichung erfüllt

a(u, v) = f(v) ∀v ∈ C_{0}^{1}[0, 1].  (3.3.3)

Für p > 0, q ≥ 0 ist die Bilinearform a(u, v) definit, stellt auf C_{0}^{1}[0, 1] also ein Innenprodukt dar und \sqrt{a(v, v)} eine Norm. Aufgrund dieser Tatsache ist die Lösung u auch Minimalstelle des konvexen Funktionals

J(v) := \frac{1}{2}a(v, v) − f(v), d.h., J(u) = inf\{J(v) : v ∈ C_{0}^{1}[0, 1]\}.  (3.3.4)

Mit beliebigem 0 ≠ v ∈ C_{0}^{1}[0, 1], 0 ≠ t ∈ R, gilt für das Element u + tv nämlich

J(u + tv) = \frac{1}{2}a(u + tv, u + tv) − f(u + tv)  (3.3.5)
= \frac{1}{2}a(u, u) − f(u) + t[a(u, v) − f(v)] + \frac{t^{2}}{2} a(v, v) = J(u) + \frac{t^{2}}{2} a(v, v)
> J(u).

Umgekehrt wäre für a(u, v)−f(v) ≠ 0 auch u keine Minimalstelle. Daher stimmen die Lösungen des RWPs (3.3.1), von (3.3.3) und (3.3.4) überein, wenn sie in C_{0}^{2} liegen. Es muß aber betont werden, daß bei weniger glatten Koeffizienten p, q, g das Minimum bei (3.3.4) in einer Funktion u ∉ C^{2}[0, 1] angenommen werden kann. Daher heißen Funktionen, die (3.3.3) oder (3.3.4) erfüllen, verallgemeinerte oder schwache Lösungen. Für eine korrekte Formulierung müssen dann andere Räume zugrundegelegt werden, vgl. (3.3.10).

Auch die Lösung eines partiellen, elliptischen Randwertproblems ist Minimalstelle eines konvexen Funktionals. Es sei wieder Ω ⊆ R^{2} ein beschränktes Gebiet mit einem Rand Γ aus stückweise stetig differenzierbaren Kurven. Nach der Greenschen Formel gilt hier die Beziehung

∫ ∫_{Ω}[u_{x}v_{x} + u_{y}v_{y}] dx dy = − ∫ ∫_{Ω}[u_{xx} + u_{yy}]v dx dy + ∮_{Γ} \frac{∂u}{∂n}v ds  (3.3.6)

für u ∈ C^{2}(Ω̄) und v ∈ C^{1}(Ω̄). Diese Relation spielt wieder eine Rolle, wenn man Störungen u + tv einer Lösung u betrachtet. Speziell beim Dirichletproblem müssen u und u + tv die gleichen Randwerte besitzen, also v auf dem Rand verschwinden. In diesem Fall entfällt daher das Randintegral in (3.3.6). Zur Vereinfachung wird wieder das Dirichletproblem bei der (etwas erweitereten) Poissongleichung mit q ≥ 0 diskutiert,

− u_{xx} − u_{yy} + qu = g in Ω, u = 0 auf Γ.  (3.3.7)

Zu diesem Randwertproblem kann schon für Funktionen aus dem Raum C_{0}^{1}(Ω), wobei

C_{0}^{k}(Ω) := {u ∈ C^{k}(Ω) : u|_{Γ} ≡ 0}, k = 0, 1,

die Bilinearform a : C_{0}^{1}(Ω) × C_{0}^{1}(Ω) → R und die Linearform f : C_{0}^{1}(Ω) → R gebildet werden,

\begin{aligned}a(u, v) & := & ∫ ∫_{Ω}[u_{x}v_{x} + u_{y}v_{y} + quv] dx dy, \\ f(v) & := & (g, v)_{2} = ∫ ∫_{Ω} uv dx dy, \end{aligned} u, v ∈ C_{0}^{1}(Ω).  (3.3.8)

Aufgrund der Greenschen Formel gilt also wieder:

u ∈ C_{0}^{2}(Ω) löst (3.3.7) ⇒ a(u, v) = f(v) ∀v ∈ C_{0}^{1}(Ω).

Auch Neumann- oder Cauchy-Randbedingungen können in einer ähnlichen Formulierung berücksichtigt werden. Dazu werden die Restriktionen an u in der Definition des Raums C_{0}^{1} nur auf Γ_{d} vorgenommen, für Γ_{n}, Γ_{c} werden dafür Randintegrale mit uv, u in die Funktionale a bzw. f aufgenommen. Im folgenden wird aber weiter nur das Dirichletproblem behandelt. Es läßt sich zeigen, daß

∥v∥_{1,2} := ( ∫ ∫_{Ω}[v_{x}^{2} + v_{y}^{2}] dx dy)^{1/2}  (3.3.9)

auf C_{0}^{1}(Ω) eine Norm darstellt. Für q ≥ 0 ist daher auch a(v, v) wieder definit. Dies war die entscheidende Eigenschaft in der Entwicklung (3.3.5), die zeigt, daß auch die Lösung u von (3.3.7) als Mimimalstelle (3.3.4) von J(v) = \frac{1}{2}a(v, v) − f(v) charakterisiert werden kann mit den Bi-Linearformen aus (3.3.8). Die Tatsache, daß dabei nur einmalige Differenzierbarkeit von u erforderlich ist, ist nach den Bemerkung zu Beginn über die Herkunft der Dgl zweiter Ordnung aus einem System 1. Ordnung, nicht so überraschend.

In diesem Minimalproblem dürfen im Prinzip alle Funktionen v betrachtet werden, für die das Funktional J sinnvoll definiert ist, also auch solche, bei der die ersten Ableitungen nur quadrat-integrierbar sind. Dazu gehören bei (3.3.8) die Elemente des Sobolev-Raums

W_{0}^{1}(Ω) := {u ∈ L_{2}(Ω) : u_{x}, u_{y} ∈ L_{2}(Ω), u|_{Γ} = 0}.  (3.3.10)

Jede klassische Lösung des RWPs (vgl. Def. 3.2.1) ist auch globale Minimalstelle von J in W_{0}^{1}(Ω). Diese Erweiterung des Grundraums von C_{0}^{1} auf W_{0}^{1} hat große Bedeutung für die Konstruktion von Verfahren, da W_{0}^{1} auch Funktionen enthält, die (nur) stückweise stetig differenzierbar sind.

3.3.2 Rayleigh-Ritz-Galerkin-Verfahren

Mit der im letzten Abschnitt hergeleiteten Minimaleigenschaft der Lösung können in elementarer Weise Näherungsverfahren konstruiert werden dadurch, daß das Minimum des Energie-Funktionals J nicht mehr im Gesamtraum V (= W_{0}^{1}[0, 1] bzw. = W_{0}^{1}(Ω)) gesucht wird, sondern nur noch in einem endlichdimensionalen linearen Unterraum V_{Δ} ⊆ V. Man kann bei V_{Δ} zunächst an die Menge aller Polynome eines bestimmten Maximalgrads denken, flexibler sind aber die im nächsten Abschnitt besprochenen Spline-Funktionen, die nur noch stückweise stetig differenzierbar sind. Es sei nun also

V_{Δ} ⊆ V, dim V_{Δ} = n ∈ N,

und {B_{i}}_{i=1}^{n} eine Basis von V_{Δ}. Die Rayleigh-Ritz-Galerkin-Lösung (RRG-Lösung) u_{Δ} ∈ V_{Δ} des Minimalproblems (3.3.4) ist definiert durch

J(u_{Δ}) = min_{v∈V_{Δ}} J(v) = min_{v∈V_{Δ}} \frac{1}{2}a(v, v) − f(v).  (3.3.11)

Für sie Analyse dieses Verfahrens gibt es ein allgemeines Prinzip. Als wesentliche Voraussetzung wird dabei die Vergleichbarkeit der von der Bilinearform a induzierten Norm und der Standard- (Sobolev-) Normen benötigt. Daher wird jetzt angenommen, daß Konstanten μ_{0}, μ_{1}, M_{1} > 0 existieren so, daß

μ_{0}∥v∥_{2}^{2} ≤ μ_{1}∥v∥_{1,2}^{2} ≤ a(v, v) ≤ M_{1}∥v∥_{1,2}^{2} ∀v ∈ V.  (3.3.12)

Diese Voraussetzung ist in den betrachteten Beispielen (3.3.2) bzw. (3.3.8) erfüllt, wenn die Koeffizienten beschränkt sind und p ≥ p^{∗} > 0, q ≥ 0 gilt.

Satz 3.3.1 Unter der Voraussetzung (3.3.12) existiert eine eindeutige RRG-Lösung u_{Δ} von (3.3.11). Der Vektor η = (η_{j})_{j=1}^{n} der Koeffizienten von u_{Δ} = ∑_{j=1}^{n} η_{j}B_{j} löst das lineare System

Aη = γ, A = (a(B_{i}, B_{j}))_{i,j=1}^{n}, γ = (f(B_{i}))_{i=1}^{n}.  (3.3.13)

Die Matrix A = (a_{ij}) ist symmetrisch, positiv definit. Äquivalent zu (3.3.13) ist die Aussage

a(u_{Δ}, v) = f(v) ∀v ∈ V_{Δ}.  (3.3.14)

Bemerkung: Die Bedingung (3.3.14) heißt Galerkin-Bedingung. Sie definiert auch dann brauchbare Lösungen, wenn die Bilinearform a nicht mehr definit ist, eine Minimallösung (3.3.11) also nicht mehr existiert.

Beweis Es sei v = ∑_{j=1}^{n} ξ_{j}B_{j} ∈ V_{Δ}, ξ := (ξ_{i}), beliebig. Dann gilt wegen der Bilinearität von a(., .) und wegen (3.3.12) zunächst

ξ^{T} Aξ = ∑_{i,j=1}^{n} ξ_{i}ξ_{j}a(B_{i}, B_{j}) = a( ∑_{i=1}^{n} ξ_{i}B_{i}, ∑_{j=1}^{n} ξ_{j}B_{j}) = a(v, v) ≥ μ_{0}∥v∥_{2}^{2} > 0,

für v ≠ 0, d.h., ξ ≠ 0. Also ist A definit und somit invertierbar. Daraus folgt

J(v) = \frac{1}{2}a( ∑_{i=1}^{n} ξ_{i}B_{i}, ∑_{j=1}^{n} ξ_{j}B_{j}) − f(∑_{i=1}^{n} ξ_{i}B_{i}) = \frac{1}{2}ξ^{T} Aξ − γ^{T} ξ
= \frac{1}{2}(ξ − A^{−1}γ)^{T} A(ξ − A^{−1}γ) − \frac{1}{2}γ^{T} A^{−1}γ.

Dieser Ausdruck wird minimal genau dann, wenn ξ = A^{−1}γ = η ist, also v = u_{Δ}. Nach Definition entspricht (3.3.13) der Bedingung a(u_{Δ}, B_{i}) = f(B_{i}), i = 1, .., n. Da {B_{i}} eine Basis von V_{Δ} darstellt (n.V.) ergibt sich (3.3.14).

Die Lösung u erfüllt (3.3.3), also f(v) = a(u, v) auch mit v = u. Für beliebige v gilt daher

J(v) − J(u) = \frac{1}{2}a(v, v) − f(v) − \frac{1}{2}a(u, u) + f(u) = \frac{1}{2}a(v, v) − a(u, v) + \frac{1}{2}a(u, u)
= \frac{1}{2}a(v − u, v − u).

Da die Minimalstellen von J(v) − J(u) und J(v) übereinstimmen (auch in V_{Δ}), bekommt man daraus eine andere Charakterisierung der RRG-Lösung, die direkt auf Fehlerschranken führen. u_{Δ} ist nämlich die Bestapproximierende an u im Raum V_{Δ}, gemessen in der ’a-Norm’.

Satz 3.3.2 Für die RRG-Lösung u_{Δ} aus (3.3.11) gilt

μ_{1}∥u_{Δ} − u∥_{1,2}^{2} ≤ a(u_{Δ} − u, u_{Δ} − u) = inf_{v∈V_{Δ}} a(v − u, v − u) ≤ M_{1} inf_{v∈V_{Δ}} ∥v − u∥_{1,2}^{2}.

Durch diesen Satz erhält man sofort eine globale Fehlerschranke in der 1, 2-Norm, also eine Schranke für den Fehler in der Ableitung von u_{Δ}, wenn man rechts eine beliebige, geeignete Funktion v ∈ V_{Δ} einsetzt, z.B., eine Interpolierende der Lösung u. Schranken für den Abstand dist_{1,2}(u, V_{Δ}) = inf_{v∈V_{Δ}} ∥v − u∥_{1,2} werden später für konkrete Räume V_{Δ} hergeleitet.

3.3.3 Spline-Räume, Finite Elemente

Die Effizienz und Flexibilität des RRG-Verfahrens hängt wesentlich von der Wahl des endlichdimensionalen Raums V_{Δ} ab. Man sieht schnell, daß ein Raum von Polynomen bei partiellen Problemen sehr unhandlich wird. Gerade bei der Behandlung unregelmäßiger Gebiete ist dagegen der folgende Zugang sehr anpassungsfähig, bei dem das Grundgebiet in kleine, einfache Elemente zerlegt wird. Hierbei ist es sehr hilfreich, daß das Minimalproblem (3.3.4) schon für stückweise differenzierbare Funktionen erklärt ist.

Zunächst wird das gewöhnliche Randwertproblem (3.3.1) auf [0, 1] betrachtet. Eine Unterteilung des Intervalls führt auf ein Gitter Δ : 0 = x_{0} < x_{1} < .. < x_{m} = 1, mit Schrittweiten h_{i} = x_{i+1} − x_{i}. Die einfachste Funktionenklasse in W_{0}^{1}[0, 1] sind die stückweise linearen Funktionen (lineare Splines, vgl. Numerik I, §4). Unter Berücksichtigung der Randbedingungen sei

V_{Δ} := {s ∈ C_{0}^{0}[0, 1] : s|_{[x_{i},x_{i+1})} ∈ Π_{1}, i = 0, .., m − 1}  (3.3.15)

der Raum aller stetigen Funktionen mit Randbedingung null, die in jedem Teilintervall Polynome vom Grad 1 (∈ Π_{1}) sind. Eine einfache Basis dieses Raums ist durch die Dach-Funktionen

B_{i}(x) := \{ \begin{matrix}(x − x_{i−1})/h_{i−1}, & x ∈ [x_{i−1}, x_{i}) \\ (x_{i+1} − x)/h_{i}, & x ∈ [x_{i}, x_{i+1}) \\ 0, & sonst \end{matrix} ,  (3.3.16)

i = 1, …, m − 1, gegeben, vgl. die Skizze.

[figure: 0; 1; B_{i}; v_{Δ}]

Diese B_{i} bilden eine Kardinalbasis, da für x = x_{k} gilt u_{Δ}(x_{k}) = ∑ η_{j}B_{j}(x_{k}) = η_{k}. Der wichtigste praktische Vorteil der genannten Basis ist aber ihre Lokalität, die dazu führt, daß nur wenige der Matrixelemente a(B_{i}, B_{j}) in (3.3.13) von Null verschieden sind. Dies reduziert den Rechenaufwand sowohl bei der Aufstellung des Gleichungssystems, als auch bei seiner Auflösung.

Als Beispiel werde die Dgl −u′′ = g mit p ≡ 1, q ≡ 0, also die Bilinearform a(u, v) = ∫_{0}^{1} u′(x)v′(x)dx betrachtet. Für die Basisfunktionen gilt nach (3.3.16)

B′_{i}(x) := \{ \begin{matrix}1/h_{i−1}, & x ∈ [x_{i−1}, x_{i}) \\ −1/h_{i}, & x ∈ [x_{i}, x_{i+1}) \\ 0, & sonst \end{matrix} .

Daher ist a(B_{i}, B_{j}) = 0 für |i − j| > 1. Für j ∈ {i − 1, i, i − 1} gilt

a(B_{i}, B_{j}) = ∫_{x_{i−1}}^{x_{i}} \frac{1}{h_{i−1}} B′_{j}dx − ∫_{x_{i}}^{x_{i+1}} \frac{1}{h_{i}} B′_{j}dx = \{ \begin{matrix}−\frac{1}{h_{i−1}} , & j = i − 1 \\ \frac{1}{h_{i−1}} + \frac{1}{h_{i}} , & j = i \\ −\frac{1}{h_{i}} , & j = i + 1 \end{matrix} .

Das Verfahren arbeitet daher, zumindestens bei der Approximation der zweiten Ableitung, ähnlich wie das Differenzenverfahren, denn es führt auf das Gleichungssystem (η_{i} = u_{Δ}(x_{i}), s.o.)

−\frac{1}{h_{i}} (η_{i+1} − η_{i}) + \frac{1}{h_{i−1}} (η_{i} − η_{i−1}) = ∫_{x_{i−1}}^{x_{i+1}} g(x)B_{i}(x)dx, i = 1, .., m − 1.

Für das RRG-Verfahren im Raum (3.3.15) sollte man daher auch einen Fehler der Größenordnung h^{2} erwarten. Dieser läßt sich mit Hilfe von Satz 3.3.2 auf den Approximationsfehler inf{∥u − v_{Δ}∥ : v_{Δ} ∈ V_{Δ}} zurückführen. Bei der Interpolierenden in V_{Δ} kann der Fehler unabhängig auf den einzelnen Teilintervallen untersucht werden. Die Interpolierende einer Funktion v ∈ V ist eine Restriktion R_{Δ}v im Sinne von §2.

Lemma 3.3.3 Zu einer Funktion v ∈ C_{0}^{2}[0, 1] sei R_{Δ}v := ∑_{i=1}^{m−1} v(x_{i})B_{i} die Interpolierende aus V_{Δ}. Dann gilt mit H := max_{i} h_{i} die Fehlerschranke

∥(R_{Δ}v − v)^{(k)}∥_{2} ≤ (\frac{H}{π} )^{2−k}∥v′′∥_{2}, k = 0, 1.

Auf diese Aussage kommt man direkt über die Standarddarstellung des Interpolationsfehlers für Polynome (Numerik I, 3.3.11), zur Verifikation der angegebenen Konstanten ist allerdings eine diffizilere Argumentation nötig. Wird das Infimum in Satz 3.3.2 durch Einsetzen von v = R_{Δ}u abgeschätzt, ergibt sich mit dem letzten Lemma die Fehlerschranke,

∥(u_{Δ} − u)′∥_{2} ≤ \sqrt{\frac{M_{1}}{μ_{1}}} \frac{H}{π} ∥u′′∥_{2}.

Mit einem zusätzlichen Argument (’Nitsche-Trick’) kann aus dieser O(H)-Schranke für die Ableitung bei vielen Randwertproblemen die globale Fehlerschranke

∥u_{Δ} − u∥_{2} ≤ CH^{2}∥u′′∥_{2} ≤ C̃H^{2}∥g∥_{2}

für die Funktionswerte hergeleitet werden. Im Vergleich zu den entsprechenden Aussagen bei Differenzenverfahren ist dabei hervorzuheben, daß hier nur die L_{2}-Integrierbarkeit von u′′ ausgenutzt wurde im Vergleich zur Voraussetzung u ∈ C^{4} bei den Differenzenverfahren. RRG-Verfahren höherer Ordnung erhält man durch Verwendung von Splines höherer Ordnung, wobei die Funktionen in V_{Δ} stückweise aus Polynomen höheren Grads zusammengesetzt sind.

Von größter praktischer Bedeutung ist die Möglichkeit, beim RRG-Verfahren mit stückweise linearen Ansatzfunktionen arbeiten zu können, bei den partiellen Problemen. Sehr einfach sind Zerlegungen bzw. Approximationen des Gebiets Ω durch Dreieckgitter. Daher ist das RRG-Verfahren wesentlich flexibler als Differenzenverfahren bei der Approximation komplizierterer Gebiete und der lokalen Gitterverfeinerung. Beispiele für regelmäßige bzw. lokal verfeinernde Triangulierungen beim L-Gebiet sind:

[figure]

Dabei sind nur Triangulierungen aus offenen, paarweise disjunkten Dreiecken T_{i} zulässig,

Δ = {T_{i}}_{i=1}^{m}, Ω̄ = ⋃_{i=1}^{m} T̄_{i},

bei denen alle Ecken eines Dreiecks auf Ecken seiner Nachbardreiecke oder den Rand fallen. Der Raum aller stetigen, stückweise linearen Funktionen auf Δ,

V_{Δ} := {s ∈ C_{0}(Ω) : s|_{T_{i}} linear},  (3.3.17)

kann wieder durch eine Basis von Dachfunktionen aufgespannt werden. Es sei {P_{i}}_{i=1}^{n} ⊆ Ω die Menge aller Eckpunkte von Dreiecken (Knoten) im Innern des Gebiets in einer geeigneten Numerierung. Dann wird durch

B_{i} ∈ V_{Δ}, B_{i}(P_{j}) = δ_{ij}, j = 1, …, n,

i = 1, …, n, eine Basis von V_{Δ} definiert. Ein Beispiel einer solchen Basisfunktion bei sechs in einer Ecke zusammenstoßenden Dreiecken ist

[figure]

Wesentliche Merkmale dieser Basis sind wieder die Kardinalität, u_{Δ} = ∑_{i=1}^{n} u_{Δ}(P_{i})B_{i}, und die Lokalität der Basis, da B_{i} nur auf den Dreiecken T_{j} mit P_{i} ∈ T̄_{j} nicht verschwindet.

Diese Eigenschaft ist gerade bei mehrdimensionalen Problemen von großer Bedeutung, da dann nur wenige der Integrale a(B_{i}, B_{j}) in der Matrix A von Null verschieden sind, und bei den restlichen für i ≠ j nur über zwei Dreiecke integriert werden muß.

a(B_{i}, B_{j}) = ∑_{P_{i}∈T_{k}∧P_{j}∈T_{k}} ∫ ∫_{T_{k}}[(B_{i})_{x}(B_{j})_{x} + (B_{i})_{y}(B_{j})_{y} + qB_{i}B_{j}] dxdy

Bei der praktischen Berechnung der Integrale a_{ij} = a(B_{i}, B_{j}), γ_{i} = f(B_{i}), ist es am einfachsten, alle Dreiecke T_{i} abzuarbeiten und die einzelnen Teilintegrale in a_{ij}, γ_{i} aufzuaddieren. Die Matrix A ist hier wieder dünn besetzt und symmetrisch, pos. definit. Daher können sehr effiziente (Iterations-) Verfahren zur Lösung des Gleichungssystems (3.3.13) eingesetzt werden.

Zur Fehlerabschätzung der zugehörigen RRG-Lösung nach Satz 3.3.2 läßt sich eine Schranke für den Interpolationsfehler wieder lokal herleiten. Dabei ergibt sich eine verwertbare Aussage allerdings nur, wenn die Dreiecke T_{k} nicht zu spitz werden.

Lemma 3.3.4 Zu einer Funktion v ∈ C_{0}^{2}(Ω) sei R_{Δ}v := ∑_{i=1}^{n} v(P_{i})B_{i} die Interpolierende aus V_{Δ}, (3.3.17). In der Dreieckzerlegung des Gebiets sei h die längste Dreieckseite und α der kleinste Innenwinkel eines Dreiecks T_{i}, i = 1, .., m. Dann gilt die Fehlerschranke

∥R_{Δ}v − v∥_{1,2} ≤ C \frac{h}{sin α}∥v∥_{2,2},

mit ∥v∥_{2,2}^{2} := ∑_{j+k=2} ∥\frac{∂^{2}}{∂x^{j}∂y^{k}} v∥_{2}^{2}.

Damit kann man wie im gewöhnlichen Fall in zwei Schritten einen Fehler O(h/ sin α)∥u∥_{2,2} in den ersten Ableitungen der RRG-Näherung u_{Δ} und eine Schranke O(h^{2}/ sin α)∥u∥_{2,2} für den Fehler ihrer Funktionswerte herleiten. Zusätzlich kann unter einfachen Annahmen an die Koeffizienten zumindest für konvexe Gebiete die Abschätzung ∥u∥_{2,2} ≤ C∥g∥_{2} gezeigt werden. Daher hat man auch im partiellen Fall die h^{2}-Konvergenz unter schwächeren Voraussetzungen als bei den Differenzenverfahren.

Allerdings ist bei diesen starken Aussagen zu beachten, daß sie nur für das exakte RRG-Verfahren zutreffen, bei dem alle Integrale a(B_{i}, B_{j}), f(B_{i}) = (g, B_{i})_{2} ohne Fehler bestimmt werden. In der Praxis berechnet man jedoch zumindestens die Integralanteile, in denen allgemeine Koeffizienten (hier q, g) auftreten, mit Hilfe von Quadraturformeln. Um auch für diese approximierten RRG-Verfahren einen O(h^{2})-Fehler garantieren zu können, müssen doch wieder stärkere Annahmen an die Koeffizienten, also auch an die Regularität der Lösung gemacht werden. Verfahren höherer Konvergenzordnung ergeben sich bei Verwendung von stückweise stetigen Polynomen höheren Grades auf den Dreieckzerlegungen.

3.4 Partielle Anfangs-Randwert-Probleme

Die nicht-elliptischen partiellen Differentialgleichungen 2. Ordnung, die parabolischen und hyperbolischen, besitzen eine gemeinsame Eigenschaft. Schreibt man den jeweiligen einfachsten Vertreter (Wärmeleitungs-, Wellengleichung) in der Form

u_{xx} − u_{y} = 0, bzw. u_{xx} − u_{yy} = 0,

so lassen sich damit sinnvolle Problemstellungen formulieren, wenn Randbedingungen nur in x-Richtung vorgegeben werden, in y-Richtung dagegen nur einseitige, also Anfangsbedingungen. Im physikalischen Zusammenhang hat die erste Variable meist die Bedeutung einer Ortskoordinate, während die zweite einen Zeitablauf beschreibt. Daher wird sie im folgenden mit t bezeichnet. Da die Geometrie des Problems sich i.a. nicht mit der Zeit ändert, ist das Grundgebiet meist ein Zylinder, z. B., [0, 1] × [0, ∞). Mit stetigen Koeffizienten p, r, g, p > 0, p_{x} stetig werden also Anfangs-Randwertprobleme betrachtet mit einer der Dgln

u_{t} = [p(x, t)u_{x}]_{x} + q(x, t)u + g(x, t), x ∈ (0, 1), t > 0, (parabolisch) (3.4.1)
u_{tt} = [p(x, t)u_{x}]_{x} + q(x, t)u + g(x, t), x ∈ (0, 1), t > 0, (hyperbolisch) (3.4.2)

und den Rand- und Anfangsbedingungen

\begin{matrix}u(x, 0) = ϕ(x), x ∈ (0, 1), u(0, t) = u(1, t) = 0, t ≥ 0, \\ u_{t}(x, 0) = ψ(x), x ∈ (0, 1) nur bei (3.4.2). \end{matrix}  (3.4.3)

Die Probleme (3.4.1,3.4.2),(3.4.3) verhalten sich also in Ortsrichtung wie (elliptische) Randwertprobleme, in Zeitrichtung wie Anfangswertprobleme erster bzw. zweiter Ordnung. Diese Tatsache kann bei der numerischen Approximation auch ausgenutzt werden durch Semidiskretisierung in der sogenannten Linienmethode. Dabei wird jeweils eines der von früher bekannten Verfahren zur Diskretisierung von Anfangs- oder Randwertproblemen getrennt auf die Zeit- und Ortsrichtung angewendet. Allerdings muß man sich dabei über die Eigenschaften der nach dem ersten Diskretisierungsschritt entstehenden gewöhnlichen Probleme klar werden. Bei der Zeitlinienmethode, z.B., wird das Problem nur in der Ortsvariablen diskretisiert, das AWP in Zeitrichtung (zunächst) beibehalten. Dabei werde zunächst das Differenzenverfahren betrachtet mit einem äquidistanten Gitter Δ = {x_{i} = ih : i = 0, .., n + 1}, h = 1/(n + 1), wie in §1.7. Der Ausdruck [pu_{x}]_{x} wird dabei durch zweifache Differenzenquotienten approximiert,

[p(x, t)u_{x}(x, t)]_{x} = \frac{1}{h^{2}} [p(x_{i+1/2}, t)(u(x_{i+1}, t)−u(x_{i}, t))−p(x_{i−1/2}, t)(u(x_{i}, t)−u(x_{i−1}, t))]+O(h^{2}),

x_{i±1/2} := x_{i} ±h/2. Unter Vernachlässigung des O(h^{2})-Fehlerterms ergibt sich für die Gitterwerte u_{Δ}(x_{i}, t) einer Näherung u_{Δ}, die immer noch Funktionen der Zeit sind, ein System von gewöhnlichen Anfangswertproblemen. Faßt man nämlich diese Funktionswerte als Komponenten einer Vektorfunktion auf, w(t) = (w_{i}(t))_{i=1}^{n}, w_{i}(t) := u_{Δ}(x_{i}, t) (und w_{0} = w_{n+1} ≡ 0), dann führt die obige Differenzenapproximation im parabolischen Fall (3.4.1) auf das AWP erster Ordnung

\begin{aligned}w′(t) & = & f(t, w(t)), t > 0, mit \\ f_{i}(t, y) & = & \frac{1}{h^{2}} [p_{i+1/2}(y_{i+1} − y_{i}) − p_{i−1/2}(y_{i} − y_{i−1})] + q_{i}y_{i} + g_{i}, \\ w_{i}(0) & = & ϕ(x_{i}), i = 1, …, n. \end{aligned}  (3.4.4)

Die Koeffizientenwerte sind, wie früher, q_{i} = q(x_{i}, t), g_{i} = g(x_{i}, t) und p_{i±1/2} = p((i±\frac{1}{2})h, t). Beim hyperbolischen Problem (3.4.2) ergibt sich das entsprechende AWP mit der Dgl zweiter Ordnung w′′ = f(t, w(t)) und der zusätzlichen Anfangsbedingung w′_{i}(0) = ψ(x_{i}). Im betrachteten Beispiel ist die Funktion f affin linear, f(t, y) = A(t)y + g⃗(t). Die Matrix A ist tridiagonal (vgl. (1.7.5)) und sogar symmetrisch mit den Elementen

a_{ii}(t) = −\frac{1}{h^{2}} [p_{i+1/2} + p_{i−1/2}] + q_{i}, a_{i,i+1}(t) = a_{i+1,i}(t) = \frac{1}{h^{2}} p_{i+1/2}.  (3.4.5)

Wie in Lemma 3.2.4 kann für p > 0 bei dieser Matrix die Aussage y^{T} A(t)y ≤ q^{∗}y^{T} y, q^{∗} = max_{i} q_{i}, gezeigt werden. Daher sind die Eigenwerte von A reell und durch q^{∗} nach oben beschränkt, also insbesondere negativ für q^{∗} ≤ 0. Allerdings wird die Norm ∥A(t)∥ ≅ \frac{4}{h^{2}} ∥p∥_{∞} für kleine Schrittweiten h sehr groß. Tatsächlich läßt sich zeigen, daß sehr kleine Eigenwerte von A (mit großem Betrag) λ ≪ 0 auftreten können. Bei (3.4.4) handelt es sich somit um (eines der Standardbeispiele für) ein steifes Anfangswertproblem. Die Diskussion aus §1.5 ist also hier anwendbar. Die wichtigste Erkenntnis in diesem Zusammenhang war, daß zur Integration des AWPs (3.4.4) nur implizite Verfahren sinnvoll sind. Das läßt sich auch anhand der einfachsten Verfahren bei der Wärmeleitungsgleichung nachvollziehen. Bei dieser (p ≡ 1, q = g ≡ 0) wird das einfache Einschrittverfahren (w(t_{k}) ≅ w^{[k]} ∈ R^{n})

\frac{1}{τ} (w^{[j+1]} − w^{[j]}) = (1 − ϑ)Aw^{[j]} + ϑAw^{[j+1]}, j = 0, 1, …,  (3.4.6)

mit Zeitschrittweite τ und einem Parameter ϑ ∈ [0, 1] betrachtet. Die Matrix A (3.4.5) ist konstant. Das Verfahren stimmt für ϑ = 0 mit dem expliziten und für ϑ = 1 mit dem impliziten Eulerverfahren, überein, für ϑ = \frac{1}{2} ist es die implizite Mittelpunktregel. Es führt auf eine Serie von Gleichungssystemen (σ := τ/h^{2}, ϑ′ := 1 − ϑ),

− ϑσw_{i−1}^{[j+1]} + (1 + 2ϑσ)w_{i}^{[j+1]} − ϑσw_{i+1}^{[j+1]} = ϑ′σw_{i−1}^{[j]} + (1 − 2ϑ′σ)w_{i}^{[j]} + ϑ′σw_{i+1}^{[j]}, (3.4.7)

i = 1, …, n. Beim expliziten Euler-Verfahren, ϑ = 0, reduziert sich dies auf

w_{i}^{[j+1]} = σw_{i−1}^{[j]} + (1 − 2σ)w_{i}^{[j]} + σw_{i+1}^{[j]}, i = 1, …, n.  (3.4.8)

Wie im letzten Paragraphen können diese Beziehungen kompakt durch Differenzensterne beschrieben werden. Zu den verschiedenen Parameterwerten ϑ gehören die Sterne

[figure: −σ; 2σ − 1; −σ; 1; ϑ = 0; −σ/2; σ − 1; −σ/2; −σ/2; σ + 1; −σ/2; ϑ = 1/2; −σ; 2σ + 1; −σ; 1; ϑ = 1]

wobei die Zeitachse t nach rechts gezeichnet ist. Für σ ≤ \frac{1}{2} sind alle Gewichte in (3.4.8) nichtnegativ, es findet also eine fortgesetzte Mittelwertbildung in den Zeitschritten statt. Dies ist ein gutes Modell für die durch die Wärmeleitungsgleichung beschriebenen Ausgleichsvorgänge. Allerdings muß dabei die Zeitschrittweite τ aus Stabilitätsgründen inakzeptabel klein gewählt werden, τ = h^{2}σ ≤ h^{2}/2. Für τ > h^{2}/2, d.h. σ > 1/2 wird das explizite Verfahren (3.4.8) nämlich instabil, denn dann tritt in der (scharfen) Abschätzung

∥w^{[j+1]}∥_{∞} = max_{i} |w_{i}^{[j+1]}| ≤ (2σ + |1 − 2σ|) ∥w^{[j]}∥_{∞},

der Vorfaktor 4σ − 1 = 4τ/h^{2} − 1 > 1 auf. Dies paßt zur Diskussion von §1.5. Denn, da für die Norm ∥A∥ = 4/h^{2} gilt und die Stabilitätsabszisse des expliziten Eulerverfahren 2 ist, entspricht dieser Fall der instabilen Situation τ∥A∥ > 2.

Beim impliziten Eulerverfahren (ϑ = 1) dagegen ist das Tridiagonalsystem (3.4.7) in jedem Zeitschritt diagonaldominant (mit einer M-Matrix) und hier kann ohne jede Schrittweitenbeschränkung die Schranke ∥w^{[j+1]}∥_{∞} ≤ ∥w^{[j]}∥_{∞} gezeigt werden. Das allgemeine Verfahren (3.4.6) ist für alle ϑ ∈ [\frac{1}{2}, 1] A-stabil. Bei diesen Verfahren muß sich die Schrittweitenwahl nur nach den Genauigkeitsanforderungen richten. Dabei ist zu beachten, daß im AWP (3.4.4) auch schon ein Diskretisierungsfehler steckt, der Fehler beim Gesamtverfahren (3.4.6) setzt sich daher aus einem Orts- und einem Zeitfehler zusammen, die jeweils getrennt durch Wahl der Ortsschrittweite h und der Zeitschrittweite τ beeinflußt werden können. Für die Mittelpunktregel (ϑ = 1/2) hat dieser globale Fehler die Form C_{x}h^{2} + C_{t}τ^{2}, für alle anderen ϑ ≠ 0 nur C_{x}h^{2} + C_{t}τ. Das Verfahren arbeitet dann am effizientestens, wenn beide Fehler ungefähr gleich groß sind. Bei der Zeitlinienmethode, mit festem Ortsgitter, kann nur der Zeitfehler beeinflußt werden. Die Zeitschrittweite τ ist daher so zu steuern, daß dieser ungefähr die Größe des Ortsfehlers besitzt.

Eine größere Anpassungsfähigkeit ist mit der Orts-Linienmethode möglich. Dabei werden die beiden Diskretisierungsschritte vertauscht. Die erste Diskretisierung erfolgt in Zeitrichtung mit einem impliziten Verfahren (z.B. dem impliziten Euler-Verfahren). Dies führt auf eine Serie von gewöhnlichen RWPen in jeder Zeitstufe. Das folgende Diagramm veranschaulicht die auftretenden Gitter für beide Möglichkeiten:

\overline{}Zeit- Linienmeth.
Orts- ↓ Linienmeth.  AWP-Lösung ↓

\overline{}RWP- Lösung

Wird bei den RWPen das Differenzenverfahren mit einer festen, konstanten Schrittweite h verwendet, ergibt sich wieder das Verfahren (3.4.6) und das im Bild gezeigte Endgitter. In dieser zweiten Interpretation ist es aber auch möglich, bei der Lösung der RWPe in jedem einzelnen Zeitschritt das Ortsgitter neu (adaptiv) zu wählen. Die im Verfahren verwendeten Werte der vorhergehenden Zeitstufe können dabei interpoliert werden. Hiermit läßt sich sowohl die Zeitals auch die Ortsunterteilung leichter an den Lösungsverlauf anpassen, die Gitter sind dann irregulär.

Auch bei den hyperbolischen Problemen (3.4.2) kann, im Prinzip, analog zum parabolischen Fall vorgegangen werden. Zur Diskretisierung der zweiten Zeitableitung u_{tt} bzw. w_{tt} (in der Zeitlinienmethode) müssen dabei allerdings drei Zeitstufen herangezogen werden. Ein explizites Verfahren ergibt sich direkt aus der Standardapproximation von u_{tt} mit Hilfe des zweiten Differenzenquotienten, für die Gitterwerte (wie oben) w_{i}^{[j]} = u_{Δ}(ih, jτ), i = 1, .., n, lautet das Gesamtverfahren

w_{i}^{[j+1]} = s^{2}w_{i−1}^{[j]} + 2(1 − s^{2})w_{i}^{[j]} + s^{2}w_{i+1}^{[j]} − w_{i}^{[j−1]}, i = 1, …, n,  (3.4.9)

mit s = τ/h. Diese Rekursion kann im Zeitschritt j = 0 gestartet werden, wenn die zweite Anfangsbedingung aus (3.4.3) mit der symmetrischen ersten Differenz approximiert wird, w_{i}^{[1]} − w_{i}^{[−1]} = 2τψ(x_{i}), die eine künstliche Zeitebene −τ einbezieht. Durch eine etwas aufwendigere Analyse als bei den parabolischen Problemen läßt sich zeigen, daß das Verfahren ebenfalls nur stabil ist unter der (schwächeren) Schrittweiteneinschränkung s = τ/h ≤ 1. Es können auch einfache implizite Verfahren angegeben werden, die ohne Schrittweitenbeschränkung stabil sind. Der Vergleich mit den expliziten Verfahren fällt dennoch anders aus als bei den parabolischen Problemen aufgrund der unterschiedlichen Eigenschaften der Lösung, die jetzt diskutiert werden.

Bei der parabolischen Gleichung (3.4.1) werden die Lösungen für wachsendes t > 0 immer glatter. Bei der Wärmeleitungsgleichung u_{t} = u_{xx} mit u(0, t) = u(1, t) = 0 besteht die Lösung aus einer Linearkombination von Funktionen der Form

e^{−k^{2}π^{2}t} sin kπx, k ∈ N,

die in t-Richtung umso schneller verschwinden, je stärker sie in x-Richtung variieren. Außerdem hängt der Lösungswert u(x, t) von allen Daten im Bereich [0, 1] × [0, t) ab. Beim hyperbolischen Problem ergibt sich dagegen ein anderes Bild. Die allgemeine Lösung der Wellengleichung kann explizit konstruiert werden. Jede Funktion der Form

u(x, t) = a(x + t) + b(x − t) löst u_{tt} = u_{xx}  (3.4.10)

mit a, b ∈ C^{2}(R). Wenn man die Randbedingungen in x-Richtung zur Vereinfachung außer Acht läßt, hat diese Lösung die Anfangsbedingungen

a(x) + b(x) = u(x, 0) = ϕ(x), a′(x) − b′(x) = u_{t}(x, 0) = ψ(x)

zu erfüllen. Durch Differentiation der ersten Gleichung können a und b bis auf Konstanten bestimmt werden. Daraus folgt die explizite Darstellung des AWPs bei der Wellengleichung,

u(x, t) = \frac{1}{2}(ϕ(x + t) + ϕ(x − t) + ∫_{x−t}^{x+t} ψ(ξ) dξ).  (3.4.11)

Zwei wichtige Eigenschaften können an dieser Lösungsdarstellung abgelesen werden. Erstens zeigt schon (3.4.10), daß die Eigenschaften der Lösung vom Startzeitpunkt praktisch unverändert fortgepflanzt werden, zweitens, daß die Lösung an einem Punkt (x, t) nur von den Startwerten im Intervall [x − t, x + t], dem sogenannten Abhängigkeitsintervall, abhängt. Daher pflanzen sich Störungen (Sprünge, Schockwellen) längs der Wellenfronten x±t = const ungedämpft fort. Diese Linien nennt man Charakteristiken. Die aufgezählten Eigenschaften sollten bei den numerischen Verfahren berücksichtigt werden. Dies bedeutet, daß implizite Verfahren (mit großen Schrittweiten) nur angebracht sind, wenn die Lösung von Anfang an sehr glatt ist und daß andernfalls bei expliziten Verfahren die diskreten Abhängigkeitsgebiete an die der Dgl angepaßt werden sollten. In folgender Skizze

[figure: h; τ; t; x]

ist zu erkennen, daß für τ < h die ’diskreten’ Charakteristiken (dicke Linien) des Verfahrens (3.4.9) eine größere Steigung besitzen als die exakten (Steigung ±1), Störungen sich in der Näherung in x-Richtung also schneller fortpflanzen als in der Realität. Im Verfahren (3.4.9) ist daher das Schrittweitenverhältnis s = 1 am günstigsten. Bei allgemeinen Dgln (3.4.2) hängen die Charakteristiken von den Koeffizienten ab, sind im allgemeinen also keine geraden Linien mehr. Hier wird die Diskussion sehr viel schwieriger.

Index

Abhängigkeitsintervall, 81 Adams-

Bashforth-Verfahren, 21 Moulton-Verfahren, 22

Anfangswert, 36 Anfangswertproblem, 5, 9

steif, 30, 78 Anlaufrechnung, 23

BDF-Verfahren, 34 Bilinearform, 70, 71

Charakteristiken, 81

Dahlquist, 26, 35 Diagonaldominanz, 62 Differenzen

-Gleichung, 9, 28

-Quotienten, 34, 42, 55, 78 Dirichlet-Problem, 54, 56, 63, 70 Dreieck-Gitter, 75

Einschritt-Verfahren, 9, 13, 78 Entwicklung

asymptotische, 19, 27 Euler-Verfahren, 9, 79

Fehler -Schätzung, 16, 19, 50 globaler, 15, 26, 44, 49, 63, 68, 73, 75,

79 lokaler, 12, 22, 23, 29, 43, 48, 63, 65, 67

Fundamentalsystem, 7, 36, 37 Funktional

konvexes, 70 lineares, 70

Galerkin-Bedingung, 72 Gauß-

Algorithmus, 40, 44

Quadraturformel, 33 Gragg-Bulirsch-Stoer-Verfahren, 29 Greensche Formel, 70 Gronwall-Lemma, 7 diskret, 13

Kollokationsverfahren, 49

Kondensation, 41

Konsistenz, 12, 22, 25, 48, 63 -Ordnung, 12

Linienmethode, 78

Lipschitzbedingung, 6, 8, 14, 26

M-Matrix, 61, 67, 69, 79 Mehrschrittverfahren, 21 implizit, 34 lineare, 24

Mehrstellenverfahren, 45, 65

Mehrzielmethode, 39 Mittelpunktregel explizite, 27 implizite, 79

Neumann-Randbedingung, 57 Neumann-Reihe, 32, 44, 61 Neville-Algorithmus, 29 Newton-Verfahren, 37, 40 Numerierung

lexikographisch, 56 Schachbrett-, 57

Ordnungs-Barriere, 26, 35 Bedingungen, 12, 25 Steuerung, 24, 29

Part.Differentialgleichung elliptische, 52–54, 70 hyperbolische, 52, 53, 77

parabolische, 53, 77 Poissongleichung, 54, 71 Polynom

charakteristisches, 35 Lagrange-, 23 Legendre-, 33 Newton-, 24

Prädiktor-Korrektor-Verfahren, 23 Prolongation, 47

Randwertproblem, 5, 36, 42, 54 Restriktion, 47, 63, 74 Richardson-Extrapolation, 19, 29 Rosenbrock-Wanner-Verfahren, 34 Rundungsfehler, 20 Runge-Kutta-Verfahren, 50

allgemeine, 11 eingebettete, 17 explizite, 11 implizite, 33 klassisches, 10 linear-implizite, 34 stetige, 19

Schrittweiten-Steuerung, 16, 17, 24, 29, 34 Spline-Funktionen, 50, 75 Stabilität, 13, 14, 26, 48, 62, 65, 79, 80

A-, 33 absolute, 32, 35

Stabilitäts-Bereich, 32 Funktion, 32, 35

Stern Fünf-Punkte-, 55, 65 Neun-Punkte-, 65, 66

Treppenmatrix, 40 Triangulierung, 75

Verfahren explizite, 21, 79 implizite, 22, 31, 79 Wärmeleitungsgleichung, 53, 78, 80 Wellengleichung, 52, 81