%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From: MATLAB An Introduction with Applications
%       Amos Gilat
%
Chapter 3   Mathematical Operations with Arrays
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preamble:

% One of the main design goals of Matlab was to make 
% calculations in linear algebra straightforward.
% Among other things, this means that operations central 
% to that branch of mathematics have a simple and natural 
% expression in the language. 

% A key example of this is that, as we shall see, 
% the operations of 

%   matrix-vector multiplication
%   matrix-matrix multiplication
%   matrix-vector division
%   matrix-matrix division

% are implemented using the '*' and '/' operators 
% that we are familar with in the context of arithmetic
% expressions using scalar quantities.  We refer to these 
% operations as "matrix" or "whole-array" or "array" 
% multiplication and division.

% However, there is a distinct type of multiplication/division
% operation we will frequently wish to use with arrays, which 
% is to apply the operators on an "element-by-element" basis.  

% Matlab thus has a separate notation/syntax for element-wise 
% multiplication, division (and exponentiation) operations,
% which is to prepend the usual operator with a '.' (dot,
% period).  The syntax can also be used with addition and 
% subtraction, but in those instances the '.' is, as we shall 
% see redundant.

% Thus, the operators for element-wise arithmetic 
% operations on arrays are 
%
%    .*     Element-by-element multiplication
%    ./     Element-by-element (right) division
%    .\     Element-by-element left division
%    .^     Element-by-element exponentiation.
%

% For both matrix and element-wise caclulations there are 
% constraints that the operands must satisfy so that the 
% operation is well defined.

% 1) For matrix operations, the dimensions (which we will
%    also refer to as sizes) of the operands must be 
%    be conformant according to the usual rules  of linear 
%    algebra.

% 2) For element-wise operations, one of the following must 
%    be true
%    
%    i)   Both of the operands are scalars.
%
%    ii)  One of the operands is a scalar and the other
%         is an array.
%
%    iii) Both operands are arrays and have the same
%         dimensions.

%    The result of an element-by-element array operation
%    has the same dimensions of the operand array(s), with
%    values that, as the nomenclature suggests, are 
%    calculated from the action of the particular 
%    binary operator (*, /, \, ^) on corresponding 
%    elements of the operand arrays.
%
% If all of this is confusing, the examples below should
% clear things up.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Arithmetic operations between a scalar and an array
%
% If one of the operands of an elementary operation
% involving +, -, *, or / is a scalar and the other 
% is an array, A, then the operation is applied 
% (distributed) to each element of A, yielding a result 
% that is also an array with the same dimensions as 
% A.  

% In short, "what you get is what you expect"!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  v = 1:5

v =

   1   2   3   4   5


>>  v1 = v + pi

v1 =

    4.1416    5.1416    6.1416    7.1416    8.1416


>>  v2 = 3.1 * v

v2 =

     3.1000     6.2000     9.3000    12.4000    15.5000


>>  v3 = v / 5

v3 =

   0.20000   0.40000   0.60000   0.80000   1.00000


>>  v4 = 2 * v + 6

v4 =

    8   10   12   14   16


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

m =

   1   2   3
   4   5   6


>>  m1 = 3 * m

m1 =

    3    6    9
   12   15   18


>>  offset = -12

offset = -12

>>  m2 = m / 2 + offset

m2 =

  -11.5000  -11.0000  -10.5000
  -10.0000   -9.5000   -9.0000


>>  -3 * ones(1,7)

ans =

  -3  -3  -3  -3  -3  -3  -3


>>  6 * ones(3)

ans =

   6   6   6
   6   6   6
   6   6   6


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NOTE: The exponentiation operator ^ can NOT be used
% in this way, e.g. v^2 will NOT return the vector 
% [1 4 9 15 25].  The element-wise syntax discussed 
% below must be used instead.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.1 Addition and Subtraction
%
NOTE: For addition and subtraction, matrix/array
% operation is always an element-by-element operation.
%
% If arrays are added or subtracted then:

%   1) They must have identical sizes.
%
%   2) The arrays are added or subtracted element-wise 
%      to produce the result.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  VectA = [8 5 4]

