>> X1=[-1.2 -0.8 -0.9 -0.2 0 0.5 1.1 1.9
2.3 3.4];
>> Y1=[ 2.2 2.3 2.4 1.9 2.1 2.2 2.4 1.7 1.1 0.8];
>> xp1=linspace(-1.2, 3.4, 30);
>> lsp(X1,Y1,5,xp1)
Lab visit October 17, 2001 - Assignment 4 due October 24, 2001
Bent E. Petersen
The purpose of this assignment and lab visit is to introduce you to using Matlab interactively and also to do a bit of Matlab programming. Assignment 4 consists of the 5 problems that you will find below.
Many features of Matlab, even commonly used ones, will be omitted. Be sure to look at Bob Higdon's notes on Matlab. They give you a much more detailed introduction than what I present here.
Higdon notes on Matlab 239 KB 37 pages PDF
When you start Matlab a window, the command window, opens up. You can type interactive commands in the command window. Editing is a bit primitive, but functional. Cursor up may be used to recall previous commands (in the current session).
Matlab commands are terminated by a newline (that is, pressing the Enter key if you are working at the Matlab command line). Several comamnds may be stacked up on one line by separating them by commas. Matlab's line continuation character is ... This is useful spreading long commands over several lines.
Matlab's assignment operator is the equals sign =. The test for equality is the double equals sign ==.
Matlab commands are contained in m files, one per file (until recently), or are entered interactively at the Matlab prompt.
A list of the directories where Matlab searches for m files can be obtained by the command
>> path
One of the directories in which Matlab always looks is the current working directory; the command
>> cd
identifies the current working directory. You can also use
>> pwd
to get the current working directory. The command
>> cd z:\
changes the working directory to z:\ which is probably a good idea if you are working in the MLC lab. Files saved on z:\ will be visible no matter what lab computer you are using.
The command
>> dir
returns a list of files in the current working directory. You can also use the command
>> ls
Many commands in Matlab have both MSDOS and UNIX style versions. To get a rapid introduction to Matlab and some feeling for what it can do, try the command
>> demo
The demos will keep you busy for quite some time.
Any discussion of Matlab data types is bound to be short. Matlab has essentially only one data type - the array. Let's create a 2-by-2 matrix
>> A = [1 2; 3 4]
A =
1 2
3 4>> A^2
ans =
7 10
15 22
Sure enough this is the square of the matrix A. Now try
>> A.^2
ans =
1 4
9 16
What happened? This is a Matlab idiom - an important one. The dot before a command (let's try to refrain from calling it a dot com) causes Matlab to carry out the command componentwise. That is, each entry in A was squared, rather than A being squared. The dot commands are called array commands, whereas the standard linear algebra commands are called matrix commands.
Here's another example:
B = [3 7; 1 -2]
B =
3 7
1 -2>> A*B, A.*B, B*A, B.*A
ans =
5 3
13 13ans =
3 14
3 -8ans =
24 34
-5 -6ans =
3 14
3 -8
Problem 1: Explain what happened, that is, explain the results above.
>> A.^B
ans =
1.0000 128.0000
3.0000 0.0625
Problem 2: Again explain what happened.
We can recall the value of A by simply entering A. Thus
>> A, exp(A), expm(A)
A =
1 2
3 4ans =
2.7183 7.3891
20.0855 54.5982ans =
51.9690 74.7366
112.1048 164.0738
Here the first matrix returned is A. The second is exp(A) which computes the exponential of each entry in A, whereas expm(A) computes the exponential of the matrix A. Thus exp() is an array operation while expm() is the familiar matrix operation traditionally denoted by exp(). Be careful!
Try
>> help exp expm
to get information on both commands.
While it is possible to do quite a bit with Matlab by just entering commands at the Matlab prompt it is troublesome if you discover that you made an error several steps back. You can use the up and down arrow keys to step through your commands, but this is hardly satisfactory. The best solution is to use m-files.
There are basically two kinds of m-files: scripts and functions.
A script m-file consists of a sequence of commands which are carried out when the script is loaded by typing its name at the command prompt. Be very careful! The main purpose of a script is to add a number of definitions to your workspace, or to carry out a sequence of commands. Everything defined in the script becomes part of your workspace. Variables may be overwritten silently.
A function m-file is an m-file which accepts parameters (including the empty set). Functions generally communicate with the command window through input and output parameters. Global variables are also available, but should be used sparingly. All other variables are local to the m-file and are not "visible" in the command window. In particular, they do not overwrite variables of the same name in the command window.
To create a new m-file open the file menu and select new. An editor window will open up. When you are done, save the file (with the file extension .m). The editor will save the file to your current default directory as set in Matlab. If you already set your default directory that's fine, but you should still check carefully to make sure you are saving to an appropriate personal file location.
A function m-file always starts with the statement declaring it to be a function m-file. Typically the statement looks like this
function[output parameters] = fname(input parameters)
The parameters are separated by commas. The square brackets may be omitted if there is only one output parameter. The square brackets and the equals sign may be omitted if there are no output parameters (for example, a special plot routine).
Here is an example (with an excess of comments) The lines with multiple equal signs are not part of the m-file. You can download this m-file as imagic.m
========== this file must be saved as imagic.m ==========
function [ M, d ] = imagic ( n )
% This totally useless function computes a magic matrix of the size
% requested and adds the identity matrix. It then returns the inverse
% and the determinant of the computed matrix.
%
% Ver. 0.01 Copyright (c) 2001 by The Mad Hatter
% Mth 355 Fall 2001
% All the lines above, but not this one, will be returned in response
% to the command help imagic
% The disp() commands below illustrate how to add text to the
% output of the command. You shouldn't overdo this!
disp('imagic.m version 0.01 Copyright (c) 2001 by The Mad Hatter')
disp('Guaranteed to be totally useless.')
A = magic(n) + eye(n);
M = inv( A );
d = det( A );
% Note the semicolon suppresses output from the statement.
% Note the call [C,a] = imagic(4) assigns the computed matrix
% to C and the computed determinant to a.
% The call C = imagic(4) will assign the computed matrix to
% C and will discard the computed determinant.
% The call imagic(4) will assign the computed matrix to ans
% and will discard the computed determinant.
========== end of imagic.m ==========
Note the first consecutive group of comment lines is returned by the help command. This is very useful - use it! To see how it works, after downloading this m-file try
>> help imagic
Matlab 4 and earlier allows only one function per m-file. This restriction often results in a plethora of small special purpose m-files. Matlab 5 and later allow subfunctions in m-files. This results in cleaner and better organized code (well, it can). Only the main function (that is, the first one) can be called from the command prompt or other m-files. The subfunctions are strictly local.
Let's write a (slightly) more useful m-file. This one computes and displays the polynomial least squares fit for a given data set. It is a very simple combination of commands already available in Matlab. Be sure to use help to look up each of them. You can download this m-file as lsp.m
========== this file must be saved as lsp.m ==========
function lsp(Xdata,Ydata,deg,xplot)
% Given data points Xdata and Ydata this routine plots
% the least squares fit polynomial of degree deg. The
% plot of the polynomial is performed at the abscissas
% given by xplot.
% Mth 355 Fall 2001
% The declaration
%
% function lsp(Xdata,Ydata,deg,xplot)
%
% shows this m-file returns no value.
p = polyfit(Xdata,Ydata,deg);
q = polyval(p,xplot);
plot(xplot,q,'c',Xdata,Ydata,'oy')
========== end of lsp.m ==========
Note how the plot() command can handle several plots. The 'c' specifies the color cyan, the 'oy' specifies little circles in yellow. Sometimes though you do not want to combine all plots in one command, or perhaps one of your plot commands is stem() (try help stem for more information) and you don't see how to combine it with regular plot() commands. In that case you can use hold (hold on and hold off - check help hold) to combine plots.
Let's check our lsp command:
>> X1=[-1.2 -0.8 -0.9 -0.2 0 0.5 1.1 1.9 2.3 3.4];
>> Y1=[ 2.2 2.3 2.4 1.9 2.1 2.2 2.4 1.7 1.1 0.8];
>> xp1=linspace(-1.2, 3.4, 30);
>> lsp(X1,Y1,5,xp1)
The linespace command above creates a list of 30 points equally spaced from -1.2 to 3.4. If you are very attentive you may have noted that the colors in the image do not match what was requested. That is because I edited the image to make it easier to see.
There is still a lot to learn about 2D graphics, but let's jump ahead to 3D graphics. Matlab has a reputation for being able to produce acurate and beautiful 3D graphics.
Here is a sample space curve
>> t=pi/30:pi/60:20*pi ;
>> plot3(t.*cos(t),t.*sin(t),log(t))
>> title('Conical Helix with log rise'), xlabel('x=t*cos(t)')
>> ylabel('y=t*sin(t)'), zlabel('z=log(t)')
The first command creates a list from pi/30 to 20*pi with stepsize pi/60. In the plot3() comamnd note that cos(t), etc., means create a list by computing the cosine at each point in the list t. That is why there is a dot in the expression t.*cos(t) - we want to multiply the list (or array) t with the list cos(t) componentwise. The next four commands, two per line (separated by commas), annotate the image in the current plot window.
Once you have a nice graphic you may wish to save it for use in another application. The command
>> print -depsc2 filename.eps
will produce a full-color encapsulated Postscript file image of your current plot window (you select a name to use in place of "filename" of course). The command
>> print -dmfile filename.m
will create an m-file which can be run to re-create your plot. This is handy if you are deep into an interactive session. The command
>> print -dmeta
will export the image in the current plot window to the Windows Clipboard (on Windows machines) as a Windows Metafile. It can then be pasted into almost any graphics application (I use Paint Shop Pro), touched up and converted as you see fit. You can use -dbitmap instead if your graphics application does not use Metafiles.
The plot3() command plots (intelligently) lists of points in 3-space. To plot the graph of a function of two variables we need to specify a (rectangular) plot domain and points at which to evaluate the function.
>> [X,Y] = meshgrid(0:0.2:4, -2:0.2:2) ;
Here X runs from 0 to 4 in 0.2 increments and Y runs from -2 to 2 in 0.2 increments. Next we compute the values of our function ( (x-2)^2 (1-y^2) ) at the grid points
>> Z = (X-2).^2 .* (1-Y.^2) ;
Note the dots here to specify array operations. We are ready to plot
>> mesh(X,Y,Z)
>> xlabel('x'), ylabel('y')
>> surf(X,Y,Z)
>> xlabel('x'), ylabel('y')
>> colormap(copper)
>> shading interp
>> surfl(X,Y,Z)
>> shading interp
This last image has a light source (which you can move). It works better if you use colormap(gray). When done with this image, to avoid surprises later do
>> colormap('default')
Try also colormap(hot) and colormap(cool). Check colormap() in help.
Try the commands
>> contour(X,Y,Z,30)
>> contour3(X,Y,Z,40)
>> pcolor(Z)
Plot a few graphs of you own choosing.
>> clear x
removes the variable or function x from your workspace, wheras
>> clear
or
>> clear variables
removes all variables and
>> clear functions
removes all compiled functions. (The first time you use an m-file it is compiled and becomes part of your current work space.)
>> clear all
removes all functions and variables from your workspace and gives you a more or less fresh start.
>> who
or better
>> whos
will list all your current variables.
The functions ones(), zeros() and eye(n) are used to create special matrices. Check help and save me some typing!
>> clear x
>> for n = 1:6
x(n) = 3^n ;
end
>> x
x =
3 9 27 81 243 729
Many traditional for loops may be replaced by Matlab's array operations. Thus
>> x = 3.^[1:6]
x =
3 9 27 81 243 729
There is a difference between these two commands. In the for loop we are assigning values to components of x. Thus if x already was defined 1 by 10 array only the first 6 components would be changed. The rest would remain. That is why I cleared x firs. In the second command we assign x as a 1 by 6 array. If it was already defined as something else, it is now redefine.
Note - try to get used to the Matlab idiom. Thus
>> [1:6].^3
ans =
1 8 27 64 125 216
That is, we cube each entry in [1:6] == [ 1 2 3 4 5 6 ].
Let see how many terms of the harmonic series we need to add up to get a sum of 6 or more.
>> s = 1; n = 1;
>> while s < 6
n = n +1 ;
s = s+1/n ;
end
>> n
n =
227
If structures follow the same pattern as in other programming languages. All you need to know really is the keywords are if elseif else end and there is no then. To test for equality use == and to test for inequality use ~=.
A polynomial is represented by its coefficients. Thus 3 x^4 - 2 x^3 + x - 5 is represented by
>> p = [3 -2 0 1 -5 ]
p =
3 -2 0 1 -5
Matlab is very good at finding roots of polynomials
>> roots(p)
ans =
0.2249 + 1.0920i
0.2249 - 1.0920i
1.2714
-1.0545
To multiply two polynomials (as polynomials) use the conv() command. Thus
>> q = [ 3 1 4]
q =
3 1 4
>> r = conv(p,q)
r =
9 -3 10 -5 -14 -1 -20
As we saw above, in lsp.m, polyval() is used to evaluate polynomials and to plot polynomials. For example,
>> x = linspace(-2,2,60) ; % 60 points
>> y = polyval(p,x) ; % evaluation at points in x
>> ax = zeros(1,60) ; % 60 points for x-axis
>> plot(x,ax,'g-.') % plot x-axis
>> hold on
>> plot(x,y,'r') % plot polynomial in red
>> hold off
>> title('polynomial plot')
>> xlabel('x axis')
>> ylabel('y axis')
Be sure to check fplot() if you want to get good plots as automatically as possible.
There is lots more to discover about Matlab but I'm getting tired of typing. Besides, I have to finish before we go to the lab - so I'll end with a few problems for you to play with.
The Engineering school favors Matlab for all sorts of serious engineering problems. As a counter-balance I will ask you to do a few whimsical problems.
Problem 3. Write a function m-file, call it stir(), to return the Stirling approximation to n factorial
sqrt(2*n*pi) * exp(-n) * nn * exp(1/(12*n))
Your routine should work correctly on vectors. Thus
>> stir([2 3 5])
ans =
2.0007 6.0006 120.0026
Problem 4. Write a function m-file, call it logstir(), to return the base 10 logarithm of Stirling's approximation to the factorial of n. Do not just do log10(stir(n)) since that will bomb for n > 143 or so.
Check your work against
>> logstir([100 200 300 400])
ans =
157.9700 374.8969 614.4858 868.8064
Problem 5. Estimate log10(10!!). What about log10(50!!).?