prepared by Mark E. Johnson, maejohns@math.princeton.edu

MAE 305

Mathematics in Engineering I


MATLAB Tutorial

An outline of this document follows:
Introduction
Getting Started
Fundamentals
Numerical Integration
Printing
Help
Nice Demo
Plotting Direction Fields

Introduction

MATLAB is an interactive computing environment for the numerical study of matrices, or more generally, systems of equations ; it can be used for numerical analysis, matrix computation and graphics. The name MATLAB comes from MATrix LABoratory and the main point of emphasis is that the basic element used in MATLAB is a matrix. Along the way, we will learn how to do some very useful things with MATLAB (invert matrices, solve systems of equations, compute eigenvalues and eigenvectors, etc.), but let's say that our goal for this tutorial will be to use MATLAB to (a) numerically solve differential equations and to (b) plot the solutions and produce hardcopies of the results.

Getting Started

Before we can do anything on the computer, we need to know how to use the computer. The beauty of common sense. In order for you to be reading this, you must have known how to log into a computer and how to "surf the web" to this page. So I will begin by assuming that you are sitting at a Silicon Graphics computer in Engineering Quad, Room E-423 and that you have logged in.

If you try to run MATLAB (by typing "matlab" in a shell window), you will fail. This is because the University has purchased a MATLAB software license for only a selection of computers. More information on the availability of MATLAB and further details of the license at Princeton can be found here.

That's OK because we can remotely log into one of the licensed machines. First you need to create a new window. To create a new window, go to the ``Tools'' bar in the toolbox and find the ``Shell'' option (for some of you, you may need to look for the ``Unix Shell'' option unser the ``Desktop'' button), OR hold down the right mouse button while the pointer is inside your console widow, drag the mouse to the ``clone'' menu and choose a size for your new shell ; release the mouse button.

You should now have two windows. In one of the two windows, log into the machine called phoenix by typing

rlogin phoenix
After entering your password, set the DISPLAY environment variable to the name of the machine at which you are working --- this tells phoenix where to open new windows. (You can determine the name of the machine you are working on my typing the ``hostname'' command in the other window.) While we are at it, set the printer to the name of the printer in Room E423, gutenberg
setenv DISPLAY YourMachineName:0
setenv PRINTER gutenberg
Don't forget to type the ":0"

Now, some of you (particularly those who have taken CivE 201 in the past) might have trouble with the previous commands. If you have problems, then try first issuing the command ``csh'', then trying the setenv commands. If this works, then you will have to type this each time you plan on setting the DISPLAY variable. If you are still having problems, then please ask me for further help!

After doing all of this, go to the other window and type ``xhost + phoenix'' to give permission to phoenix to open windows on your machine.

xhost + phoenix
phoenix being added to access control list

FYI, phoenix is not the only machine with the MATLAB license at Princeton. In fact, any of the arizona and any of the sesamest machines can be used in place of phoenix.

We are now at a point where we can begin using MATLAB. Before we do so, I would like you to create a subdirectory on your computer account called "matlab". Do this by issuing the following command in either window:

mkdir matlab
Now copy two files from my matlab directory to your matlab directory. Note that my username on phoenix is "mej" and not "maejohns"
cp ~mej/matlab/de1.m ~/matlab/de1.m
cp ~mej/matlab/de2.m ~/matlab/de2.m
cp ~mej/matlab/sol1.m ~/matlab/sol1.m
cp ~mej/matlab/sol2.m ~/matlab/sol2.m
Now we are ready to run MATLAB from the phoenix window. This is done by typing the command ``matlab'' at the prompt
matlab
You should see the following output in your window:

                            < M A T L A B (R) >
                (c) Copyright 1984-94 The MathWorks, Inc.
                            All Rights Reserved
                               Version 4.2c
                                Nov 23 1994

Commands to get started: intro, demo, help help
Commands for more information: help, whatsnew, info, subscribe
 
>>
A matlab graphics window

