Trapez-Methode

Das implizite Trapez-Verfahren ist ein Verfahren zur numerischen Lösung eines Anfangswertproblems

y^{'} (t) = f (t, y (t)), y (t_{0}) = y_{0}

Es lässt sich sowohl den Runge-Kutta-Verfahren als auch den Adams-Moulton-Verfahren zuordnen. Das Trapezverfahren ist A-stabil mit der Besonderheit, dass für die Schwingungsgleichung $y^{'} = i α y$ kein Amplitudenfehler auftritt^[1]. Das Verfahren lässt sich aus der Trapezregel herleiten:

y_{n + 1} = y_{n} + \frac{h}{2} (f_{n + 1} + f_{n})

mit

f_{n} := f (t_{n}, y_{n}) .

Herleitung

Für die Herleitung von Einschrittverfahren wird das Anfangswertproblem meist in der zu ihr äquivalenten Integralgleichung umgeformt^[2]

\begin{matrix} \dot{y} & = f (t, y), y (t_{0}) = y_{0} \\ ⟺ y (t) & = y_{0} + \int_{t_{0}}^{t} f (s, y (s)) d s . \end{matrix}

Nun besteht die Idee bei der impliziten Trapez-Methode eine simple Quadraturformel für das Integral zu benutzen: die Trapezregel. Man approximiert in jedem $k$ -ten Schritt den Integranden wie folgt

\int_{t_{k}}^{t_{k + 1}} f (s, y (s)) d s \approx \frac{h}{2} (f (t_{k}, y_{k}) + f (t_{k + 1}, y_{k + 1})) .

Zusammen ergibt dies die Trapez-Methode^[3]

y (t_{k + 1}) = y (t_{k}) + \int_{t_{k}}^{t_{k + 1}} f (s, y (s)) d s \approx y_{k} + \frac{h}{2} (f (t_{k}, y_{k}) + f (t_{k + 1}, y_{k + 1})) = : y_{k + 1} .

Lösungsmethode

Zur Lösung dieses, in der Regel nichtlinearen, Gleichungssystems können verschiedene numerische Verfahren genutzt werden. Für das quadratisch konvergente Newton-Verfahren ergibt sich konkret:

y_{n + 1}^{(k + 1)} = y_{n + 1}^{(k)} - {(I - \frac{h}{2} \frac{\partial f_{n + 1}^{(k)}}{\partial y_{n + 1}^{(k)}})}^{- 1} (y_{n + 1}^{(k)} - y_{n} - \frac{h}{2} (f_{n + 1}^{(k)} + f_{n})) .

Man erhält also ein lineares Gleichungssystem

(I - \frac{h}{2} J^{(k)}) y_{n + 1}^{(k + 1)} = - \frac{h}{2} J^{(k)} y_{n + 1}^{(k)} + y_{n} + \frac{h}{2} (f_{n + 1}^{(k)} + f_{n}),

wobei J die Jacobi-Matrix

J^{(k)} := {(\frac{\partial f}{\partial y})}_{n + 1}^{(k)}

,

$I$ die Einheitsmatrix und $k$ der Iterationsschritt ist.

Stabilität

Mit der Testgleichung $y^{'} (t) = λ y (t)$ bekommt man die Stabilitätsfunktion

R (z) = \frac{2 + z}{2 - z}, z = h λ \in ℂ .

Auf der imaginären Achse $z = i η$ gilt $| R (i η) | = 1$ , daher ist die Trapezmethode A-stabil.

Schrittweite h

Die (variable) Schrittweite kann aus folgender Beziehung berechnet werden:

| \frac{R (h λ)}{e^{h λ}} - 1 | = δ

;

$δ$ bezeichnet den zugelassenen lokalen Diskretisierungsfehler. Der Ansatz $y_{n + 1} = y_{n} + \frac{h}{2} (f_{n + 1} + f_{n}) = : R (h λ) y_{n}$ liefert für die implizite Trapez-Methode

R (h λ) = \frac{2 + h λ}{2 - h λ}

.

Dabei ist $λ := max_{j} | λ_{j} |$ der Betrag des betragsmäßig größten Eigenwerts der Jacobi-Matrix (Spektralradius). Die numerische Bestimmung der Eigenwerte ist sehr zeitaufwendig; für den Zweck der Schrittweitenberechnung ist es im Allgemeinen ausreichend die Gesamtnorm $λ = N \cdot max_{i, j} | a_{i j} |$ heranzuziehen, die immer größer oder gleich der Spektralnorm ist. N ist der Rang der Jacobi-Matrix und $a_{i j}$ deren Elemente.

Literatur

Hans R. Schwarz, Norbert Köckler: Numerische Mathematik. 5. Auflage, Teubner, Stuttgart 2004, ISBN 3-519-42960-8, S. 343.

Einzelnachweise

↑ M. Kloker: Numerische Löser (Zeitintegrationsverfahren) für die Gewöhnliche Modelldifferentialgleichung y'=αy (PDF; 2,2 MB), Universität Stuttgart, 1996
↑ Reusken, Arnold.: Numerik für Ingenieure und Naturwissenschaftler. Springer, Berlin 2006, ISBN 3-540-25544-3, S. 378.
↑ Reusken, Arnold.: Numerik für Ingenieure und Naturwissenschaftler. Springer, Berlin 2006, ISBN 3-540-25544-3, S. 383.

[1] M. Kloker: Numerische Löser (Zeitintegrationsverfahren) für die Gewöhnliche Modelldifferentialgleichung y'=αy (PDF; 2,2 MB), Universität Stuttgart, 1996

[2] Reusken, Arnold.: Numerik für Ingenieure und Naturwissenschaftler. Springer, Berlin 2006, ISBN 3-540-25544-3, S. 378.

[3] Reusken, Arnold.: Numerik für Ingenieure und Naturwissenschaftler. Springer, Berlin 2006, ISBN 3-540-25544-3, S. 383.

[1]

[2]

[3]