VectA =

   8   5   4


>>  VectB = [10 2 7]

VectB =

   10    2    7


>>  VectC = VectA + VectB

VectC =

   18    7   11


>>  VectA + 13

ans =

   21   18   17


>>  A = [5 -3 8; 9 2 10]

A =

    5   -3    8
    9    2   10


>>  B = [10 7 4; -11 15 1]

B =

   10    7    4
  -11   15    1


>>  A - B

ans =

   -5  -10    4
   20  -13    9


>>  C = A + B

C =

   15    4   12
   -2   17   11


>>  C - 8

ans =

    7   -4    4
  -10    9    3


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.2 Array (Matrix) Multiplication>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix multiplication in MATLAB is denoted by *
% and is defined according to the standard rule of 
% linear algebra: namely 
%
% If 
%     A is an m x n matrix with elements a(i,j)
%
% and 
%     B is an n x p matrix with elements b(j,k)
%
% then
%     C = A * B is an m x p matrix with elements c(i,k)
%
% defined by  
%               n
%     c(i,k) = Sum [ a(i,j) x b(j,k) ]
%              j=1
%
NOTE: For this operation to be well defined,
% the number of columns of A must equal the number of 
% rows of B.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

A =

   1   4   2
   5   7   3
   9   1   6
   4   2   8


>>  B = [6 1; 2 5; 7 3]

B =

   6   1
   2   5
   7   3


>>  C = A * B

C =

   28   27
   65   49
   98   32
   84   38


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D = B * A  

% would generate an error message since the number of 
% columns of B is not equal to the number of rows of A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  F = [1 3; 5 7]

F =

   1   3
   5   7


>>  G = [4 2; 1 6]

G =

   4   2
   1   6


>>  F * G

ans =

    7   20
   27   52


>>  G * F

ans =

   14   26
   31   45


>>  AV = [2 5 1]

AV =

   2   5   1


>>  BV = [3; 1; 4]

BV =

   3
   1
   4


>>  AV * BV

ans =  15

>>  BV * AV

ans =

    6   15    3
    2    5    1
    8   20    4


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.4 Element-By-Element Operations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% As already noted above, the following binary operators
%
%   *      multiplication
%   \      right division
%   /      left division
%   ^      exponentiation
%
% can be applied to arrays in an element-by-element
% fashion by prefixing the operator with a '.'
%
%   .*     element-by-element multiplication
%   .\     element-by-element right division
%   ./     element-by-element left division
%   .^     element-by-element exponentiation

NOTE: When using element-by-element operations
% between arrays, the arrays must be the same size.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

A =

   2   6   3
   5   8   4


>>  B = [1 4 10; 3 2 7]

B =

    1    4   10
    3    2    7


>>  A .* B

ans =

    2   24   30
   15   16   28


>>  C = A ./ B

C =

   2.00000   1.50000   0.30000
   1.66667   4.00000   0.57143


>>  D = A .\ B

D =

   0.50000   0.66667   3.33333
   0.60000   0.25000   1.75000


>>  C .* D

ans =

   1   1   1
   1   1   1


>>  B .^ 3

ans =

      1     64   1000
     27      8    343


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using element-by-element operations to compute a
% function at many values of its argument.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  x = [1:8]

x =

   1   2   3   4   5   6   7   8


>>  y = x.^2 - 4*x

y =

   -3   -4   -3    0    5   12   21   32


>>  z = [1:2:11]

z =

    1    3    5    7    9   11


>>  y = (z.^3 + 5*z) ./ (4*z.^2 - 10)

y =

  -1.0000   1.6154   1.6667   2.0323   2.4650   2.9241


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We can also use the element-wise notation with the
% addition and subtraction operators, but it is 
unnecessary.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  [ 1 2 3 4] .+ [5 6 7 8]

ans =

    6    8   10   12


>>  [ 1 2 3 4] + [5 6 7 8]