should also briefly appear then disappear. If you do not see this, then you probably have not set your DISPLAY variable correctly. Try to set it again or ask me for some help.


Fundamentals

The essential object used by MATLAB are matrices. There are a number of ways to enter a matrix, let's start by defining a matrix by entering this statement:

>> A = [1 2 3; 4 5 6; 7 8 0]

A =

     1     2     3
     4     5     6
     7     8     0


Alternatively, we could have entered the matrix by typing
>> A = [ 1 2 3
         4 5 6
         7 8 0 ]

which would produce the same output as above. Note that MATLAB is case-sensitive in that the name A is different from the name a. Before continuing, let's enter a vector b for future use
>> b = [1; 2; 0]
We can suppress the output by following the command with a semicolon (Mathematica has this same feature)
>> b =  [1; 2; 0] ;
Here is as good a place as any to let you know that the arrow keys can be used to edit previous commands. Try pushing the "up" arrow to see what happens, then try pushing the left and right arrows, then the down arrow, etc.... You get the picture.

Having defined the matrix A and the vector b, we can now perform some elementary operations.

We start with finding the transpose B=A' by using the special character ('):

>> B = A'

B =

     1     4     7
     2     5     8
     3     6     0
To find the inverse of A, call it C=1/A, we issue the ``inv'' command.
 
>> C = inv(A)
C =

   -1.7778    0.8889   -0.1111
    1.5556   -0.7778    0.2222
   -0.1111    0.2222   -0.1111
To confirm that this is indeed the inverse, multiply C and A to (hopefully) obtain the identity matrix.
>> C*A

ans =

    1.0000    0.0000         0
   -0.0000    1.0000         0
    0.0000    0.0000    1.0000
The ``det'' command is used to compute the determinant of a square matrix:
>> det(A)

ans =

    27
Keep in mind that in order to perform the following elementary operations, the dimensions of the matrices must be consistent with the definition of the operation (e.g. A+B is well-defined only if A and B are of the same dimension and AB is defined only when the number of columns of A equals the number of rows of B).
>> A + B

ans =

     2     6    10
     6    10    14
    10    14     0

>> 5*A

ans =

     5    10    15
    20    25    30
    35    40     0

>> A * B

ans =

    14    32    23
    32    77    68
    23    68   113

>> A * b

ans =

     5
    14
    23
Now suppose we would like to solve the linear system Ax=b, where x is the unknown vector, b is a known vector and A is a known matrix. If the inverse of A exists, then the solution can be written as x=(1/A)b. In MATLAB, this operation (left multiplication by the inverse) is done using the ``backslash'' character.
>> x = A\b

x =

         0
         0
    0.3333
We can doublecheck by multiplying Ax and comparing with b, defined above.

>> A*x

ans =

     1
     2
     0
To compute the eigenvalues of a matrix A, the command ``eig(A)'' does the trick, displaying all eigenvalues as output. The ``eig'' command can also be used to compute the corresponding eigenvectors by using the following syntax:

>> [X,D]=eig(A)

X =

    0.7471   -0.2998   -0.2763
   -0.6582   -0.7075   -0.3884
    0.0931   -0.6400    0.8791


D =

   -0.3884         0         0
         0   12.1229         0
         0         0   -5.7345
X and D are thus matrices : D is a diagonal matrix with the eigenvalues along the diagonal and X contains the eigenvectors in each of its columns. Mathematically, this can be written as
AX = XD.

Numerical Integration

In addition to matrix operations, MATLAB can be used to numerically determine solutions to differential equations. For MAE305, this is its most important use. The MATLAB command ``ode23'' is a 2nd/3rd order Runge-Kutta method, which is described in Section 8.4 of the text by Boyce and DiPrima.

Example 1

As our first example, we consider the first order differential equation
dy/dx = A * (sin(x)-y), y(0)=B.
For some fixed value of the parameter A and for some initial condition B, we would like to plot the solution y(x) in some interval of the independent variable x.

