MATH 353
Matlab 2

MATLAB is a very powerful and interactive software package for numerical computations. It is based upon matrices and m-files.

First, log on to your generic account {student} & {m4th4UD} . If you are unable to do that, please let me know.

Open MATLAB. First, go to the home directory by typing (in the Command Window)

>> cd ~/

Then type

>> mkdir mywork (or whatever you wish to name this directory)

and

>> cd mywork

M-files (following ``Matlab Guide'', D. Higham & N. Higham).

An M-file is a text file that has a .m extension and contains MATLAB commands. There are two types of M-files:

Script M-files have no input or input arguments and operate only on the variables in the current workspace. A script file enables you to store a sequence of commands that are to be used repeatedly or will be used at a future time. One example of script file is my_plot.m from Lab #1. Other interesting example can be found at http://www.math.udel.edu/~bacuta/M353/MatlabCodes.html

under `` Warming up''.

Function M-files contain a function definition line and can accept input arguments and return output arguments. Their internal variables are local to the function. Function M-files enable you to extend the MATLAB language by writing your own functions.

Here is an example of Function file. Type the next lines in the MATLAB editor:

------------------------------------------------------------------

% Input : r - a positive number (or a vector with positive entries)

%Output1: A -area of the disc of radius r

%Output2: C -the length of the circumference of the circle of radius r

function [A, C] = circle(r);

A= pi *r.^2;

C=2*pi*r;

----------------------------------------------------------------------

For this function, ``r'' is the input argument, and ``A,C'' are output arguments.

Save the file as circle.m in mywork. In the Command Window, type:

>> help circle

>> circle(1)

>> [a, c] = circle(1)

What is the difference in output for the above two commands?

>> r=linspace(0,4,5)

>> [a, c] = circle(r)

The dot from r.^2 operation allows r to be a vector. If r is a vector, then r.^2 is a vector of the same dimension whose components are the 2nd power of the corresponding components of r.

Here is another example of MATLAB implementation of a simple real function of two variables (f(x,y)=3x + 2y^2 ). The name of the function is fun2.m.

----------------------------------------------------------------------

% function file for the two variables function z=fun2(x,y)=3x + 2y^2

%Input: two vectors x and y (of the same length)

%Output: a vector z of the same length with x and y.

function z = fun2(x, y);

u=3*x; w=2*y.^2;

z=u+ w;

----------------------------------------------------------------------

Save the file as fun2.m in mywork. In the Command Window, type:
>> zz=fun2(1,2)
>> fun2(1,2)
>> x=ones(1,4)
>> y=2*ones(1,4)
>> zz=fun2(x, y)
>> u

You should get ``??? Undefined function or variable 'u'.'' Why?

The uncommented part of the function can be simplified to:

----------------------------------------------------------------------

function z = fun2(x, y);

z=3*x + 2*y.^2;

----------------------------------------------------------------------

In many situations one function needs to be passed as argument to another function. One way to do that is through a function handle. A function handle is a MATLAB data type that contains all the information necessary to evaluate the function to be passed as an argument. It is created by putting the @ character before the function name. Try:
>>ezplot(@sin)
>>feval(@fun2,1,2)
>> f=@fun2
>>feval(f,1,2)

