Description
MACM 316 – Computing Assignment 1 CA1 – Floating Point Errors and Matlab Plotting
This assignment is a more extensive investigation of the rounding error example we studied in class using
the Matlab code roundex.m (posted on Canvas). The polynomial function (x−a)
n with a any real number
and n a positive integer can be written in expanded form as
(x − a)
n =
Xn
k=0
n
k
x
n−k
(−a)
k = x
n −
n
1
x
n−1
a +
n
2
x
n−2
a
2 −
n
3
x
n−3
a
3 + … +
n
n
(−a)
n
(∗)
where
n
k
=
n!
(n−k)! k!
are binomial coefficients.
1. Plot f(x) = (x − 2)n
for n = 1, 2, 3, 4, 5, 6 on the domain x ∈ [0, 4], using the factored (unexpanded)
form. Display your 6 curves together on a single plot, and use different colors, line styles and/or
a legend to clearly identify the various curves. Choose appropriate limits for the y-axis that ensure
important features of the function are visible. Consider this your “exact result”.
2. Write a Matlab function in an m-file named fexpand.m that computes the polynomial in the expanded form (∗). The first line of your code should be this function definition statement:
function fx = fexpand(a, n, x)
The function has 3 input parameters (a, n and x) and returns a single output argument (the computed
value of f(x)). You can make use of the built-in Matlab function nchoosek to compute the binomial
coefficients. Test your code completely, making sure that it is able to exit gracefully with a suitable
error/warning message for any values of the input arguments that might generate an invalid result.
3. Produce 6 plots depicting your “expanded” f(x) curves that “zoom in” near the point x = 2 using a
series of successively smaller x-intervals, x ∈ [2 − δ, 2 + δ] for δ = 0.5, 0.1, 0.05, 0.025, 0.01, 0.005†
.
In
your report, you should only include the plots for δ = 0.5, 0.05, 0.005, with two select values of the
polynomial degree n, for a total of 3 × 2 = 6 plots – and choose the two n values that you think best
illustrate your conclusions! Describe any interesting behaviour or differences that arise between the
various n, δ values.
In addition, be sure to comment on the following:
(a) Identify the smallest value of the exponent n for which your plots of the expanded polynomial
(∗) differ significantly from the “exact” plots (in part 1).
(b) Describe what happens as n increases (it might help to try a few n even larger than 6).
(c) Explain why performing these computations in double-precision floating point arithmetic gives
rise to the errors you observe.
Be specific!
†Note that you should choose your plotting points to be consistent with the size of each x-interval! In other words, as δ gets
smaller, you should also decrease the spacing between plotting points so that you can clearly resolve any interesting features
MACM 316 – Computing Assignment 2 CA2 – Explorations in Root-Finding
CA2 – Explorations in Root-Finding
In this assignment, you will apply a combination of analytical techniques and root-finding algorithms to
study the roots of the nonlinear equation
f(x) = cos(x) + 1
1 + e−2x
.
(a) Start by plotting f(x) on the interval x ∈ [−6π, 6π] and describe the overall behaviour of the function,
as well as the number and rough location of its roots. Use the “zoom” feature of Matlab’s plotting
window (or change the axis limits with set or axis commands) to make sure that you identify all
roots – you may have to increase your plotting point density in order to view f in sufficient detail!
(b) Investigate analytically what happens to the function for large |x| by computing the following two
limits of the exponential term in f(x):
lim x→−∞
1
1 + e−2x
and lim x→+∞
1
1 + e−2x
.
Use these results to determine two simpler “limit functions” that approximate the negative and positive “halves” of f(x):
• f−(x) for x < 0, which approximates f for large negative values of x,
• f+(x) for x > 0, which approximates f for large positive values of x.
Plot f±(x) on their respective half-intervals along with the original function f(x). Then derive analytically the exact values for all roots of f−(x) (for x < 0) and f+(x) (for x > 0). Add these roots to
your plot and comment on how well they seem to approximate the actual zeros of f(x).
(c) Next, apply the bisection method to compute to within an absolute error tolerance of 10−6
the smallest† negative root of f(x) (call it x
∗
).
Use the bisect2.m code from lectures and choose an initial
bracket that is motivated by your plot from part (b). Using an appropriate error measure, compare
your computed root x
∗
to the smallest root of f−(x) (call it x−) that you determined analytically in
part (b). How well does x− approximate your bisection result x
∗
?
‡
Use bisection to compute the next two negative roots of f(x) and compare them with the corresponding roots of f−(x). What do you notice about the accuracy of your approximations?
†By “smallest” I mean the smallest in magnitude or closest to zero.
‡Here I want you to think of the bisection result as your “exact” root, and the analytical result as the “approximate” root.
(d) Repeat part (c) to compute the four smallest positive roots of f(x) and compare them with roots
of f+(x). Explain how you choose an appropriate initial bracket for each root, and compare your
computed roots with the corresponding estimates x+ from f+(x).
(e) Consider the following fixed point iteration§
:
xk+1 = g(xk) = arccos
−1
1 + e−2xk
. (∗)
Show that finding a fixed point of g(x) is equivalent to finding a root of f(x) = 0. Use the code
fixedpt.m from lectures to approximate two roots:
• the smallest negative root from part (c) using an initial guess of x0 = −1.5, and
• the smallest positive root from part (d) using an initial guess of x0 = 3.0.
Compare with the results you obtained using the bisection method, making sure that you apply the
same stopping criterion for x – note that this may require modifying the fixedpt code!
Describe any unexpected behaviour you observe. Can you explain the convergence of your fixed
point iterations (perhaps with the help of a sketch)? What is the real problem with using the fixed
point function in (∗) to try to compute all roots of f(x)?
§The Matlab built-in function for the “arccos” or inverse cosine function is called acos.
MACM 316 – Computing Assignment 3 CA3 – Gaussian Elimination for Random Matrices
In this computing assignment, you will study the behaviour of floating-point errors when the Gaussian elimination (GE) algorithm with partial pivoting (as implemented in Matlab’s “backslash” command) is applied to random matrices of widely differing sizes.
You have already seen in class that
when solving systems with small matrices, the errors remain small. However, even for state-of-theart implementations of GE, these errors can become significant when the matrix is very large. The
purpose of CA3 is to answer the question: How large is “large”?
To construct a random matrix system for which the exact solution (or at least a pretty good approximation of it) is known, we’ll make use of a simple trick:
• Construct the N × N matrix A with random entries chosen from the interval [−1, 1].
• Let x = [1, 1, . . . , 1]T be an N-vector of ones – this is our exact solution.
• Compute the right hand side corresponding to any particular choice of A to be b = Ax.
Then, we compute the approximate solution y = A \ b using Matlab’s backslash operator. To measure
the error between x and y, let
δ = max
i=1,…,N
|xi − yi
|
which is the maximum component-wise difference between the two N-vectors†
.
Because A is a random matrix, no single A gives a good measure of the actual error in GE; instead,
we must run this calculation for a number of trials with different realizations of A in order to get a
sense for the representative error.
For this purpose, let M be the number of trials and suppose that the
solution error for the k
th trial is δ
(k)
. Then we can define the mean or average error as
EN =
1
M
δ
(1) + δ
(2) + · · · + δ
(M)
which should be a better measure of error for an “average” random matrix A.
To help you get started, I have posted the simple code GERandom.m on Canvas, which performs
the computations just described for given values of M and N. The code determines the mean error
EN over M random trials and then generates a scatter plot showing the error in each trial.
Note that
this plot uses a log scale for error (using Matlab’s semilogy command) so that large changes in the
error magnitude are easier to see. Experiment with this code and observe how increasing N leads to
an increase in the mean error EN .
†You’ll see later in section 3b that this is just the vector max-norm or inf-norm, denoted δ = kx − yk∞.
This leads us to the goal of CA3 which is . . .
GOAL: To determine the size of the matrix N∗
for which the mean error in Gaussian elimination reaches EN∗ ≈ 1. In other words, how large does a random
matrix A have to be before the round-off errors in GE are the same order of magnitude as the solution x? Or, when does relative error exceed 100%!
In practice, your computer will not have the processing power to find N∗ directly because these
matrices are simply too huge. Instead, you should compute the mean error EN for a sequence of more
“moderate” sized N values and plot the data. You will then observe that your data points lie along
a predictable curve, which you may then extrapolate to even larger values of N to estimate roughly
where the curve hits the value EN ≈ 1.
Here is a detailed list of what is required in your 1-page report:
(a) Modify the code GERandom.m to turn it into a function with this calling sequence
function EN = GERandom(N, M)
which should make your later calculations much easier. You can also remove any unneeded
output or plotting statements.
(b) Use your modified code to experiment with different N and M, choosing values that are as large
as possible without forcing your code to take too long to execute (limit yourself to 5–10 minutes
of clock time at most). You need to strike a balance: you want your number of trials M to be
large enough to get an accurate estimate of the average error, but you also want N as large as
possible to obtain an accurate extrapolation. Carefully justify your final choice of N and M.
(c) Plot your EN versus N for different values of N ranging between 5 and the maximum value you
chose in (b). Use a log-log scale, which in Matlab can be done in two ways:
loglog(N, EN) or plot(log10(N), log10(EN))
where N and EN are vectors containing your computed data points.
(d) Observe that your data lies roughly along a straight line! This suggests seeking a linear function
that fits your log-log data, or in other words to find constants p1 and p2 so that log10(EN ) ≈
p1 log10(N) + p2. You can do this manually or else run the following Matlab command
p = polyfit(log10(N), log10(EN), 1)
which fits a polynomial of degree 1 to your log-log data, returning a 2-vector p that contains the
coefficients of the linear fit‡
. Add this line to your data plot and comment on how accurately it
approximates the data.
(e) Estimate N∗ by extrapolating your linear function to the point where it hits the value EN = 1 (or
log10 EN = 0).
‡You’ll learn more in section 4c about the algorithm implemented in polyfit to perform this curve fitting.
MACM 316 – Computing Assignment 4 CA4 – Approximating matrix exponentials
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?”
MACM 316 – Computing Assignment 5 CA5 – Parametric splines for general curves
In this assignment, your aim is to extend the definition of a cubic spline to interpolate data that cannot
be described by a single-valued function.
(a) As a first example, consider the data points listed in the table (below, left) which are sampled
from an underlying smooth curve pictured in the plot (below, right):
t x y
0 0.0 0.0
1 1.0 3.0
2 2.0 3.0
3 2.0 4.0
4 3.0 5.0
multivalued!
Clearly, the curve passes through multiple y-values at x = 2, meaning that the underlying curve
is multi-valued and doesn’t pass the “vertical line test” – so it’s not a function!! This means that
the usual cubic spline approach based on trying to interpolate a function of the form y = f(x)
will not work.
The trick to dealing with data like this is to use a parametric spline that is based
on a parametric description of the curve: x = R(t) and y = S(t) for some parameter t. For
this example, it’s easiest to assign parameter values t = 0, 1, 2, 3, 4 corresponding to the index of
the points in the table. The particular choice of t values is actually not important as long as t is
monotonically increasing. Because R(t) and S(t) are now both functions of t, we can interpolate
them with two separate cubic splines.
Based on the t-x-y data from the table, generate two cubic splines x = R(t) and y = S(t) using
the Matlab spline function with its default not-a-knot end-point conditions. Provide three
separate plots: of R versus t, S versus t, and S versus R (the last of which should reproduce the
plot above).
(b) Next use a parametric spline to interpolate the more interesting “four-leaf” curve pictured below.
The table lists coordinates (xi
, yi), i = 0, 1, . . . , 12, for 13 points lying along this curve.
t x y
0 2.75 -1.0
1 1.3 -0.75
2 -0.25 0.8
3 0.0 2.1
4 0.25 0.8
5 -1.3 -0.25
6 -2.5 0.0
7 -1.3 0.25
8 0.25 -1.3
9 0.0 -2.5
10 -0.25 -1.3
11 1.3 -0.25
12 2.75 -1.0
Determine the parametric spline x = R(t), y = S(t) that interpolates this set of points (x, y),
again using Matlab’s spline function. Plot your parametric spline curve and verify (by zooming in on your plot) that the right-most leaf is different from the other three leaves in that it is
not smooth – that is, the spline endpoints meet at a cusp.
(c) For a periodic curve like the one in part (b), it is more appropriate to use periodic end-point
conditions instead of not-a-knot conditions. That is, we should take R0
0
(t0) = R0
n−1
(tn) and
R00
0
(t0) = R00
n−1
(tn), and similarly for S(t). Use the Matlab code perspline.m posted on Canvas
to generate two periodic cubic spline approximations for R(t) and S(t), plot your parametric
curve, and compare to what you obtained in part (b). Verify that the cusp is eliminated.
(d) You can now get creative! Start by drawing your own periodic parametric curve, which can be
any smooth curve of your own design as long as the endpoints meet at the same location†
. Your
curve should be both . . .
• interesting: it should cross itself at least once, such as the 4 “leaf crossings” in part (b), and
• not too complicated: so that it can be approximated by a cubic spline with 20–40 points.
To generate your list of 20–40 data points, you may find it helpful to save your drawing as an
image file, and then use Matlab’s ginput (graphical input from mouse). The built-in function
ginput allows you to select points by clicking with the mouse at a sequence of along your parametric curve in the plotting window, and then outputs all x and y coordinates of the points you
selected.
If you type help ginput, you will see that it records mouse clicks within the plotting
window until the “enter” key is pressed, after which it returns two vectors of coordinates.
You
may find the following sequence of Matlab commands helpful:
figure(’position’, get(0,’screensize’)) % largest window possible
axes(’position’, [0 0 1 1]), axis square % make x,y-axes equal
imshow(’myimage.png’) % display your drawing on-screen
[x,y] = ginput; % record mouse clicks until ’Enter’
save mydatafile.mat x y % save x,y data points to a file
close % delete the huge window (if desired)
The save command saves your data points to a file that can later on be read back in using the
load command. You may then use your (x, y) data as input to perspline.m in the same way
you did in part (c) to construct your parametric spline and plot your results.
†
If you are looking for artistic inspiration, you can find a host of examples by entering “one line drawings” in a Google
Images search. Beware that many of these images are too complex for this assignment, so choose wisely! (or simplify).
