Our class computational lab will be located in Geology 310, which offers both MATLAB and Maple. To access MATLAB in Geology 310, double-click on the `MATLAB 5.3' icon on the left-hand-side of the screen. This and future labs will be available online through Dr. Powell's homepage at
The program MATLAB , which is a general numerical analysis/numerical linear algebra tool, is available for student use in the Engineering public access lab, EL105. This lab is in the Engineering Lab building. Dr. Powell also has an (aged) educational license for MATLAB , so if you would like to install this program on your personal computer you are welcome to check out the CD at his office.
Information about the open access engineering lab is available on the Web at
In addition, MATLAB Help documentation is available on campus at
as well as at the MathWorks Web site
Open a Web browser and link to the documentation in the Engineering Lab to support this Lab exercise. The `Getting Started' guide (handed out in class) should be there in the left hand column; click on this and then click on `Manipulating Matrices' in the right hand frame.
The MATLAB program was based originally on a set of computational routines for
doing linear algebra, and retains a linear algebra bias in its operations. To
define the matrix
A = [ 1 2 3; 4 5 6; 7 8 9]The column vector and row vector can be defined
c = [1; 2; 7] r = [2 4 6]To multiply the matrix A on the left by the row vector use the MATLAB command r*A; to multiply the column vector on the left by A use the MATLAB command A*c. Now, try these backward and see what happens!
MATLAB is a great program for doing matrix operations. For example, to find the eigenvalues of A use the command
eig(A)Type the help command, help eig, to see what else the eig command is capable of in MATLAB . In particular notice that the command
[V,D]=eig(A)returns the eigenvalues of A as columns in the matrix V. You should have noticed that one of the eigenvalues of A is zero; this is because the matrix A is singular, that is, the rows are dependent. (In this case, if you take twice the second row and subtract the first, you get exactly the third.) You could determine this using by-hand row-reduction, but for giggles try the following command to get the reduced echelon form of A:
rref(A)
Command Line Help MATLAB has a nice command-line help facility which you can use to determine what commands do or what commands might be capable of doing certain things. At any time during a MATLAB session you may type help COMMAND for MATLAB's internal information on how COMMAND may be used, including examples. If you do not know what command you are looking for, you may type lookfor SUBJECT to get MATLAB to search for commands which refer to your SUBJECT. To see what other basic matrix manipulations MATLAB is capable of, type help matfun. To see a general list of help topics, try just help.
Take some time now to work through the Matrix Manipulation section of the Getting Started with MATLAB guide. In particular, familiarize yourself with the following types of operations:
Shortcuts There are a couple of nice shortcuts in MATLAB that make things go much quicker. One of them is the key. This brings back commands that you have typed previously, which can then be editted on the command line. Try using and see what happens. Also, you can use to search for the last command that begins with certain letters. So, for example, to see the last command you typed in which begins with A = type `A =' on the command line and then hit the key.
Another shortcut: Often you do not want to see the output of a command, particularly when it is a long vector or big matrix. To suppress output, follow the MATLAB command with a semicolon.
EXERCISE 1: In Edelstein-Keshet an annual seed-plant model is presented,
with
Here
|
||
Often we are interested not in matrices but scalar functions of an independent variable. MATLAB retains its linear algebra bias; we define dependent and independent variables as row vectors in MATLAB . So, to create an `x' variable defined in ten increments from zero to one, we would use one of the equivalent commands
x=1/10*[0 1 2 3 4 5 6 7 8 9 10] x=[0:10]/10; x=linspace(0,1,11);In the latter two commands a semicolon is used after the input to suppress the output; try it without the semicolons to see what that means and to test that the three definitions of the vector `x' are equivalent.
Now let's define some functions of ; sin() and can be defined as follows:
f1 = sin(x); f2 = x.^2;Note the operation . 2, which means take the vector `x' and create a new vector by squaring each element of `x.' In MATLAB the operation 2 (without the period in front) is reserved for squaring matrices! Now suppose we want to write the function
f3 = ( x.*f1 )./(1 + 3*f2);Some things to note here. Firstly, the operation .* means element-by-element multiplication of vectors (or matrices) in MATLAB , and is distinct from the operation *, which means real matrix multiplication. To see what I mean, compare the results of A*A and A.*A. Note also the operation ./, which means element-by-element division, and is distinct from the / operation, which means `matrix quotient' in MATLAB . Finally, in creating the denominator of notice that addition and multiplication by scalars requires no special notation (breathe a sigh of relief).
EXERCISE 2: Use the commands above to define the function
on the interval , using 100 equally-spaced samples of . |
||
Now let's see about plotting some functions. After defining the function in the previous exercise you might want to see what it looks like. Use the command
plot(x,g)A new window should pop up, with your graphic output. You may want to plot against its `envelope' or amplitude, . One way to do this is to use the hold on command, which tells MATLAB not to erase the current picture when drawing a new one. This would generate the positive and negative envelope functions on the same graph as :
hold on, plot(x,1./(1+x.^2),'r') plot(x,-1./(1+x.^2),'y'), hold offWhat colors did the new functions come out in, and what does that have to do with components of the commands above?
Another way to plot several functions at once in different colors is to just do it. To make your display more pretty, type whitebg('k') first, and then type the sequence:
plot(x,g,'m',x,1./(1+x.^2),'g',x,-1./(1+x.^2),'b')Which should generate output in magenta, green and blue. Which function is which color? Use help legend to see how to use put a legend on a graph; labels can be added using xlabel and ylabel.
One may also plot data using discrete points as opposed to connected lines, as in this example which plots a polynomial approximation to sin(x) in little magenta circles:
x=linspace(-2,2,100); plot(x,sin(x),'y',x,x-x^3/6+x^5/120,'mo');To see some of the other printing options try help plot. At this point you should play with the plotting software a bit; in the Getting Started guide go to Graphics and look through Basic Plotting.
One can write simple programs, or `scripts,' in MATLAB , which are just lists of MATLAB commands written in a text file (called m-file for the `.m' file extension MATLAB uses) that can be called from within the MATLAB command window. To create one of these m-files, click on the File menu in the upper left corner of the command window, then `New' and `M-file.' A text-editor window should pop up, in which you can type a sequence of commands. For example, you could type in the following commands to compare the sin(x) with several of its Taylor polynomial approximations:
x=linspace(-2,2,100); s1=x; s2=s1-x.^3/6; s3=s2+x.^5/120; s4=s3-x.^7/(7*6*120); plot(x,sin(x),'m',x,s1,'r',x,s2,'y',x,s3,'g',x,s4,'b')To save your m-file go to the `File' menu and click `Save As.' You must get on the C:/Temp disk, or use your own floppy, then save the m-file as something, say yourname1.m (I will use jim1 below). Within the MATLAB command window, you can now type
cd c:\temp % this changes directories to the location of your m-file jim1and the sequence of commands you put in the m-file will be executed. Alternatively, on the script editor window at the top there is a little button which looks like a piece of paper with a down-arrow beside it; clicking on this button will also run your script.
One of the reasons to use an m-file is the need for loops which will iterate a
certain activity (like generations of the plant model). Recall that our
simple plant reproduction model (considering adults only) was given by
p=zeros(1,10)Now let's set initial data for the first two generations:
p(1)=0; p(2)=1;Here is a for loop which will calculate the next eight generations:
for n = 2:9 p(n+1) = 2 * p(n) + 3/10 * p(n-1) endTry putting this into an m-file along with a plot command to visualize the growth of the population.
EXERCISE 3: Find the solution to the equation
with the initial data and , where is the initial number of adult plants. Write an m-file which numerically calculates solutions for a few generations to the same equation, and which plots the numerical solution on the same set of axes as the analytic solution. |
||
As mentioned earlier, the annual plant model also has a matrix form:
p=zeros(2,10) % set up a vector to hold ten generations of plants % adults are top row, seeds are bottom row % define parameters sigma=1/2; alpha=2/5; beta=1/5; gamma=10; % pick the initial condition, one adult, no seeds p(:,1)=[1;0]; % the ':' means 'both rows' % set up the reproduction matrix M=[ alpha*sigma*gamma beta*sigma*(1-alpha); sigma*gamma 0 ]; for j=2:10 % begin loop for calculating next generations p(:,j) = M*p(:,j-1); end % end of loop % plot results with seeds on x, adults on y: plot(p(1,:), p(2,:), 'b*'), xlabel('Seeds'), ylabel('Adults')Save and run this script and observe the results. Try different initial conditions and use the hold on and hold off commands to plot the results of different initial conditions on the same graph - what do you oberve? How do the eigenvalues and eigenvectors you calculated earlier relate to what you are observing now?