Introduction to MATLAB for Biology/Math 4230

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

http://www.math.usu.edu/ powell
In the right-hand frame follow the link for Bio/Math 4230.

Access Outside Geology 310

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

http://www.engineering.usu.edu/pclab/

In addition, MATLAB Help documentation is available on campus at

http://www.engineering.usu.edu/pclab/matlab-help/helpdesk.html

as well as at the MathWorks Web site

http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml

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.

Matrices and Basic Operations

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

\begin{displaymath}
{\bf A} = \left( \begin{array}{rrr}
1 & 2 & 3 \\
4 & 5 & 6\\
7 & 8 & 9 \end{array} \right)
\end{displaymath}

in MATLAB use the command
       A = [ 1 2 3; 4 5 6; 7 8 9]
The column vector $\vec{c} = (1, 2, 7)^T$ and row vector $\vec{r}=(2,4,6)$ can be defined
       c = [1; 2; 7]
       r = [2 4 6]
To multiply the matrix A on the left by the row vector $\vec{r}$ use the MATLAB command r*A; to multiply the column vector $\vec{c}$ 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:

You may also want to check out some of the MATLAB demos. Click on the Help menu and then click on Demos. In the left hand menu of the windo that appears click on Matrices and then on Basic matrix operations in the right hand menu. This will generate a slide show which will lead you through several of the basic MATLAB matrix and plotting commands.


Shortcuts There are a couple of nice shortcuts in MATLAB that make things go much quicker. One of them is the $\uparrow$ key. This brings back commands that you have typed previously, which can then be editted on the command line. Try using $\uparrow$ and see what happens. Also, you can use $\uparrow$ 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 $\uparrow$ 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 $(s_n,p_n)$ model is presented, with

\begin{displaymath}
\left(\begin{array}{c} p_{n+1} \\ s_{n+1} \end{array}\right)...
...ight)
\left(\begin{array}{c} p_{n} \\ s_{n} \end{array}\right)
\end{displaymath}

Here
$\gamma$
is the number of seeds produced per plant in a year,
$\alpha$
is the fraction of seeds that germinate after one year,
$\beta$
is the fraction of seeds that germinate after two years, and
$\sigma$
is the fraction of seeds that survive the winter.
Let $\sigma =\frac12 , \alpha=\frac25, \beta=\frac15, \gamma=10$. Use the MATLAB command eig to find eigenvalues and eigenvectors. Test these eigenvalues and eigenvectors using matrix*vector multiplication (that is, check that the matrix times the eigenvector gives the eigenvector back, multiplied by the eigenvalue). The asymptotic (long-time) distribution of plants and seeds is given by the eigenvector corresponding to the largest eigenvalue. In this case, what is the predicted distribution of plants and seeds? What is the expected growth rate (i.e. eigenvalue) of this population?
 
     
 

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 $x$; sin($x$) and $x^2$ can be defined as follows:

       f1 = sin(x); 
       f2 = x.^2;
Note the operation . $\hat{\hspace{.1in}}$2, which means take the vector `x' and create a new vector by squaring each element of `x.' In MATLAB the operation $\hat{\hspace{.1in}}$2 (without the period in front) is reserved for squaring matrices! Now suppose we want to write the function

\begin{displaymath}
f_3 = \frac{x \sin(x)}{1 + 3 x^2}
\end{displaymath}

in its MATLAB form. We would then write
       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 $f_3$ 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

\begin{displaymath}
g = \frac{\sin(x^2)}{1 + x^2}
\end{displaymath}

on the interval $x\in [-7,7]$, using 100 equally-spaced samples of $x$.
 
     
 

Plotting Functions

Now let's see about plotting some functions. After defining the function $g(x)$ 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 $g$ against its `envelope' or amplitude, $\pm (1+x^2)^{-1}$. 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 $g$:
       hold  on, plot(x,1./(1+x.^2),'r')
       plot(x,-1./(1+x.^2),'y'), hold off
What 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.

Basic M-Files

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
       jim1
and 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.

Loops

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

\begin{displaymath}
p_{n+1} = \alpha \sigma \gamma p_n+\beta \sigma^2 (1-\alpha) \gamma p_{n-1}.
\end{displaymath}

When $\sigma =\frac12 , \alpha=\frac25, \beta=\frac15, \gamma=10$ this becomes

\begin{displaymath}
p_{n+1} = 2 p_n+ \frac3{10} p_{n-1}.
\end{displaymath}

Let's define a vector which contains the adult plant numbers for ten generations, using the zeros command to set the numbers initially to zero:
       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)
       end
Try 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

\begin{displaymath}
p_{n+1} = 2 p_n+ \frac3{10} p_{n-1}.
\end{displaymath}

with the initial data $p_1=0$ and $p_2=N$, where $N$ 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:

\begin{displaymath}
\left(\begin{array}{c} p_{n+1} \\ s_{n+1} \end{array}\right)...
...ight)
\left(\begin{array}{c} p_{n} \\ s_{n} \end{array}\right)
\end{displaymath}

Let's write a small m-file to predict the growth of this population when $\sigma =\frac12 , \alpha=\frac25, \beta=\frac15, \gamma=10$. Create a new m-file and in the editor type
       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?



James Powell
2002-01-13