## ENGG1811 Lab 8: Using Matlab

### Objectives

After completing this lab, students should be able to

• Use Matlab as a calculator working with arrays in one and two dimensions
• Understand how to plot variables in 2 dimensions
• Perform array operations
• Use the comet function to aid in visualising function movement

### Assessment

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 any time.

### 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 result.
1. Create a 2x3 array m whose first row contains the values 1, 2 and 4 and whose second row contains 3, -7 and 16.
2. Confirm how big m is with the size and length functions.
3. 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 transpose only.
4. 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
```
5. (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:
 Component: MTBF (hours): Antenna Transmitter Receiver Computer 2000 800 3000 5000
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 the documentation.

### Part B: Plotting

1. 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.
2. Repeat the exercise using linspace, assign to x2. (Hint: You can use x to see how many points you need.)
3. 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 .
4. 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).
5. 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).
6. 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 function.
7. 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 equations 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;
```
8. 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 size functions.
9. 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 gas.

In derived SI units, for oxygen, a = 138.2 kPa (L/mol)2 and b = 0.0319 L/mol.

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 form.

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 preceding comment.
• 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:

1. 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!).
2. 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.
3. Assign vectors x and y according to the formula (remember to re-execute this when you change the other variables).
4. 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;
```
5. Make sure the figure window is visible, then issue the actual display command:
```comet(x, y)
```

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.

a b δ Figure
1 1 0 Line
1 1 π/2 Circle (same as in part B)
2 1 π/4 2-lobed parabolic
5 4 π/2 5x4 lobes
1 3 π/2 You-know-who's logo

Wikipedia has 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.