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

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

exp ⁡ ( A x ) = ∑ k = 1 n p k ( x ) M k − 1 {\displaystyle \exp(Ax)=\sum _{k=1}^{n}p_{k}(x)M_{k-1}} .

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

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

Die p k {\displaystyle p_{k}} sind rekursiv durch folgende Anfangswertprobleme definiert:

{ p 1 ′ ( x ) = λ 1 p 1 ( x ) p 1 ( 0 ) = 1 { p k ′ ( x ) = λ k p k ( x ) + p k − 1 p k ( 0 ) = 0 {\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 n = 2 {\displaystyle n=2} : Sei folgende Matrix A = ( 3 − 1 1 1 ) {\displaystyle A={\begin{pmatrix}3&{-1}\\1&{1}\end{pmatrix}}} gegeben. Wir werden nun mit dem Putzer-Algorithmus das Matrixexponential exp ⁡ ( A x ) {\displaystyle \exp(Ax)} für beliebiges x ∈ R {\displaystyle x\in \mathbb {R} } berechnen. Als erstes bestimmen wir die Eigenwerte der Matrix A {\displaystyle A} unter Berücksichtigung der algebraischen Vielfachheit. Dazu berechnen wir das charakteristische Polynom und setzen es gleich 0 {\displaystyle 0} :

χ A ( λ ) := det ( A − λ I ) = det ( 3 − λ − 1 1 1 − λ ) = λ 2 − 4 λ + 4 = ( λ − 2 ) 2 = ! 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 λ = 2 {\displaystyle \lambda =2} der einzige Eigenwert der Matrix A {\displaystyle A} ist, allerdings mit algebraischer Vielfachheit 2 {\displaystyle 2} . Somit handelt es sich bei ( λ 1 , λ 2 ) = ( 2 , 2 ) {\displaystyle (\lambda _{1},\lambda _{2})=(2,2)} um eine Abzählung der Eigenwerte von A {\displaystyle A} .

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

{ p 1 ′ ( x ) = 2 p 1 ( x ) p 1 ( 0 ) = 1 { p 2 ′ ( x ) = 2 p 2 ( x ) + p 1 p 2 ( 0 ) = 0 {\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 p 1 ( x ) = exp ⁡ ( 2 x ) {\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 p 2 ( x ) = x exp ⁡ ( 2 x ) {\displaystyle p_{2}(x)=x\exp(2x)} .

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

exp ⁡ ( A x ) = ∑ k = 1 2 p k ( x ) M k − 1 = exp ⁡ ( 2 x ) ( 1 0 0 1 ) + x exp ⁡ ( 2 x ) ( 1 − 1 1 − 1 ) = exp ⁡ ( 2 x ) ( 1 + x − x x 1 − x ) {\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 k ∈ { 2 , … , n } {\displaystyle 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 λ 1 t {\displaystyle p_{1}(t)=e^{\lambda _{1}t}}

p k ( t ) = e λ k t ∫ x = 0 t p k − 1 ( x ) e − λ k x d x {\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 ( λ = λ 1 = λ 2 = … = λ n ) {\displaystyle (\lambda =\lambda _{1}=\lambda _{2}=\dotsc =\lambda _{n})} , folgt aus der integralen Darstellung für p k {\displaystyle p_{k}} mit k ≥ 2 {\displaystyle k\geq 2} , dass die Funktion f = 1 {\displaystyle f=1} gerade ( k − 1 ) {\displaystyle (k-1)} -mal integriert werden muss. Für k ≥ 1 {\displaystyle k\geq 1} gilt somit:

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

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

exp ⁡ ( A t ) = e λ t E + t e λ t ( A − λ E ) + … = e λ t ∑ i = 0 n − 1 t i i ! ( A − λ E ) i {\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 ( λ i ≠ λ j {\displaystyle \lambda _{i}\neq \lambda _{j}} für i ≠ j {\displaystyle i\neq j} ). In diesem Fall gilt für die Diskriminante des charakteristischen Polynoms D n ≠ 0 {\displaystyle D_{n}\neq 0} . Die Lösung für p 1 ( t ) {\displaystyle p_{1}(t)} bleibt davon unverändert. Ab k ≥ 2 {\displaystyle k\geq 2} folgt nach Integration:

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

p 3 ( t ) = e λ 1 t − e λ 3 t ( λ 1 − λ 2 ) ( λ 1 − λ 3 ) + e λ 2 t − e λ 3 t ( λ 2 − λ 1 ) ( λ 2 − λ 3 ) {\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

p k ( t ) = ∑ i = 1 k − 1 e λ i t − e λ k t ∏ j = 1 k ( λ i − λ j )       {\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 j ≠ i {\displaystyle j\neq i} .

Für das Exponential einer ( n × n ) {\displaystyle (n\times n)} -Matrix mit verschiedenen Eigenwerten folgt damit die Newton-Interpolation^[2]

exp ⁡ ( A t ) = e λ 1 t E + ∑ i = 2 n [ λ 1 , … , λ i ] ∏ j = 1 i − 1 ( A − λ j E ) {\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 x {\displaystyle x} abhängigen Funktionen können dabei rekursiv berechnet werden mit

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

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

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

exp ⁡ ( A t ) = V diag ⁡ ( e λ i t ) V − 1 = ∑ j = 1 n e λ i t v j y j T {\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 y j T {\displaystyle y_{j}^{T}} die j {\displaystyle j} -te Zeile von V − 1 {\displaystyle V^{-1}} ist, mit der Lagrange-Interpolation

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

wobei

A j = ∏ k = 1 n A − λ k E λ j − λ k ( k ≠ j ) {\displaystyle 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 {\displaystyle A_{j}=v_{j}y_{j}^{T}}

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