Trapez-Methode

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

<math>y'(t) = f\left(t, y(t)\right), \quad

      y(t_0) = y_0

</math>

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 <math>y'= \mathrm{i}\alpha y</math> kein Amplitudenfehler auftritt<ref>M. Kloker: Numerische Löser (Zeitintegrationsverfahren) für die Gewöhnliche Modelldifferentialgleichung y'=αy (PDF; 2,2 MB), Universität Stuttgart, 1996</ref>. Das Verfahren lässt sich aus der Trapezregel herleiten:

mit

Herleitung

Für die Herleitung von Einschrittverfahren wird das Anfangswertproblem meist in der zu ihr äquivalenten Integralgleichung umgeformt<ref>{{#invoke:Vorlage:Literatur|f}}</ref>

<math>

\begin{align}

 \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 ds \,.

\end{align} </math> Nun besteht die Idee bei der impliziten Trapez-Methode eine simple Quadraturformel für das Integral zu benutzen: die Trapezregel. Man approximiert in jedem <math> k </math>-ten Schritt den Integranden wie folgt

<math>

\int_{t_k}^{t_{k+1}} f(s,y(s)) \, \mathrm ds

 \approx \frac{h}{2} \Big( f(t_k,y_k) + f(t_{k+1},y_{k+1}) \Big) \,.

</math> Zusammen ergibt dies die Trapez-Methode<ref>{{#invoke:Vorlage:Literatur|f}}</ref>

<math>

y(t_{k+1})

 = y(t_k) + \int_{t_k}^{t_{k+1}} f(s,y(s)) \, \mathrm ds
 \approx y_k + \frac{h}{2} \Big( f(t_k,y_k) + f(t_{k+1},y_{k+1}) \Big)
 =: y_{k+1} \,.

</math>

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:

<math>y^{(k + 1)}_{n+1} =

y^{(k)}_{n+1} - \left(I - \frac{h}{2}\frac{\partial f^{(k)}_{n+1}}{\partial y^{(k)}_{n+1}}\right)^{-1}\left( y^{(k)}_{n+1} - y_n - \frac{h}{2}(f^{(k)}_{n+1} + f_n) \right). </math>

Man erhält also ein lineares Gleichungssystem

<math>(I-\frac{h}{2} J^{(k)})y^{(k + 1)}_{n+1} =

-\frac{h}{2} J^{(k)} y^{(k)}_{n+1} +y_n + \frac{h}{2}(f^{(k)}_{n+1} + f_n),</math>

wobei J die Jacobi-Matrix

<math>J^{(k)} := \left(\frac{\partial f}{\partial y}\right)^{(k)}_{n+1}</math>,

<math>I</math> die Einheitsmatrix und <math>k</math> der Iterationsschritt ist.

Stabilität

Mit der Testgleichung <math>y'(t)=\lambda y(t)</math> bekommt man die Stabilitätsfunktion

<math> R(z)=\frac{2+z}{2-z},\quad z=h\lambda\in\Complex.</math>

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

Schrittweite h

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

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

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

<math>R(h\lambda)=\frac{2+h\lambda}{2-h\lambda}</math>.

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

Literatur

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

Einzelnachweise