Note: This page has been translated by MathWorks. Please click here

To view all translated materials including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materials including this page, select Japan from the country navigator on the bottom of this page.

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:

http://web.comlab.ox.ac.uk/projects/pseudospectra/

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

```
A =
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 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

```
err1 =
1.0e-14 *
0.3553 0.1776 0.0888
0.0888 0.1332 -0.0444
0 0 -0.2665
```

err2 = E - E2

```
err2 =
1.0e-14 *
0 0 -0.1776
-0.0444 0 -0.0888
0.1776 0 0.0888
```

err3 = E - E3

```
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)

```
E1 =
-0.0996 0.0747
-0.1991 0.1494
```

E2 = expmdemo2(A)

```
E2 =
1.0e+06 *
-1.1985 -0.5908
-2.7438 -2.0442
```

E3 = expmdemo3(A)

```
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)

```
E1 =
0.3679 0.3679
0 0.3679
```

E2 = expmdemo2(A)

```
E2 =
0.3679 0.3679
0 0.3679
```

E3 = expmdemo3(A)

```
E3 =
0.3679 0
0 0.3679
```

Was this topic helpful?