Accelerating the pace of engineering and science

# Documentation Center

• Trials
• Product Updates

## Matrix Exponentials

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/

Start with the Matrix A

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

```

Scaling and Squaring

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

```

Matrix Exponential via Taylor Series

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

```

Matrix Exponential via Eigenvalues and Eigenvectors

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

```

Compare the Results

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

```

Taylor Series Fails

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

```

Eigenvalues and Vectors Fails

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?