The MATLAB function fzero.m (>>X = fzero(FUN,x0) tries to find a zero of the function FUN near x0. We can apply this function to find the roots of a simple polynomial function like f(x) =x^2-4. The function f(x) =x^2-4 is not one of the elementary functions already implemented in MATLAB. Thus, you have to create a .m function named, say fun4u.m . The body of the function could be as follows:

----------------------------------------------------------------------

%nice comments here.

function y=fun4u(x)

y=x.^2 -4;

----------------------------------------------------------------------

(type the above three commands in the MATLAB editor and save them as fun4u.m

Next, in the Command Window type
>>fzero(@fun4u,1)
>>fzero(@fun4u,-1)

To avoid writing a .m file implementing the function f(x) =x^2-4, another alternative is to call the function as an ``inline" expression:
>> fun4me='x.^2-4';

OR >>fun4me= inline( 'x.^2-4', 'x')
>> root=fzero(fun4me,8)
> ezplot(fun4me)

----------------------------------------------------------------------

Another way to work with functions is to define them as anonymous functions. Here is one example (finding the root of f(x)=x^-4 close to x=3 using MATLAB's fzero.m :


>> f= @(x) x.^2-4;
>> root=fzero(f,8)

----------------------------------------------------------------------

MATLAB supports the usual programming structures ``for'', ``if'' , ``while'' and ``switch.'' The syntax is a bit different from C syntax, but the use is much the same. Information about any of them may be obtained from the help command. The following ``for'' statement sums 0.1 n times (assuming that n is given).

sum=0;
for k=1:n
    sum=sum+0.1;
end;

Incorporate the above lines in a function named, say mysum:

----------------------------------------------------------------------

function S=mysum(n)

S=0;

for k=1:n

    S=S+0.1;

end

----------------------------------------------------------------------

What should be the name of the function? Now type:
>>format long
>>S=mysum(10);
>>S=mysum(1000);

Why is the answer NOT S=100?

Other examples of functions which contain the ``for'' structure are Table1.m and fixpt.m (see the class home-page). The function ``bisect2.m'' contains ``if'' and ``while'' statements.

Here is a simple example of how to use the ``while'' structure. Type:
>> n=101;
>> while n>0   %hit enter
disp(sprintf('at this step n=%d and sqrt(n)=%2.5f ', n, sqrt(n))) %hit enter
n=n-4;   %hit enter
end   %hit enter

Comments? Make sure that you understand what each of the above command lines do.

To illustrate the ``if'' structure, type:
>> help if

Copy the statements given by MATLAB help for ``if'' in the Example and create the following function named stiff.m.

----------------------------------------------------------------------

function A=stiff(n)

%%Building a square matrix which has 2 on the main diagonal,

%-1 on the diagonals above and bellow the main diagonal, and

% 0 off the first three diagonals.

for I=1:n

    for J=1:n

          if I==J

                A(I,J) = 2;

          elseif abs(I-J) == 1

                A(I,J) = -1;

          else

                A(I,J) = 0;

          end

    end

end

---------------------------------------------------------------------- Now type
>>A=stiff(8);

Do you get a tridiagonal matrix?

Make sure that you understand what each of the above command lines do.

Feel free to ask (many) questions.        

----------------------------------------------------------------------

Finding roots by using a fix-point theorem

Let us try to use the fpi.m function discussed in class, and fixptpt.m to find the fixed point of the function g=exp(x/4) in the interval (0, 2). First you need to write or download the function fpi.m from the class home page. Second, the function g(x) = exp(x/4) is not one of the elementary functions already implemented in MATLAB. Thus, you have to create a .m function named, say myexp.m . The body of the function could be as follows:

----------------------------------------------------------------------

%nice comments here.

function y=myexp(x)

y=exp(x/4);

----------------------------------------------------------------------

Third, you will need to pass some values for the rest of the arguments of the function fpi.m. Type >>help fpi . For the initial guess, you should plot y=g(x) and y=x, on (0,2) in the same system of coordinates. One fast way to do so is:
>>x=linspace(0,2);
>>y=myexp(x);
>> plot(x,y,x,x)

From the plot one can easily guess that a good initial guess is x0= 1.

To find the fixed point of myexp.m using fpi.m, type
>> xc=fpi(@myexp, 1,10) IF you do not to create a .m file for y=exp(x/4), then you can define the function in the Command Window as an anonymous function:
>> sameexp =@(x) exp(x/4)
>> >> fpi(sameexp,1,10)

OR inline function:
>> samexp = inline('exp(x/4)')
>> >> fpi(sameexp,1,10)

Note the difference in the way the first argument is called for the function myexp.m vs. the anonymous or inline function sameexp.m

Finally, to find the fixed point using fixptpt.m type the following:
>>tol=1e-8
>> maxIter=100
[k, xc, err, x]=fixpt(@myexp,1, tol, maxIter)

How many iterations were performed?

----------------------------------------------------------------------