Matlab Tutorial #3

Contents:
1.  For loops
2.  While loops
3.  IF statements
4.  integration (trapz)
5.  fmin
6.  fmins
7.  leastsq
8.  examples
	linear regression
9.  differential equations

1.  For loops

A "for loop" in any programming language is used to perform a set 
of operations a finite number of times.  For example, assume
 you want to start with a number, say x=2, and double this number 
eight times.  This means that the "for loop" must perform a 
simple multiplication, x*2, eight times.

x = 2;
for i=1:8
	x=x*2;
end

This will give us one final 'x' value.  Notice that i is not used in the 
operation, but only as an iteration variable. 

Notice the "end" statement.  The "for" and the "end" encase the repeated 
operations.  Therefore, if you have six hundred statements in between
 "for" and "end", those 600 statements will be repeated eight times. 
 With that in mind, I will show you a generic "for" format

for iteration variable = start:stepsize:finish
	operations...
end


Now, suppose we want to store an array of the doubled values.  
'i' becomes more useful.  We'll use 'y' as the storage for 2*x.


x=2;
y=[];	% this defines an empty array, that we will fill in the for loop
for i=1:8
	y(i) = x*2*i;
end

Not that much more complicated.  i is used as the counter variable. 
 It is also used in the expression.  We could have done this the first time if we had wanted.  What's the difference?  The value of 'x' is maintained.

2.  While loops

While loops also perform a set of operations repeatedly.  However, 
unlike the 'for loop' which functions a finite number of times, 'while' 
performs until a certain condition is met.  This condition can be greater
than, less than, or equal to.  Here is the generic format.  Notice the 'end'
statement again.

while (some condition)
	operations
end

Let's consider our "doubling" example above.  Suppose we want to double 
x until we reach 30.

x=2;
while (x<30)
	x=x*2;
end

Now, let's add in the complexity of making this 'y' array.

x=2;
y = [0];
i = 2;
while (y(i-1)<30)
	y(i) = x*2*i;
	i=i+1;
end

Why did I start at i=2?  To avoid two problems.  
Problem 1.  I iterate i after performing the operation.  If I check y(i),
 it's value really doesn't exist yet, so I have to check y(i-1).  Simple
Problem 2.  If I start at i=1, the first value I check will be y(0).  This
 value doesn't exist and matlab will give you an error message.  My solution,
 was to start at i=2.  This will give me a dummy value of 0 as the first 
element in the y array.  As long as I'm aware of that, it really doesn't 
matter.  I have still solved the problem.

3.  The IF statement

The 'if' statement is a one-time conditional statement.  (Of course, 
you can use it inside loops).  The generic structure is as follows.

if (some condition)
	operations...
end

Is this starting to look familiar?  Again, these "some conditions" can 
be greater than, less than, or equal to.  You can also have if-else 
statements

if (some condition)
	operations
else
	other operations
end

And finally, you can have if, if-else, ... else

if (some condition)
	operations
elseif (some other condition)
	operations
else
	operations

For a simple example,  let's suppose we want to perform different 
operations on a number 'x' depending on it's value.  If it is less
 than 100, it is too small so we'll add 50, if it is between 100 and
 500 we only want to add 5, and if it's greater than 500, it's too big,
 so we'll subtract 100.

if (x<100)
	x = x+50;
elseif (x<500)
	x = x+5;
else
	x = x-100;
end




4.  Numerical integration.

Right now, you really know all you need to know to do numerical integration.  You can pull out your Schaum's outline or some other book on numerical methods and write the algorithms required to perform numerical integration using the 'for loop'.  However, Matlab does have some integration functions.  The most popular is trapz.m.  This is the trapezoidal rule and it's very simple to use.

integrated value = trapz(independent variable array, dep var array)
	- or -
integrated value = step size*trapz(array to integrate)

So, let's integrate a simple function, y = x^3-2x^2+7, from x = 0 to x=10.

x = [0:10];	% This makes an array from 0 to 10 with a step size of 1
y = x.^3 - 2.*x.^2 +7;	% this creates the array to integrate.
intval = trapz(y);

Easy right?  How good is this integration?  Would it be better if 
we used a smaller step size?  Probably

x = [0:0.1:10];
y = x.^3 - 2.*x.^2 +7;
intval = 0.1*trapz(y)

How do the values compare?  What should I do to make sure I have the 
best trapezoidal integration possible?  Simple make the step size smaller
and smaller until it makes no difference.  0.1 is probably sufficient, but
try 0.01 and see if 'intval' is different.

5.  fmin

Very often, you will want to find the minimum or maximum value of a 
function, for optimization purposes or solving equations and unknowns.  
'fmin' will find the value of ONE variable to minimize a function.

For example, suppose you want to the x value that minimizes the value of a 
parabolic function, y = 2*x^2 + 4.  This is very simple to solve 
analytically, but we'll do it in Matlab.

In order to use fmin, you will need to create a function file.  This is
 the file that stores your expressions, in this case just the one equation
 for y.  

