CS101 - Lab activity #5

Objectives

By completing this Lab exercise you should become proficient in the following learning objectives:

Use global variables to write functions that modify the value of variables in the current environment.
Use the Matlab audioplayer function to produce sound from .wav files.
Use Matlab to solve simple dynamical problems.

References

Lecture 5 slides 11-18 on Global Variables.
Matlab Help: Using Matlab Mathematics Functions.
Getting Started with Matlab: Section 4.3.3 (about Global variables).
Getting Started with Matlab: Section 5.5 (ordinary differential equations)

Instructions

In this lab, you will work in groups.

Part 1: Global variables in Matlab

Use the Matlab editor to create the following function named dummy . Make sure to include the ; (semicolon).

function result = dummy(parameter)
% notice that a function may not have inputs and outputs
x = 4;
result = x + parameter;

In the Matlab command window type:

>> clear
>> x = 5;
>> y = dummy(-4);

Answer Question 1 on the answer sheet.
#1 Write the value that Matlab returns when you type the following at the command prompt:
>> x
       x =
                  _____________________________
>> y
       y =
                  _____________________________
>> result
        result =
                      ___________________________

You have already learned that all variables in a function are local . Local means that those variables aren't in your workspace. Up to this point, the only way a function can change a variable in your workspace is if you assign the value the function returns to that variable.

Answer Question 2 on the answer sheet.
#2 Write the value that Matlab returns when you type the following at the command prompt:
>> result = dummy(-4)
       result =
                  _____________________________
>> y = dummy(0)
       y =
                  _____________________________
>> result
        result =
                      ___________________________

Using a global variable is a second way to have a function change a variable in your workspace. The global command is like a hand-shake, to share a variable between a function and a workspace you must execute the global command in both. For example,

Answer Question 3 on the answer sheet.
#3 In the editor modify the code for the dummy function by adding the global command as shown below,

function result = dummy(parameter)
global x
  x = 4;
result = x + parameter;

Write the value that Matlab returns when you type the following at the command prompt:
>> global x
>> x
       x =
                  _____________________________
>> y = dummy(-4)
       y =
                  _____________________________
>> x
        x =
                      ___________________________

The handshake that the global commands make works both ways.

Answer Question 4 on the answer sheet.
#4 In the editor modify the code for the dummy function by removing the statement x =4; so that your function should now look like the following:

function result = dummy(parameter)
global x
result = x + parameter;

Write the value that Matlab returns when you type the following at the command prompt:
>> x = 10;
>> y = dummy(-4)
       y =
                  _____________________________

The global command works with more than one variable. If in some other function you wanted to share variables named a , b , c then you would type global a b c   .

Part 2: Discretization of sound

SETI has recorded a message from deep space. They are sure that the message comes from a truly superior life form that we just can't comprehend. Your goal is to decode the message that will enlighten mankind. Download the file containing the message here to your computer by right-clicking and selecting "Save Link Target As" and keep the default file name 'backward.wav' .
Digital sound is really discretized sound. For example, the analog-to-digital sound card on your (home) computer has a clock that ticks thousands of times a second. Every tick of the clock a voltage measurement is made of the analog signal coming into your sound card, (for instance you can plug a microphone into your A/D card). The individual discrete voltage readings are saved to a file.
We will use the Matlab function wavread to read the voltage values from the file 'backward.wav'. Actually these are not the true voltage values but scaled values so that that max is 1 and min is -1. To use the wavread function type,
>> [volts, frequency, bits] = wavread('backward.wav');
wavread is an example of a Matlab function that returns more than one value. In this case volts is a vector of voltage values, frequency is the frequency of the clock on the A/D card (samples per second). We won't worry about bits for now.
You can play a file in Matlab by using the play function. First ask your TA for a set of headphones. Then type,

>> message = audioplayer(volts,frequency);
>>play(message)

The message is unintelligible but perhaps this superior life form has mastered traveling back in time and perhaps this message has been temporally displaced (recorded backwards).

