Trapez-Methode

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

${\displaystyle 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 ${\displaystyle y'=\mathrm {i} \alpha y}$ kein Amplitudenfehler auftritt^[1]. Das Verfahren lässt sich aus der Trapezregel herleiten:

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

mit

${\displaystyle 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]

${\displaystyle {\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 ${\displaystyle k}$-ten Schritt den Integranden wie folgt

${\displaystyle \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^[3]

${\displaystyle 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:

${\displaystyle 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

${\displaystyle (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

${\displaystyle J^{(k)}:=\left({\frac {\partial f}{\partial y}}\right)_{n+1}^{(k)}}$,

${\displaystyle I}$ die Einheitsmatrix und ${\displaystyle k}$ der Iterationsschritt ist.

Stabilität

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

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

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

Schrittweite h

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

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

${\displaystyle \delta }$ bezeichnet den zugelassenen lokalen Diskretisierungsfehler. Der Ansatz ${\displaystyle 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

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

Dabei ist ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle 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.