Trapez-Methode

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

$y'(t)=f\left(t,y(t)\right),\quad 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'=\mathrm {i} \alpha y$ kein Amplitudenfehler auftritt. 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

${\begin{aligned}{\dot {y}}&=f(t,y)\,,\qquad y(t_{0})=y_{0}\\\Longleftrightarrow \quad y(t)&=y_{0}+\int _{t_{0}}^{t}f(s,y(s))\,\mathrm {d} s\,.\end{aligned}}$

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))\,\mathrm {d} s\approx {\frac {h}{2}}{\Big (}f(t_{k},y_{k})+f(t_{k+1},y_{k+1}){\Big )}\,.$

Zusammen ergibt dies die Trapez-Methode

$y(t_{k+1})=y(t_{k})+\int _{t_{k}}^{t_{k+1}}f(s,y(s))\,\mathrm {d} s\approx y_{k}+{\frac {h}{2}}{\Big (}f(t_{k},y_{k})+f(t_{k+1},y_{k+1}){\Big )}=: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)}-\left(I-{\frac {h}{2}}{\frac {\partial f_{n+1}^{(k)}}{\partial y_{n+1}^{(k)}}}\right)^{-1}\left(y_{n+1}^{(k)}-y_{n}-{\frac {h}{2}}(f_{n+1}^{(k)}+f_{n})\right).$

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)}:=\left({\frac {\partial f}{\partial y}}\right)_{n+1}^{(k)}$,

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

Stabilität

Mit der Testgleichung $y'(t)=\lambda y(t)$ bekommt man die Stabilitätsfunktion

$R(z)={\frac {2+z}{2-z}},\quad z=h\lambda \in \mathbb {C} .$

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

Schrittweite h

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

$\left\vert {\frac {R(h\lambda )}{\mathrm {e} ^{h\lambda }}}-1\right\vert =\delta$;

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

$R(h\lambda )={\frac {2+h\lambda }{2-h\lambda }}$.

Dabei ist $\lambda \,:=\max _{j}{|\lambda _{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 $\lambda =N\cdot \max _{i,j}|a_{ij}|$ heranzuziehen, die immer größer oder gleich der Spektralnorm ist. N ist der Rang der Jacobi-Matrix und $a_{ij}$ deren Elemente.

Literatur

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