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

Für $k\in \{1,\dots ,n\}$ definieren wir nun rekursiv stetig differenzierbare Funktionen $p_k\in C^1(ℝ,ℂ)$ und $(n\times n)$-Matrizen $M_{k-1}\in {ℂ}^{n\times n}$, so dass für $x\in ℝ$ 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:

$\{\begin{array}{cc}{p_1}^{\prime }(x)& ={\lambda }_1{p}_1(x)\\ p_1(0)& =1\end{array}\{\begin{array}{cc}{p_k}^{\prime }(x)& ={\lambda }_k{p}_k(x)+p_{k-1}\\ p_k(0)& =0\end{array}$

Beispiel

Illustration des Verfahrens für $n=2$: Sei folgende Matrix $A=\left(\begin{array}{cc}3& -1\\ 1& 1\end{array}\right)$ gegeben. Wir werden nun mit dem Putzer-Algorithmus das Matrixexponential $\exp (Ax)$ für beliebiges $x\in ℝ$ 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 \left(\begin{array}{cc}3-\lambda & -1\\ 1& 1-\lambda \end{array}\right)={\lambda }^2-4\lambda +4=(\lambda -2{)}^2\stackrel{!}{=}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=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$ und entsprechend $M_1=(A-{\lambda }_1{I}_2)M_0=(A-{\lambda }_1{I}_2)=\left(\begin{array}{cc}1& -1\\ 1& -1\end{array}\right)$.

Zuletzt berechnen wir für $k\in \{1,2\}$ die Funktionen $p_k$, die über folgende Anfangswertprobleme definiert sind:

$\{\begin{array}{c}{p_1}^{\prime }(x)=2p_1(x)\\ p_1(0)=1\end{array}\{\begin{array}{c}{p_2}^{\prime }(x)=2p_2(x)+p_1\\ p_2(0)=0\end{array}$

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 ℝ$ angeben:

$\exp (Ax)=\sum _{k=1}^2{p}_k(x)M_{k-1}=\exp (2x)\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+x\exp (2x)\left(\begin{array}{cc}1& -1\\ 1& -1\end{array}\right)=\exp (2x)\left(\begin{array}{cc}1+x& -x\\ x& 1-x\end{array}\right)$

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,\dots ,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}\underset{x=0}{\overset{t}{\int }}{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=\dots ={\lambda }_n)$, folgt aus der integralen Darstellung für $p_k$ mit $k\ge 2$, dass die Funktion $f=1$ gerade $(k-1)$-mal integriert werden muss. Für $k\ge 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)+\dots =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\ne {\lambda }_j$ für $i\ne j$). In diesem Fall gilt für die Diskriminante des charakteristischen Polynoms $D_n\ne 0$. Die Lösung für $p_1(t)$ bleibt davon unverändert. Ab $k\ge 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\ne 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,\dots ,{\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,\dots ,{\lambda }_{k+1}]=\frac{[{\lambda }_1,\dots ,{\lambda }_k]-[{\lambda }_2,\dots ,{\lambda }_{k+1}]}{{\lambda }_1-{\lambda }_{k+1}}$ $(k\ge 2)$.

Aufwand von Polynom-Lösungsmethoden

Der Putzer-Algorithmus hat einen Aufwand der Größenordnung $𝒪(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 $𝒪(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\mathrm{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}(k\ne j)$,

dann folgt daraus

$A_j=v_j{y}_j^T$

und der Aufwand der Größenordnung $𝒪(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.
• Vorlage:YouTube