MS554: Principles of Numerical Computing.
Lecture 1: Introduction to matrix calculations, matlab, and FORTRAN.

August 28, 2003. In this course you will be performing a variety of numerical computations, using both FORTRAN and matlab. FORTRAN is a standard for scientific computing, and you will need to write, compile, run, and debug FORTRAN code during your career as a numerically minded scientist or engineer. Matlab is a higher-level computational engine that is useful for completing numerical calculations, visualizing the results, and if you use the symbolic-math toolbox - can sometimes be used to compare numerical approximations to analytical results.

In this lecture I will present some of the basics of representing data in matlab, running matlab macros and functions, and using matlab to plot results. I have copied a lot of these notes from a class that I audited at U.Va.; George Hornberger and Pat Wiberg's Numerical Methods in Hydrology.

1. Computation with Matlab

In this, the first lecture of the class we introduce some of the capabilities of MATLAB and ask you to do a simple project using the program.

The basics: matrices and vectors

MATLAB operates on matrices, which are arrays of numbers. That is, a matrix is a "table" of numbers with, say, N rows and M columns. For example, if N=3, and M=3, a matrix A has dimension 3x2. Assignment of such a matrix in MATLAB can be accomplished as follows.

Type in a "[" to let MATLAB know that you are going to input a matrix. Type the numbers for the first row of the matrix (say "1" and "2') and then type a ";" to let MATLAB know that you have reached the end of a row. Type the next two rows in the same way.
>>A = [1 2; 3 4; 5 6]

After typing this in; MATLAB echoes with a matrix A:

A =
    1    2
    3    4
    5    6

A vector of numbers is simply a matrix with one of the dimensions equal to one. For example, you might have data on air temperatures taken at noon every day over a week. These can be represented as a vector in MATLAB:

>>T=[12.1 13.6 9.5 8.2 10.4 11.7 11.9]

Typing this line into MATLAB results in:

T=
    12.1000 13.6000 9.5000 8.2000 11.7000 11.9000

As defined, T is a "row vector". It can be converted to a column vector by taking its transpose, an action that is accomplished in MATLAB by placing a single quote after the matrix or vector. That is, to define the transpose of T, we type

>>T'

ans=
    12.1000
    13.6000
    9.5000
    8.2000
    11.7000
    11.9000

Setting Vectors and matrices with internal MATLAB commands.

Special vectors and matrices can be constructed with built-in MATLAB statements. For example, it often is useful to construct a vector that covers a certain range and has evenly-spaced entries. This would be the case where we want to examine a sequence of regularly spaced data where we did not explicitly have the independent variable recorded. If, for example, we had water levels recorded every 5 minutes, and wanted to plot the data versus a time (hour) variable, we can form such a variable as follows:

>>t=1:5/60:24;

Note that we put a semicolon at the end of this assignment statement. This tells MATLAB not to echo the 365-element vector of times.

Other MATLAB functions that are useful in setting vectors and matrices are "zeros" (sets a matrix with all zero elements), "ones" (sets a matrix with elements equal to 1), and "rand" and "randn" (sets a matrix with elements chosen using a pseudo-random-number generator for a uniform or normal distribution, respectively, on the interval [0,1]). Try typing "a=rand(3,4)" and "b=ones(1,10)" to get a feel of these utilities.
 

Simple functions and plots


MATLAB also has many built-in functions. Thus, we can form a sine function for the M2-portion of a site having a 4m tidal range using:

>>M2_PER=12.42;   % period of M2 Tide
>>M2_AMP=2.0; % approximate range of M2 Tide
>>M2_PHA=0.0;  % phase of tide; assuming zero here
>>MSL=10; % mean sea level in meters.
>>level=M2_AMP*sin(2*pi*t/M2_PER + M2_PHA) + MSL;

We can look at the result by using MATLAB's graphics capabilities:

>>plot(t, level,'-');
>>xlabel('Time Hours')
>>ylabel('Tidal Level in meters')
>>set(gca,'ylim',[0 25])

Try this and see the results.

{figure 1 here}

One of the most important things to note about MATLAB as you learn to use it is the extensive on-line HELP facility. Typing "help" by itself gives a series of choices from which you can refine your search. Typing "help item" gives help on the item chosen. For example, to get help with the MATLAB plot utility:

>>help plot

 PLOT   Linear plot.
    PLOT(X,Y) plots vector Y versus vector X. If X or Y is a matrix,
    then the vector is plotted versus the rows or columns of the matrix,
    whichever line up.  If X is a scalar and Y is a vector, length(Y)
    disconnected points are plotted.

    PLOT(Y) plots the columns of Y versus their index.
    If Y is complex, PLOT(Y) is equivalent to PLOT(real(Y),imag(Y)).
    In all other uses of PLOT, the imaginary part is ignored.

    Various line types, plot symbols and colors may be obtained with
    PLOT(X,Y,S) where S is a character string made from one element
    from any or all the following 3 columns:

           y     yellow        .     point              -     solid
           m     magenta       o     circle             :     dotted
           c     cyan          x     x-mark             -.    dashdot
           r     red           +     plus               --    dashed
           g     green         *     star
           b     blue          s     square
           w     white         d     diamond
           k     black         v     triangle (down)
                               ^     triangle (up)
                               <     triangle (left)
                               >     triangle (right)
                               p     pentagram
                               h     hexagram

    For example, PLOT(X,Y,'c+:') plots a cyan dotted line with a plus
    at each data point; PLOT(X,Y,'bd') plots blue diamond at each data
    point but does not draw any line.

    PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,...) combines the plots defined by
    the (X,Y,S) triples, where the X's and Y's are vectors or matrices
    and the S's are strings.

    For example, PLOT(X,Y,'y-',X,Y,'go') plots the data twice, with a
    solid yellow line interpolating green circles at the data points.

    The PLOT command, if no color is specified, makes automatic use of
    the colors specified by the axes ColorOrder property.  The default
    ColorOrder is listed in the table above for color systems where the
    default is yellow for one line, and for multiple lines, to cycle
    through the first six colors in the table.  For monochrome systems,
    PLOT cycles over the axes LineStyleOrder property.

    PLOT returns a column vector of handles to LINE objects, one
    handle per line.

    The X,Y pairs, or X,Y,S triples, can be followed by
    parameter/value pairs to specify additional properties
    of the lines.

    See also SEMILOGX, SEMILOGY, LOGLOG, GRID, CLF, CLC, TITLE, XLABEL,
    YLABEL, AXIS, AXES, HOLD, COLORDEF, LEGEND, SUBPLOT, and STEM.