ans =

    6    8   10   12


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NOTE: It is very important to understand the
% distinction between "whole array" operations and 
% element by element operations for *, / and ^ 
%
% It is a very common error for a Matlab user to 
% type something like this
%
>> v1 = [1 2 3]
>> v2 = [4 5 6]
>> v3 = v1 * v2
%
% expecting v3 to be assigned [4 10 18].  
%
% Instead, the following error message results:
%
Error using *
Inner matrix dimensions must agree.
%
% since * in this context denotes matrix
% multiplication.  
%
% Similarly, attempting to square each element of v1
% using 
%
>> v4 = v1^2

% generates the error message:
%
Error using ^
Inputs must be a scalar and a square matrix.
To compute elementwise POWER, use POWER (.^) instead.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.5 Using Arrays In MATLAB Built-in Math Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% When any of MATLAB's built in math functions is 
% invoked with an array argument, the output is an
% array with identical size of the input, and with 
% elements whose values are the results of applying
% the function to the corresponding elements of 
% the input array.

% Again this is a case where "what you get is what you 
% expect"!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  x = [0:pi/6:pi]

x =

   0.00000   0.52360   1.04720   1.57080   2.09440   2.61799   3.14159


>>  y = cos(x)

y =

 Columns 1 through 6:

   1.0000e+00   8.6603e-01   5.0000e-01   6.1232e-17  -5.0000e-01  -8.6603e-01

 Column 7:

  -1.0000e+00


>>  d = [1 4 9; 16 25 36; 49 64 81]

d =

    1    4    9
   16   25   36
   49   64   81


>>  h = sqrt(d)

h =

   1   2   3
   4   5   6
   7   8   9


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.6 Built In Functions for Analyzing Arrays
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Reduction Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Many of the functions described below (those marked
% with an (*)) operate on arrays of general dimensions,
% as well as on vectors.  Furthermore, when operating 
% on a vector, many of these functions return a 
single value.
%
% For higher dimensional arrays, and for matrices in
% particular, such functions generally operate along
% columns
 (i.e. along the first dimension of the 
% array), returning a row vector of values: each value
% in the row vector is the result of applying the 
% function to the elements of the corresponding 
% column of the original matrix.

% One net effect of these functions is thus to reduce 
% the dimensionality of the array by one (or, more 
% properly, to collapse the length of the first 
% dimension to 1): such functions are therefore often 
% termed reductions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% max(<A>) (*)

%   If <A> is a vector, returns the largest element in 
%   <A>.
%
%   If <A> is a matrix, returns a row vector containing 
%   the largest element of each column of <A>.
%
% [<d><n>] = max(<A>) (*)

%   If <A> is a vector, <d> is the largest element
%   in <A> and <n> is the (first) position of the 
%   largest element.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  A = [5 9 2 4 11 6 11 1]

A =

    5    9    2    4   11    6   11    1


>>  C = max(A)

C =  11

>>  [d, n] = max(A)

d =  11
n =  5

>>  AA = [1 3 2; 4 1 3; 8 11 12; -2, 6, 100]

AA =

     1     3     2
     4     1     3
     8    11    12
    -2     6   100


>>  CC = max(AA)

CC =

     8    11   100


>>  [dd, nn] = max(AA)

dd =

     8    11   100

nn =

   3   3   4


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% min(<A>) (*)

%   Returns the smallest element in <A>.
%
% [<d><n>] = min(<A>) (*)

%   <d> is the smallest element in <A> and <n> 
%   is the (first) position of the smallest element.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  A = [5 9 2 4 11 6 11 1]

A =

    5    9    2    4   11    6   11    1


>>  C = min(A)

C =  1

>>  [d, n] = min(A)

d =  1
n =  8

>>  AA = [1 3 2; 4 1 3; 8 11 12; -2, 6, 100]

AA =

     1     3     2
     4     1     3
     8    11    12
    -2     6   100


>>  CC = min(AA)

CC =

  -2   1   2


>>  [dd, nn] = min(AA)

dd =

  -2   1   2

nn =

   4   2   1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sum(<A>) (*)

%   Returns the sum of elements of <A>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  A = [5 9 2 4]

