Beim Putzer-Algorithmus (nach Eugene James Putzer) handelt es sich um eine rekursive Methode zur Berechnung des Matrixexponentials. Ziel des Verfahrens ist es also, für gegebenes
und
das Matrixexponential
zu bestimmen.
Der Algorithmus spielt wie auch das Matrixexponential vor allem in der Theorie gewöhnlicher Differentialgleichungen eine Rolle, da durch ihn Lösungen linearer Differentialgleichungssysteme bestimmt werden können.[1]
Sei
eine quadratische
-Matrix. Sei weiter
eine Abzählung der Eigenwerte der Matrix
unter Berücksichtigung der algebraischen Vielfachheit. Diese Abzählung existiert, da
algebraisch abgeschlossen ist.
Für
definieren wir nun rekursiv stetig differenzierbare Funktionen
und
-Matrizen
, so dass für
folgende Beziehung gilt:
.
Die Rekursionsvorschrift für die Matrizen
lautet
und
.
Die
sind rekursiv durch folgende Anfangswertprobleme definiert:
![{\displaystyle \left\{{\begin{aligned}p_{1}'(x)&=\lambda _{1}p_{1}(x)\\p_{1}(0)&=1\end{aligned}}\right.\quad \left\{{\begin{aligned}p_{k}'(x)&=\lambda _{k}p_{k}(x)+p_{k-1}\\p_{k}(0)&=0\end{aligned}}\right.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c8908df377749ba1fb1f9aea3b409ec52d09ee07)
Illustration des Verfahrens für
: Sei folgende Matrix
gegeben. Wir werden nun mit dem Putzer-Algorithmus das Matrixexponential
für beliebiges
berechnen. Als erstes bestimmen wir die Eigenwerte der Matrix
unter Berücksichtigung der algebraischen Vielfachheit. Dazu berechnen wir das charakteristische Polynom und setzen es gleich
:
![{\displaystyle \chi _{A}(\lambda ):=\det(A-\lambda I)=\det {\begin{pmatrix}3-\lambda &{-1}\\1&{1-\lambda }\end{pmatrix}}=\lambda ^{2}-4\lambda +4=(\lambda -2)^{2}\,{\overset {!}{=}}\,0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/69aa11a13ed1f5c673c778e8bd69bd862b66e8ff)
Es folgt, dass
der einzige Eigenwert der Matrix
ist, allerdings mit algebraischer Vielfachheit
. Somit handelt es sich bei
um eine Abzählung der Eigenwerte von
.
Als nächstes bestimmen wir mit den Rekursionsformeln von oben die Matrizen
. Es folgt
und entsprechend
.
Zuletzt berechnen wir für
die Funktionen
, die über folgende Anfangswertprobleme definiert sind:
![{\displaystyle \left\{{\begin{array}{ll}p_{1}'(x)=2p_{1}(x)\,\\p_{1}(0)=1\end{array}}\right.\quad \left\{{\begin{array}{ll}p_{2}'(x)=2p_{2}(x)+p_{1}\,\\p_{2}(0)=0\end{array}}\right.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c0c614a3fbac8cc4b119e27aae56bbf743b35eb1)
Das erste Anfangswertproblem lässt sich mittels der Lösungsmethode Trennung der Veränderlichen für gewöhnliche Differentialgleichungen leicht als
bestimmen. Das zweite Anfangswertproblem lässt sich ebenfalls sehr einfach über die Methode Variation der Konstanten berechnen und wir erhalten als Lösung
.
Da wir nun alles beisammen haben, können wir
für ein
angeben:
![{\displaystyle \exp(Ax)=\sum _{k=1}^{2}p_{k}(x)M_{k-1}=\exp(2x){\begin{pmatrix}1&{0}\\0&{1}\end{pmatrix}}+x\exp(2x){\begin{pmatrix}1&{-1}\\1&{-1}\end{pmatrix}}=\exp(2x){\begin{pmatrix}1+x&{-x}\\x&{1-x}\end{pmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6612b6cb2b5b4d2c2dd089242d0ed37a1682038a)
Wie im Beispiel gezeigt, lässt sich das erste Anfangswertproblem mittels der Lösungsmethode Trennung der Veränderlichen für gewöhnliche Differentialgleichungen bestimmen. Die weiteren Anfangswertprobleme mit
lassen sich ebenfalls einfach über die Methode Variation der Konstanten berechnen. Dementsprechend folgen zunächst die Lösungen:
![{\displaystyle p_{1}(t)=e^{\lambda _{1}t}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/16a9a8f2c9afb0cdbcb04cfa797a8c6050e2fe64)
![{\displaystyle p_{k}(t)=e^{\lambda _{k}t}\int _{x=0}^{t}p_{k-1}(x)\,e^{-\lambda _{k}x}\,\mathrm {d} x}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ed06e5ddaf5142f091844487ff783e4673523111)
Hieraus lassen sich wiederum weitere spezielle Lösungen abhängig von den Eigenwerten ableiten.
Sind alle Eigenwerte gleich
, folgt aus der integralen Darstellung für
mit
, dass die Funktion
gerade
-mal integriert werden muss. Für
gilt somit:
![{\displaystyle p_{k}(t)={\frac {t^{k-1}}{(k-1)!}}e^{\lambda t}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/749b697928ea0ae90f80e03cec618fd0926e8114)
Das Matrix-Exponential einer
-Matrix mit gleichen Eigenwerten
kann schließlich explizit mit folgender Formel bestimmt werden:
![{\displaystyle \exp(At)=e^{\lambda t}E+t\,e^{\lambda t}(A-\lambda \,E)+\dotsc =e^{\lambda t}\,\sum _{i=0}^{n-1}{\frac {t^{i}}{i!}}(A-\lambda \,E)^{i}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8fa58c2c086ddb2c3a93a76f804a5b2d7d7b019d)
Häufig sind alle Eigenwerte der Matrix verschieden (
für
). In diesem Fall gilt für die Diskriminante des charakteristischen Polynoms
. Die Lösung für
bleibt davon unverändert. Ab
folgt nach Integration:
,
![{\displaystyle p_{3}(t)={\frac {e^{\lambda _{1}t}-e^{\lambda _{3}t}}{(\lambda _{1}-\lambda _{2})(\lambda _{1}-\lambda _{3})}}+{\frac {e^{\lambda _{2}t}-e^{\lambda _{3}t}}{(\lambda _{2}-\lambda _{1})(\lambda _{2}-\lambda _{3})}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6c4a3cb8f83ec78da9c26ddc530bfba31c9cf840)
- usw.
Die Lösung folgt hierbei offensichtlich der Systematik
mit
.
Für das Exponential einer
-Matrix mit verschiedenen Eigenwerten folgt damit die Newton-Interpolation[2]
![{\displaystyle \exp(At)=e^{\lambda _{1}t}E+\sum _{i=2}^{n}[\lambda _{1},\dotsc ,\lambda _{i}]\prod _{j=1}^{i-1}(A-\lambda _{j}\,E)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4551ca4eaa42787b4a354ba732b384f4136fb35d)
für die Darstellung der Matrix-Exponentialfunktion als Polynom. Die von
abhängigen Funktionen können dabei rekursiv berechnet werden mit
,
.
Der Putzer-Algorithmus hat einen Aufwand der Größenordnung
(im Wesentlichen
Matrizenmultiplikationen). Damit ist der Aufwand deutlich höher als bei Verwendung von Methoden auf Basis der Taylorreihe oder auf Basis von Matrix-Zerlegungsmethoden (wie Diagonalisierung oder QR-Algorithmus), mit jeweils einem Aufwand der Größenordnung
. Der Algorithmus ist also eher nur für kleinere Matrizen geeignet, jedoch erhält man hier vorteilhaft eine vollständig symbolische Lösung.
Vergleicht man zum Aufwand die Lösung für die Matrix-Exponentialfunktion mittels Diagonalisieren der Matrix
,
wobei
die
-te Zeile von
ist, mit der Lagrange-Interpolation
,
wobei
,
dann folgt daraus
![{\displaystyle A_{j}=v_{j}y_{j}^{T}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/daa6515dfa935b8d063f3790e091c45b5d727235)
und der Aufwand der Größenordnung
zur Berechnung von
ist völlig unnötig.[3]
- Paul Waltman: Differential Equations and Dynamical Systems, Dover Publications, 2004, ISBN 978-0486434780, S. 49–60
- Eugene J. Putzer: Avoiding the Jordan Canonical Form in the Discussion of Linear Systems with Constant Coefficients, The American Mathematical Monthly, Vol. 73(1), 1966, S. 2–7, DOI:10.1080/00029890.1966.11970714.
- DorFuchs: Der Putzer-Algorithmus, den kaum jemand kennt, zur Bestimmung der Matrixexponentialfunktion auf YouTube, 13. Februar 2018, abgerufen am 9. Februar 2023.
- ↑ Lebovitz 2016 (.pdf)
- ↑ Moler: Cleve Moler, Charles F. Van Loan: Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. In: SIAM Review. Band 45, Nr. 1, 2003, ISSN 1095-7200, S. 16 - 18, doi:10.1137/S00361445024180 (cornell.edu [PDF]).
- ↑ Moler2: Cleve Moler, Charles F. Van Loan: Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. In: SIAM Review. Band 45, Nr. 1, 2003, ISSN 1095-7200, S. 22 - 23, doi:10.1137/S00361445024180 (cornell.edu [PDF]).