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 -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)}}$ ,

die Einheitsmatrix und 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

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

$\left\vert {\frac {R(h\lambda )}{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.

Basierend auf einem Artikel in:

Wikipedia.de