A =

   5   9   2   4


>>  sum(A)

ans =  20

>>  AA = [1 3 2; 4 1 3; 8 11 12; -2, 6, 100]

AA =

     1     3     2
     4     1     3
     8    11    12
    -2     6   100


>>  sum(AA)

ans =

    11    21   117


>>  sum(sum(AA))

ans =  149

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Statistical Functions (assuming <A> is a vector)
%
% (These are all reductions.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% mean(<A>) (*)
%
%   Returns the mean value of elements of <A>.
%
% median(<A>) (*)
%
%   Returns the median value of elements of <A>.
%
% std(<A>) (*)
%
%   Returns the standard deviation of elements of <A>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  vals = [1 2 3 3 4 4 4 5 6 6 7 7 7 7 7 8 9]

vals =

   1   2   3   3   4   4   4   5   6   6   7   7   7   7   7   8   9


>>  mean(vals)

ans =  5.2941

>>  median(vals)

ans =  6

>>  std(vals)

ans =  2.2573

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miscellaneous Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sort(<A>) (*)
%
% sort(<A>, 'descend') (*)

%   Returns the elements of <A> sorted in ascending 
%   order.
%
NOTE: sort is not a reduction.  It 
% returns an array of the same size as its input.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  A = [5 9 2 4]

A =

   5   9   2   4


>>  sort(A)

ans =

   2   4   5   9


>>  sort(A, 'descend')

ans =

   9   5   4   2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For a matrix, columns are sorted
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  AA = [1 3 2; 4 1 3; 8 11 12; -2, 6, 100]

AA =

     1     3     2
     4     1     3
     8    11    12
    -2     6   100


>>  sort(AA)

ans =

    -2     1     2
     1     3     3
     4     6    12
     8    11   100


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sort rows using transpose operator (twice).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  sort(AA')'

ans =

     1     2     3
     1     3     4
     8    11    12
    -2     6   100


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dot(<a>,<b>

%   Returns the scalar (dot) product of two vectors <a>
%   and <b>, which can be either row or column 
%   vectors, but which must be the same length.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  v1 = 1:10

v1 =

    1    2    3    4    5    6    7    8    9   10


>>  v2 = 11:20

v2 =

   11   12   13   14   15   16   17   18   19   20


>>  dot(v1, v2)

ans =  935

>>  dot(v1, v2')

ans =  935

>>  dot(v1', v2)

ans =  935

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cross(<a>,<b>)
%
%   Calculates the cross product of two vectors <a> and 
%   <b>, each of which must be of length 3.  Returns
%   a vector which also has 3 elements.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  u1 = [1 0 0]

u1 =

   1   0   0


>>  u2 = [0 1 0]

u2 =

   0   1   0


>>  cross(u1, u2)

ans =

   0   0   1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.7 Generation of (Pseudo)-Random Numbers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
NOTE: As with many languages and software 
% systems, the "random" numbers generated by MATLAB
% are actually computed using specific, deterministic
% algorithms, and are thus not truly random: hence
% the nomenclature "pseudo-random".
% However, "good" random number generators, such as 
% those used by MATLAB / octave are designed so that 
% the sequences of numbers produced share many 
% statistical properties with true random values.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rand
%
%   Generates a single pseudo-random number, drawn from 
%   a uniform distribution between 0 and 1 (i.e. the 
%   probability of getting any number in the range is 
%   the same).  Also note that the generator will never
%   produce exactly 0 or exactly 1.  We will refer to
%   such a value as a "standard random number" 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  rand

ans =  0.73942

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note the use of commas to separate multiple stmnts
% on a single line in the following.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  for ir = 1:10, rand, end

ans =  0.20548
ans =  0.62264
ans =  0.53670
ans =  0.14572
ans =  0.24181
ans =  0.67596
ans =  0.070774
ans =  0.89507
ans =  5.7565e-04
ans =  0.33630

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rand(1,n)
%
%   Generates an n element row vector of standard 
%   random numbers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  a = rand(1,4)

a =

   0.79375   0.10076   0.16928   0.57121


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rand(n)
%
%   Generates an n x n matrix of standard random 
%   numbers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  b = rand(3)

b =

   0.452238   0.143796   0.341164
   0.182840   0.945288   0.356127
   0.647507   0.028950   0.110610


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rand(m, n)
%
%   Generates an m x n matrix of standard random 
%   numbers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  c = rand(2, 4)

c =

   0.88551   0.44796   0.60787   0.98743
   0.92627   0.77347   0.12980   0.77255


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% randi(imax)
%
%   Generates a uniformly distributed random integer 
%   in the range 1 to imax inclusive

%   NOTE! May be useful in homework.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  for ir = 1:10, randi(3), end

ans =  1
ans =  3
ans =  2
ans =  1
ans =  1
ans =  1
ans =  2
ans =  2
ans =  1
ans =  2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Generating Uniformly Distributed Pseudo-Random 
% Numbers on an Arbitrary Range

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To generate pseudo-random numbers uniformly 
% distributed on the interval (a,b) use

% a + (b - a) * rand 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: Generate 10 random numbers on the interval
% (-1,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  vr = -1 + 2 * rand(1,10)

vr =

 Columns 1 through 8:

  -0.28236   0.89734  -0.14936   0.22884  -0.97204   0.14496  -0.90659   0.14713

 Columns 9 and 10:

   0.47531  -0.46155


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Normally distributed pseudo-random numbers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Both MATLAB and octave have a function randn which
% operates analogously to rand, but which generates 
% pseudo-random numbers that are normally distributed 
% with zero mean and unit standard deviation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate row vector with 20 pseudo-random numbers

% randn(20) would retutn a 20 x 20 matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  vrn = randn(1,20)

vrn =

 Columns 1 through 7:

   1.406586   0.236024   0.763061   1.468343   1.116320   0.062045   0.593305

 Columns 8 through 14:

   0.762610  -1.093805  -0.950297   0.341393  -0.151328   0.469343   1.011659

 Columns 15 through 20:

   0.283399  -0.518756  -0.382473  -0.243508  -0.925791   1.638263


>>  for i = [10 100 1000 10000 100000], vrn = randn(1,i);, ...
mean(vrn), std(vrn), disp(''), end

ans = -0.73079
ans =  0.82566

ans =  0.039684
ans =  0.97673

ans =  0.012221
ans =  1.0324

ans = -0.0024515
ans =  0.99253

ans =  4.3502e-04
ans =  1.0010


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.3 Array (Matrix) Division
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix division in MATLAB is also associated with
% the rules of linear algebra; and specifically with 
% the solution of matrix equations.  We thus first 
% review some pertinent definitions and concepts.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Identity Matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% The identity matrix is a square matrix in which the
% elements along the main diagonal are 1, and all other
% elements are 0.  As we have already seen, the n x n
% identity matrix can be generated in MATLAB using 
% the eye command.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

>>  I = eye(3)

I =

Diagonal Matrix

   1   0   0
   0   1   0
   0   0   1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% In linear algebra, the identity matrix is to general 
% matrices and vectors as the constant 1 is to real
% numbers, i.e. for any n x m array, A, multiplication
% on the left by the n x n identity matrix I yields
% the original array, A, that is 
%
%  I A = A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

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

A =

   1   2
   3   4
   5   6


>>  I * A

ans =

   1   2
   3   4
   5   6


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% In the special case that A is n x n (i.e. square),
% we have

I A = A I = A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

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

A =

   1   2   3
   4   5   6
   7   8   9


>>  I * A

ans =

   1   2   3
   4   5   6
   7   8   9


>>  A * I

ans =

   1   2   3
   4   5   6
   7   8   9


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Inverse of a Matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
% If A is a n x n (square) matrix, then the n x n (square)
% matrix B is the inverse of A if and only if
%
A * B = B * A = I

% where I is the n x n identity matrix.

%                                    -1
% The inverse of A is often denoted A   and can be 
% computed in MATLAB either by raising A to the -1
% power, or by using the inv function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

>>  A = [2 1 4; 4 1 8; 2 -1 3]

A =

   2   1   4
   4   1   8
   2  -1   3


>>  B = inv(A)

B =

   5.50000  -3.50000   2.00000
   2.00000  -1.00000   0.00000
  -3.00000   2.00000  -1.00000


>>  A * B

ans =

   1   0   0
   0   1   0
   0   0   1


>>  B * A

ans =

   1   0   0
   0   1   0
   0   0   1


>>  A^-1

ans =

   5.50000  -3.50000   2.00000
   2.00000  -1.00000   0.00000
  -3.00000   2.00000  -1.00000


>>  A * A^-1

ans =

   1   0   0
   0   1   0
   0   0   1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Not all square matrices have inverses.  Indeed, 
% the necessary and sufficient condition for the n x n
% matrix A to have an inverse is that its determinant
% be non-vanishing.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Determinants
%
% I assume that you have seen the definition of the 
% determinant of a matrix, if not, refer to any text
% on linear algebra (or wikipedia)
%
% In MATLAB, the determinant of a (square) matrix can 
% be computed using the det function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  A = [1 2; 3 4]

A =

   1   2
   3   4


>>  det(A)

ans = -2

>>  A = [3 4; 1 2]

A =

   3   4
   1   2


>>  det(A)

ans =  2

>>  A = [4 3; 2 1]

A =

   4   3
   2   1


>>  det(A)

ans = -2

>>  A = [2 1; 4 3]

A =

   2   1
   4   3


>>  det(A)

ans =  2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Array (Matrix) Division

% MATLAB has 2 types of matrix division, namely
% right and left division.  Both are best defined
% in the context of solving matrix equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Left Division \
%
% Left division is used to solve the matrix equation

A X = B
%
% where A is an n x n matrix, while X and B are column
vectors of length n.

% This equation can be solved by left-multiplying both 
% sides with the inverse of A:
%
%  -1       -1
% A  A X = A  B

%           -1 
%    I X = A  B
%
%           -1
%      X = A  B
%
% This last expression can be defined in terms of the 
left division operator \ as
%
X = A \ B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Right Division /
%
% Right division is used to solve the matrix equation

%  X C = D
%
% where C is an n x n matrix, while X and D are row
vectors of length n.

% This equation can be solved by right-multiplying both 
% sides with the inverse of C
%
%       -1     -1
%  X C C  = D C
%
%              -1 
%     X I = D C
%
%              -1
%       X = D C
%
% This last expression can be defined in terms of the 
right division operator / as
%
X = D / C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Example: Solving three linear equations
%
% Use matrix operations to solve the following system
% of linear equations
%
%    4x -  2y + 6z = 8
%    2x +  8y + 2z = 4
%    6x + 10y + 3z = 0
%
%  This can be solved in either of the following forms

%  -            -  -   -   -   -
%  |  4  -2   6 |  | x |   | 8 |
%  |  2   8   2 |  | y | = | 4 |
%  |  6  10   3 |  | z |   | 0 |
%  -            -  -   -   -   -
%
%               -            -         
%  -         -  |  4   2   6 |    -         -
%  | x  y  z |  | -2   8  10 |  = | 8  4  0 |
%  -         -  |  6   2   3 |    -         - 
%               -            -          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>  A = [4 -2 6; 2 8 2; 6 10 3]

A =

    4   -2    6
    2    8    2
    6   10    3


>>  B = [8; 4; 0]

B =

   8
   4
   0


>>  X = A \ B

X =

  -1.80488
   0.29268
   2.63415


>>  Xb = inv(A) * B

Xb =

  -1.80488
   0.29268
   2.63415


>>  C = [4 2 6; -2 8 10; 6 2 3]

C =

    4    2    6
   -2    8   10
    6    2    3


>>  D = [8 4 0]

D =

   8   4   0


>>  Xc = D / C

Xc =

  -1.80488   0.29268   2.63415


>>  Xd = D * inv(C)

Xd =

  -1.80488   0.29268   2.63415