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 $A\in \mathbb {C} ^{n\times n}$ und $x\in \mathbb {R}$ das Matrixexponential $\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.

Das Verfahren

Sei $A\in \mathbb {C} ^{n\times n}$ eine quadratische $(n\times n)$-Matrix. Sei weiter $\{\lambda _{j}\}_{j=1,\ldots ,n}$ eine Abzählung der Eigenwerte der Matrix $A$ unter Berücksichtigung der algebraischen Vielfachheit. Diese Abzählung existiert, da $\mathbb {C}$ algebraisch abgeschlossen ist.

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

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

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

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

Die $p_{k}$ sind rekursiv durch folgende Anfangswertprobleme definiert:

$\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.$

Beispiel

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

$\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 $\lambda =2$ der einzige Eigenwert der Matrix $A$ ist, allerdings mit algebraischer Vielfachheit $2$. Somit handelt es sich bei $(\lambda _{1},\lambda _{2})=(2,2)$ um eine Abzählung der Eigenwerte von $A$.

Als nächstes bestimmen wir mit den Rekursionsformeln von oben die Matrizen $M_{k},k\in \{0,1\}$. Es folgt $M_{0}={\begin{pmatrix}1&{0}\\0&{1}\end{pmatrix}}$ und entsprechend $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 $k\in \{1,2\}$ die Funktionen $p_{k}$, die über folgende Anfangswertprobleme definiert sind:

$\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 $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 $p_{2}(x)=x\exp(2x)$.

Da wir nun alles beisammen haben, können wir $\exp(Ax)$ für ein $x\in \mathbb {R}$ angeben:

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

Spezielle Lösungen

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 $k\in \{2,\ldots ,n\}$ lassen sich ebenfalls einfach über die Methode Variation der Konstanten berechnen. Dementsprechend folgen zunächst die Lösungen:

$p_{1}(t)=e^{\lambda _{1}t}$

$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.

Gleiche Eigenwerte

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

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

Das Matrix-Exponential einer $(n\times n)$-Matrix mit gleichen Eigenwerten $\lambda$ kann schließlich explizit mit folgender Formel bestimmt werden:

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

Ungleiche Eigenwerte

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

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

$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

$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 $j\neq i$.

Für das Exponential einer $(n\times n)$-Matrix mit verschiedenen Eigenwerten folgt damit die Newton-Interpolation

$\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 $x$ abhängigen Funktionen können dabei rekursiv berechnet werden mit

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

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

Aufwand von Polynom-Lösungsmethoden

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

$\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 $y_{j}^{T}$ die $j$-te Zeile von $V^{-1}$ ist, mit der Lagrange-Interpolation

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

wobei

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

dann folgt daraus

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

und der Aufwand der Größenordnung ${\mathcal {O}}(n^{4})$ zur Berechnung von $A_{j}$ ist völlig unnötig.

Literatur

• 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.