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


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)