Colorado State University, Pueblo; Spring 2014
Math 242 — MATLAB Computation
Course Schedule & Homework Assignments

Here is a link back to the course syllabus/policy page.

In the following, the image means that class was videoed that day and you can get a link by e-mailing me. Note that if you know ahead of time that you will miss a class, you should tell me and I will be sure to video that day.

This schedule is will be changing very frequently, please check it at least every class day, and before starting work on any assignment (in case the content of the assignment has changed).

• M:
  • installing Linux Mint on thumbdrives distributed to the class
  • discussion of what is an operating system, and what is open-source software
• T:
  • Read:
    1. Optional but fun: here is an article written by the science fiction author Neal Stephensen about the history of operating systems. It is a tiny bit dated, having been written 15 years ago, but quite amusing and well-written. (This version taken from Stephensen's web site http://www.cryptonomicon.com/beginning.html.)
    2. Here is an entire site devoted to teaching people how to use "the command line", that way of doing things in Linux where you type commands rather than just pointing and clicking. Read at least the pages numbered 1 to 6 in the list near the bottom of the page http://linuxcommand.org/lc3_learning_the_shell.php.
  • booting from the Linux Mint thumbdrive
  • starting syaptic and installing new programs: in particular, GNU Octave
  • opening a terminal window, basic commands there:
  • cd directoryname [change to a given directory]
  • cd [go to your home directory]
  • cd .. [go to a direction up step up the hierarchy]
  • ls [list the contents of the current directory]
  • rmdir [remove a given directory]
  • pwd [print working (current) directory]
  • starting Octave in a terminal window; basics:
  • arithmetic: just type it. be careful of operator precedence (or use lots of parentheses)
  • it's an interpreted language so there is no need to write a program and compile it... instead, it compiles and executes everything immediately
  • variables can be set -- carefule of the "=" sign, it represents "set this variable to have that value", not the usual thing with "=" in mathematics
  • legal variable names
  • values of variables (or other expressions) can be integer or floating point ... or matrices!
  • expression [$a_{11}$ $a_{12}$; $a_{21}$ $a_{22}$] to define a $2\times 2$ matrix
  • Octave knows to use matrix multiplication when you type $a*b$, if $a$ and $b$ are matrices
  • Do HW0: Send me e-mail (to jonathan.poritz@gmail.com) telling me:
    1. Your name.
    2. Your e-mail address. (Please give me one that you actually check fairly frequently, since I may use it to contact you during the term.)
    3. Your year/program/major at CSUP.
    4. The reason you are taking this course.
    5. What you intend to do after CSUP, in so far as you have an idea.
    6. Past math classes you've had.
    7. Other math and science classes you are taking this term, and others you intend to take in coming terms.
    8. Your favorite mathematical subject.
    9. Your favorite mathematical result/theorem/technique/example/problem.
    10. Anything else you think I should know (disabilities, employment or other things that take a lot of time, etc.)
    17. [Optional:] If you could meet one historical figure, who would you choose and why?
• W:
  • Read:
    1. Here is a useful quick reference ["cheatsheet"] for many basic Linux commands
    2. This is an entire book entitled Octave Programming Tutorial, which we will use a fair bit of. Please read the "getting started" and "vectors and matrices" sections thereof.
  • a few more useful commands in a terminal window:
  • rm [remove a file]
  • cat [display the contents of a file on the screen]
  • less [display the contents of a long file on the screen, one screenful at a time: press space to get to the next screenful]
  • gimp [huge, complicated wonderful program which can display -- and modify in many ways -- images ... it's a free equivalent of Photoshop, only better]
  • students worked in Octave through the two sections described in today's Read #2, above, for which it was necessary first to:
  • install with the Package Manager (synaptic) the package gnuplot-x11
  • in Octave, use the diary command to take a transcript of your commands and its responses; see this manual page
  • diary won't save the plots, though ... you can use the print command, though -- see this page. For example, "print -dpng <filename.png>" will save the current plot in the image format PNG.
• $\Theta$:
  • here is the Octave diary from today's class
  • clc
  • who
  • whos
  • format short
  • format long
  • bank
  • format +
  • format bit
  • format native-bit
  • format hex
  • length(x)
  • k:l:m
  • rand(n)
  • A(1,1:5)
  • A(2:4,2:4)
  • roots([1,0,0,1])
  • A\v
  • det(A)
  • B=diag([1 2 3]))
  • eye(n)
  • diag(A)
• F:
  • Since we did a long class on Tuesday, we will not have a Friday class this week.
  • Today [Friday] is the last day to add classes.

• Today [Monday] is the last day to drop classes without a grade being recorded.
• T:
  • some Linux customization, UI: background, key repeat rates, screen layout, focus on mouseover
  • shell: choosing, using, defining shortcuts
  • emacs: start-up file (see model here)
  • error handling: error (with and without "\n")
  • argument sanity checking: how to tell if a number is an integer? (one way is to use if n==floor(n) ... )
  • printf
  • return
• W:
  • Today we play Sudoku!
  • We will represent a sudoku board by a $9\times 9$ matrix whose entries can be any integer from the set $\{0, 1, 2, \dots, 9\}$, where $0$ means "has not been filled yet", and the other values are themselves.
  • Write a function showparts which takes as input a sudoku board and prints out the 9 columms, 9 rows, and 9 $3\times3$ submatrices which must be examined to see if it is a valid [partial] Sudoku solution.
  • Now write a function isvalidpart which takes as input a matrix of size $1\times9$, $9\times1$, or $3\times3$ and return false if any number $\{1,\dots,9\}$ appears more than once as an entry of the input, or returns true otherwise. Here is one way you might do that:
    • In isvalidpart, make a vector seen of length 9 which starts out filled with zeros (try the zeros command). Think of this vector as Booleans whose ith entry contains the value of the condition "have we seen the number i in the input yet".
    • Now loop through the input and look at each element, call it i:
    • If i is 0, skip over it (you might want to use the continue command).
    • If not, check that i is a whole number in the range 1 to 9 (inclusive), otherwise throw an error.
    • If i is such a whole number, look at the value seen(i).
    • If this value is true, return from isvalidpart immediately with return value false.
    • Otherwise, keep looping.
    • If the loop finishes, return with value true.
  • Now make a version of your showparts function which also prints out if the parts are valid or not.
  • Make a function isvalid which takes an entire sudoku board as input and tells if it is valid or not, returning true or false
  • Go back to isvalidpart and make a version which returns more information: if the part is valid, it returns a vector of all the values which were not used in that part; if invalid, it returns a value which occurred more than once in that part.
  • Using the new isvalidpart, can you make a function which prints a version of the board which in each location prints out the value in that location, if one is defined, and prints a little matrix of the values which would be valid in that location?
  • Can you use the above pieces to make a program which does a "brute-force" solution of a given, partially filled-in Sudoku board?
• $\Theta$:
  • Remember: many Octave functions and operations that usually work on single values (scalars) can also be applied to vectors (or matrices), to give a vector (or matrix) of answers. These include conditionals
    • vec==17 makes a vector of the same size as vec where each component is a 0 if that component of vec is not 17, and a 1 if it is.
    • vec(i j k) makes a vector of values [vec(i) vec(j) vec(k)]
    • ~vec makes a vector which is 0 where vec is nonzero, and 1 where vec is zero.
  • The command find(vec) which finds the locations of vec which are non-zero
  • Therefore vec(find(vec)) makes a vector of the non-zero components (in order) out of vec
  • The command unique(vec) which makes a vector of the elements of vec sorted in order and with repeats removed
  • Ff you're done with the assignments from yesterday, make a fancier version of isvalidpart which returns, for valid parts, a list of all the values which could be filled in for whichever 0's may be in that board part and still make it a valid sudoku board
• F:
  • Cleaning up some of the sudoku code we'be been writing.
  • To finish the complete, combined isvalid and have it return a single true or false: boolean operations &, |, &&, and || — the difference between short-circuiting boolean operations (&& and ||) and the non-short-circuiting ones & and |.

• T:
  • A little about variables in Octave being created on demand and of the right type, dynamically
  • Creating a variable as a matrix element implicitly builds the variable as a matrix of the smallest rectangular size which would contain that element: a(2,3)=17 for a previously undefined variable a would automatically create a matrix a of size $2\times3$ whose $(2,3)$ entry was 17 and all other entries were 0s.
  • Doing the same thing with strings builds a rectangular matrix padded with spaces: so b=["a";"string"] makes b to be a $2\times6$ matrix of characters whose first row has the character "a" followed by 5 spaces, and whose second row is the vector of characters which forms the string "string"
  • We've seen sort of a vector: it just puts the entries in increasing order. Also works on a string, putting the characters in "increasing" order (first capitals in alphabetic order, then lower case letters; some symbols mixed in, in a strange order).
  • sort on a vector sorts each column separately in increasing order. This is quite weird for matrices of strings, as the above example b
  • Starting towards a program which will put a matrix of strings into alphabetic order, let's start with putting a column vector of numbers into increasing order, as follows:
  • Given a (numerical) column vector v, loop through the elements of the vector, with index i. we'll assume that the vector from location 1 to i-1 is already sorted, so the idea is to see where element v(i) goes.
  • Loop through the entries 1 to i-1 with a loop index j. if v(j)<v(i), do nothing in this j loop.
  • If v(j)>v(i), we need to insert v(i) at that location. do this by saving the value of v(i), say in a variable x. move the part of v in locations j to i-1 up one step, then put x into location j of v
  • To "move a part of a vector up one step": do something like v((m+1):(n+1))=v(m:n)
  • Write this code, test it, and submit it by e-mail
• W:
  • The end goal of this week we're going to simulate Pachinko to see the connection with a famous distribution in statistics — see this video for a real-world version of this.
  • but first, we need two tools: plotting histograms and generating random numbers
  • rand
  • nargin
  • Now, students please make a function coin which optionally takes an argument p (given value of .5 if no argument is given) and returns true with probability p, otherwise returns false
  • Test the theoretical reliability of your function in as many ways as you can think of: make graphs, compute averages, etc.
  • Now write code which takes an input n and flips the coin n times, then makes a list of the lengths of runs of heads (well, of trues). make a histogram of this result, using hist
  • Let's make our own hist program!
  • Write a function counts which takes a data vector v and a vector bins which mark of the bin boundaries for a histogram, and produces a vector of counts of the data in the bins.
  • Use bar to plot the data from counts, so making your own version of hist
• $\Theta$:
  • In applied mathematics, techniques which use randomness where there is none inherent in the original problem are called Monte Carlo techniques. Today we'll try one such method, as practice using rand
  • Let's compute $\pi$. A very direct way to do so is to compute the area of a unit circle. Write a program pibyRS ["pi by Riemann Sums"] which takes an input n and divides the interval from -1 to 1 into n subintervals, picks a point $x_j$ in each subinterval (left endpoint, right endpoint, midpoint, whichever you like), and finds the area of the unit circle by computing $$\int_{-1}^1 \sqrt{1-x^2}\,dx \approx \sum_{j=1}^n (2/n)\cdot2\sqrt{1-x_j^2}$$
  • Try running your pibyRS program several times to see how much accuracy you get by different values of n: make a graph of the approximate $\pi$ values versus n.
  • Now do this with a Monte Carlo approach: write a program pibyMC which takes an input n, picks n random points in the unit square in the $x-y$ place, counts them as inside the unit circle if $x^2+y^2<1$, and reports its guess of $\pi$ as $$\pi\approx4\cdot\left(\frac{\text{number of points inside unit circle}}{\text{total number of random points in unit square}}\right)$$
  • Likewise make a graph of approximate $\pi$ values agains n, now for pibyMC
• F:
  • Now, really do Pachinko. So write a program which takes two inputs n, the number of steel balls to be dropped into the Pachinko machine, and r, the number of rows it will have. Make a for loop over the n balls, and for each of them, loop over the r row. The ball should start out at location x=0 and as it goes past each row, it should either bounce left or right, so either do x-- or x++, based on the output of a fair coin flip.
  • At the last row, figure out the totals of the numbers of balls in each bin, and then make a histogram of them.
  • A variation: try doing the above with an unfair coin. Or with a complex obstacle on each row, which leaves the ball alone with probability $p_0$, moves it right one step with probability $p_1$, left one step with probability $p_{-1}$, right two steps with probability $p_2$, left with probability $p_{-2}$, etc., where the probabilities are chosen so that $\sum_{i=-N}^N p_i = 1$.
  • If you finish all of the above, make a graphical version of the Pachinko program! Use the commands sleep, set(a,'visibile','off') (where a is a return value from a plot command), and plot(x,y,'o')

• T:
  • A basic introduction to polynomial fitting through data:
  • Get the data in this file.
  • This data comes from Galileo Galilei's experiments in rolling balls down inclined planes. An explanation can be found here.
  • Now use polyfit to fit a best quadratic polynomial to the data (it's in two data sets).
  • Plot the data and the fitted curve on the same axes.
• W:
  • Here is a worksheet to follow to get some basic practice with polynomial fitting to data. The data files can be found at this page and this page.
  • Next we talked in general about fitting polynomials of degree $n$ to datasets with $n+1$ points, as follows:
  • Take the coordinates of the data to be $(x_1,y_1),\dots,(x_{n_+1},y_{n+1})$. Make an Octave vector y=[$y_1;\dots;y_{n+1}$]. Think of the vector a as representing the coefficients of the degree-$n+1$ polynomial we are looking for. Put the $x$ values together into an $(n+1)\times(n+1)$ matrix V (this is called the Vandermonde matrix) as follows: $$V=\begin{pmatrix}x_1^n&x^{n-1}_1&\dots&x_1&1\\x_2^n&x^{n-1}_2&\dots&x_2&1\\\vdots&\vdots&\ddots&\vdots&\vdots\\x_{n+1}^n&x^{n-1}_{n+1}&\dots&x_{n+1}&1\end{pmatrix}$$ Now the matrix equation V*a=y represents exactly the fact that the polynomial function with coefficients a has y values $y_1,\dots,y_{n+1}$ at the $x$ coordinates $x_1,\dots,x_{n+1}$ (think about it!).
  • Fact: the Vandemonde matrix is always invertible if the $x$ coordinates of the data points are all distinct.
  • We can thus solve for the polynomial coefficients by a=V^-1y, or (and this is more efficient and stable) using the Octave left division: a=V\y.
  • If we have more data points — say $\ell$ — than the degree of the polynomial we are seeking — still $n$ — then the system $V \vec{a}=\vec{y}$, which is an $\ell\times n$-matrix times an $n\times1$ vector equalling an $\ell\times1$ vector,is overdetermined. It may (most likely does) not have a solution. The closest thing to a solution would be a vector $\vec{a}$ which minimizes the norm-squared of the distance from $V\vec{a}$ to $\vec{y}$, $||V\vec{a}-\vec{y}||^2$ &mdash this is the least squares solution.
  • We find the least squares solution by solving the normal equations, which are $(V^T V)\vec{a}=V^T\vec{y}$. [Note $V^T V$ is $n\times n$ and $\vec{a}$ is $n\times1$, as it $V^T\vec{y}$.] We can solve the normal equations as $\vec{a}=(V^T V)^{-1} V^T\vec{y}$ or, what is better in Octave, a=(V'*V)\(V'*y).
  • Write some code! For the first part, write a function which takes two input vectors x and y, assumed to be of the same length, representing the $x$ and $y$ coordinates of the data points you are going to study. The function should return the coefficients of the polynomial of degree length(x)-1 which fits that data, and graph the data points and the polynomial on the same plot (which plot should have $x$ range min(x) to max(x) and $y$ range min(y) to max(y).
  • For the second, continue to accept the same x and y inputs, but now have another, n, specifying the degree of the polynomial to fit to those data points. Use the normal equations method to do this. Again, the function should return the coefficients of the polynomial of degree n which fits that data, and graph the data points and the polynomial on the same plot which plot should have $x$ range min(x) to max(x) and $y$ range min(y) to max(y).
  • Please comment you code, particularly make sure the comments which come up when you type help functionname are clear and completely describe what kind of inputs your functions take, any assumptions you are making about those inputs, and what kind of outputs your function returns.
• $\Theta$:
  • Continue the above project. Don't forget those comments!
  • Here is a nice source of important data, if you have time please download some of the files and put nice curves through them.
• F:
  • Continue; comments. If you are finished, feel free to go back and finish some previous week's work which might need a little final polishing. Or else work on the extra credit project which follows.
  • Extra credit: Do as much as you have time and interest in:
    1. Make a version of your previous program, the one which takes the degree of the polynomial to use as the approximating function, which calculates and prints a measure of how good the approximation is. The best measure of this type would be the sum of the squares of the residuals. That is, if the data values are $(x_1,y_1),\dots,(x_n,y_n)$ and the function (at first a polynomial of some degree $n$, although we will change that in a moment) your code has found as the best approximator is $f(x)$, the error measure would be $$\sum_{i=1}^n \left(y_i-f(x_i)\right)^2\ .$$ Your program should therefore do all of what is described above and also print "Least squares error is: yadda-yadda."
    2. Make a version of your curve fitting program which does a polynomial fit of $y$ against $x$, but also one of $y$ against $\log(x)$, $\log(y)$ against $x$, and $\log(y)$ against $\log(x)$. For each, calculate the least squares error, and pick the one with the smallest such error to report and to graph.

• T:
  • This week's project will also count as our midterm! So you may not search aimlessly on the 'net for code which might be useful -- although if you need to research some specific Octave feature, that is fine: feel free to read online documentation on some function like hist.
  • We shall work on the Caesar cipher. In the original version, the origial message (the cleartext or plaintex) was written without spaces or punctuation and all in lowercase. To make the coded version (to encrypt the plaintext, getting the ciphertext), every letter was moved three letters further down the alphabet, wrapping around to A if necessary.
  • In our version, the magic number 3 will be replaced by any number the user chooses — the particular value is called the key.
  • If we are given a message which was encrypted with an unknown key, but we want to recover the orignal message (so, we want to get the plaintext from a given ciphertext with unknown key), one approach is to try all possible decryptions with different keys (there are only 26 possible!) and pick the one which gives something that looks like English text.
  • Our goal this week is write an entirely automatic Caesar ciper-breaking program. We can still still try all possible keys, but we want our program to pick out which key is most likely to be giving the correct decryption ... without human intervention.
  • One way to do that might be to look for the most frequent letter in the ciphertext we were given and to assume that that was the encryption of the letter E, which is the most common letter in English text. However, some odd bit of English text, particularly if it is fairly short, might not have E as the most frequent letter, so this method of decryption will fail.
  • Instead, we will compare the entire histogram of frequencies of the ciphertext letters to the histogram of frequencies of English letters. We will try all possible (all 26) shifts of the ciphertext histograms and choose which is the "most similar" to the standard English frequencies. That will be the one which we will report as the most likely shift (Caesar cipher key value).
  • The last thing we need to agree upon is how to measure how similar are two histograms, or tables with frequencies of all 26 letters. We shall use what is called the $L^2$ distance: given functions $f$ and $g$ defined on the same domain, being the whole numbers from $0$ to $25$, we shall say the $L^2$ distance between them is the number $$||f-g||^2_{L^2}=\sum_{j=0}^{25} (f(j)-g(j))^2$$ (actually, this is the square of the $L^2$ distance, but it is easier to minimize the square than to have an extra square root to take all the time and therefor to be working with the real $L^2$ distance; and it gets the same answer).
  • here are some Octave fragments we tried today which were useful:
    • To take a string m and make it Caesar encryption ready, it was good to do m=lower(m(isletter(m))).
    • Characters in Octave are stored internally as integers taking up one byte. The correspondence between characters and integers is by the famous ASCII (American Standard Code for Information Interchange, dating from the early 1960's; see this page for more information about ASCII) code. All lowercase letters are in a sequential block starting, however, with 'a' not having the value 0. In particular, to transform a string of lowercase letters to an integer vector with values 0 to 25 (corresponding to the letters 'a' to 'z'), it suffices to do str-'a'.
    • Another way to convert a string to a vector of integers str+0; or, more directly by int32(m); these methods do not force A to correspond to 0.
    • A vector of integers can be converted to the corresponding characters (so, to a string) by str=char(vec). (When doing this, it is important that the integers are the ASCII codes.)
    • We can Caesar encrypt with key $3$ by the command e=char(mod((m-'a'+3),26)+'a'). (Think through why this works!)
• W:
  • Start by writing a program which takes a string and a key value and returns the Caeser-type encryption of that string with that key. You need to make sure you remove all non-letters from the string and convert it to lower case. Test your program. Exchange an encrypted message with a classmate and decrypt -- ask your classmate politely for the key.
  • Note you can make a special decryption program which would simply subtract the key (rather than add it) to the letters of the ciphertext. Or you can use exactly the encryption program with the decryption key which is 26 minus the encryption key.
  • Now work on the automatic decryption program: it should take as input a string and try all possible decryption keys on it. For each, it should compute for each ciphertext letter what fraction of the whole ciphertet was comprised of that latter. Then it should compute the $L^2$ distance between that frequency table and the standard English one -- which can be found here. It should return the string whose frequency table has smallest $L^2$ distance to the model English letter frequency table.
  • The Octave command readline can be used to read content from a file; try help readline. So can be textscan. fprintf can write to files.
• $\Theta$:
  • If you are ambitious, do as much of the above as you can for the Vigenère cipher, which is like a whole bunch of Caeser's done next to each other. That is, one takes as a key a whole word, which is repeated as many times as necessary to make as many letters as are in the message. Then one adds the two strings together (mod $26$, of course) to get the ciphertext. So you could write some code which does Vigenère encryption and decryption.
  • To do automatic decrypting of a Vigenère cipher text, one breaks the text into as many pieces as there are letters in the key and does a separate attach on it like the one we did above for Caesar. The difficulty is what to do if we don't know how many letters are in the key (answer: try a bunch of different possibilities).
  • As a historical note: Vigenère was considered unbreakable for a least a couple hundred years.
  • If you want to make your ciphertexts and cleartexts prettier, you can do the following: with a string s, save where it has letters or non-letters by where_lets=isletter(s), where_nonlets=!isletter(s). Then do t=s(isletter(s)), transform t in some way (encryption, probably), then at the end do s(where_lets)=t to put the transformed letter back mixed in among the original non-letters.
• F:
  • Finish your code — your midterm! Some things to keep in mind:
    1. Your code should run correctly -- make sure you check it completely!
    2. Make sure your code sanity checks all its input.
    3. Make sure you hand in (by e-mail) all the necessary programs (and other files, if necessary) so that your project runs, and which enable it to be tested -- e.g., include a Vigère encryption function if you want your Vigère decryption function to be tested.
    4. Make extensive comments in your code. In particular, every major control structure (sub-function, loop, condition, etc.) should be explained, and every variable should be described. Really. Always. (You don't have to explain somthing twice: so if you have a loop for i=1:n, you can explain what the loop is iterating through and then you don't need to explain the variable i.

• T:
  • a useful shell (general Linux) command: wget URL, downloads the page at that URL into your current local directory.
  • Some useful Octave commands relating to matrix operations:
    1. eig
    2. rref
    3. null
    4. norm, such as norm(A,"fro")
• W:
  • We're doing a worksheet on finding eigenvalues, which can be found here.
• $\Theta$:
  • Here's a supplemental project you can do to make an optional, extra credit part of this week's project..
  • This supplemental part is to make a program which creates a graphical tool to find eigenvectors of $2\times2$ matrices.
  • Write a function which takes two inputs, a $2\times2$ matrix A and a positive integer n. First, you program should draw n equally spaces line segments from the origin to points on the unit circle.
  • Now, think of each radial line segment as a (unit!) vector $\vec{v}$, compute the vector $\vec{w}=A\vec{v}$, then normalize it to a vector $\vec{u}=\vec{w}/||\vec{w}||$. Draw another line segment which is a vector from $\vec{v}$ to $\vec{v}+\vec{u}$. (Again, this will look nicer if you draw
  • Actually, it looks a bit nicer if you
    1. draw only half of the radial line segments, the half which ends on the unit circle;
    2. draw only half of the linear transformed vectors, again the half which touches the unit circle;
    3. make the two sets of vectors you draw have different colors (by e.g., plot([v(1) w(1)],[v(2) w(2)],'g'); to get green, 'm' to get magenta, etc.)
    4. make the plot have a square aspect ratio, (by e.g., axis([-2 2 -2 2],"square");)
  • This process should make a picture something like: The idea of such images is that a vector $\vec{v}$ is an eigenvector of $A$ if $A\vec{v}$ is parallel to $\vec{v}$ — it will look like it is in the same (or the opposite) direction — the two segments will line up!
• F:
  • Today [Friday] is the last day to withdraw (with a W) from classes

• T:
  • talked about some ideas of basics of Fourier Transforms as a kind of appoximation scheme for periodic functions which another famous approximation scheme, Taylor Series, does not handle well.
  • Some useful relevant Octave:
    1. fft — [Fast] Fourier Transform
    2. ifft — Inverse [Fast] Fourier Transform
    3. stem — plotting which doesn't attempt to connect the dots into a continuous function
  • Mentioned relationship to AM and FM broadcasting, smoothing noisy signals, etc.
• W:
  • We're doing a worksheet on Fourier (and cosine) transformations, which can be found here.
• $\Theta$:
  • Continuing the weekly project. A few notes:
  • There are few small errors in the worksheet. This is (in part) intentional: the hope is that you will work them through and it will be a useful part of the exercise.
  • For example, there is a place where imagesc will complain that it cannot display image matrices of complex numbers. So perhaps you have to extract a real number corresponding to a complex one (coming from an ifft — as we have said, fft, and therefore ifft are complex-valued). There are several ways to do this: real, imag, and abs all work. Perhaps you should compare these three — can you say why one works very poorly, and the other two are fine?
  • Another example: the code in the last section of the worksheet does not work exactly as it is written. In particular, there is an issue with the use of end -- the worksheet is based on some code which was taken from a book that was designed for Matlab, which apparently has a different approach to the keyword end. In your code, however, you can figure out what value to use in place of end by looking at the length of one of the vectors involved.
  • A final problem may occur if you use your own image file which is in jpg format, but is a specialized grayscale jpg. This will work fine at first, it will be read in in the usual way... but you will not have to sum the R, G, and B values, as there is only one value, the grayscale. Worse, the values will be uint8's, which do not work well with other numbers. If this happens to you, just conver the input image matrix into a matrix of doubles with the double command.
• F:
  • Please finish up as much of the week's main project as you can. When it is done, send me the parts as specified at the beginning of the worksheet itself. If you are done, here is an extra credit project you can tackle for this last class of the week:
  • Let's use the fft to do smoothing and to seek underlying patters in time-varying signals, unlike the signals which varied in space (actually in two spacial dimensions) which we did in the worksheet, the images.
  • So start by getting some interesting time series data. I leave this to you to find -- try to get around a hundred numbers spaced regularly with time. You can use stock market data every day for a year, or temperature data (e.g., here is some historical temperature data for Pueblo, CO from this NASA website, your favorite other data.
  • Your data should be equally spaced along the time (we'll use the $x$) axis; if not, make up some values.
  • Take the (one-dimensional) fft of this data.
  • Use your polynomial interpolation from earlier this term and subtract from each data point the best linear approximation to that data (the least squares regression line), because this will make a new data set whose mean $y$ is 0 and which has no overall trend. But keep that regression line around to use later!
  • Compute the fft. Plot it for fun.
  • Now filter off the higher frequencies by multiplying by a nice function for this purpose, called a Butterworth function of order $n$: $$h(\omega)=\sqrt{\frac{1}{1-(\omega/\omega_c)^{2n}}}$$ This function is 1 for $\omega=0$ and goes smoothly to 0 for large $\omega$. If we multiply the Fourier transform of our data, it will smoothly ramp down the high frequencies to zero, which is what we want to do in order to make the curve smoother.
    The $\omega_c$ is called the critical frequency and tells where the function is decaying fastest.
  • Graph some different Butterworth functions (often we use $n=5$ or $10$ or something of that size). In Octave, the function might be

                          w = sqrt(1./(1 + (freq./wc).^(2*b_order)));
    

  • Multiply a well-chosen Bullterworth filter with your data's fft, take the inverse tranform, and plot both the original data and the smoothed version on the same axes.

• Spring Break! No classes, of course, but ...
• ... please do start thinking about what you would like to do for your final project. Here is a page describing what a project should be, and giving some suggestions that might help you get started with finding a topic.

• T:
  • By this date (but much earlier — e.g., during Spring Break! — is better!) you must have communicated the general area of your final project to your instructor.
  • This week we're going to work on more elaborate plotting, particularly of data which is more than two-dimensional.
    • mseshgrid
    • mesh
    • meshc
    • meshz
    • hidden
    • surf
    • surfc
    • surfl
    • contour
    • contour3
    • quiver and quiver3
    • compass
    • feather
    • rose
    • stairs
    • stem
    • stem3
    • scatter
    • plotmatrix
    • stem3
    • errorbar
    • polar
    • pie
    • pie3
    • pcolor
    • comet
• W:
  • Start this worksheet.
• $\Theta$:
  • By this date you must have emailed to your instructor the specific topic of your final project, described in about a paragraph of text.
  • After you finish the parametric worksheet, please do the following problem:
    For a Calc 3 class, the instructor wants to do as examples these triple integrals: $$\int_0^2\int_{-2}^1\int_0^3 2x^2+y^2z\,dxdydz\qquad\text{and}\qquad\int_0^1\int_0^{1-x^2}\int_0^{1-x^2-y} x\,dzdydx\ \ .$$ Help out this instructor by making 3-D plots of the regions in space over which the integration will take place.
    Notes:
    • When plotting a region in space, the way to do it is to plot its surface, maybe (usually) broken into several pieces which you can plot one after the other with hold on to get then all together.
    • If you are unfamiliar with this notation, here is an explanation: $$\int_a^b\int_{p(x)}^{q(x)}\int_{r(x,y)}^{s(x,y)} f(x,y,z)\,dzdydx$$ means to integrate the function $f(x,y,z)$ over the region described as: $$\begin{eqnarray*}&\text{$x$ ranges over the interval }[a,b]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\\text{for a particular }x,\ \ \ &\text{$y$ ranges over the interval }[p(x),q(x)]\ \ \ \ \ \ \ \\\text{for a particular }(x,y)\ \ \ &\text{$z$ ranges over the interval }[r(x,y),s(x,y)]\end{eqnarray*}$$
    • Be careful of the order of the variables (and hence of the $dx$, $dy$, and $dz$) in the two above integrals.
    • Please label your axes for clarity.
• F:

  This weeks work should consist of your graphs and code for the parametric surfaces and either the Thursday extra problem on pictures of domains for two particular triple integrals, or the following problem: Write an Octave program called deathmask which takes a single string input specifying a filename. It then reads in that file — it will be an image file — processess it the way we did on the FFT worksheet (combining the RGB values to make it a single value of intensity), and then does a surf or mesh plot of this data (maybe both? could be another input to your function chooses which type to make). Try running it on the Mona Lisa image, or at least the part which contains her face, to make something like an ancient deathmask ... should be something like the famous deathmask of Agamemnon, as the image to the right. (Try to get your colormap to make this nice effect of beaten gold.) Submit to your instructor the code/script/diary and images to all parts of this week's work.

• T: This class has been taken over by the following special event:
  Off to Infinity in Finite Time
  In 1994, Xia solved a 100-year-old problem on the existence of singularities in the N-body problem, i.e. the solution of Newton's equations of motion for N "planets" acting under the mutual influence gravitation via Newton's inverse square law. Graham Harper will give a presentation based on Xia and Saari's 1995 article in the Notices and Xia's 1994 thesis that appeared in the Annals of Math. Join us for some interesting dynamics.
• W, $\Theta$, and F: Finishing up final projects. Please be in frequent contact with your instructor, so there are no surprises about how your project is moving forward (or — hopefully not — getting stuck). We may we have some presentations during class time this week, for those students who are ready. Here are the final project presentations which are scheduled so far:
  1. $\Theta$:
     • Jamen
  2. F:
     • Correy
     • Graham
     • Joel
• There are some notes on the project description web page about what is expected for your report and presentation.