% 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');