In order to do this, we first must create a matlab file (a ".m", said "dot m") which contains the definition of the differential equation. Remember a little while back you copied the files called "de1.m" and "de2.m" to your matlab directory? Well this differential equation is defined in the file called "de1.m", so let's take a look at it. The easiest way to view and edit files on a Silicon Graphics computer is to use the jot command. So, in your other window, type

jot ~/matlab/de1.m
A window containing the following should appear:
function yprime = de1(x,y)
global A;
yprime = A * (sin(x)-y);
We do not need to do anything to this file - I just wanted you to look at it. There are three lines: the first one defines the name of the function, the second declares the variable A to be global (which means we can change it in matlab and it will be changed in this function, too -- if you don't understand, don't worry about it now) and the third line defines the actual differential equation as seen above.

Now we want to solve this differential equation. This is what we do from within matlab: First we need to declare A a global variable in MATLAB by typing

global A;
Then initialize A as 4 by typing
A=4;
Next we define the limits of integration (where is x going to start and where is it going to end?)
x0 = 0; xf = 20;
Now what is our initial value of y (i.e. what is B) ?
y0 = 12 ;
Now we can integrate using the MATLAB command ode23
[x,y] = ode23('de1',x0,xf,y0);
The ode23 command fills the two arrays "x" and "y" on the left hand side of the above expression with values along the solution. If you just type "x", you should see a whole bunch of numbers which range from x0 to xf. If you type "y", then you should see a list of numbers from the solution curve In order to make sense out of all of these numbers, we can finally plot the solution curve with the command
plot(x,y)
Now on your own try new values of A and B to see what happens. To save yourself a little effort, look at the file called sol1.m by using the jot editor. All of the above commands are already typed into this file for you ; all you need to do is change the various values, save the file and type sol1 from within matlab. A new plot should appear.

Example 2

As our second example, we use a system of two first order differential equations
dy/dx = z
dz/dx = -A*z - C*C*y+Sin(x), y(0)=a, z(0)=b.
with two parameters. This equation is defined in the de2.m and looks as follows:
function yprime = de2(x,y)
global A;global C
yprime=zeros(2,1);
yprime(1) = y(2);
yprime(2) = -A*y(2) - C*C*y(1)+sin(x);
The only real difference is that is Example 1, we were solving one equation, whereas here we are solving two coupled equations. In order to do this, we use arrays which are 2-by-1. The third line initializes yprime to be a 2-by-1 array filled with zeros. The equation is then defined in the following two lines.

To integrate, we do pretty much the same thing as above ; the only difference is that the initial condition must now by a 2-by-1 array instead of just a number.

global A;global C;
A=4;C=3;
x0 = 0; xf = 20;
y0 = [12; -5] ;
[x,y] = ode23('de2',x0,xf,y0);
Now we have two options for plotting. The first option is performed just as in Example 1
plot(x,y)
and produces a plot with two curves : one for y1(x) and one for y2(x) (which is the same as y(x) and z(x) from the original equation).
To plot the phase portrait (y2 versus y1), try a different version of the plot command
plot(y(:,1),y(:,2))
and you should see the curve spiralling into a point.
Just as in Example 1, you can edit the file sol2.m to make some changes to the values of A, C, a, and b. Give it a try and see what happens. Also experiment with different limits of integration. You might also try to modify the equation in de2.m to see what happens....

Example 3

We will not cover this example during the tutorial session, but you might want to try it on your own. The file vdpol.m can be found in my directory ~mej/matlab/vdpol.m just as the others were. Enjoy!

As our last example, we consider the Van der Pol equation, which is the second order equation

 x'' + (x*x-1) x' + x = 0
which can be rewritten as a system of two first order differential equations
 x1' = x1(1-x2*x2)-x2
 x2' = x1
