This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

Matrix Exponentials

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;

Method 1: Scaling and Squaring

expmdemo1 is an implementation of algorithm 11.3.1 in the book:

Golub, Gene H. and Charles Van Loan. Matrix Computations, 3rd edition. Balitmore, 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

Method 2: Taylor Series

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

is the identity matrix with the same dimensions as . 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

Method 3: Eigenvalues and Eigenvectors

expmdemo3 assumes that the matrix has a full set of eigenvectors such that . The matrix exponential can be calculated by exponentiating the diagonal matrix of eigenvalues:

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

Compare Results

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

E = expm(Asave);
err1 = E - E1
err1 = 3×3
10-14 ×

    0.3553    0.1776    0.0888
    0.0888    0.1332   -0.0444
         0         0   -0.2665

err2 = E - E2
err2 = 3×3
10-14 ×

         0         0   -0.1776
   -0.0444         0   -0.0888
    0.1776         0    0.0888

err3 = E - E3
err3 = 3×3
10-13 ×

   -0.0711   -0.0444   -0.0799
   -0.0622   -0.0488   -0.0933
   -0.0711   -0.0533   -0.1066

Taylor Series Failure

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×2
106 ×

   -1.1985   -0.5908
   -2.7438   -2.0442

E3 = expmdemo3(A)
E3 = 2×2

   -0.0996    0.0747
   -0.1991    0.1494

Eigenvalues and Eigenvectors Failure

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

See Also