Matlab tutorial


What follows below is a short and very basic practical introduction to
matlab/octave. You can find a great deal more information about
matlab or octave on the web. Octave is a free/openware (GNU) version of
matlab: it is not completely identical with matlab and its graphing
facilities are not as fancy, but it works well and is free (both
price wise and following the free/openware philosophy).

This document was created by Sohrab Ismail-Beigi in November, 2003,
primarily for use in the ENAS 856a/PHYS 650 course.

This web page was last modified on

var lm=document.lastModified;
document.write(” ” + lm + “.”);


The basics

Matlab is a program that is geared for performing numerical
calculations in an easy and straightforward way. It is text-based so
you type in commands, hit return, and this creates output or plots.
For example, if you type 2+2 and hit return, you get 4. If you
type cos(pi) you get -1.

Operators

The usual mathematical operations of interest are +, -, *
(multiply), / (divide), and ^ (exponentiation). Whether these
act on numbers or vectors or matrices depends on the case (see below).
So, just to be clear, 2^3=8, (3+4)/14=0.5, etc.

Variables

You will have noticed that matlab spits out an ans= before giving
you the answer in the two examples above. This means it is printing
the value of the variable ans which holds the answer to the
calculation. If you type the name of a variable by itself and hit
return, matlab prints its value (e.g. type ans and see its value).

We can create and assign variables ourselves. For example, typing
x=2*3+4 will result in x=10. From now on, the variable x is
defined and currently has value 10. If you type x by itself, matlab
will print its value.

To see a list of variables currently defined, use the who command.
The whos command is a longer version of the same thing showing the
type of variable and its size.

Vectors and matrices

We are not limited to simple numbers. We can define vectors or
matrices easily. For example, v=[1 2 3] sets the variable v to
be a row vector of size 3 whereas w=[1 2 ; 3 4] creates a 2x2
matrix. z=[5 6 7]' creates a 3D column-vector (the prime
denotes Hermitian transpose). The creation of w and z could also have
been accomplished with the longer but more clear:

  w=[1 2
     3 4]

and

  z=[5
     6
     7]

We can now take a dot product of v and z by typing
v*z. We can do matrix-vector multiplication with w*[1 1]' or
[1 1]*w (they are different!)

Complex numbers and complex conjugation

If you do not define a variable named i, then it is by default the
unit imaginary number, sqrt(-1). Creating complex numbers is easy:

  z = 2+3*i

does the trick. Complex conjugation can be done in two ways: either
by using the function conj, or by using the Hermitian conjugation
operator (single quote):

   z = 2+3*i;
   conj(z)
   
     ans = 
            2.000 - 3.000i
   z'
     ans = 
            2.000 - 3.000i

Transposition and Hermitian conjugation

For real matrices these are the same and most easily accomplished
with the the operator (single quote): if M is a real matrix,
then M' is its transpose.

For complex matrices, the transpose is given by the transpose
function; Hermitian conjugation by the operator, as above. Thus
we have that M' = transpose(conj(M)).

Reducing and paging matlab screen output

If your variables are large, you probably do not want matlab to print
the values to the screen. You can silence matlab by putting a
semicolon ; at the end of a command: e.g. v=[1 2 3]; will set
the variable v quietly. A more useful example would be
a=linspace(0,1,100); which sets a to be a linearly spaced vector
of 100 numbers from 0 to 1 (inclusive).

If you want to see the output but it has scrolled by too fast, then
you need to turn on paging with the more command: more on. When
paging is on, each full page of text shows a --more-- and waits for
you to press a key: press space for another page, return for another
line, and q to quit the output. You can turn paging off with
more off.

Aborting a command

You can always abort any matlab command by pressing control-c (the
control or “Ctrl” key and the “c” key at the same time).

Creating a program file: .m files

When you type a name that is not a built-in command or variable,
matlab looks at the current directory for any program files that have
the name you typed plus a “.m” added at the end and will execute them
if found. You can find the current directory by using the pwd
command and change directories with the cd.

A matlab .m file is nothing more that a list of commands put into a
file, one line at a time. When you type the file’s name, the commands
get executed. As an example, let us say we have created a file called
myfirst.m containing the four lines

  clear
  a=2;
  b=3;
  c=a+b;

The clear command clears out any existing variables in memory.
When you type myfirst at the matlab prompt, nothing seems to
happen, but you can check that the variable c now contains the
value 5 and that a and b are defined (with the who command).