At this point, make sure that you have a file called vdpol.m in your matlab directory which contains the following text:
function xdot = vdpol(t,x)
xdot = zeros(2,1);
xdot(1) = x(1).*(1-x(2).^2)-x(2);
xdot(2) = x(1);
To simulate the differential equation defined in vdpol over the interval 0 <= t <= 20, invoke ode23
t0 = 0; tf = 20 ;
x0 = [0 0.25]'
[t,x] = ode23('vdpol', t0, tf, x0);
plot(t,x)
You should see a plot of x1 versus t superimposed with a plot of x2 versus t. To plot the phase portrait (x2 versus x1), try
plot(x(:,1),x(:,2))

Printing

Now that you have learned how to generate output in a window on your screen, it is time to learn how to print out these results on the printer (so you can save them forever! -- or hand them in for a grade! -- or whatever). This is actually quite simple in MATLAB. Try typing the command
help print
You should see a whole mess of stuff appear on your screen about the print option. When you type the command print output.ps, the contents of the current display window are written to the PostScript file called output.ps in your directory. There are a number of options to the print command, which I will not go into now, but feel free to ask about them or just to experiment with them on your own.

To print the file output.ps, go to you non-MATLAB window and type the ls command, which lists all of your files and directories in the current directory. If you see the file called output.ps, then go ahead and print it with the command

lpr -Pgutenberg output.ps
or if you have set the printer environment variable, you should be able to just use the command
lpr output.ps

There are many alternatives to this method of printing, so if you know of another, feel free to use it. Remember, the purpose of this tutorial is go get you started. You are now responsible to keep it going....


Help

Two commands which will come in handy when using MATLAB are the help and the lookfor commands. Try typing ``help eig'' to get a full description of what the eig command does. Now try typing ``help inverse''. You probably see a reply which resembles ``inverse not found''. A more general information search is provided by the lookfor command. Try typing
lookfor inverse

Nice Demo

Finally, everything described here plus a whole lot more (including cool graphics!) is explained in the MATLAB demonstration. To run it, type "demo" and enjoy!

Some Exercises to Try


MAE305ers: Return to Mathematica Tutorial


For your browsing pleasure, the MATLAB webpage can be found here.

Plotting Direction Fields

This section was added for the April 10th computer session. Since so many people were interested on how to plot direction fields using Matlab, I thought I would add this short little section describing how to go about it. The command you need to know is called quiver. The very first thing you should do is read the information provided to you through the help utility by typing the command
help quiver
Although this information may seem a little bit confusing and technical at first glance, I wanted to make sure that you all knew how to access it again in the future.

In order to use the quiver command, you need to supply it with four objects, which are labeled X,Y,U,V in the help page. For the purposes here, the variable X and Y will be vectors with the grid points in the x and y-directions, respectively. What I mean is that if we want to ultimately plot a direction field by drawing arrows at a bunch of points in a grid, then these vectors will contain their x and y coordinates. To clarify what I am saying, let's just continue. First we set

x = -2:.2:2 ;
y = -1:.15:2 ;
Now the two vectors contain sequences of numbers. The x sequence looks something like: -2, -1.8, -1.6, ... You can experiment with the stepsize and the limits of these sequences at your leisure. The point is that these two vectors define a collection of "mesh points" in the plane, where the coordinates are given by (x_i,y_j).

Now, assuming that you have defined a function (which comes from a 2nd order inhomogeneous problem!) (like vdpol from above), then you can run a command (which you will need to copy over from my directory by issuing the UNIX command "cp ~mej/matlab/dirfld.m ~/matlab/.") from my directory and executing it as follows: (MAKE SURE YOU HAVE ALSO COPIED THE VDPOL EXAMPLE WITH "cp ~mej/matlab/vdpol.m ~/matlab/."):

[px,py] = dirfld('vdpol',x,y);
This generates two matrices called px and py, which we can now send to quiver by doing
quiver(x,y,px,py)
For your information, you can make the arrows longer by scaling them (multiplying their length) by a number you add to the end of the argument list you pass to quiver. See below
quiver(x,y,px,py,1.7)
Voila! And that is how to make a direction field.

As a reminder, we discussed how to print in a previous section.