Trapez-Methode

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

$y^{\prime }(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^{\prime }=\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{array}{cc}\dot{y}& =f(t,y),y(t_0)=y_0\\ ⟺y(t)& =y_0+\underset{t_0}{\overset{t}{\int }}f(s,y(s))\mathrm{d}s.\end{array}$

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

$\underset{t_k}{\overset{t_{k+1}}{\int }}f(s,y(s))\mathrm{d}s\approx \frac{h}{2}(f(t_k,y_k)+f(t_{k+1},y_{k+1})).$

Zusammen ergibt dies die Trapez-Methode

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

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

Stabilität

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

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

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:

$|\frac{R(h\lambda )}{{\mathrm{e}}^{h\lambda }}-1|=\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 :=\underset{j}{\max }|{\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 \underset{i,j}{\max }|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.