Getting help

The matlab help command is a gold mine of information. Typing
help alone gives you a list of matlab documentation
subjects. You get help on a subject by typing help
subject
(e.g. help general).
Typing help followed by a command name produces the
command documentation (e.g. help more).


Some useful mathematical operations

Creating a vector of equally spaced numbers

This can be done with loops or with the colon operator. For
example, to create a list (row-vector) of integers staring at 3 and
going up to 12 with spacing of 3 you coul
type [3:3:12] or [3:3:14] (the
second possibity is based on the fact that you start at 3 and go in
steps of 3 until you are above 14). Another very convenient way is
with the linspace command:

   x = linspace(x1,x2,n);

which returns a row vector with n equally spaced entries with the
first entry being x1 and the last one being x2 (so the spacing is
(x2-x1)/(n-1)).

Selecting parts of vectors and matrices: the colon operator

Sometimes we want to access or manipulate the parts of a vector or
matrix, usually a range of indices. For example, if v is a vector
with 10 elements, then v(2:5) is a vector with four elements
composed of the entries 2, 3, 4, and 5 of v (in that order). In
reality, the notation 2:5 creates a vector of numbers starting at 2
and ending at 5 (with increment 1 implied) and then we use this vector
to index v. To specify a different increment, put it between the
start and end of the list: 2:2:6 creates the vector [2 4 6] and
10:-1:1 creates [10 9 8 7 6 5 4 3 2 1].

The colon operator can be used on either side of an assignment:

   v = 1:10;
   w = v(2:5);
   v(1:4) = w;

Sections of matrices work the same way, but we specify both entries.
A : without any range means the whole row or columns. Therefore

   M(2:5,1:3)   % a 4 x 3 submatrix of M
   M(:,4)       % 4th column of M
   M(2:4,:)     % the second through fourth rows of M as a submatrix

Creating identity and zero matrices

The function eye(n) returns a n x n identity matrix. The function
zeros(n) returns a n x n matrix of zeros; zeros(m,n) returns a
m x n matrix of zeros.

Creating a vector or matrix full of ones

We may need a vector or matrix filled with the number one (1). The ones
function is what we need. ones(n,m) returns a n x m matrix filled
with ones. If m=1 then it is a column vector (n=1 is a row vector).

Creating a diagonal matrix

The diag function creates a diagonal matrix given a vector
argument, and puts the entries of the vector on the diagonal of the
resulting matrix. Thus diag([1 2 3]) returns a 3 x 3 diagonal
matrix with 1, 2, and 3 on the diagonal:

   diag([1 2 3])
     ans =
     1     0     0
     0     2     0
     0     0     3

Creating off diagonal/triangular matrices

Often, we wants to fill the non-diagonal entries in a matrix in a
rather regular way. For example, we want to fill the diagonal above
or below the main diagonal of a matrix with the same numbers. The
diag() function performs this for us when we give it an extra
option. For example, let us take a 5x5 diagonal matrix with the
numbers [10 20 30 40 50] on its diagonal and make its upper off
diagonal have the numbers [-3 -4 -5 -6]:

  M = diag([10 20 30 40 50]) + diag([-3 -4 -5 -6],1)
  M =
    10    -3     0     0     0
     0    20    -4     0     0
     0     0    30    -5     0
     0     0     0    40    -6
     0     0     0     0    50

We created a diagonal matrix and added another matrix to it with the
appropriate off diagonal entries. The trailing “1” option to
diag() means the upper diagonal; “-1” would have meant lower
diagonal. See help diag for more details.

If instead we wanted to fill the off diagonal with the same number,
say -0.5, one way to do so is to use the ones() function discussed
above: we create a vector with four unit entries with ones(4,1)
which we then scale by -0.5 and feed to diag():

  M = diag([10 20 30 40 50]) + diag(-0.5*ones(4,1),1)
  M =
   10.0000   -0.5000         0         0         0
         0   20.0000   -0.5000         0         0
         0         0   30.0000   -0.5000         0
         0         0         0   40.0000   -0.5000
         0         0         0         0   50.0000

Inverting a matrix

The inv function inverts matrices. So A = inv(B) makes A
the inverse of B so that A*B and B*A are both the identity matrix.

Diagonalizing a matrix

Given a square matrix stored in variable M, we find its
eigenvalues with the eig command:

  e = eig(M);

will put the eigenvalues of M into the column vector e. If
they are real, they are sorted in ascending order.