The "lookfor" command is often useful for finding what might be available on a given topic. For example, one might want to know what is available in MATLAB for interpolation.

>> lookfor interpolation

BILIN uses bilinear interpolation to find the values at
interp_r.m: % INTERP_R:  A no-loop 1D interpolation routine for time series data
 triterp -- Linear interpolation on triangles.
 splinep -- Spline interpolation with end-slope control.
 splineq -- Equally-spaced spline interpolation.
 splinesafe -- Spline interpolation with end-slope control.
INTERP1 1-D interpolation (table lookup).
INTERP1Q Quick 1-D linear interpolation..
INTERP2 2-D interpolation (table lookup).
INTERP3 3-D interpolation (table lookup).
INTERPFT 1-D interpolation using FFT method.
INTERPN N-D interpolation (table lookup).
SPLINE Cubic spline data interpolation.
SPLNCORE N-D Spline interpolation.
NDGRID Generation of arrays for N-D functions and interpolation.
 QHULLDEMO Tessellation and interpolation of scattered data.
INTERP Resample data at a higher rate using lowpass interpolation.
INTFILT Interpolation (and Decimation) FIR Filter Design
PDENTRP  Interpolation helper function for PDEPE.
BLKC1INTERP This block performs constant interpolation between a given set of input and
BLKL1INTERP This block performs linear interpolation between a given set of input and
BLKLINTABLE This block defines two kinds of linear table interpolation. One is interpolation of
BLKS1INTERP This block performs cubic spline interpolation between a given set of input and
MBLKMDL Driver for the multidimensional linear interpolation block.

