- The basics
- Operators
- Variables
- Vectors and matrices
- Complex numbers and complex conjugation
- Transposition and Hermitian conjugation
- Reducing and paging matlab screen output
- Aborting a command
- Creating a program file: .m files
- Getting help
- Some useful mathematical operations
- Creating a vector of equally spaced numbers
- Selecting parts of vectors and matrices: the colon operator
- Creating identity and zero matrices
- Creating a vector or matrix full of ones
- Creating a diagonal matrix
- Creating off diagonal/triangular matrices
- Inverting a matrix
- Diagonalizing a matrix
- Some useful matlab/octave programming
- Visualizing and plotting
- Examples

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.*

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`

.

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.

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.

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!)

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

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))`

.

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`

.

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

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).

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`

).

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)).

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

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.

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).

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

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

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.

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.

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

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

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.

The **grid** command draws some grid lines with uniform spacing.
`grid`

turns it on, `grid off`

turns it off.

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.

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.

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.

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).

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.

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.

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.

The best way to learn is to see some examples.

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')

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

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.

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.

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)