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:***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.)- 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*:- Your name.
- 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.)
- Your year/program/major at CSUP.
- The reason you are taking this course.
- What you intend to do after CSUP, in so far as you have an idea.
- Past math classes you've had.
- Other math and science classes you are taking this term, and others you intend to take in coming terms.
- Your favorite mathematical subject.
- Your favorite mathematical result/theorem/technique/example/problem.
- Anything else you think I should know (disabilities, employment
or other things that take a lot of time,
*etc.*) - [Optional:] If you could meet one historical figure, who would you choose and why?

*W:***Read:**- 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)**

- here is the
*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.**

*T:*- installed
**emacs**— which is an*editor*, unlike the*word processors***Microsoft Word**and**LibreOffice Writer**, makes straight files of characters, without lots of extra cruft that blows the file up - the linux commend
**ls -l**, file size, owner, group of owner, RWX permissions. - some new
*Octave:***save****load***functions*— local variables in a function*do not*persist in the*Octave*workspace*scripts*— local variables in a scipt*do*persist in the*Octave*workspace**edit**

- installed
*W:*- some new
*Octave:* **deal**... switching variables**disp***strings*are just row vectors;*concatenating strings**scripts:*names, calling**clear**(companion to*save*and*load*)**save**for specific variables**clock**- don't forget
**help** - elementwise operations on vectors....
**a.*b**and**a.^b** **num2str**- Exercises: please do these; for each, if it is a script, send
me the script; in any case, make a
*diary*of any work in*Octave*and send me that as well.- first, do the Emacs tutorial. make a few sample files
- make a script
`helloWorld`which displays the text "Hello, World!", then, on another line, the text "I just started learning Octave!" - current date and time:
- create a variable
`start`with the function**clock** - what is its size? is it a row or column?
- what does it contain? use
`help clock` - convert
`start`to a string using**datestr** - save
`start`and your string variable equivalent in a file`startTime` - modify
`helloWorld`so that it reads in those variables with the date and time and uses them as the actual times in the declaration of when you started learnign*Octave*

- create a variable
- learning programming is exponential!
- your learning time constant is 1.5 days. calculate the
number of seconds in 1.5 days and put it in a variable
`tau` - a week lasts 5 days. calculate the number of seconds in
5 days and put in a variable
`endOfWeek` - your knowledge as a function of $t$ is given by $k=1-e^{-t/\tau}$
- how well will you know
*Octave*after a week? Name this`knowAtEnd`display the phrase "after a week, I will know *hint:*to convert a number to a string, use`num2str`- in
`helloWorld`make a linear time vector`tVec`with 10,000 samples between 0 and endOfWeek - calculate the value of your knowledge at each of these time points using the same equation as above, store the result in a vector
- can you graph the results you just made?

**X**% of*Octave*" - your learning time constant is 1.5 days. calculate the
number of seconds in 1.5 days and put it in a variable

- some new
*$\Theta$:*- an emphasis on
*Octave*commands for**control flow**, meaning not simply executing each statement in a script one after the other in order, but executing a statment or block of statements only if some condition is true, or repeatedly as long as something is true, or one time after another repeatedly: **if***condition1*

*statements1*

**elseif***condition2*

*statements2*

**else**

*statements3*

**end**

Note you can have zero, one, or many of the**elseif**blocks, and the**else**is optional; the rest is required. This performs the block of statements*statements1*if*condition1*is true, otherwise, it performs the block*statements2*if*condition2*is true,*etc.,*until, if all conditions are false, it does the*statements3*.- conditional expressions:
*value1==value2**value1~=value2**value1<value2**value1<=value2**value1>value2**value1>=value2*

- note that the value of a conditional expression is
**0**if the result is**false**, while anything non-zero counts as**true**(although*Octave*will give it the value**1**if it has a choice) **while***condition*

*statements*

**end**

Keeps permforming the block of*statements*as long as*condition*is true.**for i=***vector*

*statements*

**end**

Performs the block of*statements*repeatedly, each time with the variable**i**set to the next value taken out of the*vector*. For example, one often uses:

**for i=1:n**

*statements*

**end**

to do the block of*statements***n**times.

- an emphasis on
*F:*- Discussion of
*Octave*functions: - in the file
**fun.m**, start (after optional comments — which are a good idea, and in this class should always have the students name, date, and problem) with the declaration