Look at the help information for fmin and notice that you have several 
options.  I wont go into all of the options here, but notice the first
 paragraph.  fmin will find the minimum value for a function within a 
specified range.

Let's write the function file.

function y = parabola(x)
y = 2*x^2 + 4;

on the command line type
minimum = fmin('parabola',-10,10)

This will search for the and x value between -10 and 10 that minimizes
 the value of y.  The result, of course, is zero.

Let's edit our function file and put in a slightly more complex function.
y = 2*x^2 -4*x  + 4

Now the result is 1.  Notice that if you put in different ranges for x,
 you will end up with different results.  For example, use the range 10 to 
100.  This is the difference between local, and global minimums.  This 
concept can be extremely important depending on the problem you are trying 
to solve.

6.  fmins

Let's move on to the more powerful, and slightly more complex function,
 fmins.  You can perform the same type of operation with fmins, but you can 
search for more than one variable.

Read the help file for fmins.  Notice again that you have a lot of 
flexibility.  The main differences between fmin and fmins is that you 
can search for more than one value, and instead of supplying a search 
range, you supply initial guesses.  You still have to consider the problem 
of local and global minima.  You can find different roots by supplying
 different initial guesses.

Let's look at some different applications for fmins.  Instead of minimizing
 a function value, let's focus on solving equations and unknowns.  We'll
 use the above equation for y.  Suppose we want to find the value for x,
 that will result in a value of 5 for y.  We could solve the equation 
algebraically for x, but that would require time and effort.  As graduate
 students, we'll try to avoid this as much as possible.  The alternate 
solution is to allow fmins to search for the solution.

First, we'll write the function file.
function error = search(x)
y = 2*x^2 - 4*x + 4;
error = (5-y)^2;

Why did I square the error?  To avoid the problems of negative values. 
 This will get me as close to five as the tolerance of fmins.  If you read
 more into fmins, you'll notice that this tolerance is adjustable.

Now select an initial guess.  Let's say x0 = 0.  Type in the following 
commands
x0=0
xval = fmins('search',x0)

The result I got was -0.2248.  Lets try an initial guess of 5 and see what 
happens.  Pretty much the same result.

7.  Leastsq

