% Introduction to Matlab --- Raymond Spiteri

% Matlab has a number of advantages for doing numerical computing over
% traditional methods (e.g., writing programs in FORTRAN or C and
% calling numerical libraries).

% These advantages include:

% a very high-level language for quick and easy coding

% minimal attention required for data structures (e.g., they do not
% need to be declared before they are used)

% an interactive interface for convenient experimentation and debugging

% high-quality graphics and visualization
% (try doing that easily in FORTRAN or C!)

% portability of code across platforms
% (i.e., code gives exactly the same results on different computers)

% toolboxes are available to provide specialized facilities
% (we will be using the symbolic toolbox)

% a wide range of contributed M-files (Matlab programs) are freely
% available on the web (google matlabcentral)

% NOTE: Matlab is interpreted, not compiled, so it is slower than
% FORTRAN or C, but that is not of concern to us in this course

%%%%%%%%%%%%%%%%% begin tutorial %%%%%%%%%%%%%%%%%%%%%%%

% notice that everything after a '%' is a comment
% no matter where it appears on the line

%%%%%%%%%%%%%%%% getting started %%%%%%%%%%%%%%%%%%%%%%%

% learning Matlab is no more a spectator sport than is learning how to
% swim; the best way to learn Matlab is by trying things for yourself

help      % this is a very important command to know!
doc       % so is this one

help demo % in general "help topic" when you know what you want
lookfor   % when you're not sure what you want!

demo      % for touring the capabilities of Matlab

% you can probably guess what the next two commands do
% (if you can't, you'll probably figure it out pretty quick!)
quit 
exit 

% some other useful commmands

pwd   % print working directory
cd    % on its own, this is equivalent to pwd
      % but it can be used to change directory
    
dir   % both of these show you the contents of the directory
ls    % but they are not exactly the same!

path  % shows you what's in Matlab's path
      % i.e., what directories it will look in for files
     
who   % the names of the variables in your workspace
whos  % the names and sizes

clear % clears the current workspace
      % it is often good practice to start scripts with this command!      
      
diary 'tutorial.txt'  % start recording your session in file
                      % tutorial.txt
% in the meantime, let's look at some pre-defined constants
pi
eps
Inf
NaN
diary off          % stop recording session

% scientific notation is OK

a = 1e-10

% note: there is a limit on the number of characters in an identifier
% (variable) name that matlab recognizes

namelengthmax

