Lecture 30: Intro to Matlab.
----------------------------

diary 'h:\matlec1.txt'
% the diary command produces a log file,
% this log file is only for us to see,
% not for matlab to use
a = rand(3,3)

a =

    0.9501    0.4860    0.4565
    0.2311    0.8913    0.0185
    0.6068    0.7621    0.8214

save 'h:\mydata'
% the save command saves all data in the workspace
% the mydata files is not readable, only for matlab
clear
% we just cleared the workspace
% the matrix a is gone, but we can load in a:
load 'h:\mydata'
a

a =

    0.9501    0.4860    0.4565
    0.2311    0.8913    0.0185
    0.6068    0.7621    0.8214

% by default, the type is double float,
% shown with four decimal places
format long e % long scientific notation
a

a =

  Columns 1 through 2 

    9.501292851471754e-001    4.859824687092997e-001
    2.311385135742878e-001    8.912989661489016e-001
    6.068425835417866e-001    7.620968330273947e-001

  Column 3 

    4.564676651683414e-001
    1.850364324822440e-002
    8.214071642952533e-001

format short % default
a

a =

    0.9501    0.4860    0.4565
    0.2311    0.8913    0.0185
    0.6068    0.7621    0.8214

% to solve a linear system, let us generate 
% a random right hand side vector b
b = rand(3,1)

b =

    0.4447
    0.6154
    0.7919

% we now have a and b defining a*x = b
% we can solve the system in three different ways:
x = a\b  % the operator notation

x =

   -0.0638
    0.6995
    0.3622

% let us verify whether x is a solution,
% by computing the residual r = b - a*x
r = b - a*x

r =

  1.0e-015 *

         0
   -0.1110
         0

format long e
r

r =

                         0
   -1.110223024625157e-016
                         0

% this residual is not quite zero,
% but of the magnitude the machine precision
eps

ans =

    2.220446049250313e-016

% the second way to solve a*x = b
% is via row reduction, or reduced row echelon form,
% provided in the rref command:
%   to apply rref, we first need to stack the 
% right hand side vector b to the matrix a:
ab = [a b] % constructs a new matrix

ab =

  Columns 1 through 2 

    9.501292851471754e-001    4.859824687092997e-001
    2.311385135742878e-001    8.912989661489016e-001
    6.068425835417866e-001    7.620968330273947e-001

  Columns 3 through 4 

    4.564676651683414e-001    4.447033643531942e-001
    1.850364324822440e-002    6.154323481000947e-001
    8.214071642952533e-001    7.919370374270354e-001

format short
ab

ab =

    0.9501    0.4860    0.4565    0.4447
    0.2311    0.8913    0.0185    0.6154
    0.6068    0.7621    0.8214    0.7919

rref(ab) % invoke the reduced row echelon form

ans =

    1.0000         0         0   -0.0638
         0    1.0000         0    0.6995
         0         0    1.0000    0.3622

% to verify the solution, we must select the
% last column from the answer matrix "ans":
y = ans(:,4) % select all rows, fourth column

y =

   -0.0638
    0.6995
    0.3622

x

x =

   -0.0638
    0.6995
    0.3622

norm(x-y)

ans =

  1.5944e-016

% we see that both solutions are the same,
% up to the machine precision
% The third way to compute the solution of a linear
% system is via the LU-factorization of the matrix a
[l,u,p]=lu(a)

l =

    1.0000         0         0
    0.2433    1.0000         0
    0.6387    0.5843    1.0000


u =

    0.9501    0.4860    0.4565
         0    0.7731   -0.0925
         0         0    0.5839


p =

     1     0     0
     0     1     0
     0     0     1

% we just decomposed a as a product of l (lower
% triangular) and u (upper triangular):
r = a - l*u % residual for the LU-factorization

r =

  1.0e-016 *

         0         0         0
    0.2776         0         0
         0         0         0

norm(r)

ans =

  2.7756e-017

% now we have to solve two triangular systems:
yy = l\b % solve L*y = b

yy =

    0.4447
    0.5072
    0.2115

xx = u\yy % solve U*x = y

xx =

   -0.0638
    0.6995
    0.3622

norm(x-xx)

ans =

  1.2413e-016

diary off % stops the writing to the log file