This example shows 3 of the 19 ways to compute the exponential of a matrix.

For background on the computation of matrix exponentials, see:

Moler, C. and C. Van Loan. "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later." *SIAM Review.* Vol. 45, Number 1, 2003, pp. 3-49.

Another recommended resource is the Pseudospectra Gateway website.

Start by creating a matrix `A`

.

A = [0 1 2; 0.5 0 1; 2 1 0]

`A = `*3×3*
0 1.0000 2.0000
0.5000 0 1.0000
2.0000 1.0000 0

Asave = A;

`expmdemo1`

is an implementation of algorithm 11.3.1 in the book:

Golub, Gene H. and Charles Van Loan. *Matrix Computations, 3rd edition.* Baltimore, MD: Johns Hopkins University Press, 1996.

% Scale A by power of 2 so that its norm is < 1/2 . [f,e] = log2(norm(A,'inf')); s = max(0,e+1); A = A/2^s; % Pade approximation for exp(A) X = A; c = 1/2; E = eye(size(A)) + c*A; D = eye(size(A)) - c*A; q = 6; p = 1; for k = 2:q c = c * (q-k+1) / (k*(2*q-k+1)); X = A*X; cX = c*X; E = E + cX; if p D = D + cX; else D = D - cX; end p = ~p; end E = D\E; % Undo scaling by repeated squaring for k = 1:s E = E*E; end E1 = E

`E1 = `*3×3*
5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132

`expmdemo2`

uses the classic definition for the matrix exponential given by the power series

$${e}^{A}=\sum _{k=0}^{\infty}\frac{1}{k!}{A}^{k}.$$

${\mathit{A}}^{0}$ is the identity matrix with the same dimensions as $\mathit{A}$. As a practical numerical method, this approach is slow and inaccurate if `norm(A)`

is too large.

A = Asave; % Taylor series for exp(A) E = zeros(size(A)); F = eye(size(A)); k = 1; while norm(E+F-E,1) > 0 E = E + F; F = A*F/k; k = k+1; end E2 = E

`E2 = `*3×3*
5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132

`expmdemo3`

assumes that the matrix has a full set of eigenvectors $\mathit{V}$ such that $\mathit{A}=\mathrm{VD}{\mathit{V}}^{-1}$. The matrix exponential can be calculated by exponentiating the diagonal matrix of eigenvalues:

$${e}^{A}=V{e}^{D}{V}^{-1}.$$

As a practical numerical method, the accuracy is determined by the condition of the eigenvector matrix.

A = Asave; [V,D] = eig(A); E = V * diag(exp(diag(D))) / V; E3 = E

`E3 = `*3×3*
5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132

For the matrix in this example, all three methods work equally well.

E = expm(Asave); err1 = E - E1

err1 =3×310^{-14}× 0.3553 0.1776 0.0888 0.0888 0.1332 -0.0444 0 0 -0.2665

err2 = E - E2

err2 =3×310^{-14}× 0 0 -0.1776 -0.0444 0 -0.0888 0.1776 0 0.0888

err3 = E - E3

err3 =3×310^{-14}× -0.7105 -0.5329 -0.7105 -0.6661 -0.5773 -0.8882 -0.7105 -0.7105 -0.9770

For some matrices the terms in the Taylor series become very large before they go to zero. Consequently, `expmdemo2`

fails.

A = [-147 72; -192 93]; E1 = expmdemo1(A)

`E1 = `*2×2*
-0.0996 0.0747
-0.1991 0.1494

E2 = expmdemo2(A)

E2 =2×210^{6}× -1.1985 -0.5908 -2.7438 -2.0442

E3 = expmdemo3(A)

`E3 = `*2×2*
-0.0996 0.0747
-0.1991 0.1494

Here is a matrix that does not have a full set of eigenvectors. Consequently, `expmdemo3`

fails.

A = [-1 1; 0 -1]; E1 = expmdemo1(A)

`E1 = `*2×2*
0.3679 0.3679
0 0.3679

E2 = expmdemo2(A)

`E2 = `*2×2*
0.3679 0.3679
0 0.3679

E3 = expmdemo3(A)

`E3 = `*2×2*
0.3679 0
0 0.3679