ENGG1811 Lab 8: Using Matlab
After completing this lab, students should be able to
- Use Matlab as a calculator working with arrays in one and two
- Understand how to plot variables in 2 dimensions
- Perform array operations
- Use the comet function to aid in visualising function movement
This lab consists of 4 parts. You should use a script to perform the work
in each part. They will check your Part C script fully, and will withhold
some marks if the script is not properly documented. You will also need to
show your tutor an example from Part D.
The on-line assessment is about matrix operations, and can be done at
Organising your work
After starting up matlab, you will be working in the default folder
MATLAB, We recommend that you organise your work by creating folders under
MATLAB. Create a folder called ENGG1811. Change into ENGG1811 folder and
create a folder called lab08. Change into the directory lab08 and you can
store all your files for this week in this directory.
Part A: Arrays
There are many steps in this part, but most are very simple, almost
trivial. Don't add a semicolon to the answers, we want the results to
appear on screen.
Example: typing the Matlab command
id = [1, 0; 0 1]
is one way of creating the identity matrix of order 2 and displaying the
- Create a 2x3 array m whose first row contains the values 1, 2
and 4 and whose second row contains 3, -7 and 16.
- Confirm how big m is with the size and length functions.
- Assign the transpose of m to mt. You may use either '
and .' for transpose because the array contains only real number, but
don't forget that ' is transpose and complex conjugate, while .' is
- We'll do matrix multiplication more seriously later on, but type
these expressions and interpret the results in terms of the shapes of
the two operands and that of the result (if there is any):
m * mt
mt * m
m .* mt
mt .* m
- (Chapman exercise 5.38, abridged). Component failure rates are
expressed as mean time between failures (MTBF). A system that relies on
several components will have an overall MTBF that's a function of the
individual component MTBFs:
Given the following component MTBFs, calculate the system MTBF in hours:
You should use just two Matlab statements: one to assign the four
numbers to a vector, and one to perform the calculation using suitable
array operations on the vector. You will find the sum function handy for
this question. You can learn more about the sum function by looking up
Part B: Plotting
- Create a vector x that has values from -2 to 2, with each
value differing from the previous by 0.01. Use the colon operator.
- Repeat the exercise using linspace, assign to x2.
(Hint: You can use x to see how
many points you need.)
- The vectors are not necessarily identical (the isequal
function can compare them). The default display is to show 4 decimal
places. You can turn on maximum display precision with
format long g
which shows 14 decimal places. Compute the difference between the
second-last element of x and x2.
You may see a small non-zero value which shows there is a small
difference between them. We will discuss numerical errors later on in
the course .
- Consider the following two functions:
Calculate and plot them on the same figure. Use separate calls to the plot
function, using hold to overlay them (type "help hold"
if you don't know how to use that command).
- Plot can display multiple plots with additional pairs of
arguments. Close the figure window and try it. Add titles and labels to
each axis (don't worry about fancy formatting).
- The figure window has some useful features, including a zoom (the
magnifier-with-a-plus icon). Use it to examine the bottom of the cusp
- The plots so far have been pure functions, with a single value of y
for each x. This time plot a circle of unit radius using the
Start with 200 points for theta, from 0 to 2π and create x
and y vectors from that.
Hint: to ensure you get a circle rather than an ellipse, before
the plot command type
clf; axis equal; hold on;
- Plot a circle on a polar plot using the polar function. If
you are not sure how to use it, look at the documentation. The function
polar accepts angle and
magnitude vectors, and the magnitude of every point on a unit circle is,
...? You can produce a really compact solution with the ones and
- Unit magnitude is boring. What would happen if the magnitude were a
monotonic function of θ, say a proportion or square root? Test your
hypothesis with θ ranging from 0 to 6π
Part C: Van der Waals Equation
The relationship between the volume V
of a vessel holding n moles a certain gas
and pressure P in the gas at a certain
absolute temperature T is approximated by
van der Waals equation, a modification of the ideal gas law:
where R is the universal gas constant (You
can use 8.314 J/K/mol, which is not the most accurate version, for this
lab) and the two values a
(representing attraction between particles) and b
(representing the particle volume) are characteristic of the particular
In derived SI units, for oxygen, a =
138.2 kPa (L/mol)2 and b =
The P-V curve has a different shape depending on a critical temperature
and associated volume:
Below the critical temperature the substance is in a mixed gas and liquid
Write a Matlab script that has the following characteristics:
- It is properly documented, including the definitions of constants,
problem parameters, variables, and separate code blocks with a short
- It uses problem parameters to record the oxygen a
and b values.
- It asks for the number of moles from the user and saves the result.
- It calculates the critical temperature and displays it to the user.
It also calculates the critical volume.
- It asks the user for the temperature to plot the P-V curve in terms
of a proportion of Tc (so
the input 0.5 would mean use half of Tc).
- It plots the P-V curve for the specified temperature, with a title,
axis labels and legend. The range of volume displayed should be from
0.35Vc to 1.2Vc.
A suitable format is shown at right (here room temperature is 1.94 times
Tc and the curve is close to
that of an ideal gas). Confirm that the shape changes for temperatures
less than Tc. If you've done
it correctly then for 100 moles you should see a minimum for 0.22Tc
at about V=4L.
To get numbers into a string for disp, input, title or legend, use
something like the following. num2str converts a number to a
string, with a optionally specified precision. The following example
specify that the square root of x
is to be displayed to 8 significant figures.
x = 1.23456;
disp(['The square root of ', num2str(x), ' is ', num2str(sqrt(x), 8)])
Part D: Using the comet command
You are asked to write a script to perform the operations specified
below. There are five set of parameters specified below, you only need to
use one of them. Make sure your script is properly documented.
The comet function (or command, there's no difference) is an
animated plotter that displays the trace progressively, though at a fixed,
fast speed. The newly displayed points start out in one colour (the head
of the comet), and gradually change colour as the head moves on.
A good use for comet is to trace Lissajous figures. They are described by
equations of the form
where a and b
are usually small integers and δ
represents a phase shift. They are relevant in signal processing.
For a nifty animation:
- First create a 2000-element vector for the independent variable t,
ranging from 0 to 2π (you need at least that many to slow down the
display to a convenient speed!).
- Assign to variables a, b and delta one
of the combinations in the table below. Use one line so you can easily
return to the statement and change the numbers.
- Assign vectors x and y according to the formula
(remember to re-execute this when you change the other variables).
- Copy the following series of commands, which clears the display and
rescales it so that the figures don't graze the axes.
clf; axis([-1.2, 1.2, -1.2, 1.2]); axis('equal'); hold on;
- Make sure the figure window is visible, then issue the actual display
With the up-arrow key it's easy to change the parameters, or you could
put the steps in a script and ask for a, b and delta
each time it's run.
||Circle (same as in part B)
plenty of info about Lissajous figures, including a reference to an
oscillating echanical signal light called a Mars
Light, which traces a Lissajous curve with a=1,
b=2 to make it more noticeable.
References: Chapman (2012), Exercises 5.38 and 4.11
(starting point only).
Part C background material:
Salzman, WR (2004). Critical Constants of the van der Waals Gas, http://www.chem.arizona.edu/~salzmanr/480a/480ants/vdwcrit/vdwcrit.html.
Last accessed 2014-05-03.
Rose-Petruck, C (1998). Chem 201: First-Order Phase Transitions,
Brown University. http://casey.brown.edu/research/crp/Edu/Documents/00_Chem201/7_phase_equil_1/7-phase_equilibria1-frames.htm.
Last accessed 2014-05-03.
Van der Waals Constants for Gas courses.chem.indiana.edu/c360/documents/vanderwaalsconstants.pdf.
Last accessed 2014-05-03.