**function [ ret ] = fun ()**

*statements*

**endfunction**

Or, if there are multiple inputs and outputs, the form might br

**function [ a, b ] = fun (x, y)**

*statements*

**endfunction**

where the*statements*should, somewhere along the way, give values to the variables**a**,**b**. - Note a function can call itself, usually for a simpler
set of inputs than the original one — this is called
**recursion** - Did an example of a function to compute the factorial of a
number both by
*recursion*and**iteration***[the latter meaning simply using a loop]* - discussed how to time a command: the command
**time()**returns the number of seconds, accurate to milliseconds, since the*epoch*, which in the Linux/Unix world was midnigh January 1st, 1970. - wrote a function to compute the next pari of elements of the Fibonacci sequence, giving the current and last elements as input
- then wrote another function which computed the first 10,000 elements of the Fibonacci sequence but displaye the sequence of rations of successive elements, and then plotted this seqeuce of ratios
- please send me your versions of these Fibonacci functions, the diary showing your use of them, and the plot of ratios.

- Discussion of

**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*i*th 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**.

- In
- 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?

- Today we play
*$\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

- Remember: many
*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

- A little about variables in
*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**head**s (well, of**true**s). 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**

- The end goal of this week we're going to simulate
*$\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**

- In applied mathematics, techniques which use randomness where there is
none inherent in the original problem are called
*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')**

- Now, really do

*T:*- correcting something from last Friday: rather than simply making
various parts of your graphs invisible, if you are not going to need
them again, it is much better to
**delete**them - function comments, which are visible when you say
**help***functionname* - function types: anonymous, nested, private
- function handles
- the command
**fminbnd** - global variables

- correcting something from last Friday: rather than simply making
various parts of your graphs invisible, if you are not going to need
them again, it is much better to
*W:*- a little syntax: variable numbers of return values from functions
- Project for the week: writing some function which finds the maximum of a function, by several different methods.
- First, working with a function of a single variable (which is
specified by its function handle), make a direct, brute-force
maximizer
**BFmax(fun_handle, left_endpt, right_endpt, resolution)**which simply divides the interval from**left_endpt**to**right_endpt**into equal pieces of size**resolution**and checks the values of the function with handle**fun_handle**at those points, and picks the largest one.**BFmax**should return both the maximum value of the function on the given intervale and the*x*coordinate locating the maximum. - Please make good comments in your code, and do a lot of sanity checking of the inputs. You might also consider having the program check if the run is going to be very long and warn the user, maybe by doing the first 1000 test evaluations, and then figuring out how much longer it will need to run and reporting that to the user.
- Consider how your function will work if there are several equal maxima....
- Here's another method to find a maximum: pick a starting point in
the interval (one of the endpoints, or the midpoint, or a random
point in the middle — you choose), look at the value of the
function at that point, and at two nearby points (each a distance
**resolution**away) to the left and the right. If the function is bigger in one direction, keep moving in that direction as long as the function is increasing or until you hit the end of the interval. The place where you stop will be a kind of*hilltop*, so this method is called*hill-climbing*. Write a program**HCmax**. - If your
**resolution**is very small (because you want a very accurate answer), you do hill-climbing with steps of that size, and your interval is large, this program may run a long time. One way to improve this might be as follows: start out by taking steps of a size $\Delta x$ which is much bigger than**resolution**. Then when you are at a hilltop in the sense that the function at the current $x$ value is bigger than its value at $x-\Delta x$ and $x+\Delta x$, try hill-climbing again with a $\Delta x$ which is half its previous size. - What can be a problem of the hill-climbing method -- where might it get hung up at a point which was not the maximumn of the function on the given interval?
- One way to work around the problem is to start hill-climbing from a number of random points in the inverval and to compare the resulting hill-tops which you find. Program this as well.
- Add a feature to any of your programs above which animates the steps it is performing.
- Write a whole other set of programs, modeled on the above strategies,
which works for
*functions of two variables.*Again, an animation would be very cool.

