Lab (1) -Introduction to MATLAB
Rahman Tashakkori & Darren Greene, CS - Appalachian State University, Distributed Computing Workshop, 2006
Arrays and their operations
Matrix and their operations
Some Linear Algebra in MATLAB
MATLAB is a matrix-based application.  We will utilize the tools and commands that are already built into the MATLAB and its corresponding Toolboxes to do some simulations.   In this lab, we will learn some of those commands and we also review some of the basic concepts in linear algebra.

Preparations

MATLAB can be run on either our Linux server or lab PCs.  For this lab, we will use the Windows-based version.  So, please go to start, CS Software, MATLAB, and run MATLAB. You will get the MATLAB prompt in a couple of minutes.  Be patient, the MATLAB screen will eventually appear.

Arrays (Vectors)

Creating an array (1-D matrix):
To create an array, t, with five elements (1 through 5) use:

>> t = [1 2 3 4 5]

This creates the array t with 5 members, 1, 2, 3, 4, and 5.  Note that in this case all 5 members are sequential.  We could use:

>> t = 1:5

to create the same array.

We can also create an array by giving the beginning number, the end number and an step size (increment).  In the above example, the step size is 1, thus, we could do the following:

>> t = 1:1:5

In this case the middle number is the step size. Try the following command to see the difference:

>> t = 1:0.5:5

In this case the step size is 0.5.  There are many built-in functions in MATLAB, for example: sin, cos, sqrt, .... MATLAB even allows you to add your own function.

Note that the first array member in the MATLAB is t(1).  In order to access a particular member of an array we use the index corresponding to that member.  For instance, in array t the 3rd member is t(3).  What is the value of t(3)?  To find the answer simply type:

>> t(3)

Suppose we want to multiply t by itself.  You have already generated t as an array with one row and 9 columns.  Thus, you have to write the t as a 9 row by 1 column to be able to multiply it by itself.   Therefore, we need to transpose the second t.

To compute the transpose of a matrix we will use, '.  Here is how we create the transpose of t in MATLAB:

>> tp = t'

Transpose of a matrix is the same matrix when rows and columns are switched. Now you can multiply t by itself using:
r = t*tp;

There are some powerful predefined operators and functions in MATLAB that make our life much easier.  For example, if you want to compute the sum of all members in t, you can use the sum function:

mysum = sum(t)

One powerful feature of MATLAB is its plotting capabilities. 

Matrices

A matrix is what we refer to as an array with more than one dimension.  There are many commands designed in MATLAB to do matrix manipulation.  Entering matrices into MATLAB is the same as entering a vector, except each row of elements is separated by a semicolon (;) or a return.  Let's look at an example.

Suppose we have a 3x3 matrix, A, defined as;

A =
20 22 18
12 16 12
8  10  10

and another one with the same size, B, defined as;
B =
2 4 8
1 -1 6
-2 2 4

How do we enter these matrices in MATLAB.
First method (Note that by placing a ";" at the end of your command you can avoid displaying the outcome):
>> A = [20 22 18; 12 16 12; 8  10  10]

The output will looks like this:
A =
    20    22    18
    12    16    12
     8    10    10

Second method:
>> B = [2 4 8
1 -1 6
-2 2 4];   % this time no display

But, you can view B at any time:

>> B

The output looks like this:
B =
     2     4     8
     1    -1     6
     -2     2     4

Let's try a few of the operations in linear algebra here.  Suppose, you want to find the transpose of matrix A, rows are interchanged with columns, and call the new matrix C.
>> C = A'
This generates the following output;
C =
    20    12     8
    22    16    10
    18    12    10

To compute D = A * B.

>> D = A*B
D =
    26    94   364
    16    56   240
    6    42   164

So far, you have tried several of the MATLAB commands individually. You could put all of the commands that you want to execute together in an m-file, and just run the file. If you put all of your m-files in the same directory from which you run MATLAB, then MATLAB will always find them. To create an m file, you can use the MATLAB Editor (or any other word processor or editor) to type the commands and once you are done, save the file with .m extension.  Then at the MATLAB prompt type the name of the file without the .m to run the commands in that file.   Next we will create a .m file and will run it. But before that we need some theoretical background on what we plan to compute.

Computing pi using Monte Carlo Simulation

