lab2

Math 181
Computational Mathematics
Lab 2

R. Robinson
Due 1/31/01

Welcome to Lab 2. In this lab we will learn very basic MATLAB.

Basic Linear Algebra in MATLAB

Matlab was originally conceived as a general purpose matrix manipulationtool for teaching purposes. It was based on the public domain linear algebrasubroutine package LINPAK (written in FORTRAN). However it quikly caughton with numerical analysts, scientists and engineers. Since then it has been vastly extended, and has evolved into the computing environment of choice in many scientific and engineeringapplications.

MATLAB predates the "symbolic algebra" packages like Maple and Mathematica. Like those packages it is interpreted rather than compiled. Also it has many built in high level mathematical functions. It can be used in a "super calculator" mode and also has a good built in programming language. This saves a lot of devbelopmenttime compared to FORTRAN or C.

Unlike Maple and Mathematica, MATLABis based on machine floating point numbers, which makes it FAST. In myexperience it can frequently solve problems that Maple or Mathematica istoo slow to solve.

To start MATLAB, click on the MATLAB icon.

Please make a record of your session. This is what you will pass it in (i.e. post on your web-site). Like the Congressional Record, you may edit it after the fact!

>> diary lab2.txt

Note that the Matlab comment character is a "%". Please add your name the date and the name of the lab to the file as follows:

>> %Y. Name
>> % Lab 2
>> % 1/24/01

In Matlab everything is a matrix. A number is a 1x1 matrixand a vector is a nx1 or 1xn matrix. All numbers are floating point .

Entering Matrices.

>> A=[1 1 1; 2 0 -1; 0 1 1 ]

(I've supressed the outputs for the sake of suspense!)

Note that row elements are not delimited (optionally you can put in commas). Columns are separated by semicolons.

Here's the RHS in Ax=b

>> b=[6 3 4]

Note that we get a row vector! I can fix this by transposing.

>> b=b'

The lab assignment for today is as follows. Create a diary of today'ssession that, in addition to the lab lecture, contains answers to all thequestions. You need not answer the questions correctly the firsttime or even at the same time. You can edit your diary files into a "amooth"version to pass in (they are ASCII files).

Question 1: How could you enter column b without taking the transpose?

Now we look at several ways of solving Ax=b.

>> inv(A)

>> x=inv(A)*b

An alternative to this would have been x=ans*b.

Several comments:
"inv" is matrix inverse (you can also type A^(-1)). "*" is matrix multiplication. Of course shape restrictions apply.

Question 2: How would you multiply A by its inverse? Try it.

Note: You can multiply two matrices of the same size entry by entry using ",*". Compare the two commands A*A and A.*A (the latter has a "."before the *). How are they different?

Here is an alternative way to solve this system

>> x=A\b

...this is matrix division.

The way you were probably taught to solve systems in your linear algebraclass was to construct the augmented matrix and perform row reduction (the Gauss-Jordan algorithm).Let's try this here.

Matlab notation makes the augmentation process easy (here we augmentA with b).

>> Aaug=[A b]

Then we perform row reduction. Notice the solution to Ax=b.

>> rref(Aaug)

Do you remember how to extract the solution form this? What if the matrix A were singular?

By the way, MATLAB has a good help file. Type

>> help rref

Timing.

If you type

>> flops

you will get a count of the total number of floating point operations used so far (since starting matlab). Typing this before and after a calculation and subtracting gives a good way to "time" a computation.

Question 3: Which of the 3 ways to solve Ax=b above is faster?

Again type

>> inv(A);

Make an identity matrix

>> B=eye(3)

and let

>> AB=[A B]

>> RAB=rref(AB)

What do you see on the right hand side?

Suppose we want to strip off the left half? MATLAB has a good way of doing this.

To understand it we first learn some basic MATLAB commands.

Type

>> AB(1,1)
>> AB(1,2)
>> AB(1,3)

Do you see how this works?

Now type

>> 1:3
>> 1:.5:3

Do you see how this works? How about this?

>> AB(1,1:6)

Question 4: How do you get the inverse from AB

Question 5: Solve the problem x+Cx=d, where

>> C=[.5 .4 .2; .2 .3 .1; .1 .1 .3]

Use the demand vector

>> d=[50; 30; 20]

This is an example of an "Leontief Input-Output" model from Economics. C is called the consumption matrix, d is the demand and x is the production. See any elementary linear algebra book for a description.

Now we move on to singular systems. The easiest way to guarentee
singularity is to have more variables than equations (underdetrermined)

>> A=[1 1 1; 2 0 -1]

>> b=[6; 3]

>> Aaug=[A b]

>> rref(Aaug)

To make this easier to read try

>> format rat

To restore the default format (short)

>> format

You can now solve Ax=b by hand (do you remember how?). There are infinitely many solutions. Somewhat surprisingly, MATLAB has no built-in routine to do this (Maple does). Perhaps this is because Gauss-Jordan is not viewed as the optimal
numerical algorithm to solve this problem.

Here is another example:

>> A=[1 2 8 7; 2 4 16 14; 0 1 3 4]

>> b = [18; 36; 8]

Find the general solution to Ax=b using

>> rref([A b])

Now we want to concentrate on the homogeneous equation Ax=0.

The solution set to this is called the nullspace. We can find it via

>> b=zeros(3,1)
>> rref([A b])

But note that the last row is redundant. It's enough to look at just

>> rref(A).

By the way the number of nonzero rows is the rank of A.

Question 7: Can you find the general solution to Ax=0 from this?(Use paper and pencil). It will have two parameters and will be a linear combination of two vectors. These two vectors are a basis for the nullspace.

Please read the description of the singular value decomposition. This is a PDF file and you'll need to install the free Adobe Acrobat Reader to view it (if it isn't already installed on the computer ypou are using). By the way, I produced this document in "LaTeX" and you can view my TeX source file (as ASCII text).

Question 8: (After class) Solve the problems from questions 1 2 & 5 using svd.

Question 9: Find the best linear and quadratic polynomial (least squares sense) through the points. Use a graph (see below) to see which is best:

(a) (2,3), (3,2), (5,1), (6,0), (7,0)

(b) (1,1.8), (2,2.7), (3,3.4), (4,3.8), (5,3.9)

A few other Matlab concepts.

Graphing

>> x=0:.1:20
>> y = (x.^2).*sin(x)

(put a colon at the end to avoid the printout)

>> plot(x,y)
>> xlabel('time')
>> ylabel('money')
>> title('Sample plot')
>> print plot1.eps -deps

The last command creates a postscript file withe the plot in it. You can put this on your web site (but it is large). It can be compressed by winzip or by converting to pdf.

Here is an example of a two dimensional plot z=f(x,y)

>> u=-5:.2:5;
>> [X,Y]=meshgrid(u,u);
>> Z=X.^2-Y.^2;
>> surf(X,Y,Z)
>> print plot2.eps -deps