CS101 - Lab activity #5

Name: _________________________

Netid: ___________________ Section: ___________

Objectives

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

References

Instructions



Part 1: Prelab

You should use the material found in Lecture 5,6 and 7 of the course notes and/or Matlab by Gilat Chapters 5.2, 9.2, 10.1, 10.3 and 10.4  . As with all prelabs/labs, if you are not sure of the answer you may go to any EWS lab on campus and check your results by using Matlab, or you can ask anyone of the CS101 staff for help.

1. Using the built-in function fzero write the Matlab command(s) to find a root of the function

y=x7- ln4(x)
with an initial guess of x = 5. (Don't write the root 0.6359 approx.)

_____________________________________________________________

2. Using the built-in function fplot write the Matlab command to plot the function

y=x4 - ln4(x) on the interval [1,2].


____________________________________________________________

3. Using the built-in function quadl write the Matlab command to compute the definite integral of the function




on the interval [1,5]. (Don't write the answer 15.3662, Just the Matlab code that gives the answer.)

____________________________________________________________

4. Using the built-in function trapz write the Matlab command to compute the definite integral of the function



the interval [1,5] using 100 points. (Don't write the answer 15.3660. Just the Matlab code that gives the answer.) You may use multiple lines.

x = __________________________________________________________


y = __________________________________________________________


trapz( ____________________________ , ________________________________)



5. Complete the following Matlab code, by filling in the blanks, to plot the surface described by the equation,

z = cos( sqrt(x2 + y2)) / sqrt(x2 + y2)

over the interval [-10,10] x [-10,10] in the x-y plane.

>> x = linspace(-10,10,400);

>> y = x;

>> [XMAT, YMAT ] = meshgrid(x,y);

>> ZMAT = _________________________________________________________________________

>> surf( XMAT , YMAT , ZMAT )

6. We want to model the motion of a free falling object using the first order ordinary differential equation: http://en.wikipedia.org/wiki/Free_fall



(see also Classical Mechanics: A Modern Perspective 2nd ed. Barger and Olsson p. 8) where the mass of the object m = 75 Kg ,the acceleration due to gravity g = 9.8 m/s2 and


and

drag coefficient CD = 0.51, http://en.wikipedia.org/wiki/Drag_coefficient

cross-sectional area of the object, perpendicular to air flow A = 1.1 m2 ,

air density ρ = 1 Kg/ m3 ,

We want to solve for the vertical velocity v (which is a function of t written v = v(t) ) using the Matlab ODE solver ode45. Assume that at time t = 0, v = 0

( v(0) = 0).

To solve this equation we would type,

>> [ t,v] = ode45(@vprime, [0,100] , 0);

But we first need to complete the function named vprime.

a) Complete the function vprime by filling in the blank below.

function dv = vprime(t,v)

m = 75;

CD = 0.51;

A = 1.1;

rho = 1;

g = 9.8;

dv = ____________________________________________________________________________

b) Using this model what is the terminal velocity ? _____________________

(Hint: try plot(t,v) or v(end) )


Part 2: 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(-3);


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. Note that when you called the function dummy you saw x = 4 display in the command window. Since we omitted the ; to suppress output we then saw the execution of that command although it wasn't executed in our workspace.


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:

>> clear
>> global x

>> x = 5;
>> 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;

Save your modified code. 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 3: Discretization of sound

SETI has recorded a message from deep space. They are sure that the message comes from a truly superior life form whose intelligence is beyond our comprehension. 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 (scaled) 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. (If you are logged in remotely then you will not hear any sound since this program is runing on an EWS computer.) 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 4: Dynamics

Intro : You may skip to Free Fall: without air resistance 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/s2. 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 a first order 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) * t2 + 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) * t2 + 5000 (equation 2)

Free Fall : without air resistance

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

dy/dt = v (definition of velocity)


dv/dt = -g (from Newton's second law, read previous section)



Since there are two ordinary differential equations we will need two initial conditions to obtain a unique solutionx, so let 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. The function ode45 will call this function which we will name Derivatives.
The function ode45 will pass two values to our function Derivatives, a scalar t (time) and a column vector yv (containg two scalars, a value for y and for v).

function dydv = Derivatives(t , yv)
y = yv(1); % yv = [ some_y_value some_v_value] so use subscripting to obtain y
v = yv(2);
g = 9.8;
dydv = [ v ; -g]; % since dy = v and dv = -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', [0 100], [ 5000 ; 0]) ; % note that we are telling ode45 that at time =0 that y = 5000 and v = 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 : with air resistance

From prelab 5 we can describe the equation of motion system of first order ordinary differential equations:

dy/dt = v (definition of velocity)


dv/dt = (c/m) * v2 -g (from prelab 5)


Answer question #9.

9. Fill in the blank to complete the code for the function named Derivatives_resistance .

function dydv = Derivatives_resistance(t, yv)
m = 70;
CD = .5;
A = 1;
rho = 1;
g = 9.8;
y = yv(1);
v = yv(2);

dydv = ____________________________________________________________________________

Now solve the system of ordinary differential equations using ode45 with initial conditions at t = 0 , y = 5000, v = 0.

[t , yv] = ode45(@Derivatives_resistance, [ 0,120], [5000 ; 0]);

Answer questions 10, 11, 12.

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 #13.

13. 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 dxdydz = Derivatives_xyz(t , xyz)
x = xyz(1);
y = xyz(2);
z = xyz(3);
dxdydz = [ ___________________________; ________________________________ ; _____________________________________];


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 #14.

14. The actual data points (x,y,z) are plotted using
A) squares
B) stars
C) circles
D) diamonds


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!