This tutorial is designed to reacquaint you with the MATLAB ideas initially presented by Tom Huber at the Envision It! workshop held at Arlington High School in late January. My hope is that this material will prove to be good preparation for Tom's follow up workshop on April 12th.
Most of the material in this tutorial is derived from the worksheets presented at the January workshops. These documents can be accessed directly by visiting Tom's homepage at Gustavus Adolphus College. Check under the Envision It! item.
A matrix is simply a rectangular list of numbers. The "size" of a matrix is given as two numbers, the first is traditionally the number of rows in the matrix while the second is the number of columns in the matrix. Matrices are usually written in tabular form contained between two large parentheses or square brackets.
For example, the matrix below is a 2x3 (this size is read "two by
three") matrix without the enclosing parenteses or brackets.
34 56 31 -45 6 43
Arrays, sometimes called vectors, are special cases of matrices, namely arrays have a single row or a single column of numbers. Arrays with a single row of numbers are called row vectors to contrast them with column vectors. Only one subscript is needed to reference a particular entry in an array.
To create a matrix and assign it a name within the MATLAB program simply fire up MATLAB and at the prompt enter the matrix's name, an equal sign (=), and an open bracket ([). Next enter the numbers, separated by either commas or blank spaces, which make up the first row of the matrix. If your matrix has only one row, type in the close bracket (]) and hit RETURN. If you have several rows to enter, signify the end of each row with a semicolon before proceeding to enter the values for the next row. (No semicolon is necessary after the last row.) When you have finished entering all the entries of the matrix type a close bracket and hit RETURN. Your matrix should appear on the screen.
As an example, consider the 2x3 matrix given above. To give this matrix the name A we would enter the following at the MATLAB prompt:
A = [ 34 56 31; -45 6 43 ]
MATLAB includes a command which returns the size of an array or matrix. The command is the size command. It returns a row vector whose first entry is the number of rows in the argument and whose second entry is the number of columns.
As an example the command
size(A);would return the vector [2 3] if A is the matrix defined in the example above.
Many special matrices appear so frequently that the designers of MATLAB have created shortcuts for creating them.
temp = zeros(10)will do the trick.
To create a non-square matrix consisting of all zeros it is necessary to use a second version of the command. This version takes two arguments, the number of rows and the number of columns in the matrix being created. For example, to create a 3x5 matrix of all zeros, the command
many0 = zeros(3,5)could be used.
lonely = ones(1,99)
ident = eye(5)creates the 5x5 identity matrix.
test = 0:4:20would create the row vector [0 4 8 12 16 20] and names it test. The first entry in the vector is given by the first entry in the command, the second entry in the command gives the interval between entries (a.k.a. the stepsize), 4 in this case. The last entry tells the process when to stop, at 20 in this example.
If only two numbers are given, e.g. 1:10, MATLAB assumes the stepsize is 1.
In some of these examples you may have been overwhelmed with output after you entered the command. One nice MATLAB trick is to suppress the printed output by ending any command with a semicolon.
a1 = zeros(10) a2 = zeros(10);
MATLAB's power is in its ability to perform matrix arithmetic very efficiently. MATLAB divides its matrix arithmetic into two varieties, componentwise operations in which the indicated operation is performed on two matrices of exactly the same size and the resulting matrix is also of this common size. The entries in the resulting matrix are calculated by performing the operation on entries occupied the same position within each matrix.
The second type of matrix arithemetic is the more traditional type of matrix multiplication and exponentiation taught in high school and college algebra classes.
Aside from operations on pairs of matrices, MATLAB also allows for operations on a matrix and a number (also known as a scalar).
In the sections below we cover each of these types of arithmetic operations.
There are four scalar operations, addition, subtraction, multiplication and division, denoted by the symbols + - * and / respectively. In each case the indicated operation is performed on each entry in the matrix.
As an example, consider the two MATLAB commands shown below:
fdj = [ 1 2 3;5 4 3; 6 5 8 ]; abc = 3*fdj;The matrix abc contains entries created by multiplying the corresponding entries of fdj.
abcwhich should show you the contents of the matrix abc.
fdj+3 fdj-6 fdj/2.in each case the indicated operation should be performed on each entry of the matrix fdj.
MATLAB's componentwise matrix addition and subtraction matches the traditional form of the same operations. For two matrices to be added or subtracted they must be of the same size. When two matrices are added the resulting matrix is of the same size as the original matrices. The entries are computed by adding or subtracting the corresponding entries in the two original matrices. This addition or subtraction is indicated to MATLAB by placing a + or - sign between two matrices. For example, the MATLAB commands:
abc = [1 2;3 4] def = abc - [6 -1; 2 0]would print out the square matrix with -5 and 3 in the first row and 1 and 4 in the second row. This new matrix would be given the name def.
abc = 1:10; def = 5:14; ghi = 3*abc + def
abc = [1 2 3 4;5 6 7 8]; def = [4 3 2 1;0 -1 -3 3]; abc + def
For example, the three MATLAB commands:
abc = [1 2 3; 4 5 6]; def = [1 2 1; -1 3 -2]; abc.*defcreate a matrix whose first row has entries 1, 4 and 3. The second row has entries -4, 15 and -12.
Note that MATLAB does NOT use the asterisk to denote this componentwise multiplication. The asterisk, by itself, represents the more traditional form of matrix multiplication. To indicate componentwise multiplication MATLAB uses the dot-star notation, namely a period followed immediately by an asterisk.
The same idea holds for componentwise division. The MATLAB symbol for this is dot-slash (./).
MATLAB uses the dot-carat (.^) to signify componentwise exponentiation. The exponent is a scalar and the .^ operation simply means to raise each entry in the matrix to the indicated power. For example,
abc = [1 2;3 4]; abc.^2would print out a 2x2 matrix whose entries were 1, 4, 9 and 16.
Traditional matrix multiplication is a very useful operation which is somewhat tedious to define. It is beyond the scope of this MATLAB tutorial, but I encourage you to discuss the matter with a mathematics teacher at your school if you are interested. For this tutorial we will stick to a very useful special case of matrix multiplication, namely the multiplication of a row vector by a column vector.
For the product of a row vector by a column vector to be defined, the two vectors must contain the same number of entries. If this is true, then the product is a scalar computed by first multiplying the corresponding entries in the two vectors and then summing up the resulting products. This idea is almost identical to the idea of a dot product or inner product often discussed in vector algebra sections.
The symbol MATLAB uses to indicate traditional matrix multiplication is the asterisk without a dot.
This rather strange definition of matrix multiplication can be extremely useful. If the product turns out to be zero, this means the vectors, considered as arrows on a graph, are perpendicular to each other. If the two vectors being multiplied have the same entries in the same order, then the square root of the product is the vector's length. Lastly, multiplying a vector by an appropriately sized vectors of all ones is a quick way to add up the entries in the original vector.
A common task facing the MATLAB user is to change a column or row vector into the other. This can be accomplished using the MATLAB transpose operator, the apostrophe ('). If A is a row vector, the notation A' signifies a column vector whose entries are identical to those of A. A similar transformation occurs if you begin with a column vector.
abc = [1 2 3 4]; def = [2 5 4 3]; abc * def
abc = [1 2 3 4]; def = [2 5 4 3]; abc * def'
A*ones(1,123)'which multiplies the vector A with the column (after the transpose) vector made up of 123 ones. The result of this product is the total of the 123 entries in A.
To explore this technique create a row vector with 10 entries, and use the idea above to direct MATLAB to compute the sum of the entries.
MATLAB supports a large suite of functions, including the usual trigonometry functions as well as functions designed specifically for matrices, such as the matrix inverse function.
Most of the functions we will use are functions which are often applied to single numbers when used outside the context of MATLAB. An example would be the sine trigonometric function. When applied to a matrix in MATLAB these functions perform in a componentwise fashion. As an example, consider the matrix
nums=[0 5*pi/6 pi/2; pi/6 pi 2*pi]The command
sin(nums)would return the 2x3 matrix whose first row is 0, 0.5 and 1 while the second row is 0.5, 0 and 0.
To create two dimensional graphs of functions or other data, it is necessary to create two row vectors which contain the x and y coordinates, respectively, of the points to be plotted. MATLAB plot command connects the points indicated by these coordinates with a series of straight lines.
For example, the following three MATLAB commands generate a unit square whose lower left hand corner is the origin.
x = [ 0 1 1 0 0]; y = [ 0 0 1 1 0]; plot(x,y)While the physical drawing of the square occurs too quickly to watch, MATLAB begins at the origin and draws the square in the direction indicated by the points, in this case counterclockwise.
More traditional mathematical graphs can be created by first defining
the range of x values over which you want to graph and then
creating a vector which contains a good number of value in this
x range. For example, to create a graph of the function
y = sin(x)*cos(x)^2
over the range -2*pi to 2*pi, we could enter the following MATLAB commands:
x = -2*pi:0.1:2*pi; y = sin(x).*cos(x).^2; plot(x,y)Note the componentwise arithmetic here. This is critical to the proper functioning of the plot command.
It is possible to place more than one graph on a single plot. The MATLAB commands below plot both the sine and cosine curve on the same plot.
x=-2*pi:0.1:2*pi; y=sin(x); z=cos(x); plot(x,y,x,z);
theta = 0 : 0.05 : 2*pi;
hold on
axis('square')
plot(cos(theta),sin(theta))
When you are finished using the image, type hold off. To
see the effect of the axis('square') statement, enter the
MATLAB command axis('normal') and redraw the graph.
theta = 0 : 0.1 : 2*pi; r = sin(3*t); plot(r .* cos(theta), r .* sin(theta))Note the componentwise multiplication used in this problem, indicated by the dot-asterisk (.*) notation.
As an optional exercise (I'm a mathematician after all) guess what happens if the sin(3*t) term is changed to sin(4*t). Then test your guess by creating the graph.
To create a 3-D graph we need to create three matrices. The first gives the x coordinates, the second gives the y coordinates and the third gives the z coordinates.
These matrices can be created in a series of steps. As an example,
let's consider the graph of the function
z = (2 sin(3x))/( (x^2+2)(y^2+1) ).
for values of x and y between -3 and 3.
The steps for achieving this are:
The following MATLAB commands carry out these three steps for the example we are considering.
x = -3:0.2:3; y = -3:0.2:3; [X Y] = meshgrid(x,y); Z = 2*sin(3*X)./( (X.^2 + 2).*(Y.^2 + 1) ); surf(X,Y,Z)Try these out for yourself and see the results. Note how the formula for Z is related to the mathematics with which we started. See if you can explain why each period in the Z formula must be there.
If you prefer contour maps to surface maps, try the contour command:
contour(X,Y,Z)and a graph of the contour lines for the function will eventually appear.
One of the current social issues receiving a fair amount of publicity in the media is the question of the age distribution of the human population in the United States. The usual context for this problem is the question of the future solvency of the social security system, but there are many other important questions dealing with future populations and the age structure of that population. (The need for new elementary, middle and high schools, for example.)
In this section of the tutorial you will build a simple model which forecasts an age structured population many years into the future. The assumptions are based on the premise that certain demographic parameters remain constant over the period of the time being studied. As the tutorial progresses, I will develop an example of this type of model called the Leslie Matrix model while encouraging you to develop a similar model on your own.
The main idea of these models is to split a population into several distinct age classes and model how individuals are created and subsequently move from one age class to the next. In the Leslie model every member of one age class either dies or moves into the next age class at every iteration of the model. In the model you will be constructing individuals will experience one of three possibilities: death, stasis (remaining in the same age class) and maturation (moving on to the next age class).
Focusing for awhile on the Leslie model, there are two parameters associated with each age class. The percent survival rate for the age class represents the proportion of the individuals in an age class which survive to move into the next age class during the next iteration of the model. The second parameter is the fecundity parameter which represents how many offspring (age class zero) are produced, on average, by each member of the age class. (Some argue convincingly that fecundity is most closely related to the number of females in the population and thus create Leslie models which count only the females in a population.)
To continue our development of the Leslie Model we require an
introduction of some mathematical notation. Suppose P(i)
represents the population of age class i. If s(i)
represents the survival rate of age class i, then the
population in age class i+1 in the next iteration can be
computed using the current population in age class i through
the formula
next iteration's P(i+1) = s(i)*P(i).
If we let N be the last age class in our model (all
individuals in age class N die before the next iteration)
then the formula above works for values of i from 0 up to
N-1, inclusive. This still leaves us without a population
for the next iteration's youngest age class, age class 0.
The fecundity parameters are used to determine how many newborns are
created during the next iteration of the model. If we let
f(i) be the fecundity (average birth rate) per individual in
age class i, then the total number of newborns will be:
f(0)*P(0) + f(1)*P(1) + f(2)*P(2) + ... + f(N)*P(N)
If P is a column vector of population by age class, then it is possible to compute the P vector for the next iteration through a simple matrix (non-componentwise) multiplication. In particular, if A is the (N+1)x(N+1) matrix whose first row is the fecundity rates f(i) and whose diagonal below the main diagonal is made up of the survival rates s(i), then the next iteration's population column vector can be found by performing the matrix multiplication A*P.
As an example, suppose we have a population of critters who never grow past five years old. Suppose, further, that 80% of the critters survive each year to move onto the next age class, with the exception of the five year old critters which always die. If newborns and yearlings have no offspring, and if middle aged critters (2, 3 and 4 year olds) have, on average, 0.35 newborns each year, while the oldest critters (5 year olds) have only 0.1 newborns on average, how can we use current populations to forecast future populations and the age structure of that population?
In this model the parameters are given in the table below:
Age Class Survival Rate Birthrate
i s(i) f(i)
0 0.8 0.00
1 0.8 0.00
2 0.8 0.35
3 0.8 0.35
4 0.8 0.35
5 --- 0.10
To create the matrix A mentioned above for this particular set of parameters, the MATLAB command below would do the trick. (Note the elipses at the end of some lines. This is how we enter a MATLAB command which takes more than one line to complete.)
A = [ 0.00 0.00 0.35 0.35 0.35 0.10; 0.8 0 0 0 0 0; ... 0 0.8 0 0 0 0; 0 0 0.8 0 0 0; 0 0 0 0.8 0 0; ... 0 0 0 0 0.8 0]which creates the matrix A shown below:
A = 0.0000 0.0000 0.3500 0.3500 0.3500 0.1000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000
To forecast the future age structured population of our colony of critters we need only determine the current age structured population, represent this population as a column vector P and then move forward one year at a time in our model by performing the multiplication A*P.
For example, if our initial population consists of sixty critters evenly spread throughout the six age classes, our initial P vector would be:
P = [10;10;10;10;10;10]Note that the semicolons cause P to be a column vector. We could have achieved the same end by separating the 10's with spaces or commas, creating a row vector, and then using the matrix transform operator (an apostrophe).
To determine the population distribution next year we multiply A and P giving the vector
11.5000 8.0000 8.0000 8.0000 8.0000 8.0000If we want the population two years in the future there are a number of ways to proceed. First, we can take the vector above and multiply it by A to get the population for the next iteration. Proceeding in this fashion we get a series of vectors, each of which gives the population for the succeeding iteration.
If we only want populations for specific future years a second way may make more sense. This method is based on noting that for each year into the future we go, we arrive at the population by multiplying by the matrix A. Thus, if we want to go ten years into the future we get there by multiplying by A ten times. We can simplify this process by using the (noncomponentwise) matrix exponentiation operator carat (^).
To move ten years into the future we start with the current population (held in the column vector P) and multiply it by A ten times. This could be accomplished with the MATLAB command
A*A*A*A*A*A*A*A*A*A*PThere is an easier way however. Note that in the computation above, we multiply the matrix A by itself 10 times. (This is traditional matrix multiplication, not componentwise multiplication.) We can use the exponention operator (^) as a way to indicate repeated matrix multiplications without writing them all out. Using matrix exponentiation, the MATLAB command given above can be recast as:
A^10*PIt is important to note that the carat (^) without a leading dot indicates repeated traditional matrix multiplication. This contrasts with the dot-carat (.^) version of the exponentiation operator which would work in a componentwise fashion, raising each entry in a matrix to the indicated power.
In our example, the command A^10*P yields the result:
3.0288 2.7583 2.4885 2.3087 2.1601 1.9545a result which leads us to be very concerned about the long term viability of this population.
ones(1,6)*(a^10*p)
To investigate this question we must update the entries in the A matrix to reflect the hypothetical introduction of the new drug. Only three entries inthe matrix A need be changed, so rather than recreate A from scratch we will make the three individual changes.
The first change is that the birth rate for two year olds must be changed to 2.3. Since the birth rate for two year olds is found in the first row, third column of the matrix A, this change can be effected by the MATLAB command:
A(1,3)=2.3;
Make the remaining two changes in A, changes which reflected the lowered survival rates for age classes 3 and 4.
Begin with the same initial population (all 10's) as before. Does the population seem healthier (whatever that means) under the drug regime? How does the drug effect the age distribution of the population?
If fecundity rates remain unchanged, what survival rates are needed for our critter population to remain stable? What proportion of a stable population falls into each age class?
It is fun to see how populations evolve over time. Toward this end it is possible to create a three dimensional graph showing the evolution of a Leslie population over time by creating a matrix whose columns are "snapshots" of the population at various times under study.
Suppose P is the initial population (all 10's in our example). We will use the matrix TimePop to be the matrix whose columns will contain the population at given points in time. Since the first population we are interested in is the initial population, we will set the first column of TimePop to the initial population using the MATLAB command
TimePop=P;We will use a second matrix called PP to hold the most recent population computed. Thus, PP should start off holding the initial population
PP=P;The task now facing us is to generate the population vector for the next time period, add this column vector to the matrix TimePop and then do the same again, moving ahead yet another time step. This can be achieved by executing the following MATLAB commands repeatedly:
PP = A * PP; TimePop = [TimePop PP];The first command moves the population vector PP forward one time step. The second command adds a new column to the TimePop matrix. This new column contains the most recently computed population.
Using the up-arrow key, repeat the line above 15 to 20 times, creating a TimePop matrix with many population snapshots.
The evolution of the critter population over time can be viewed by creating a three dimensional surface plot of the array TimePop. To create this surface plot enter the following MATLAB commands and note the plot which is generated. Pay special attention to the creation of labels and the colorbar.
surf(TimePop)
view(0,90)
colormap(jet)
colorbar
xlabel('Year')
ylabel('Age Class')
The colon (:) notation we discussed earlier is a convenient way to create large arrays where the entries are evenly spaced. A second use of the notation is to extract parts of existing matrices. In particular, a single colon is a reference to every row or every column in a matrix, depending on its position.
For example, the notation A(2,:) refers to every column of the second row of the matrix A, or, said another way, A(2,:) is the entire second row of A.
Similarly, B(:,5) refers to the entire fifth column of the matrix B.
plot(TimePop(3,:))
Using calculators it is difficult to efficiently deal with Leslie models consisting of more than about ten age classes. Using MATLAB, however, it is relatively easy to deal with populations with dozens of age classes. Human populations are generally modeled with 101 age classes, allowing for ages between 0 and 100. It is also much easier to move far into the future, spotting long term trends.
In this section we discuss modifications to the Leslie Population Model which describe a slightly different situation than that modeled by Leslie. This section introduces no new MATLAB skills and is present for those who want to undertake the challenge of creating a new model. I will describe the mathematics of the new model, but will leave it entirely up to you to implement the model in MATLAB. A knowledge of traditional matrix multplication, the non-componentwise variety, will be helpful.
In the Leslie model only two things can happen to individuals in an age class; they can die or they can move onto the next age class. In the model presented here we allow for a third option, remaining alive and in the same age class.
Given that it is possible to remain in the same class for several time steps, it is probably best to leave off the word age when describing these classes.
As far as I know, this version of the Leslie model has no well-known name, so we get to name the model for the purposes of this tutorial. Feeling singularly uncreative at the moment, I choose to call it Model X.
Model X requires three sets of parameters. Denote the fecundity of class i with the symbol f(i), the proportion of class i moving on to the next class at each iteration with the symbol s(i) and the proportion of each class which survives but remains behind with the symbol r(i).
Proceeding in a fashion reminiscent of the Leslie model development, if we create a model with classes numbered from zero to N, and let P(i) denote the number of individuals in class i, then we can derive the following equations which give the values of the various P(i) for the next time step in terms of the P(i) values for the current time step.
Considering, first, the values of P(i) in the next time step
for classes numbered from one to N we see:
Next period's P(i+1) = s(i)*P(i) + r(i+1)*P(i+1)
This equation tells us that next year's class i+1 is made up
of those individuals from this year's class i who move up a
class as well as those individuals from this year's class i+1
who survive but do not advance in class. (a.k.a. flunking?) This
equation holds for values of i from zero to N-1
inclusive.
A careful look at this situation reveals that we have not yet specified the population of class zero in next year's population. Recall that members of class zero are either true newborns or previous members of class zero who do not advance for some reason.
The mathematical expression which gives the number of class zero
residents in next year's population is given by:
Next period's P(0)=
(r(0)+f(0))*P(0) + f(1)*P(1) + f(2)*P(2) + ... + f(N)*P(N)
Pay special attention to the first term of this equation. It does not
follow the pattern of the other terms.
Model X has a structure similar to the Leslie model. In fact, if P is a column vector whose N+1 entries are the populations of each class at some point in time, then there is an (N+1)x(N+1) matrix A so that the population next year can be determined by calculating the traditional matrix product
A*P