If the eigenvectors are also desired, then we call eig but change
the format of the output:

  [v,e] = eig(M);

will put the eigenvectors into the columns of v and the
eigenvalues into the diagonal of e: note that unlike the first use
of eig for eigenvalues where e was a vector, here e is a
square matrix which is diagonal with the eigenvalues on the diagonal.
To get a vector with the eigenvalues, use diag(e).

To be very clear, the relation of v, e, and M is
M = v*e*inv(v). If M is Hermitian, then inv(v)=v' so
v'*v and v*v' are both the identity matrix.


Some useful matlab/octave programming

Program or .m files

Read the above section on Creating a program file: .m files.
Programs are always longer, so it is good to put all the commands
into a file rather than having to type them over and over and…

Loops

At some point, everyone has to program a loop. To loop the variable
j over the numbers 1 to 5 (inclusive), and to print 2*j+1 for each
value of j, we use the for and end commands:

  for j=1:5
    2*j+1
  end

You can go backwards, but you have to give an extra argument: for
j=5:-1:-2
will loop j starting from 5 and ending at -2 by adding
-1 each time (thus in for j=1:5, the increment of 1 is implied).

For a more serious example, let us create a vector v of n elements
which contain the squares of 1 to n, and let us create a n x n matrix
M with some entries:

  for j=1:n
    v(j) = j^2;
  end
  for j=1:n
    for k=1:n
      M(j,k) = 2*j^2 + cos(k*j*pi/50);
    end
  end

Visualizing and plotting

Making a two-dimensional plot

At some point, you want to look at the data you’ve been working on so
hard. The most basic thing you can do is to create a standard two
dimensional plot: y as a function of x.

What follows below is mostly
specialized for matlab as octave has a slightly different notation.

If you have a vector of length n of numbers v, then issuing the plot
command

  plot(v)

will create a new graphical window and plot the values of v using
connected line segments. The values are the vertical axis (y), while
since no x data was given, the number 1 through n are used. If you
have your own x values, then

  plot(x,v)

will plot v versus x. This assumes the data are in corresponding
order: x(1) and v(1) go together, x(2) and v(2), etc.

Getting some grid lines

The grid command draws some grid lines with uniform spacing.
grid turns it on, grid off turns it off.

Plotting the rows/columns of a matrix

It is quite easy to plot columns or rows of a matrix M:

  plot(M(2,:))
  plot(M(:,7))

The first command plots the second row of M and the second command
plots the seventh column.

Superimposing multiple plots

You will notice that each time you issue the plot command, it
clears the screen and starts fresh. Sometimes you want to
superimpose two or more plots. This is accomplished with the hold
on
command: after you plot the first set of data, issue hold on.
This tells the program that future plots are to be superimposed:

  plot(x,y1)
  hold on
  plot(x,y2)

will plot y1 and y2 versus x on the same plot.

Clearing the plotting screen

If you make a mistake or want to start fresh, the clf command
will clear the screen. It will also turn hold off, so you have to
hold on again when needed.

Choosing symbols and colors for plots

The plot command takes many possible arguments. To plot y versus
x with red circle symbols, do

  plot(x,y,'ro')

The r stands for red while the o stands for an open circle. Run help plot for all the possibilities. The color comes first then the
symbol. If one of them is not specified, some default is used.

Some useful symbols are: o (circle), x (x-marks), +, *, s
(square), d (diamond), ^ (up triangle), - (solid line, the default), :
(dotted line), – (dashed line), -. (dash-dotted line). Colors are
one letter y(ellow), m(agenta), c(yan), r(ed), g(reen), b(lue),
w(hite), k(black).

Logarithmic plots

In principle, this is useless, since one could just use
plot(x,log10(y)) to plot the log base 10 of y versus x. In
practice, it is nice to have a log plot with the right type of grid
lines. To make the y axis logarithmic, use the semilogy command:
it is otherwise identical with the plot command as far as
arguments and usage are concerned. There are also the semilogx
and loglog commands.

Putting titles and labels on the plot

The title command puts a title on the plot:

  title('My favorite plot')

The xlabel and ylabel commands have the same format as title
and put labels on the x and y axes.

Printing

The print command alone will send the current plot to the default
printer. print -Pprintername will send the plot to the named
printer. To print to a file, supply print with a filename:

  print plotfilename

