• Please read the Guidelines for Assignments first.

• Submit a one-page PDF report to Canvas and upload you Matlab scripts (as m-files). Do not use any other file formats.

• Keep in mind that Canvas discussions are open forums.

• Acknowledge any collaborations and assistance from colleagues, TAs, instructors etc. Computing the matrix exponential

Recall from calculus that the exponential of a real number x has the Taylor series expansion

exp(A) = 1 +x+x²1/2! + x³1/3! + x⁴1/4!+.....

Let A be an 7/ x n matrix. The exponential of A is defined as

exp(A) = 1 +A +A²1/2! + A³1/3! + A⁴1/4!+.....

where I is the n x n identity matrix. Matrix exponentials are useful in a number of problems in computational mathematics, including the solution of systems of differential equations.

Note: the exponential of a matrix is not equal to the exponential of its entries. For example:

A simple algorithm for computing (an approximation to) the matrix exponential involves summing the first k terms of the infinite series (1). Your goal in this assignment to investigate the accuracy, efficiency and robustness of this algorithm. You will be testing on a particular matrix A which is contained in the file CA5matrix.mat. Download this file from Canvas and load it into Matlab using the following command:

1. load('CA3matrix.mat');

This loads the matrix as a 1000 x 1000 array A.

(a) Let k = 50 and k = 150. Run your algorithm to get the output, which you should call expAk. Plot the result using the following commands:

1. imagesc (real (expAk) ) ;

2 colormap gray

(b) Run your algorithm for k = 10, 20, 30, ... ,150 and use the tic and toc commands (recall the in-class demo TicToe.m) to plot the computational time versus k. How does the computational time appear to depend on k? Does this agree with what you would expect from counting the number of flops in your algorithm? Explain.

(c) Matlab has a built-in function for computing the matrix exponential exp(A):

1 expA = expm(A);

Using this, compute the error of your algorithm:

1 err = norm(epr - eprk);

Note: if you don't know What a norm is, don't worry. It's just a way of measuring the error.

Now plot the error against 1: for the same range as in part (b) (you might want to use a log scale for the y-axis). How does the error behave? "What can you say about the accuracy of your algorithm?

What can you say about the robustness of your algorithm?