Cleaner algorithm

The Putzer algorithm is a recursive method for calculating the matrix exponential . The aim of the procedure is to determine the matrix exponential for a given and . ${\displaystyle A\in \mathbb {C} ^{n\times n}}$ ${\displaystyle x\in \mathbb {R} }$ ${\displaystyle \textstyle \exp(Ax):=\sum _{k=0}^{\infty }{\frac {(Ax)^{k}}{k!}}}$

Like the matrix exponential, the algorithm plays a role primarily in the theory of ordinary differential equations , since it can be used to determine solutions to systems of linear differential equations .

The procedure

Let be a square matrix. Let us also count the eigenvalues ​​of the matrix taking into account the algebraic multiplicity. This enumeration exists because it is algebraically closed . ${\displaystyle A\in \mathbb {C} ^{n\times n}}$ ${\displaystyle (n\times n)}$ ${\displaystyle \{\lambda _{j}\}_{j=1,\ldots ,n}}$ ${\displaystyle A}$ ${\displaystyle \mathbb {C} }$

For we now define recursively continuously differentiable functions and quadratic matrices , so that the following relationship applies: ${\displaystyle k\in \{1,\ldots ,n\}}$ ${\displaystyle p_{k}\in C^{1}(\mathbb {R} ,\mathbb {C} )}$ ${\displaystyle (n\times n)}$ ${\displaystyle M_{k-1}\in \mathbb {C} ^{n\times n}}$ ${\displaystyle x\in \mathbb {R} }$

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

The recursion rule for the matrices is: ${\displaystyle M_{k}}$

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

They are defined recursively by the following initial value problems : ${\displaystyle p_{k}}$

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

example

Illustration of the procedure for : Let the following matrix be given. We will now calculate the matrix exponential for anything using the Putzer algorithm . First we determine the eigenvalues ​​of the matrix , taking into account the algebraic multiplicity. To do this, we calculate the characteristic polynomial and equate it : ${\displaystyle n=2}$ ${\displaystyle A={\begin{pmatrix}3&{-1}\\1&{1}\end{pmatrix}}}$ ${\displaystyle \exp(Ax)}$ ${\displaystyle x\in \mathbb {R} }$ ${\displaystyle A}$ ${\displaystyle 0}$

${\displaystyle \chi _{A}(\lambda ):=\det(A-\lambda I)=\det {\begin{pmatrix}3&{-1}\\1&{1}\end{pmatrix}}=\lambda ^{2}-4\lambda +4=(\lambda -2)^{2}\,{\overset {!}{=}}\,0}$

It follows that the only eigenvalue of the matrix is, but with algebraic multiplicity . Thus it is a matter of a counting of the eigenvalues ​​of . ${\displaystyle \lambda =2}$ ${\displaystyle A}$ ${\displaystyle 2}$ ${\displaystyle (\lambda _{1},\lambda _{2})=(2,2)}$ ${\displaystyle A}$

Next, we determine the matrices using the recursion formulas from above . It follows and accordingly . ${\displaystyle M_{k},k\in \{0,1\}}$ ${\displaystyle M_{0}={\begin{pmatrix}1&{0}\\0&{1}\end{pmatrix}}}$ ${\displaystyle M_{1}=(A-\lambda _{1}I_{2})M_{0}=(A-\lambda _{1}I_{2})={\begin{pmatrix}1&{-1}\\1&{-1}\end{pmatrix}}}$

Finally we calculate for the functions that are defined by the following initial value problems: ${\displaystyle k\in \{1,2\}}$ ${\displaystyle p_{k}}$

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

The first initial value problem can easily be determined as using the solution method separation of the variables for ordinary differential equations . The second initial value problem can also be calculated very easily using the method variation of the constants and obtained as a solution . ${\displaystyle p_{1}(x)=\exp(2x)}$ ${\displaystyle p_{2}(x)=x\exp(2x)}$

Now that we have everything together, we can state for one : ${\displaystyle \exp(Ax)}$ ${\displaystyle x\in \mathbb {R} }$

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

Special solutions

As shown in the example, the first initial value problem can be determined using the solution method Separation of Variables for ordinary differential equations. The other initial value problems with can also be easily calculated using the method variation of the constants. Accordingly, the solutions follow first: ${\displaystyle k\in \{2,\ldots ,n\}}$

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

From this, in turn, other special solutions can be derived depending on the eigenvalues.

Same eigenvalues

If all eigenvalues equal , follows from the integral representation for having the function just be integrated even have. The following applies to: ${\displaystyle (\lambda =\lambda _{1}=\lambda _{2}=\dotsc =\lambda _{n})}$ ${\displaystyle p_{k}}$ ${\displaystyle k\geq 2}$ ${\displaystyle f=1}$ ${\displaystyle (k-1)}$ ${\displaystyle k\geq 1}$

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

The matrix exponential of a matrix with the same eigenvalues can finally be determined explicitly with the following formula: ${\displaystyle (n\times n)}$ ${\displaystyle \lambda }$

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

Unequal eigenvalues

Often all eigenvalues ​​of the matrix are not the same . In this case applies to the discriminant of the characteristic polynomial . The solution for remains unchanged. From after integration follows: ${\displaystyle (\lambda _{1}\neq \lambda _{2}\neq \dotsc \neq \lambda _{n})}$ ${\displaystyle D_{n}\neq 0}$ ${\displaystyle p_{1}(t)}$ ${\displaystyle k\geq 2}$

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

etc.

The solution here obviously follows the system:

${\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})}}~~~}$ With ${\displaystyle j\neq i}$

The Newton interpolation follows for the exponential of a matrix with unequal eigenvalues ${\displaystyle (n\times n)}$

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

for the representation of the matrix exponential function as a polynomial. The functions that depend on can be calculated recursively with: ${\displaystyle x}$

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

Effort of polynomial solution methods

The Putzer algorithm has an expense of the order of magnitude (essentially matrix multiplications ). This means that the effort is significantly higher than when using methods based on the Taylor series or based on matrix decomposition methods (such as diagonalization or QR algorithm ), each with an effort of the order of magnitude . The algorithm is therefore only suitable for smaller matrices, but a completely symbolic solution is advantageously obtained here . ${\displaystyle {\mathcal {O}}(n^{4})}$ ${\displaystyle n}$ ${\displaystyle {\mathcal {O}}(n^{3})}$

If one compares the complexity of the solution for the matrix exponential function by diagonalizing the 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}}$,

here is the -th line of , with the Lagrange interpolation ${\displaystyle y_{j}^{T}}$ ${\displaystyle j}$ ${\displaystyle V^{-1}}$

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

in which

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

is, then it follows

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

and the effort of the order of magnitude to calculate is completely unnecessary. ${\displaystyle {\mathcal {O}}(n^{4})}$ ${\displaystyle A_{j}}$

literature

• Paul Waltman: Differential Equations and Dynamical Systems , Dover Publications, 2004, ISBN 978-0486434780 , pp. 49-60

Individual proof

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 . tape 45 , no. 1 , 2003, ISSN  1095-7200 , p. 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 . tape 45 , no. 1 , 2003, ISSN  1095-7200 , p. 22-23 , doi : 10.1137 / S00361445024180 ( cornell.edu [PDF]).