The default format for the file is probably postscript. You can
control the format with some options. Here are some self explanatory
examples:

  print -dps  file.ps
  print -deps file.eps
  print -dtiff file.tiff
  print -djpeg90 file.jpg    % creates jpeg with quality level 90
  print -dpng file.png

Run help print to see what the details are for the system you are
running on.


Examples

The best way to learn is to see some examples.

Plotting a function of x

Plot sin(x) from 0 to 10*pi with 100 points. Draw the plot and put
red circles where the data resides:

   clear
   x = linspace(0,10*pi,100);
   y = sin(x);
   clf
   plot(x,y)
   hold on
   plot(x,y,'ro')

Creating a tridiagonal matrix

Tridiagonal matrices occur often in math: they have the diagonal and
the elements immediately above and below the diagonal nonzero, and
the rest are zero. Here we create a 100 x 100 tridiagonal matrix:

   clear
   n = 100;
   M = zeros(n);
   v = 1:n;
   for j=1:n
     M(j,j) = v(j);
   end
   for j=1:n-1
     M(j,j+1) = v(j)^2;
     M(j+1,j) = 2*v(j)+1;
   end

Actually, this can be accomplished more easily with a special form of
the diag function which takes an extra argument and fills up a
particular diagonal (see help diag):

   clear
   n = 100;
   v = 1:n;
   M = diag(v) + diag(v(1:n-1).^2,1) + diag(2*v(1:n-1)+1,-1);

Here we choose the parts of v we want for the off diagonals, and the
.^2 squares every element of its argument individually.

To make sure we understand what is going on, here is what the 5x5
upper section of M should look like: after running the above code,
typing M(1:5,1:5) produces

   1   1   0   0   0
   3   2   4   0   0
   0   5   3   9   0
   0   0   7   4  16
   0   0   0   9   5

Diagonalizing a matrix and testing orthogonality of eigenvectors

Let us create a 10x10 Hermetian matrix with random complex entries.
The randn() function returns random numbers distributed according
to the normal distribution with mean zero and standard deviation
one. So we do

   M = randn(10,10)+i*randn(10,10);

This gives us a 10x10 matrix of random complex numbers which is *NOT*
Hermetian. To make it so, we can symmetrize it by averaging it with
its Hermetian conjugate M':

   H = (M+M')/2;

H is Hermetian. We diagonalize it to get the eigenvectors v
and eigenvalues e with the command

  [v,e] = eig(H);

Both v and e are 10x10 matrices. v has the eigenvectors in
the columns, and e has the eigenvalues on the diagonal of the
matrix. You can check that e is diagonal with real entries. You
can extract the diagonal with diag(e).

Let us check that the eigenvectors of the Hermetian matrix are
orthonormal. For example, let us check that the first and second
vectors are orthogonal:

   v(:,1)'*v(:,2)
   ans=
      1.0786e-16 - 8.8048e-17i

The way to read this is that we are taking the first and second
columns of the matrix v and taking the complex scalar product: the
prime ' operator will take the column vector v(:,1) and
conjugate-transpose it into a row vector which when multiplied by
v(:,2) results in a number (scalar), which is quite close to zero.

We can check the orthogonality of all vectors at once by computing
all mutual dot products by the condensed formula

  v'*v

which returns a 10x10 matrix of dot products, which should be close
to the identity matrix.

Calculating an outer product

Let us say someone has given you a set of vectors in the columns of a
matrix v and you want to compute the “outer product” of the first
n columns. That is, you want to create a matrix O so that the
i,j entry of O will be the following sum


          n
O(i,j) = Sum { v(i,k)*conjugate(v(j,k)) }
         k=1

(Note that this is not a matlab line above.) We could accomplish
this by writing three nested for loops over i,j,k, but that is
rather inefficient in matlab. Rather, consider the following short
matlab expression that computes the same thing:

  O = v(:,1:n)*v(:,1:n)'

We take v and form a new matrix out of its first n columns by
writing v(:,1:n) (the first colon means we take all rows). We
then multiply this result by its own Hermetian conjugate on the
right. You can check that the entries of O are given as per the
above formula.

Take home message: Most sums and products can be written in
matrix and evaluated in a very efficient and compact way in matlab.
The most common ingredients are inner and outer products.

Plotting a vector on a log plot

We have a vector v and want to plot its entries on a natural log plot.
One way to do so is to do

  plot(log(v))

This is fine. However, the semilogy function is specialized for
log plots and makes the axes and grids work in a nice way:

  semilogy(v)