% Math22 Matlab session example commands. Alex Barnett, 10/4/16

% Anything on a line after this percent symbol is a comment, ignored by matlab
% Use comments to aid human comprehension!

% Use Matlab's edit command to open this and edit, and/or paste lines into
% the main matlab command window.

A = [1 2 3; 4 5 6]          % fill a matrix, call it A
x = [1;0;-1]                % 3x1 column vector
A*x                         % matrix-vector product
A*A               % cannot do!
B = transpose(A)  % fill B with A^T
B = A'            % abbreviation for transpose
B*A               % matrix-matrix mult
A*B
5*A               % scalar mult
C = [-1 1 -1;1 -1 1];    % same shape as A
A.*C          % elementwise mult (note .* not *).  Ie, ans_ij = a_ij*c_ij

whos         % what's in your workspace?
save temp.mat   % save your workspace to a file
clear          % clear the workspace
load temp.mat    % load in your saves workspace

A = [1 2 3;4 5 6;7 8 9];
rref(A)           % 2 pivots   (as expect)
b = A*[1;1;-1]         % make a RHS we know is a lin comb of cols of A
x = A\b;       % "\" is called "backslash": solves a (possibly
               % rectangular) linear system. Aka mldivide, linsolve
A*x - b       % check if x *is* a solution: should give (close to) zero vector
% Note: not exactly zero since all arithmetic is done to around 16 digits
% accuracy. This is called "round-off error".
format long g             % displays all decimal digits

det(A)      % close to zero since in exact arithmetic it's singular.
inv(A)      % scary large values - note it warns you too!

% Round-off error can bite you: doing x = inv(A)*b is disaster if singular...
y = inv(A)*b;   % mathematically, y should be a solution, but...
A*y-b           % oops! is not the zero vector, not even close. Fail.
% This is because of round-off error in doing the mat-vec inv(A)*b

doc inv       % find all about any command
help inv      % similar but appears in the command window

rref([1 1; 0 1])               % gives I_2 as expect
rref([1e-10 1; 0 1])           % ditto
rref([5e-16 1; 0 1])           % ditto
rref([2e-16 1; 0 1])           % due to round-off error, have missing pivot!

% show how fast things are for large matrices (n will set the square n*n size):
n = 100; A = rand(n,n); b = rand(n,1); tic; x = A\b; toc
n = 1000; A = rand(n,n); b = rand(n,1); tic; x = A\b; toc  % fast
n = 10000; A = rand(n,n); b = rand(n,1); tic; x = A\b; toc  % slow (10 sec)
% in general, time is proportional to number of arithmetic operations (flops),
% which grows like n^3 for

A = rand(40,60);               % rectangular matrix of random numbers in [0,1]
imagesc(A); colorbar;          % show color-scale (density) plot of matrix
U = rref(A);                   % compute its REF form
spy(U);                        % show diagram of nonzero entries

% LISTS: the colon notation gives a row vector
y = 1:10
y = 1:2:10          % in steps of 2
y = (1:10)'         % this gives a column vector

% LOOPS

% simple loop doing something many times. Eg, evolving a difference eqn:
A = [.9 .1; .1 .9];     % migration matrix
x = [1;0];             % initial populations
for i=1:20
  x = A*x                % will print out vector each year
end
% problem: we lost all the intermediate values. The next version stores them.
% Better version, store the complete evolution of x...
x = zeros(2,20); % 20 columns since we want x_0, x_1, ..., x_20
x(:,1) = [1;0];  % set the initial values
for i=1:19
  x(:,i+1) = A*x(:,i);
end
% now x has complete history, useful, can plot it...
figure; plot(1:20, x(1,:), '+-');

% make a professional labelled plot with both populations as a func of time...
figure; plot(1:20, x, '+-');
xlabel('time (years)'); ylabel('populations (millions)'); title('evolution');

