Putzer-Algorithmus

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 ${\displaystyle A\in \mathbb {C} ^{n\times n}}$ und ${\displaystyle x\in \mathbb {R} }$ das Matrixexponential ${\displaystyle \textstyle \exp(Ax):=\sum _{k=0}^{\infty }{\frac {(Ax)^{k}}{k!}}}$ 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 ${\displaystyle A\in \mathbb {C} ^{n\times n}}$ eine quadratische ${\displaystyle (n\times n)}$-Matrix. Sei weiter ${\displaystyle \{\lambda _{j}\}_{j=1,\ldots ,n}}$ eine Abzählung der Eigenwerte der Matrix ${\displaystyle A}$ unter Berücksichtigung der algebraischen Vielfachheit. Diese Abzählung existiert, da ${\displaystyle \mathbb {C} }$ algebraisch abgeschlossen ist.

Für ${\displaystyle k\in \{1,\ldots ,n\}}$ definieren wir nun rekursiv stetig differenzierbare Funktionen ${\displaystyle p_{k}\in C^{1}(\mathbb {R} ,\mathbb {C} )}$ und ${\displaystyle (n\times n)}$-Matrizen ${\displaystyle M_{k-1}\in \mathbb {C} ^{n\times n}}$, so dass für ${\displaystyle x\in \mathbb {R} }$ folgende Beziehung gilt:

${\displaystyle \exp(Ax)=\sum _{k=1}^{n}p_{k}(x)M_{k-1}}$.

Die Rekursionsvorschrift für die Matrizen ${\displaystyle M_{k}}$ lautet

${\displaystyle M_{0}:=I_{n}}$ und ${\displaystyle M_{k}:=(A-\lambda _{k}\cdot I)M_{k-1}}$.

Die ${\displaystyle p_{k}}$ 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.}$

Illustration des Verfahrens für ${\displaystyle n=2}$: Sei folgende Matrix ${\displaystyle A={\begin{pmatrix}3&{-1}\\1&{1}\end{pmatrix}}}$ gegeben. Wir werden nun mit dem Putzer-Algorithmus das Matrixexponential ${\displaystyle \exp(Ax)}$ für beliebiges ${\displaystyle x\in \mathbb {R} }$ berechnen. Als erstes bestimmen wir die Eigenwerte der Matrix ${\displaystyle A}$ unter Berücksichtigung der algebraischen Vielfachheit. Dazu berechnen wir das charakteristische Polynom und setzen es gleich ${\displaystyle 0}$:

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

Es folgt, dass ${\displaystyle \lambda =2}$ der einzige Eigenwert der Matrix ${\displaystyle A}$ ist, allerdings mit algebraischer Vielfachheit ${\displaystyle 2}$. Somit handelt es sich bei ${\displaystyle (\lambda _{1},\lambda _{2})=(2,2)}$ um eine Abzählung der Eigenwerte von ${\displaystyle A}$.

Als nächstes bestimmen wir mit den Rekursionsformeln von oben die Matrizen ${\displaystyle M_{k},k\in \{0,1\}}$. Es folgt ${\displaystyle M_{0}={\begin{pmatrix}1&{0}\\0&{1}\end{pmatrix}}}$ und entsprechend ${\displaystyle M_{1}=(A-\lambda _{1}I_{2})M_{0}=(A-\lambda _{1}I_{2})={\begin{pmatrix}1&{-1}\\1&{-1}\end{pmatrix}}}$.

Zuletzt berechnen wir für ${\displaystyle k\in \{1,2\}}$ die Funktionen ${\displaystyle p_{k}}$, 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.}$

Das erste Anfangswertproblem lässt sich mittels der Lösungsmethode Trennung der Veränderlichen für gewöhnliche Differentialgleichungen leicht als ${\displaystyle p_{1}(x)=\exp(2x)}$ bestimmen. Das zweite Anfangswertproblem lässt sich ebenfalls sehr einfach über die Methode Variation der Konstanten berechnen und wir erhalten als Lösung ${\displaystyle p_{2}(x)=x\exp(2x)}$.

Da wir nun alles beisammen haben, können wir ${\displaystyle \exp(Ax)}$ für ein ${\displaystyle x\in \mathbb {R} }$ 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}}}$

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 ${\displaystyle k\in \{2,\ldots ,n\}}$ 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}}$