Let's solve the exact same problem using leastsq.  What's the difference? 
 The optimization algorithm is different.  Leastsq uses a Marquark (maybe
 that's spelled right) method to find the solution.  fmins uses a different
 method.  Optimization methods are beyond the scope of this tutorial, but
 it is important to keep these in mind when solving real problems, so that
 you can use the most effective method.

Look at the help information on leastsq.  Notice that you should return
only the 'error' value, not the square of the 'error' value.  leastsq will
 take care of the squaring.

Write the function file:
funtion error = search2(x)
y = 2*x^2 -4*x + 4;
error = 5-y;

You still have to supply an initial guess, so enter the following commands
x0=0
xval = leastsq('search2',x0)

My result was -0.2247.  This is the same as fmins gave us.  At this point, 
we probably feel pretty confident about this result.

8.  Practical example - Regression

This is going to be a rather lengthy example of non-linear regression.  But
 it will show you how to solve for more than one variable using fmins and 
leastsq.  First, we will have to create some data and then add error to our
 measurements.  We will write a program file to solve this problem.  I will 
walk you through it step by step with comments.  Remember each comment is
 preceded by a '%'.  You do not have to type my comments into your program,
 but it might help you if you're looking at it later.

%Here's our program file
% demonstration file for regression
% first make some data  ( I chose this data totally arbitrarily
x1 = [1:10]';	% I'll use the ' because I prefer to work with columns
x2 = [10:-1:1]';
x3 = [1:0.1:1.9]';

% now I'll make up some parameters
b1 = 2;
b2 = 2.5;
b3 = 1.5;

% now i'll make up an equation
Y = b1*x1 - b2*x2 + b3./x3;

% assume that Y is some measured value in an experiment.  Let's add some error
% to these measurements and see if we can estimate the b values.  Remeber
% realistically we wouldn't know these values
error = rand(size(Y))*0.2 - 0.1
% this generates a vector of errors.  The vector is the same size as Y and holds values
% between -0.1 and +0.1.  This will generate a maximum of 10% error.
y = Y + Y.*error;
% ok, y is our new measurement vector.  Let's use our data and our y values to 
% estimate the b parameters

% Let's first set up an X matrix where the columns are the data sets
X = [x1 x2 x3];

% now let's make a vector of initial guesses
bg = [1 1 1];

% first, we'll solve the problem using fmins.  This will be a little more 
% advanced than the last section, but I'll explain each element in
% the command
params = fmins('freg',bg,[],[],X,y)
% for this problem, fmins requires some addition information.  
% It's going to need the data matrix, and the actual y values to compare
% with it's calculated values the []'s are just place holders.  You 
% will see that they have to be placed there by re-reading the help 
% fmins information.  The first [] takes the place of the options
% vector, (we are not changing any options), the second [] is a place holder.
% X and y can now be used in the function 'freg' which we will define later.

% we'll also estimate the parameters using leastsq
params2 = leastsq('sqreg',bg,[],[],X,y)
% the []'s serve the same purpose as explained above.  This way 'sqreg' 
% can use X and y.

ok, we're done with the main program file.  Now, we have to write two 
function files, 'freg', and 'sqreg'.  Let's start with 'freg'.

function error = freg(b,X,y)
x1 = X(:,1);	% all rows, first column
x2 = X(:,2);
x3 = X(:,3);

yhat = b(1)*x1 - b(2)*x2 + b(3)./x3;

error = sum((yhat-y).^2);
% we're minimizing the sum of squares, as we learned to do in 
% statistics (think back!).

Next, we'll write 'sqreg'
function error = sqreg(b,X,y)
x1 = X(:,1);	% all rows, first column
x2 = X(:,2);
x3 = X(:,3);

yhat = b(1)*x1 - b(2)*x2 + b(3)./x3;

error = (yhat-y);
% remember, leastsq takes care of summing and squaring.

Save the program file as regdemo.m or whatever you want. 
Run the program file. Look at:
 params
 params2

Wow, these are really lousy estimates. Hmmm!  Go into your program file and
change the bg vector to [2 2.5 1.5].  This is cheating I realize, but I'm 
trying to make a point.
Run the program again and look at the parameter estimates.  They're still
 not very good are they?  Remember we did put in 10% error.  The moral of
 the story, is your measurements should be a little better.

NOTE:  if you're doing linear regressions, there's a matlab function called 
regress.m.  Make sure you read about it.  It's quite powerful and very easy 
to use.

9.  differential equations

Hey, we're almost done.  The last function in this tutorial is ode23.  
It is a differential equation solver.  Read the help files on ode23 and 
odefile.  You will see that this differential equation solver can do a 
lot of work for you.  Let's start with a very simple example.

Suppose you are tracking some object whose position changes with time. 
 All you know about the object is its velocity and it's starting position 
y=0.  The velocity equation dy/dt= 0.4*t - 3.  

If you've read the help file, you'll notice that you'll need to write a 
function file, and define some variable called TSPAN.  TSPAN can be either 
the span of the independent variable, or an array of specific values. 
 Let's start with setting a span from 0 to 100.
TSPAN = 100

Now we'll need a function file.  Follow the format set up by matlab.  Did
 you read odefile?  
function result = position(t,a)
result = 0.4*t-3;

You must pass your function file an independent and a dependent variable,
 even if you don't use both of them.  And they must be in that order;
 independent then dependent.  In this equation dy/dt depends only on t.  
In this case, we will not use 'a'.  This is just an array of values that 
ode23 comes up with during it's solution.  We'll use it later.

Now type in the command
[T P] = ode23('position',TSPAN,0)
0 is our initial value.  If the object started at some other position,
you would put in it's coordinate.
Now type
[T P]
Look at how the position changes with respect to time.  Matlab chose time 
values within your range.  Try the same exercise only change TSPAN to a
 vector [1:100] and see what happens.

Now let's suppose that the object is moving in two dimensions, x and y. 
 You know the change in x with respect to time, dx/dt = -0.1*t + 4;
Change the function file so it looks like this
function result=position(t,a)
yval = 0.4*t - 3;
xval = -0.1*t +4;
result = [xval;yval]

IMPORTANT:  always return a COLUMN vector to ode23.  

Run the commands:
TSPAN = 100;
init = [0;0]	% we are making a column vector of initial values.  Remember
		% we must have one initial value per variable (x,y)
[T P] = ode23('position',TSPAN,init)

Look at P.  The first column represents x and the second represents y. 
 Why are they in this order?  Because this is the order we passed them 
to the ode23, (see result above).

This works quite nicely for independent equations.  We can also solve 
dependent equations.  Assume the y position was dependent on the x position
 and vice versa.  Let's make up some new equations
dy/dt = 0.4*t - 3 + x
dx/dt = -0.1*t + 4 -y

Now we'll have to use those extra inputs in our function file.  Let's
 also get a few things clear in our mind to avoid confusion.  When
 ordering inputs and outputs (result), we'll make our x coordinate the 
first value and our y coordinate the second value.

Again, change the function file
function result = position(t,a)
yval = 0.4*t - 3 + a(1);
xval = -0.1*t + 4 - a(2);
result = [xval;yval];

Hey what's this 'a' stuff?  These are the ode calculated values.  'a' is
 a vector remember.  The first value is x and the second is y.  We decided 
on this order when we defined result and we have to be consistent.  Now 
type some commands
TSPAN = 100;
init = [0;0]
[T P] = ode23('position',TSPAN,init)

You will get a P matrix of x and y positions at various times.

Congratulations on completing this rather lengthy tutorial.  There are lots
 more things to do in MATLAB, but this should get you started.

Giving Credit where credit is due:
These MATLAB tutorials are a service of Omega Chi Epsilon, the honor
 society of the chemical engineering department.
Mat3 tutorial was written by: Kimberly Rogers