Answer question #5-#8 on the Lab 5 answer sheet for Part 2.

#5 If the sound has been recorded backwards then let's solve the problem of understanding the mystery message. Specifically, the problem is that if
volts = [v1 v2 v3 ... vn] then we actually want [ vn ... v3 v2 v1]. Fill in the blanks and then type in this command to listen to the message.

>> flipvolts = volts( _________ : _________ : __________);
>> message = audioplayer(flipvolts,frequency);
>> play(message)

#6 What is the mystery message that with enlighten mankind?

_____________________________________________________________________________________________

You can modify the frequency by typing,
>> message = audioplayer(flipvolts,2*frequency);
>> play(message)

(If you want to stop the sound playing try hitting the <control>-Break or <control>-c keys.

#7 Write the command to display the number of values the vector flipvolts contains, (don't use size, that's for matrices)

>> __________________________________

Add an echo effect, we will create a new vector called shifted that is exactly like flipvolts except its values are shifted to the right.
>> shifted = [ 0 * (1:1000)' ; flipvolts(1:end-1000)]; ( first 1000 values are zero then the remaining values are flipvolts)
We don't want to over-drive the speakers so we will normalize the output peak to one so type,
>> peak = max( abs(flipvolts + shifted));

#8 Fill in the blank below with the new vector which is the old sound + new sound. You should hear the echo effect.

>> player = audioplayer( (_______________________________)/peak , frequency);
>> play(player)

Part 3: Dynamics

Intro : You may skip to Free Fall : Parachuting if you know physics with calculus

Consider a parachutist jumping from an airplane. We will consider on his/her vertical motion so this is a one dimensional problem. Ignoring air resistance, a free falling object (the parachutist) of mass m has only the force -m*g acting on it where m is the mass in kg and g = 9.8 m/s². From Newton's second law

    m*a = F

where a is the acceleration and F is the force and thus

   m*a = -m*g

Dividing both sides by m and since a = dv/dt where v is the velocity , t means time and dv/dt means the derivative of velocity with respect to time we have,

   dv/dt = -g

This equation is called an Ordinary Differential Equation.
Note that a solution to the above equation isn't a number, the solution is a function. That is, we have suppressed the fact that v is a function of time or mathematically v = v(t), so that more accurately we should write,

dv(t)/dt = -g

It's easy to see that v(t) = -g * t is a solution since the derivative of (-g * t) equals -g.
However v(t) = -g*t + 1 is also a solution and so is v(t) = -g*t + 2 and so on...
Thus an ordinary differential equation may have an infinite number of solutions.
If we knew in advance the velocity of the parachutists then there is one unique solution. For example if the initial velociy (at t = 0) v(0) = 0 then the unique solution is

v(t) = -g * t                                                                                                                     (equation 1)

Further since v = dy/dt , where y is the height of the parchutist above the ground, we can substitute and get

dy/dt = -g * t

We can guess that the solution is y(t) = -(g/2) * t² + constant . So if we knew the initial height of the parachutitst above the ground, say 5,000m then the unique solution for y(t) would be

y(t) = -(g/2) * t² + 5000                                                                                                  (equation 2)

Free Fall : Parachuting , without air resistance

Neglecting air friction, a paracutist motion is described by the following system of ordinary differential equations

dy/dt = v

dv/dt = -g

with initial conditions v(0) = 0 , y(0) = 5000 and the constant g = 9.8.

First we will show you how to use the Matlab function named ode45 to solve this system of ODE (find the functions y(t) and v(t)). Then we will have you solve the same problem except with air resistance.
Step 1. Use Matlab editor to write a function to compute the derivates of all the dependent variables, in this case y and v. ( ode45 will call this function)
     function Dyv = Derivatives_wo_resistance(t , yv)
               y = yv(1);
               v = yv(2);
               g = 9.8;
               Dyv = [ v ; -g];
Step 2. At the Matlab prompt call the ode45 function with initial values y(0) = 5000 , v(0) = 0 over the time range 0 to 100 seconds as follows:

>>   [ t , yv ] = ode45('Derivatives_wo_resistance', [0 100], [ 5000 ; 0]) ;

and to plot the solution type,

>> y = yv(: , 1); % yv is a two column matrix of values, the first column contains y values
>> v = yv(: , 2); % the second column of yv contains v values

>> plot( t , y)

>> plot(t , v)

Notice that ode45 doesn't give us the soluton to this system of ODE's as shown in equation 1 and equation 2 above, rather ode45 gives a discretized solution.

The arguments to the ode45 function are described below:
ode45( 'name of function that calculates the derivates of dependent variables' , vector containing range of values for independent variable*, vector of initital values)
* - the vector containing the range of values of the independent variable must start at the same value as the initial values are given, (in our example the vector [0 100] starts at 0 since at t= 0 we were given y(0) = 5000 and v(0) = 0).

Free Fall : Parachuting , with air resistance

Including air friction, a paracutist motion is described by the following system of ordinary differential equations

falling object

dy/dt = v
dv/dt = -g - (b/m) * v

with initial conditions v(0) = 0 , y(0) = 5000 and constants g = 9.8 , m = 100 and b = 10.

Answer question #9.

#9 Fill in the blanks to complete the code for the function named Derivatives_resistance .
     function Dyv = Derivatives_resistance(t , yv)

               y = _____________________________________

               v = _____________________________________
               g = 9.8;

             _________________________________________

               __________________________________________

               Dyv = _______________________________________

Step 2. At the Matlab prompt call the ode45 function with initial values y(0) = 5000 , v(0) = 0 over the time range 0 to 100 seconds as follows:

>>   [ t , yv ] = ode45('Derivatives_resistance', [0 100], [ 5000 ; 0]) ;

and plot just the solution for v(t),

>> v = yv(: , 2);

>> plot( t , v)

Answer question #10.

#10 What is the terminal velocity( the velocity at t = 100)? (approximately; and the units are m/s)

Chaos : Lorenz Equations -- a simple model of the weather.

Lorenz derived a simple model for weather prediction (see also poster ) . The equations are:

   dx/dt = 10* (y - x)

    dy/dt = -y + 28*x - x*z

    dz/dt = (-8/3)*z + x*y

The independent variable t represents time but the dependent variables x, y and z in this example aren't coordinates in 3D space but we will pretend that they are. As in the previous exampls x = x(t) , y=y(t) and z= z(t) but we have suppressed the t in the above equations. The problem is to find the curves x(t), y(t) and z(t). There are an infinite number of solutions to the Lorenz equations but if we say that at t=0 that x(0) = 1, y(0) = 1 and z(0) = 1 then there is a unique solution.

Answer question #11.

11.Complete the function named xyzprime that will compute the derivatives of the three dependent variables x, y and z , by filling in the blanks.

function Dxyz = Derivatives_xyz(t , xyz)
x = xyz(1);
y = xyz(2);
z = xyz(3);
Dxyz = [ ___________________________; ________________________________ ; _____________________________________];

Open the Matlab editor and copy/paste the above along with your solution.

We want to predict the weather so the question is when t=100 sec what are the values of x , y and z ? Further, we would like to see the results plotted dynamically, that is, we would like to see the point (x, y, z) plotted in 3D as ode45 is solving the problem at each time t. Do this by typing at the Matlab command window,

>> options = odeset( 'outputfcn', 'odephas3');
>> [t ,xyz] = ode45('Derivatives_xyz', [0 , 100], [1 ; 1 ; 1], options); % initial values of (x,y,z) are (1,1,1) at time = 0 , run ode45 on the time interval [0,100]

Answer question #12.

#12 Show your lab TA your dynamic plot for a signature.

Now you can solve any system of ODE's you like, but be careful, sometimes ode45 will not give you the correct results. Matlab offers a suite of ODE solvers each solver optimized for specific classes of problems. To learn how these solvers (like ode45) work you can take CS459 .

--. That's it, you're done with Lab 5!