Computing pi is not a new problem but has always been an interesting one.  There are several ways to compute pi.  We look into two methods here.  One is based on the Monte Carlo simulation and the other is based on its expansion in series form.

Computing Pi - Monte Calrlo Approach

In order to compute the pi using Monte Carlo Simulation, we draw a unit circle and place the tightest square we can around it.  That would be a one-by-one square.

full circle
1/4 circle
Figure (1) - A unit circle and a unit square tightly surrounding it
1/4 of the circle

We know the area of the square.  Therefore, we can post the board on a wall and throw darts at it in a random fashion.  We only count those darts that hit within the square so we will get a chance to re-try those that hit outside the square.  Among the darts that we randomly throw, some will hit inside the circle.  If we throw a large number of darts, at the end the ratio of the darts inside the circle to the total number of darts that hits inside the square will give us a good estimate of what the area of the circle is.  This is shown on the left hand side of Figure (1).  We can simulate this on a computer by creating random x and y values and by determining whether they produce a radius smaller than (or perhaps equal to) the  radius of the circle.  If they do, we count them inside the circle and in the mean time we keep track of the total number of x and y points we have created. 

eq1

We can do this for 1/4 of the square and multiply the result by 4 to get the total count.  This is shown on the right-hand-side of Figure (1). So the above formula can be written as:

eq2 

General algorithm looks something like this:

Initialize the inside hits to 0, insidehits = 0
Generate a random x
Generate a random y

Compute the distance between the point (x,y) and the origin: d = sqrt(x^2 + y^2)

if (d <= r)   % r is the radius of the circle
     insidehits = insidehits + 1.0

Compute the ratio.
pi = 4*ratio

Following is a MATLAB code that computes pi using Monte Carlo described above.  You can either open the MATLAB Editor, cut and paste this program, and save the file with orgpi_comp.m name in your Work Directory.  Or simply Right Click on this file orgpi_comp.m and store it in your Work Directory. 

% Orgpi_comp.m
% Monte Carlo Computation of pi.

% n is the total number histories (random darts) we will use
n = input(' Enter n: ');

% To keep the number of histories (darts) that hits
% inside the circle
inside_count = 0;

%  Generate n random points in range [0,1]
%  scale it to X[0,1]. 1/4 of the circle is defined in that range as
%  x^2+y^2 <= 1 which will be approximately pi/4.

% Compute sum of x square for variance and STD for error
x_er = 0;

start_time = cputime;       % Start timing of algorithm

for i = 1:n,
     x = rand;
     y = rand;
     if x^2 + y^2 <= 1,
          inside_count = inside_count + 1;
          x_er = x_er+1;
     end;
end;

%To get the correct pi multiply the ratio by 4
pi_approx = 4*(inside_count/n),
err = pi - pi_approx,

x_er = x_er/n;

sig = 4*sqrt(x_er)/n;

elapsed_time = cputime - start_time;     %Calculate elapsed time to perform algorithm

fprintf('1 STD Error STD is: %f \n', sig)
fprintf('pi-STD is %f: \n', pi_approx-sig)
fprintf('pi+STD is %f: \n', pi_approx+sig)

fprintf('Computing time is: %f \n',elapsed_time )

Run the above program for the following n values and record the error and the computing time.

n
computed pi
Estimated Error (1 STD)
Computing Time
10



100



500



1000



10000



50000



100000




Computing Pi - Series Expansion

In order to compute this we use the discrete form of
integral
as a sum.

The following MATLAB function (piece.m) will compute the value of pi based on the above integration.  It computes the integral using a discrete method called Trapezoidal Technique:

% piece.m
function A = piece( p, np, n )
%-------------------------------------------
f = inline( '(4/(1+x*x))', 'x' );
%-------------------------------------------
h = 1 / n;  % delta x
A = 0;  % set the initial value to 0
for i = p : np : n
     x = (i-1)*h;
     A = A + (f(x)+f(x+h))/2*h;
end

And you can run the program after you saved it using piece(xlow, step-size, xhigh):

piece(1, 1, 100) 

which runs the function from 1 to 100 with step size 1.

When I tried this, I got:
ans =

    3.1416

Homework
1) Compute the variance and standard deviation for vector t that is created using t = 1:0.25:5.