% define two vectors
b = [1 2 3]
% note: commas can be put between elements in the same row
% so the following command is equivalent to the one we just typed:
% b = [1, 2, 3]
c = [4; 5; 6]
% a semi-colon starts a new row
% there is a transpose operator (') as well, 
% so c could have been defined by
% c = [4 5 6]'

% as an aside, ": notation" is a very powerful construct in Matlab
% in general you don't have to type out all the numbers in a sequence where
% the numbers are generated in an incremental fashion
% notation: 
%                   startingValue : increment : finalValue

% so we could have written c = [4:1:6]' 
% note: not much typing is saved in this example, but you can imagine
% the savings in typing if you needed to set c = [1:2:701]

% increment = 1 is assumed unless otherwise specified
% so c = [4:6]' is equivalent to what we just did 

% this notation can be used in general
% i.e., it does not have to be associated with an array
v = 2:7:30

% non-integer increments are also OK
v = 0:0.1:1

% b is 1x3, c is 3x1, so b*c makes sense and produces a 1x1 (scalar)
% value
b*c
% this is of course really a dot product, so we could have typed
dot(b,c)

% the product c*b also makes sense; what is the result?
D = c*b

% the product c*c does not make sense; why not?
c*c

% neither does the sum b+c
b+c

% in Matlab, everything (scalars, vectors, matrices) is treated like a
% matrix

A = [1 2 3; 4 5 6; 7 8 9] % give me a matrix

% don't worry -- Matlab is case sensitive, so the old (scalar) a
% still exists!
% note that [] and () are not equivalent / interchangeable

size(A)                   % dimensions of A

A'                        % transpose

A + A                     % matrix addition

2*A                       % scalar-matrix multiplication

A - A                     % matrix subtraction

A*A                       % matrix-matrix multiplication

A^2                       % powers of matrices

A/A                       % matrix division?? this really means A*inv(A)

A^(-1)                    % same as inv(A)

det(A)                    % determinant

trace(A)                  % sum of diagonal elements

eig(A)                    % eigenvalues

[V,lambda] = eig(A)       % eigenvectors and eigenvalues

% note that you can extract different information from a function
% depending on what appears on the left-hand side of the assignment
% statement 

% access a chunk of A
% we know a little about ": notation"
A(2:3,1:2)
% when no starting or ending values are specified, 
% the values are treated as "wildcards"

A(:,1)                    % first column of A
A(2,:)                    % second row of A

% there is also a handy little construct called "end"
A(2:end,3:end)            % everything in the submatrix below A(2,3)

% you can construct matrices by gluing together other matrices
B = [b; c'; ones(1,3); [0 pi 7]]

% in general you have to be careful that you are aware as to whether an
% operation is being applied component-wise to a matrix or as a matrix
% function

% component-wise operations

A.*A                      % not the same as A^2 !!!

A.^2                      % same as A.*A

A./A                      % matrix of ones (pretty boring!)

exp(A)                    % take exponential of every element

% we have seen that numbers are displayed to 4 decimals by default
% this can be changed (see help format for all the options)
% for now, we would like long floating-point format
format long e

log(A)             % take natural logarithm (ln) of every element

% matrix functions

expm(A)             % matrix exponential! I + A + A^2/2! + A^3/3! + ...

logm(A)             % matrix logarithm (complex because of 0 eigenvalue)

% common matrices
% as an aside, you can type more than one command on a line
% separated by commas or semi-colons

m = 4; n = 2;             % the semi-colon suppresses the output
                          % very useful when writing scripts 
                          % (improves performance!)
                          
ones(m,n)                 % matrix of ones  

zeros(m)                  % matrix of zeros

eye(n)                    % identity matrix (math joke)

A = diag([1,2,3,4])           % a diagonal matrix

diag([1,2,3,4],-1)        % a matrix with a subdiagonal

% note that a function can take a different number of arguments and
% matlab figures out (usually successfully) what you mean

% matlab has "dynamic memory allocation", which means it will
% figure out how to manage memory for you

% suppose we forgot to set A(5,5) = 5 two commands ago ...
A(5,5) = 5; % no problem

% but beware! matlab fills extra elements with zeros
% this may or may not be what you want

% the previous command was OK, but ...
A      = 5*ones(4)
A(5,5) = 5

% also, dynamic array allocation is expensive
% it is better to pre-allocate arrays

m = 2000;
clear A
tic 
  for i=1:m, 
    for j=1:m, 
      A(i,j) = i+j; 
    end 
  end 
toc

clear A
A = zeros(m);
tic 
  for i=1:m, 
    for j=1:m, 
      A(i,j) = i+j; 
    end 
  end 
toc

% linear system solving: given matrix A and vector b such that Ax=b,
% find x

A = [3 1; 2 2]; 
b = [0; 1];     
x = A\b

% random numbers

x = rand            % a random number between 0 and 1 (uniformly chosen)
A = rand(m)         % uniformly distributed mxm matrix of random numbers
B = randn(m,n)      % normally distributed mxn matrix of random numbers 
                    % (mean 0, standard deviation 1)

seed = 0;           % used to 'control' randomness 
rand('state',seed)  % (make random runs reproducible)

% let's generate (and save) some random data for later
random_data = randn(1000,1);
save random_data.dat random_data -ascii
clear
load random_data.dat
random_data
more on             % let's stop it so we can look at it!

% simple plotting

ezplot('sin(x)/x')                        % default axes
ezplot('sin(x)/x',[-3*pi, 3*pi,-0.5,1.5]) % customized axes

% data plotting

plot(random_data)
hold on             % do not make new figure (set "hold" to "on")
plot(random_data,'ro')
hold off
semilogy(random_data)

% basic statistical functions

max(random_data)
min(random_data)
mean(random_data)                    % average
std(random_data)                     % standard deviation
median(random_data)
sum(random_data)
prod(random_data)
sort(random_data)                    % sort data in increasing order

hist(random_data,40)                 % create a histogram of data frequencies

a = [0.5 1 1.6 1.2 .8 2.1];          % create some other phony data
pie(a)                               % make a pie chart
pie(a,a==max(a))                     % pull out the largest slice
                                     % (dieting can wait)

a = sort(a);                         % ascending order by default
t = 1:length(a);
plot(a,'ko')
hold on
coeff = polyfit(t,a,1)               % linear regression 
                                     % (line of best fit)
fit = coeff(1)*t + coeff(2)
plot(t,fit)

coeff = polyfit(t,a,2)               % parabola of best fit
fit   = coeff(1)*t.^2 + coeff(2)*t + coeff(3) 
plot(t,fit,'g--')

% simple looping constructs

% for loop
x = 0;
for i=1:10,
  x = x + 1;
end
disp(x) % there is a more sophisticated way to print results using
        % sprintf

% while loop
x = 0;
while (x < 10),
  x = x + 1;
end
disp(x)

% now add a logical branching construct (if-statement)
x = 0;
for i=1:10,
  if (i > 5)
      x = x + 1;
  end
end
x

x = 0;
for i=1:10,
  if (i > 5)
    x = x + 1;
  else
    x = x - 1;
  end
end
x

% note: these looping and branching constructs are nice, but if you
% are using them to manipulate matrices, it is almost surely not a
% good use of them (use "data parallel" constructs instead)

m = 2e7;
v = ones(m,1);
tic
  for i=1:m,
    v(i) = exp(v(i));
  end
toc

tic
  v = exp(v);
toc

% symbolics

S1 = sym('sin(x)') % declare a symbolic function
% symbolic variables can be defined using the "syms" construct
% e.g., syms x y z

% now differentiate S1 once, the four times, then integrate it
diff(S1)
diff(S1,4)
int(S1)

% numerical quadrature: quad and quadl

a = input('Please enter left  end of integration: ');
b = input('Please enter right end of integration: ');

quad('sin(x)',a,b)

% extra long lines can be continued using ...

% here is an extra long message displaying the previous result
% note the use of the num2str function to convert numbers to
% strings so they can be printed

disp(['The value of the integral of sin(x) from x = ', ...
      num2str(a), ' to x = ', num2str(b), ' is ', ...
      num2str(ans), '.'])

% to perform a sequence of related commands, you can write a script
% scripts are easy to write, use, and debug
% they are much more convenient than having to type things in over and
% over again!
% or (gasp!) cutting and pasting into the window

% in fact this file is itself a long (and convoluted) script
% we will give a more effective illustration of a script in a
% separate file

% we have also seen the use of functions
% we will give an example of how to write your own functions
% separately as well 