*$\Theta$:*- We're heading towards
*Newton's Method*, for which we need the derivative of our function. So let's develop some code for a very simple type of function for which we will be able to work with the derivative as well —*polynomials*. Note*Octave*has built-in functions for dealing with polynomials, but we'll write our own (because fun). - Let's use vectors to stand for polynomials — the numbers in the
vector will be the coefficients of the polynomial. So if the vector is
**v**with components $v_1,\dots,v_n$, then we will think of it as representing the polynomial $$p(x)=v_1+v_2 x+\dots+v_n x^{n-1} = \sum_{j=1}^n v_j x^{j-1}\ .$$ [This is backwards to the way*Octave*stores the coefficients of polynomials.] - Write a function
**pval(p,x)**which takes a polynomial**p**(in our form as a vector) and a number**x**and computes the value of $p(x)$. - Write another function
**pderiv(p)**which takes a vector representing a polynomial and returns the vector representing the derivative of that polynomial. - Now we are in position to use Newton's method to find the zeros of a
function (we'' do a slightly different version for maximization soon).
Acutally, this works for any function whose derivative we know.

Pick a random point $x_0$. Then repeat $$x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)}$$ - Write a function which uses Newton's method to find a zero of a polynomial

- We're heading towards
*F:*- The key idea of Newton's methods is to pretend, at each iteration, that the function is equal to its Taylor polynomial of some (small) degree: for root-finding, we assume the function is linear and solve for where the line crosses the $x$-axis; for optimization, we assume the function is quadratic and solve for where it would have a local extremum.
- At a point $x_n$, the degree two Taylor polynomial is $y=f(x_n)+f^\prime(x_n)(x-x_n)+\frac12f^{\prime\prime}(x_n)(x-x_n)^2$ whose extremum happens when $x-x_n=-\frac{f^\prime(x_n)}{f^{\prime\prime}(x_n)}$, so we get the next iterative value by using $$x_{n+1}=x_n-\frac{f^\prime(x_n)}{f^{\prime\prime}(x_n)}\ .$$
- The first priority for today's in-class coding was to finish up and e-mail final versions of the programs from the last few weeks. If [and only if] you have completed that, you may proceed to writing a function which looks for maxima by using Newton's optimization method.

*T:*- a nice feature:
**tic**and**toc** - getting things out of a file, use
**fh=fopen**, then use the**fh**to read and write out of a file. - This week, we'll work on a project with strings, so today some basics
of built-in functions which handle strings:
- since a string is an array of characters, we can do simple direct
substitutions: e.g.,
**str="yadda yadda";str=str(str=='a')='o';** **repmat****isletter**and**isspace****upper**and**lower**

- since a string is an array of characters, we can do simple direct
substitutions: e.g.,
- here are some string-related tools which we did not have time to
discuss in class, but which might be useful (do
**help**for more information):*name*- stepping through ranges of letters
**strtrim**and**deblank****strfind**,**strmatch**, and**strtok****strcmp**,**strcmpi**,**strncmp**, and**strncmpi****char**and**abs****dec2hex**,**hex2dec**,**dec2bin**, and**bin2dec****num2str**[**mat2str**] and**str2num****sprintf**and**fprintf**- Problems with arrays of strings (if they have different lengths)!
But there is
**strvcat** **sortrows**and**strjust**

*cell arrays*! Like vectors (arrays) but use { } instead of [ ]:**class**- if
**c={[1 2],[1 2; 3 4],'a string'}**then look at**c(1)**compared to**c{1}** - we can do
**c{4}='another string'**or**c(4)={'another string'}**, both having the same effect - generally, we would not do
**c{4}={'another string'}**as it nests cells so we would need to access it by**c{4}{1}** - we can access groups of cells by operations like
**c{1:4}**in order to do simulataneous assignments or passing multiple arguments to a function **cellfun(fun_handle, cell_array [, 'UniformOutput', true/false])**

- again, more elaborate cell operations which we did not cover but
which could be useful:
**num2cell**,**mat2cell**, and**cell2mat****cellstr**is a kind of inverse to**strvcat****strmatch**with cell arrays

- a nice feature:
*W:***textscan(file_descriptor,**(where pattern might simply be "%s") [similar to*pattern*)**strread**]- use the string and cell array yoga discussed this week to write a
program which
- opens a file (use
*fopen*) - reads all the words of the file into a cell array (use
*textscan*) - convert the words to lower case (use
*lower*) - remove cells that have no letters at all (
*isletter*helps) - remove punctuation and other non-letter characters from all cells
(again,
*isletter*; you could do this step before the previous, then you would remove the entirely empty words) - find the unique words in the text (
*unique*works on cell arrays) - count how often each word occurs (
*histc*can do it for you, or you can do it yourself) - sort the unique words by order of decreasing frequency
(
*sort*) - print out the words in this order, each with its frequency

