Description
Recall from calculus that the exponential of a real number x has the Taylor series expansion
e
x = exp(x) = 1 + x +
1
2!x
2 +
1
3!x
3 + · · ·
If you have any experience with complex variables like z = x + iy, then you may have encountered
the complex exponential e
z = e
x
e
iy which can also be represented with this Taylor series.
But it gets
even better! If A is an n × n matrix, then the exponential of A can be defined analogously as
e
A = exp(A) = I + A +
1
2!A
2 +
1
3!A
3 + · · · (∗)
where I is the n × n identity matrix and Ak
is represents the matrix product A · A · · · A | {z }
k times
.
Matrix exponentials are useful in a number of problems in computational mathematics, especially for the solution
of systems of differential equations.
WARNING: The exponential of a matrix is not equal to the exponential of its entries!! For example:
exp 0 1
1 0 =
1.5430806348152 1.1752011936438
1.1752011936438 1.5430806348152
6=
e
0
e
1
e
1
e
0
=
1 e
e 1
.
In this CA, you will implement a simple algorithm for approximating the matrix exponential by
summing the first k terms in the infinite series (∗). The aim is to investigate the accuracy, efficiency
and robustness of this algorithm.
You will test your code on a special matrix A that is stored in one of
four files named CA4matrix#.mat†
, where #=1, 2, 3 or 4. Select one of the matrix files from Canvas
and load it into Matlab using
load(’CA4matrix#.mat’);
which assigns your chosen matrix to the variable A. Look at the entries of this matrix and note that it
is a full, 500 × 500 matrix with complex entries.
†These files are in “mat-file” format, which provides a simple and efficient method for storing and accessing variables
saved in a previous Matlab session.
(a) Write a Matlab code that implements the truncated series approximation for e
A. Test your code
on your chosen matrix A with k = 50 and k = 150 terms in the series, and assign the result to
the variable expAk.
You can visualize the size of the entries in your approximation of the matrix
exponential using the following two commands:
imagesc(real(expAk))
colormap gray
This generates a gray-scale plot with each pixel depicting the size of a matrix entry in e
A by its
colour: large entries are white, small entries are black. Describe the differences between your
two images in terms of “image quality”.
HINT: When writing your code, you should make use of the following formula for the k
th term
in the series:
1
k!
A
k =
1
k
A ·
1
(k − 1)!A
k−1
This simple trick can be used to increase the efficiency of your code by reducing the number of
matrix multiplications required in your truncated series approximation!!
(b) Run your algorithm again using a wider range of k values with k = 5, 10, 15, . . . , 150. Determine
the execution time required in each case using Matlab’s tic and toc commands, which are
invoked as follows:
tic;
{{ my code }}
cpu_time = toc;
The variable cpu time is the elapsed CPU time required to execute the code between the tic/toc
commands.Plot your CPU times versus k and describe how they depend on k. Estimate the number of floating point operations or “flops” required by your algorithm as a function of k, and then
explain whether or not your flop estimate is consistent with your plot.
(c) Matlab has a built-in command called expm for computing the matrix exponential exp(A):
expAexact = expm(A);
First, determine the CPU time required for this expm calculation and compare it with the timings
for your algorithm from part (b).
If you then treat this as your “exact result,” you can estimate the error in your truncated series
approximation using the matrix 2-norm:
err = norm(expAexact-expAk, 2);
Plot the error against k for the same range of k values as in part (b), this time using a logarithmic
scale for the error. How does the error depend on k?
(d) Complete your report with a summary discussion of how your algorithm compares to Matlab’s
expm in terms of the following three criteria:
• accuracy: measured by the values of err,
• efficiency: comparing CPU time for both methods,
• robustness: in terms of which method is most “reliable”.‡
‡Go back to the lectures notes for Section 1a, where I define robustness, which addresses questions such as “Is the
algorithm guaranteed to return a ‘good’ result?” and “How strongly do accuracy and efficiency depend on the inputs?”