${\displaystyle p_{k}(t)=e^{\lambda _{k}t}\int _{x=0}^{t}p_{k-1}(x)\,e^{-\lambda _{k}x}\,\mathrm {d} x}$

Hieraus lassen sich wiederum weitere spezielle Lösungen abhängig von den Eigenwerten ableiten.

Sind alle Eigenwerte gleich ${\displaystyle (\lambda =\lambda _{1}=\lambda _{2}=\dotsc =\lambda _{n})}$, folgt aus der integralen Darstellung für ${\displaystyle p_{k}}$ mit ${\displaystyle k\geq 2}$, dass die Funktion ${\displaystyle f=1}$ gerade ${\displaystyle (k-1)}$-mal integriert werden muss. Für ${\displaystyle k\geq 1}$ gilt somit:

${\displaystyle p_{k}(t)={\frac {t^{k-1}}{(k-1)!}}e^{\lambda t}}$

Das Matrix-Exponential einer ${\displaystyle (n\times n)}$-Matrix mit gleichen Eigenwerten ${\displaystyle \lambda }$ 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}}$

Häufig sind alle Eigenwerte der Matrix verschieden (${\displaystyle \lambda _{i}\neq \lambda _{j}}$ für ${\displaystyle i\neq j}$). In diesem Fall gilt für die Diskriminante des charakteristischen Polynoms ${\displaystyle D_{n}\neq 0}$. Die Lösung für ${\displaystyle p_{1}(t)}$ bleibt davon unverändert. Ab ${\displaystyle k\geq 2}$ folgt nach Integration:

${\displaystyle p_{2}(t)={\frac {e^{\lambda _{1}t}-e^{\lambda _{2}t}}{\lambda _{1}-\lambda _{2}}}}$,

${\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})}}}$

usw.

Die Lösung folgt hierbei offensichtlich der Systematik

${\displaystyle p_{k}(t)=\sum _{i=1}^{k-1}{\frac {e^{\lambda _{i}t}-e^{\lambda _{k}t}}{\prod _{j=1}^{k}(\lambda _{i}-\lambda _{j})}}~~~}$ mit ${\displaystyle j\neq i}$.

Für das Exponential einer ${\displaystyle (n\times n)}$-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)}$

für die Darstellung der Matrix-Exponentialfunktion als Polynom. Die von ${\displaystyle x}$ abhängigen Funktionen können dabei rekursiv berechnet werden mit

${\displaystyle [\lambda _{1},\lambda _{2}]={\frac {e^{\lambda _{1}t}-e^{\lambda _{2}t}}{\lambda _{1}-\lambda _{2}}}}$,

${\displaystyle [\lambda _{1},\dotsc ,\lambda _{k+1}]={\frac {[\lambda _{1},\dotsc ,\lambda _{k}]-[\lambda _{2},\dotsc ,\lambda _{k+1}]}{\lambda _{1}-\lambda _{k+1}}}}$ ${\displaystyle ~~~~(k\geq 2)}$.

Der Putzer-Algorithmus hat einen Aufwand der Größenordnung ${\displaystyle {\mathcal {O}}(n^{4})}$ (im Wesentlichen ${\displaystyle n}$ 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 ${\displaystyle {\mathcal {O}}(n^{3})}$. 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

${\displaystyle \exp(At)=V\operatorname {diag} (e^{\lambda _{i}t})\,V^{-1}=\sum _{j=1}^{n}e^{\lambda _{i}t}v_{j}y_{j}^{T}}$,

wobei ${\displaystyle y_{j}^{T}}$ die ${\displaystyle j}$-te Zeile von ${\displaystyle V^{-1}}$ ist, mit der Lagrange-Interpolation

${\displaystyle \exp(At)=\sum _{j=1}^{n}e^{\lambda _{i}t}A_{j}}$,

wobei

${\displaystyle A_{j}=\prod _{k=1}^{n}{\frac {A-\lambda _{k}\,E}{\lambda _{j}-\lambda _{k}}}\quad (k\neq j)}$,

dann folgt daraus

${\displaystyle A_{j}=v_{j}y_{j}^{T}}$

und der Aufwand der Größenordnung ${\displaystyle {\mathcal {O}}(n^{4})}$ zur Berechnung von ${\displaystyle A_{j}}$ 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.

1. ↑ Lebovitz 2016 (.pdf)
2. ↑ 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]). 
3. ↑ 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]).