% reviewtut.m % A review tutorial for the first lab session of Inverse Theory ESS523, % to get students brushed up on both Matlab usage and linear algebra. % Includes material from Appendix A of Aster, Borchers, & Thurber, % chapter 1 of Gelb, and the very dredges of Andy's mind... % Andy Ganse, 1 Oct 2005 % --------------------------------------------------------------------- % JUST START READING AND FOLLOW THE INSTRUCTIONS AS YOU GO... % --------------------------------------------------------------------- % Okay all you inverse theory fans out there! % Let's get started with the ultra ultra basics of Matlab: % a.) Notice that these Matlab "script" files always end in ".m" % b.) You can run them in Matlab by entering the filename without the ".m" % at the Matlab command prompt. Alternatively you can open them in % the editor of the GUI version of Matlab and click the "Run" button. % c.) If you want to "comment out" a line in the script, i.e. use it for % writing rather than for code that gets run, put a "%" as the first % character. As you'll see below, Matlab actually considers anything % to the right of a "%" as a comment (note that includes the beginning % of the line). %%% %% By the way, when "commenting out" a block of code that already has %%%% % comments in it you may end up with several %s in there - no problem. % d.) There are several help mechanisms: % - to get info about a single function like cos() for example, enter % "help cos" at the Matlab prompt. % - to get a list of commands which relate to a subject like % convolution for example), enter "lookfor convolution" at the % Matlab prompt. % - online there is general documentation at: % http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml % and a users' guide at: % http://www.mathworks.com/access/helpdesk/help/techdoc/math/ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % As you go thru the tutorial in the rest of this file, you will "comment % out" and "uncomment out" various lines, and change some of the variable % values, and then rerun the script to note the changes in computations and % output. You can do this anywhere you like, but the lines labeled "HEY!" % in the file point out particular places to do so, to emphasize key % points. Deciphering which lines are code to uncomment, and which are my % notes and explanations, is perhaps part of the learning process - but % please don't hesitate to ask for clarification and help! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SETTING VARIABLES % ----------------- a=2.718 % stored as a double precision floating point number b=4e-6 % ditto c='yeehaw' % this is a "string", ie a word rather than number. % A string could actually contain numerical digits but it can't % be used in a numerical calculation (there are a few weird % exceptions to that but please don't worry about those for now) % Notice that without a semicolon at the end of the statement, Matlab echos % the values at the command prompt. % When we do vectors and matrices below, placing that semicolon at the end % of the line is crucial, so that (for example) a vector of length 10,000 % doesn't take 5 minutes to scroll down your screen. %% HEY!: Put semicolons after the statements above now and rerun the %% script to compare what happens. % A FEW STANDARD MATHEMATICAL OPERATORS AND FUNCTIONS % -------------------------------------------------- %% HEY!: keep the above lines specifying a, b, & c uncommented, and now %% uncomment the lines below: % d=a+b % e=3*a+4/b % disp('howdy'); % a handy function for outputting one-line notes % f=strcat(c,' man') % note there are bunch of these "str___" commands % % for manipulating strings of characters % g=4*log(b)+sin(pi*a) % note: pi is built-in, and sin() expects radians, % % and log() is ln... for base10 use log10(). % disp('Note the following produces an error: why?'); % stringnumber='1.5708' % h=sin(stringnumber) % why does this line produce an error? %% HEY!: what can you do to the previous two lines to fix the error? %% Lastly, remember you can always look up the use of other math functions %% with "help" & "lookfor" % COMPLEX NUMBERS % --------------- % Matlab accepts both "i" and "j" as the imaginary number. %% HEY!: now comment out everything above, and uncomment the lines below % a=1.3+4.5i; % b=6.7+2.1j; % c=a+b % phase=angle(c) % amplitude=abs(c) % foo=real(c) % boo=imag(c) % VECTORS AND MATRICES % (when starting this section, comment out all above, and gradually % uncomment lines a few at a time in this section, running each time) % ------------------------------------------------------------------- % Matlab will both calculate and display row and column vectors in the % correct directions. % rowvector=[1,2,3,4,5] % or also simply [1 2 3 4 5] % colvector=[1;2;3;4;5] % matrix=[1,2,3; 2,4,9; 7,11,13] % Recall what a matrix transpose is. Matlab uses the ' mark as a % transpose operator. %% HEY!: note the output of the transposes below to remind yourself what %% a transpose is (uncomment the next 3 lines): % rowvector' % colvector' % matrix' % % We can also create vectors by specifying a range: % rowvector=1:5 % colvector=(1:5)' % rowvector2=2:.1:2.6 % you can specify a step increment too. % colvectorlong=(1:1000)'; % what happens if you remove the semicolon here? % We can access just one element of a vector or matrix if needed: % p=rowvector(1) % q=colvector(3) % r=matrix(2,3) % (first arg here is row, 2nd arg is column) % We can specify "all" of a row or column like this: % s=matrix(2,:) % 2nd row, all columns % Recall we can multiply vectors and matrices, but unlike scalars, here % order is important % innerproduct=rowvector*colvector % also called dot product % outerproduct=colvector*rowvector %% HEY!: compare the innerproduct and outerproduct of vectors - which is %% a scalar and which is a matrix? What are the dimensions of the matrix? %% What is the geometric interpretation of a dot product of two vectors? % % Most of the other mathematical operators and functions work on vectors & % matrices too, just with some twists: % rowvector+rowvector % rowvector+2 % rowvector*2 % rowvector.*rowvector % the "." on * and / and ^ means element-by-element % disp('the following will cause an error'); % rowvector*rowvector % HEY!: why does this cause an error? % rowvector.^2 % matrix*2 % matrix.^2 % matrix^2 %% HEY!: closely investigate the output of the previous two lines to %% notice the different between using ".^" and "^" for the matrix. %% What is matrix^0 ? And matrix.^0 ? % % We can also multiply a vector and a matrix. This is the foundation % for a linear problem such as Ax=b: % x=[1,2,3]' % A=matrix % let's reuse the earlier "matrix". (is it uncommented?) % b=A*x % The above represents 3 equations and 3 variables in a linear problem. %% HEY!: what are the 3 equations? is the solution unique or could there %% be more than one to this problem? % % To take the inverse (if it exists) of a matrix, use the inv() function. % If the inverse doesn't exist or if the matrix is even "almost" singular % Matlab will tell you in a warning message. % Ainv=inv(A) % So if you have a linear problem like the above, and as long as the matrix % A is "nice", ie invertible (square, full rank, and so on), you could % solve for x like this if you didn't have it: % xsolve=inv(A)*b %% HEY!: Now let's mess up the matrix A and see what happens. Change its %% 2nd row to equal the first and try running it again: % A(2,:)=A(1,:) % b=A*x % xsolve=inv(A)*b %% HEY!: So what happened? Why'd the computation mess up this time? % % Now say we had four equations and three variables. This is an over- % determined problem; we have more pieces of information than variables % to solve for. In this case our matrix A has more rows than columns. % There may not be a solution that exactly satisfies all the equations so % we find a least-squares (LS) solution. % x=[1,2,3]'; % A=[1,2,3; 2,4,9; 7,11,13; 8,14,2]; % 4 rows, 3 cols % b=A*x+[.1,-.1,.1,-.1]'; % (sortof artificially making up the problem here!) % We can use our class material to solve this: % xbest=inv(A'*A)*A'*b % HEY!: how close is xbest to your original x? % % HEY!: Now what if we had three equations and four variables? This is % an underdetermined problem, where you don't have enough information to % uniquely solve the problem, so you have an infinite number of solutions. % An example might be "fit the best plane to two given points". % % Anyways, many math functions work as usual on vectors and matrices, just % repeated on each element: % sin(rowvector) % exp(matrix) % log(5-matrix) %% HEY!: where do the complex numbers come from in that matrix output? % % Matlab has a function to create an identity matrix of whatever size: % I=eye(4) % And two functions to create matrices of ones or zeros: % M=ones(5,3) % N=zeros(2,4) % A function to create a vector or matrix of uniformly-dist random numbers: % P=rand(1000,1); % this is a col vector, row vector would have been (1,50) % And a vector or matrix of normally-dist (Gaussian) random numbers: % Q=randn(1000,1); % this is a col vector, row vector would have been (1,50) % We'll compare the values in P and Q in a few moments when we learn how % plot a histogram. Note the semicolons there as these vectors are so % long, in order to make sense of their histograms. Ok, back to this soon. % %% HEY!: now to remind yourself of some linear algebra identities: % let's load up the "good" matrix again (the one we didn't mess up): % A=matrix; % % samematrix=eye(3)*A % the following two results should be the same: % Ainv' % inv(A') % the following two results should be the same: % Ainv*A % is the result exactly an identity matrix? or at least close? % A*Ainv % why might it be a teeny tiny bit off? % Let's bring in another matrix: % B=[2,8,11; 1,7,12; 6,10,4] % the following two results should be the same: % (A*B)' % B'*A' % the following two results should be the same: % inv(A*B) % inv(B)*inv(A) % HEY!: (does the order matter on those multiplications?) % PLOTTING - FINALLY! % See, we had to learn how to create and manipulate vectors and matrices % before we could make sense of this section... % ------------------------------------------------------------ % (You can comment out everything in earlier sections now, except for the % P & Q random vectors - we'll finally plot those histograms.) % % myXvector=1:50; % myYvector=sin(.2*myXvector); % plot(myXvector,myYvector); % this will open a figure window % grid on; % put a grid on the plot % let's open another figure window and show how to fancy up a plot: % figure; % plot(myXvector,myYvector,'r--',... % myXvector,myYvector.*myXvector/25,'ks',... % 'LineWidth',2,... % 'MarkerEdgeColor','k',... % 'MarkerFaceColor','g',... % 'MarkerSize',10); % title('This is my Plot','FontWeight','bold','FontSize',12); % xlabel('FooFoos'); % ylabel('BooBoos'); % Check "help plot" for some info on plot options, and if you want to % really get into the rest, learn about "handle graphics". You can check % and set every option for a plot via its "handle" which can be optionally % returned like this: % h=plot(myXvector,myYvector); % set(h) % tells you all the options available for that plot % get(h) % tells you the current values of all properties of that plot % set(h,'LineWidth',1,'Marker','o'), % changes these two properties % % And here's a handy way to plot matrix contents: % figure; % C=randn(30); % 30x30 matrix of Gaussian values % imagesc(C); % plot with colors scaled between max & min elements % colorbar; % add a colorbar legend to the plot % a handy option for this type of plot, to move the origin to the bottom % left (default is top left). % axis xy % A quick digression, simply to make a cool-looking plot(!). One way to % smooth an image or surface like that random matrix is to filter it by % convolving it with a smoothing kernel. Matlab has a built-in convolution % function: % x=linspace(-1,1,15)'; % 15 points from -1 to 1 for the bell-curve below: % g=exp(-x.^2/2/(.4)^2)/sqrt(2*pi)/.4; % points on a Gaussian bell-curve % kernel=g*g'; % HEY!: what is the dimension of kernel here: matrix? scalar? % D=conv2(C,kernel); % imagesc(D); % colorbar; % HEY!: try changing the number of points in the bell-curve above - ie the % size of the smoothing kernel, and look at the effect in the rerun plot. % How's it look with 5 points, 10 points, 30 points? % % Ok, let's finally look at those histograms for P & Q, the uniformly- % distributed and normally-distributed (respecitively) random vectors we % created earlier. (Are they still uncommented?) % Matlab has a handy histogram plotting function that does it all for you: % hist(P,20); % specifying 20 bins to histogram over. default is 10. % figure; % start a 2nd plot % hist(Q,20); % specifying 20 bins to histogram over. default is 10. % HEY!: okay, see the difference between being uniformly distributed % between 0 & 1, and being normally distributed about 0 with stdev 1? % Try this: shrink and expand the 20 bins to say 4 and to 200. How's that % affect the histogram plots, and your interpretation of the distributions? % % % % % % % % % YOUR HOMEWORK FOR THIS LAB: % ------------------------------- % Use the example bits in this file to piece together a script of your % own which demonstrates the gravity problem discussed in class and the % book, where the ball's elevation over time is a parabola. (This is % example 1.1 in the Aster book.) % Do as above in this file, where we made up some x ("m" in the book, ie % the initial position, initial velocity, and acceleration of gravity). % Create a vector of times t_i, say maybe 1:.1:10 (or whatever range). % Create an A matrix ("G" in the book) with the 3 columns based on the t_i. % Calculate the trajectory y_i of elevations over times t_i, a parabola. % Now add some Gaussian noise to y using the randn() function (play with % how big its standard deviation is). % Use the normal equations to find the LS soln "mhat" for the initial pos, % initial vel, and gravitational accel. % Calculate the trajectory yhat of the best fit to the data, from mhat. % % Print out and turn in a plot showing the noisy data y and the best fit % yhat to the data, on top of each other. And show your original m and % your mhat as well for comparison. % % Don't hesitate to holler if any of this is not clear or if you're having % Matlab trouble... % % Cheers! % -Andy