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 phoenixAfter 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 gutenbergDon'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 matlabNow 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.mNow we are ready to run MATLAB from the phoenix window. This is done by typing the command ``matlab'' at the prompt
matlabYou 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
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 0Alternatively, 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 0To 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.1111To 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.0000The ``det'' command is used to compute the determinant of a square matrix:
>> det(A) ans = 27Keep 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 23Now 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.3333We can doublecheck by multiplying Ax and comparing with b, defined above.
>> A*x ans = 1 2 0To 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.7345X 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
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.mA 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
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).
plot(y(:,1),y(:,2))and you should see the curve spiralling into a point.
As our last example, we consider the Van der Pol equation, which is the second order equation
x'' + (x*x-1) x' + x = 0which can be rewritten as a system of two first order differential equations
x1' = x1(1-x2*x2)-x2 x2' = x1At 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
help printYou 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.psor if you have set the printer environment variable, you should be able to just use the command
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 quiverAlthough 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.