After finding the items available, one can use the "help" command to get more information (e.g. 'help interp1).
 

Reading data into MATLAB


How about getting data into MATLAB without typing it directly from the keyboard? Suppose you have data from a tide gauge in Sandy Hook NJ (data from http://www.co-ops.nos.noaa.gov/data_res.html) and your data is stored in a file called sandyhook_tide.dat that looks like:

alpha3: more sandyhook_tide.dat

% Verified Hourly Water Level Data (W2)
%
%       Station   --   Unique seven character identifier for the station
%       Date Time --   Date and time the data were collected by the DCP
%       WL_Value  --   Water level height
%       Sigma     --   Standard deviation of 1 second samples used to
%                      compute the water level height
%       I         --   A flag that indicates that the water level value
%                      has been inferred.
%       L         --   A flag that when set to 1 indicates that either the
%                      maximum or minimum expected water level height limit
%                      was exceeded.
%
% Data are in Meters above MSL
% Times are on UTC (GMT)

% 8531680 SANDY HOOK , NJ  from  19991201 until 19991231
% ------------------------------------------------------------------------

% Station Date       Time      WL    Sigma I L
8531680 1999 12 01 00 00  -0.394   0.013 0 0
8531680 1999 12 01 01 00  -0.555   0.009 0 0
8531680 1999 12 01 02 00  -0.567   0.009 0 0
8531680 1999 12 01 03 00  -0.326   0.013 0 0
8531680 1999 12 01 04 00   0.018   0.013 0 0
8531680 1999 12 01 05 00   0.379   0.009 0 0
8531680 1999 12 01 06 00   0.661   0.013 0 0
8531680 1999 12 01 07 00   0.834   0.017 0 0
8531680 1999 12 01 08 00   0.854   0.018 0 0
8531680 1999 12 01 09 00   0.732   0.016 0 0
8531680 1999 12 01 10 00   0.453   0.019 0 0
8531680 1999 12 01 11 00   0.208   0.013 0 0
8531680 1999 12 01 12 00  -0.031   0.020 0 0
8531680 1999 12 01 13 00  -0.229   0.014 0 0
8531680 1999 12 01 14 00  -0.276   0.015 0 0

Note that MATLAB neglects the comment lines (preceeded by "%") at the beginning of data files. You can use the "load command to get the data into MATLAB:

>>load sandyhook_tide.dat
>>whos
 Name                 Size         Bytes  Class

  sandyhook_tide    3672x10        293760  double array

Grand total is 36720 elements using 293760 bytes
>>size(sandyhook_tide)
ans =

        3672          10

Note that the files to be loaded into MATLAB MUST have an extension (here ".dat"), otherwise an extension of ".mat" will be assumed. If it is a *.mat file; MATLAB will assume that it has been stored in matlab-format; otherwise it assumes that it is nice ASCII text.  You can get a copy of shook_99.dat off the course web page (www.vims.edu/~ckharris/MS554/sandyhook_tide.dat); and try loading it into MATLAB. You should now have 3,672 records of water elevation from Sandy Hook, NJ, stored in a 3,672 X 10 matrix. The next thing to do is to assign dates and plot the water level:

>>yr=sandyhook_tide(:,2);
>>mnth=sandyhook_tide(:,3);
>>dy=sandyhook_tide(:,4);
>>hr=sandyhook_tide(:,5);
>>mnt=sandyhook_tide(:,6);
>>sec=0*mnt;    % record neglects seconds in timestamp
>>wlevel=sandyhook_tide(:,7);

>>jd=datenum(yr, mnth, dy, hr, mnt, sec); % datenum gives julian date since year 0000
>>jd2000=jd-datenum(2000,1,1);        % convert jd to julian date / 2000.

Note that in the above example, the colons are used to select entire columns of the matrix. That is "sandyhook_tide(:,2)" means "all elements in the 2nd column of sandyhook_tide". To look at what sea level was doing in Sandy Hook, NJ, you can now plot the data.

>>plot (jd2000, wlevel);

That (figure2) looks pretty nice.

>>plot(jd,wlevel)
>>datetick('x',12,'keeplimits')
>>ylabel('Water Level, m')
>>xlabel('Date')
>>title('Sea Level at Sandy Hook, NJ')

Wow - that (figure3) looks snazzy.

Use the "help" command to learn more about datenum and datetick, if you are interested.

You can load data into MATLAB as long as the data are in an array that contains only numbers, with the exception of comment lines ('% blahblahblah yada yada') at the beginning of the file. The text file must also represent a matrix- that is you must have the same number of columns of numbers on each row. Missing values can be represented with the string "nan", for "not-a-number", or by some "special value", for example we could use -999 to represent missing tide gage data, because we can be pretty sure that the data will never be near -999. You can not leave missing data blank.
 

Mathematical Operations with Matrices


MATLAB is a powerful engine for data analysis and modeling primarily because of its design for performing matrix calculations. All standard mathematical operations are supported. You just have to remember that matrices on which you are operating must be "conformable" for that particular operation; and remember that it is a matrix- not a group of scalars.
 

Elemetary operations: addition, subtraction, and so forth


Matrices (vectors as a special case) are conformable for addition (and subtraction) as long as they are the same size. Thus, if

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

A =

     1     2     3
     4     5     6

and
>> B=[7 8 9; 10 11 12]
B=

     7     8     9
    10    11    12

then the sum and difference are what you would expect:

>> A+B

ans =

     8    10    12
    14    16    18

>> A-B

ans =

    -6    -6    -6
    -6    -6    -6

Multiplication by a scalar (a 1x1 matrix if you like) is also straightforward.
>> 2*A

ans =

     2     4     6
     8    10    12

All of the functions available in MATLAB can be applied to matrices. These include sin, cos, exp, log, and others. For example:
>> sqrt(A)

ans =

    1.0000    1.4142    1.7321
    2.0000    2.2361    2.4495

Note that this functions operated on each element of the matrix. Other operations operate on the matrix as a whole.
 

Matrix Multiplication

Two matrices, say C of dimension n by m, and D of dimension k by j, are conformable for multiplication if m=k. Matrices A and B, above, are therefore not conformable for multiplication because A is 2x3 and B is 2x3. (That is, "m" is 3 and "k" is 2.) If one tries to multiply the two, MATLAB gives an error message.

>> A*B
??? Error using ==> *
Inner matrix dimensions must agree.

Now, the transpose of B is a 3x2 matrix, and the transpose is formed using a single quotation mark after the matrix name:

>> BT=B'

BT =

     7    10
     8    11
     9    12

>> BT*A

ans =

    47    64    81
    52    71    90
    57    78    99

By the rules of matrix multiplication, A and BT are conformable for multiplication.

The case of multiplication of a matrix and a vector is of particular importance as systems of equations are represented in this form. The form of a set of equations is Ax=B, where A is a matrix of (known) coefficients, x is the unknown vector, and B is a vector of known quantities. The system of equations:

3x + 3y = 9
2x - 2y = -2

is represented in matrix notation as Ax = B; where A = [ 3 3; 2 -2]; x=[x; y]; and B=[9; -2];

Our main interest in such matrix equations is solving them-- determining values for x and y that satisfy the equations. We will explore this topic in some detail later in the course.  For now; just accept that the "backslash" operator does the trick.

>> A = [ 3 3; 2 -2]; B=[9; -2];
>> xy=A\B

xy =

     1
     2

Let's check that:

>>A*xy - B
ans =
   0
    0
 

Multiplication on an element-by-element basis

Suppose you need to operate on a vector element by element. MATLAB uses a period at the end of the vector to indicate operation on an element-by-element basis. For example, suppose we want to take the vector of sea-level "wlevel" that we defined above and square each element. Note that the command "wlevel2=wlevel*wlevel" won't work because the "*" implies matrix multiplication and the only time a matrix is conformable for itself for multiplication is when it is square. And remember that, even for a square matrix, A*A is not equivalent to taking the square of each element of A. To have MATLAB square the elements of wlevel, we can use either "wlevel.*wlevel" or "wlevel.^2".

>> wlevel(1:6).*wlevel(1:6)

ans =

    0.1552
    0.3080
    0.3215
    0.1063
    0.0003
    0.1436

>> wlevel(1:6).^2

ans =

    0.1552
    0.3080
    0.3215
    0.1063
    0.0003
    0.1436

Note that we looked only at the first 6 elements of the vector by indicating the range 1:6. this is standard notation in MATLAB, A(n1:n2, m1:m2) indicates the submatrix of A with rows n1-n2, and columns m1-m2 of the original matrix A.

Another example of the use of element-by-element computation is in finding the reciprocal of the elements of the vector. Suppose we need the reciprocal (1/) of the elements in the matrix B. As it turns out, we can not use 1/B because division for matrices is not defined. Also, "inv(B)", read "inverse of B" won't work because the inverse of a matrix is not the same as inverting the matrix. (The inverse of a square matrix is another square matrix, that when multiplied by itself gives the original matrix.) So again, we use the "dot" at the end of a vector to denote the element-by-element operation.
 

>> B

B =

     9
    -2

>> 1./B

ans =

    0.1111
   -0.5000

As a final example, consider wanting to estimate (dH/dt); the change in water level per time. You can use the MATLAB "diff" function to calculate differences between consecutive elements of each vector:

>> dH=diff(wlevel);
>> dt=diff(jd2000)*24; % dt in hours
>> dHdt=dH./dt; % in m/hour
>>
>> plot(jd2000(2:end), dHdt)
>>
 

Programming in MATLAB

To use the full power of MATLAB, you will need to learn to write "m-files". These files are programs to execute a sequence of calculations and operations. For example, you might want to have a program to take a set of data and perform some statistical calculation. The files can contain any legitimate MATLAB command. The macro file can have any legal name with the extension ".m". MATLAB looks for files with a ".m" extension when the prefix is typed, interprets the commands, and executes them in sequence.

Script Files

The first type of m-file that you will use is a straight script file. Script files are, literally, a list of commands that are interpreted and executed. To illustrate, suppose you want to create a program to load the water-level data from Sandy Hook (considered earlier), and calculate some statistics on the data. Using your favorite word processor type in the following into a file named process_shook.m:

% loads sandyhook_tide.dat and plots the time-series.
% created by anerd@vims.edu, August, 2001.
load sandyhook_tide.dat
yr=sandyhook_tide(:,2);
mnth=sandyhook_tide(:,3);
dy=sandyhook_tide(:,4);
hr=sandyhook_tide(:,5);
mnt=sandyhook_tide(:,6);
sec=zeros(size(mnt));
wlevel=sandyhook_tide(:,7);

jd=datenum(yr, mnth, dy, hr, mnt, sec);
jd2000=jd-datenum(2000,1,1);

max_H=max(wlevel);
min_H=min(wlevel);
std_H=std(wlevel);   % standard deviation

clf
plot(jd, wlevel,'r')
hold on
plot(jd, std_H*ones(size(jd)),'b--')
plot(jd, -1*std_H*ones(size(jd)),'b--')
plot(jd, 0*jd, 'k-');  % put in a line for zero.

datetick('x',12,'keeplimits')

legend('water level','\pm 1 stdev')

ylabel('Water Level, m')
xlabel('Date')
title('Sandy Hook, NJ')

Note that, again, comments are preceeded by "%". Then, from within MATLAB, type "help process_shook" and
then process_shook.m
>> help process_shook

  loads sandyhook_tide.dat and plots the time-series.
  created by anerd@vims.edu, August, 2001.
>>process_shook

The results are in figure4.png.

The process_shook.m file is a pretty simple m-file, but you should be able to appreciate how easy it is to program more complicated tasks.
 

Function Files

The other type of m-file that you will use is a function file. A function file operates on vectors and matrices specified in the function call and returns variables defined in the function statement. The MATLAB help on this topic illustrates the concept nicely.

>> help function

 FUNCTION Add new function.
    New functions may be added to MATLAB's vocabulary if they
    are expressed in terms of other existing functions. The
    commands and functions that comprise the new function must
    be put in a file whose name defines the name of the new
    function, with a filename extension of '.m'. At the top of
    the file must be a line that contains the syntax definition
    for the new function. For example, the existence of a file
    on disk called STAT.M with:

            function [mean,stdev] = stat(x)
            %STAT Interesting statistics.
            n = length(x);
            mean = sum(x) / n;
            stdev = sqrt(sum((x - mean).^2)/n);

    defines a new function called STAT that calculates the
    mean and standard deviation of a vector. The variables
    within the body of the function are all local variables.
    See SCRIPT for procedures that work globally on the work-
    space.

    A subfunction that is visible to the other functions in the
    same file is created by defining a new function with the FUNCTION
    keyword after the body of the preceding function or subfunction.
    For example, avg is a subfunction within the file STAT.M:

           function [mean,stdev] = stat(x)
           %STAT Interesting statistics.
           n = length(x);
           mean = avg(x,n);
           stdev = sqrt(sum((x-avg(x,n)).^2)/n);

           %-------------------------
           function mean = avg(x,n)
           %MEAN subfunction
          mean = sum(x)/n;

    Subfunctions are not visible outside the file where they are defined.
    Normally functions return when the end of the function is reached.
    A RETURN statement can be used to force an early return.

    See also SCRIPT, RETURN, VARARGIN, VARARGOUT, NARGIN, NARGOUT,
             INPUTNAME, MFILENAME.

As an example of a computation useful in oceanography, below is a function that low-pass filters a time-series, that is, it removes the high-frequency components and returns the residual low-frequency components.  The function name is lpfilt.m; and it was written by Dr. Chris Sherwood, who now works for the USGS Coastal and Marine Geology Program.

alpha3: cat lpfilt.m

function [fdata] = lpfilt(data,delta_t,cutoff_f)
%  fdata = lpfilt(data,delta_t,cutoff_f)
%
%  Peforms low-pass filter by multiplication in frequency domain.
%  Uses three-point taper (in frequency space) between pass-band and
%  stop band.

%  Coefficients suggested by D. Coats, Battelle, Ventura.
%  Even better coefficients for the taper can be found in:
%  Rabiner, L.R., Gold,B., and McGonegal, C.A. (1980).  An Approach to
%  the approximation problem for nonrecursive digital filters.  IEEE Tran.
%  vol. AU-18(2):83-106.

%  Written by Chris Sherwood, Battelle PNL
%  Last revised: 9/03/89

n = length(data);
mn = mean(data);
data = data-mn;
P = fft(data);
N = length(P);
filt = ones(N,1);
k = floor(cutoff_f*N*delta_t);
filt(1:k) = filt(1:k)*1;
filt(k+1) = .715;
filt(k+2) = .24;
filt(k+3) = .024;
filt(k+4:N-(k+4)) = filt(k+4:N-(k+4))*0.;
filt(N-(k+3)) = .024;
filt(N-(k+2)) = .24;
filt(N-(k+1)) = .715;
P = P .* filt;
fdata = real(ifft(P));
fdata = fdata(1:n)+mn;

The function, lpfilt, takes as input the variables "data", "delta_t", and "cutoff_f", and returns the variable fdata. You will learn more about filtering time-series later in the class. For now lets see what lpfilt can do with our tidal elevation data. Suppose we want to see how frequently sea level sets up, or sets down at Sandy Hook, NJ.

>> process_shook;   % read in our tidal data.
>> dt=mean(diff(jd))*24;  % delta_t in hours; we made sure data was evenly spaced in time.
>> cutoff=1/34;  % remove components with periods < 34 hours- that should get rid of tides.
>> lpH=lpfilt(wlevel, dt, cutoff);   % use the lpfilt function to low-pass filter the data.
>> plot(jd, lpH, 'b');  % add the low-frequency values to the hourly measured values.

The resulting figure is in figure5.png on the website.

Summary

Throughout this course we will be using MATLAB for computation and for programming. You should familiarize yourself with this software. Read these notes, read the documentation (in the Student Lab), and use the on-line help facility in the program, and start practicing by writing some M-files.
 

FORTRAN

Just as important- you will need to write, compile, and debug FORTRAN files.  Make sure that you have the ability to do so by downloading the short FORTRAN file that is on the class web-site. Compile it and run it on the computer that you will be using for the rest of class.  Each compiler uses a different syntax. I tend to use fortran from within a Unix environment; and would compile the code using the command sequence

alpha3:  f77 -o ws ws.f

From within unix, the man command (alpha3:  man f77) would give information about the syntax of f77. It would tell you, for example, that the object code created, the executable code in this case, would be named ws.