- opens a file (use

*$\Theta$:*- continuing the big weekly project

*F:*- still continuing the big weekly project
- a function which could be useful (although you can also use
*isletter*) is**isstrprop(str,**where*'propertyname'*)*propertyname*is one of:*alpha:*True for characters that are alphabetic (letters).*alnum*or*alphanum:*True for characters that are alphabetic or digits.*lower:*True for lowercase letters.*upper:*True for uppercase letters.*digit:*True for decimal digits (0-9).*xdigit:*True for hexadecimal digits (a-fA-F0-9).*space:*and*wspace:*True for whitespace characters (space, formfeed, newline, carriage return, tab, vertical tab).*punct:*True for punctuation characters (printing characters except space or letter or digit).*cntrl:*True for control characters.*graph:*or*graphic:*True for printing characters except space.*print:*True for printing characters including space.*ascii:*True for characters that are in the range of ASCII encoding.

*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**, or (and this is more efficient and stable) using the^{-1}y*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**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.*functionname*

*$\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:- 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*." - 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.

- 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:

*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!)

- To take a string

- 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
*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.

- If you are ambitious, do as much of the above as you can for the
*F:*- Finish your code — your midterm! Some things to keep in mind:
- Your code should run correctly -- make sure you check it completely!
- Make sure your code sanity checks all its input.
- 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. - 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**.

- Finish your code — your midterm! Some things to keep in mind:

*T:*- a useful shell (general Linux) command:
**wget**, downloads the page at that URL into your current local directory.*URL* - Some useful
*Octave*commands relating to matrix operations:**eig****rref****null****norm**, such as**norm(A,"fro")**

- a useful shell (general Linux) command:
*W:**$\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
- draw only half of the radial line segments, the half which ends on the unit circle;
- draw only half of the linear transformed vectors, again the half which touches the unit circle;
- 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.*) - 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!

- Here's a supplemental project you can do to make an optional,
*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:***fft**— [Fast] Fourier Transform**ifft**— Inverse [Fast] Fourier Transform**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.*

- talked about some ideas of basics of
*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 bew = 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**

- By this date (but much earlier —
*W:**$\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.

- When plotting a region in space, the way to do it is to
plot its

*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.

- This weeks work should consist of your graphs and code for the
parametric surfaces

*T, W, $\Theta$, & F:*Working on final projects. In order for this to be productive time, please take 2-5 minutes at the beginning of class to send an e-mail to your instructor to say what you are going to be working on that class period, and then again at the end to say how it went.- See the assignment due next Monday.

*M:*No class meeting of course (we never have class on Monday), but please send an e-mail to your instructor with:- Current status of your final project, in an over-all sense of where you are in terms of getting data, researching theoretical background, designing code, implementing/debugging code, and so on.
- Give a fairly complete list of the references you have consulted: web
pages, books, articles,
*etc.* - Talk a bit about the over-all architecture of your project: problem-solving or calculation approach, and corresponding program architecture.
- If there is any data to be used in your work, explain where you are getting it and how you will have to (or already have) process it to be useful for your program.
- Attach code samples, if you have them.
- Please state clearly what you think is the biggest remaining obstacle to your finishing this project, and how you are planning to overcome it.
- If there is anything you would like help on, please state that clearly in this status report!

*T, W, $\Theta$, & F:*Working on final projects. In order for this to be productive time, please take 2-5 minutes at the beginning of class to send an e-mail to your instructor to say what you are going to be working on that class period, and then again at the end to say how it went.

*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:- $\Theta$:
- Jamen

- F:
- Correy
- Graham
- Joel

- $\Theta$:
- There are some notes on the project description web page about what is expected for your report and presentation.

**Exam week**, no classes.- The scheduled final exam times are
**Monday, April 28th from 1-3:20pm**and**Tuesday, April 29th from 3:30-5:50pm, both in our usual classroom** - Scheduled final project presentationa:
- Monday:
- Chris
- Matthew
- Shelby
- Stephen

- Monday:

Jonathan Poritz (jonathan.poritz@gmail.com) |
Page last modified: Monday, 21-Jul-2014 05:14:10 UTC |