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

For background on the computation of matrix exponentials, see "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later," SIAM Review 45, 3-49, 2003.

The Pseudospectra Gateway is also highly recommended. The web site is:

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

A = 0 1.0000 2.0000 0.5000 0 1.0000 2.0000 1.0000 0

`expmdemo1`

is an implementation of algorithm 11.3.1 in Golub and Van Loan, Matrix Computations, 3rd edition.

% 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 = 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. As a practical numerical method, this 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 = 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. 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 = 5.3091 4.0012 5.5778 2.8088 2.8845 3.1930 5.1737 4.0012 5.7132

For this matrix, they all do equally well.

E = expm(Asave); err1 = E - E1 err2 = E - E2 err3 = E - E3

err1 = 1.0e-14 * 0.3553 0.1776 0.0888 0.0888 0.1332 -0.0444 0 0 -0.2665 err2 = 1.0e-14 * 0 0 -0.1776 -0.0444 0 -0.0888 0.1776 0 0.0888 err3 = 1.0e-13 * -0.0799 -0.0444 -0.0622 -0.0622 -0.0488 -0.0933 -0.0711 -0.0533 -0.1066

Here is a matrix where 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) E2 = expmdemo2(A) E3 = expmdemo3(A)

E1 = -0.0996 0.0747 -0.1991 0.1494 E2 = 1.0e+06 * -1.1985 -0.5908 -2.7438 -2.0442 E3 = -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) E2 = expmdemo2(A) E3 = expmdemo3(A)

E1 = 0.3679 0.3679 0 0.3679 E2 = 0.3679 0.3679 0 0.3679 E3 = 0.3679 0 0 0.3679

Was this topic helpful?