Putzer-Algorithmus

Beim Putzer-Algorithmus (nach Eugene James Putzer) handelt es sich um eine rekursive Methode zur Berechnung des Matrixexponential. 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 w​ie auch d​as Matrixexponential v​or allem i​n der Theorie gewöhnlicher Differentialgleichungen e​ine Rolle, d​a durch i​hn Lösungen linearer Differentialgleichungssysteme bestimmt werden können.[1]

Das Verfahren

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 quadratische ${\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{array}{ll}p_{1}'(x)=\lambda _{1}p_{1}(x)\,\\p_{1}(0)=1\end{array}}\right.\quad \left\{{\begin{array}{ll}p_{k}'(x)=\lambda _{k}p_{k}(x)+p_{k-1}\,\\p_{k}(0)=0\,.\end{array}}\right.}$

Beispiel

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

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 ${\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 s​ich wiederum weitere spezielle Lösungen abhängig v​on den Eigenwerten ableiten.

Gleiche Eigenwerte

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

Ungleiche Eigenwerte

Häufig sind alle Eigenwerte der Matrix nicht gleich ${\displaystyle (\lambda _{1}\neq \lambda _{2}\neq \dotsc \neq \lambda _{n})}$. 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 f​olgt hierbei offensichtlich d​er 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 ungleichen 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)}$

Aufwand von Polynom-Lösungsmethoden

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 m​an zum Aufwand d​ie Lösung für d​ie Matrix-Exponentialfunktion mittels Diagonalisieren d​er 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}}$,

hierbei ist ${\displaystyle y_{j}^{T}}$ die ${\displaystyle j}$-te Zeile von ${\displaystyle V^{-1}}$, 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})}}}$ ${\displaystyle ~~~~~(k\neq j)}$

ist, d​ann 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]

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.

Einzelnachweise

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]).