notes from 11/2 – using colormaps, image, and imagesc; calculating the evolution of a one-dimensional array of cells under an averaging rule (see diffusionexample.m), checking the rate of convergence for a linearly convergent sequence using ratios of differences.

diffusionexample.m

Homework

Homework 1, Due 9/11 before 1:30

solutions

Homework 2, Due 9/25 before 1:30

solutions

wherezero.m whereis.m

Prob 1b: Using ‘break’ in a ‘for’ loop as a fast, straightforward, though less elegant, alternative to a while loop

Prob 2: A more straightforward version, that builds up only the factorial term.

Homework 3, Due 10/21 before class

pmult.m pdivide.m myprimes.m D.m recenter.m inorout.m

Homework 4, Due 11/1 before midnight

hwdata.mat contains two sequences X1 and X2. Download the file and then

>>load hwdata

will get the arrays into your workspace

Solutions of hw4

Homework5, Due 12/2 before class

Annotated MATLAB sessions, with topics

8/26/09

8/28/09

Mathematical operators and functions. Entering arrays, transposing. Creating random matrices. The colon operator. Array operators.

8/31/09

Using the colon operator. Array division. Accessing elements of an array – subsets of rows, subsets of columns, submatrices. Inserting and deleting rows and columns. Relational operators applied to arrays (applied elementwise). Logical addressing. Creating matrices of zeros, and matrices of all 1’s. Characters and strings in MATLAB. The size(A) function.

9/2/09 Notes on discussion of base 2

Scientific notation, base 2 representations, basic plotting. MATLAB representation of floating point numbers in base 2 using 64 bits (8 bytes)

9/4/09 script1.m Example of a very simple script file

9/9/09 mysquare.m

Several topics – use of the keyword ‘end’ in MATLAB for accessing elements of arrays, use of the ‘hold’ command in plotting, specification of a very basic MATLAB function mysquare, introduction to loops in MATLAB.

9/11/09 eulers.m taylor.m hilbert.m

Three script files illustrating the use of loops and how to program sums. eulers.m investigates the partial sums of the harmonic series and their relationship to the natural log function; taylor.m plots partial sums of the Taylor series for sin(x); hilbert.m constructs a matrix using nested loops to visit each (i,j) position of the matrix and define A(i,j) in terms of i and j. Please study these and think about how you might modify or extend these to address other situations or questions.

9/14/09 ramp.m

Use of the if statement, use of the if statement in defining a piecewise defined function. The while statement. Example of using ‘while’ to sum a series until the terms get to a certain size. The Ulam sequence – an application that uses ‘while’ and ‘if’ to generate a sequence. If n is even divide by 2, if odd the next n is 3n+1.

9/16/09 ulam.m ulam2.m ulamdata.m timedata.m

Modifying the ulam program to investigate behavior of the sequences. Tic/toc for timing computations. The function N=ulam(n) produces the Ulam sequence starting at n. The function m=ulam2(n) just gives the length m of the sequence starting at n, employing a counter. ulamdata.m is a script file for assembling an array L of sequence lengths, where L(n) is the length of the sequence with starting value n – the length of L is the variable r, which is set in the command window or another script file. timedata.m is a short script for gathering data on how long it takes to assemble L as a function of the size of L. It employs tic/toc.

9/18/09 riemann.m

Using functions as arguments (inputs) to other functions – just include an argument name in the list that will carry a function handle or function name, and use feval to evaluate function values: example: riemann(f,a,b,n) is a function that calculates Riemann sums for an arbitrary input function f. The argument f can be a function handle, such as f=@(x) x.^2+1, or a function name in quotes such as ‘exp’, or the name of a function you yourself have defined via an M-file.

9/21/09 rotations.m riemann.m

Three items:

1) How to chop the top off of a graph with vertical asymptotes (tan(x) in this case)

2) The polar curve r=cos(3Θ) graphed and rotated

3) The function riemann.m extended to include choice of point selection

9/23/09 Several programs:

sums.m Calculates the sum of the positive entries and the sum of the negative entries in an input matrix A. The algorithm is to visit each element of A using nested loops, and add A(i,j) into accumulating positive and negative sums according to the sign of A(i,j).

mysum.m A general function to sum all the elements of a matrix.

sums2.m Does the same job as sums.m, using a different method: Using B=A.*(A>0) we make all the negative entries in A equal to zero in B, while leaving the positive ones unchanged, then we use mysum.m to find the sum of all the elements in B to obtain the positive sum. Then mysum.m applied to B=A.*(A<0) will give us the negative sum. This example shows how we can use any functions we make within other functions.

mymin.m Finds the minimum value in a one-dimensional array and its index in the array. We simply start with the first element and compare one by one with subsequent entries, replacing the current minimum and the corresponding index whenever we encounter a smaller value.

mysort.m Sorts the elements of a one-dimensional array from lowest to highest by using mymin.m to successively pull out the minimum element of the remaining array, and appendng that element on the end of the sorted array that is being built up.

mysortr.m A recursive version of the sorting algorithm. Sometimes recursion is described as a function calling itself. In this case, the recursive approach is based on the fact that after we pull out the minimum element of an array, we are back in the same situation, in that we need to sort what’s left and tack it on. We simply tell MATLAB to sort what’s left, by calling the mysortr.m function again. Eventually, there is only one element left, which is already sorted, and that information propagates back up the sequence of function calls until the answer is obtained. This method is not recommended as it is very inefficient both in time and storage, but sometimes is very helpful in small problems. Imagine, for instance, that we want to sort a 5000 element array – the recursive method will stack up 4999 function calls waiting on answers!

9/25/09

The following are examples of images compressed in various formats. The can be imported into MATLAB as uncompressed images as either colormap images or “true color” RGB images.

flower.gif Example of a .gif image. The value of each element of the array corresponds to a rectangle color through the colormap when imported into MATLAB.

weather.jpg An image in the JPEG format. Each element of the array is given a 3-element color in the RGB system.

tempforecast.png An image from the National Weather Service in the .png format, also uses a colormap when imported.

9/28/09 Primarily review material

sumsquares.m Sum of squares of positive and negative entries in a matrix

usingbreak.m When does a sum first exceed a threshold? Using a for loop with a break when the threshold is first encountered is a very simple, if not so elegant, approach.

See also links associated with Homework 2 above, and review material for exam 1.

10/2/09 basep.m mypolyval.m

We developed two MATLAB functions.

1) M=basep(n,p)

creates the array M of digits of a (decimal) number n converted to base p. It employs a recursive algorithm: the first digit (from the right) is the remainder r when n is divided by p; then replace n by (n-r)/p and shift one place to the left, continuing the process. Recursively we have basep(n,p)=[basep((n-p)/r,p) r] where r is the remainder.

2) y=mypolyval(p,x)

evaluates the polynomial with coefficients in the array p (arranged from lowest to highest degrees) at the points in the array x. The algorithm is called Horner’s method, consists of successively starting with the last coefficient, multiplying by x, and adding the next coefficient going down.

10/5/09 Setting up a white-black colormap for arrays of 0’s and 1’s

Rule110 The secret of the universe in an eight-bit string.

wolfram09.m we redid this in class as

Some elementary cellular automata rules and their behavior A full discussion of the simplest nearest-neighbor rules and their behavior in one-dimensional cellular automata.

10/7/09 We completed the Wolfram one-dimensional cellular automata function cellauto.m

We note rule 22, which, starting from a single activated cell, leads to the Sierpinski triangle/gasket. Rule 110, is considered the most interesting rule as it produces for almost all inputs, discrete reproducible structures which persist, move around, and interact with each other in predictable ways.

10/9/09 Finishing up the Game of Life.

gol.m This is what we came up with. You might try modifying it to work with a general array size, rather than 300x300. This size could be an input argument.

acorn.mat A particularly long-lived initial condition

gun.mat This initial condition periodically produces (shoots out) a glider.

10/14/09 Discussion on saving data (see lecture notes above) and importing data from excel files. Saving variables as .mat files (internal) and as text files (for using elsewhere). Iteration (see lecture notes above) x(n+1)=f(x(n)). Convergence to a fixed point. Stable/unstable fixed points. Mathematical analysis: |f’(x*)|<1 is stable, |f’(x)|>1 is unstable, behavior of error near a fixed pt: e(n+1)=f’(c)e(n) where c is approximately x*. Looking at the ratio of successive errors using MATLAB.

10/16/09 More examples of the analysis of iterations. Written notes. For linearly convergent sequence close to its fixed point we have e(n)~=c*r^n for some r with |r|<1. A sequence that exactly behaves this way is the sequence generated by the partial sums of a geometric series (why?). Behavior of unstable iterations.

10/19/09 Investigation of the “chaotic map”, a very simple nonlinear map, but with generic behavior. Watch the parameter ‘a’, that’s where the action is. The iteration: x(n+1)=a*x(n)*(1-x(n)) For small values of a>1, there is a unique stable fixed point. At a=3 the fixed point “disappears”, replaced by a period 2 sequential limit, that successively “splits” into period 4, 8, 16,… over a finite interval in ‘a’, after which we get “chaos”.

chaos.m a simple script file for following the evolution of the sequence x(n) at different values of ‘a’.

10/21/09 Two more famous iterations, this time in the plane, so x(n+1)=f(x(n)) but this time x is a 2x1 array.

The Henon map is a two-dimensional nonlinear map with two parameters. The script file henon.m follows the path of a point under this map for given values of the parameters. The parameter set given results in the famous “strange attractor” but other values of the parameters can result in a single fixed point or periodic behavior. The script file fern.m is an example of a so-called iterated function system. Here, one of four linear transformations is chosen according to specified probabilities – these probabilities, as well as the entries in the linear transformations, are the parameters of the problem, so there are a lot of parameters here. We also started Newton’s method.

10/23/09 Implementing Newton’s method and numerically analyzing the results. The theoretical derivation of quadratic convergence in Newton’s method. Newton’s method for solving n nonlinear equations in n independent variables – idea is the same, linearization of the (vector) function using the Jacobian, and then solving the linear problem to give the next guess. The MATLAB commands look almost the exactly the same!

10/26/09 Newton’s method for a 2x2 system. Written notes including Newton’s method for systems and introduction to interpolation.

10/28/09 Discussion of interpolation; polynomial interpolation: assembling the interpolation matrix A, finding the coefficients via c=A\y.

10/30/09 More interpolation problems, using our polynomial interpolation function pinterp.m Graphical data input using MATLAB’s command ginput.

11/2/09 Review material on colormaps, evolution of a diffusion-type cellular automoton, and analysis of sequences for linear convergence.

11/6/09 Experiments with interpolation:

gdatainput.m A script file for inputting interpolation data using the mouse – you can input however many points you like wherever you like, right click on the mouse to end and interpolation will be carried out.

interpex1.fig A figure generated using gdatainput.m that illustrates the fact that a polynomial interpolant can have wild swings many times larger than swings in the data.

pinterp2.m Uses mypolyval.m for polynomial evaluation, the only difference with pinterp.m

lagrange.m Given a set of interpolation abscissas (x values), the lagrange fundamental polynomial l_j(x) represents the effect of datavalue y_j on the interpolant; it interpolates data which is zero at all the datapoints except at x_j where its value is 1. When the number of interpolation points gets larger, the lagrange polynomial l_j(x) can be huge, and far away from x_j.

The example of Runge is discussed. This is the function 1/(1+x^2) on the interval [-5,5]. As more interpolation points are taken, the interpolant grows worse and worse as an approximation of the underlying function.

11/9/09 Some notes on interpolation error.

We look at interpolation error. We obtain a “global upper bound” for error in the case of uniformly spaced interpolation points – this turns out to be too pessimistic in the case of our experiment with sin(pi*x). We explore, in the function cubeinterp.m, cubic interpolation in a table of function values (sometimes called a “lookup table”) that would be stored data in a computer. Given a value of x, we look at the four nearest stored values, interpolate those points with a cubic, and evaluate the interpolant at x to give our approximation of f(x). The error estimate here is much better because the operation is “local”, uses data that is localized over a relatively small interval.

11/11/09 Written class notes: Cubic interpolation in a table, and least squares.

We review the error estimate for interpolation in a table and apply the estimate to the case of tabulated values of sin(x), using our function cubeinterp.m . We give a simple example of least squares data fitting – the problem is formulated much like interpolation, except there is too much data to interpolate, and we instead choose the coefficients that result in the smallest value of the sum of the squares of the interpolation errors at the datapoints. The equations that give the least squares solution are called the normal equations: The least squares solution of Ac=y is given by the solution of A’*Ac=A’*y.

11/13/09 Least squares datafitting. Several examples are presented. You can input your own data using pointinputter.m . Approximating parametric curves: We developed a parameter T using arc length and then separately approximated x(t) and y(t) using the (T,X) and (T,Y) data respectively. Using trigonometric polynomials instead of polynomials. You should know how to set up and compute these problems yourselves without the prepackaged functions.

pdatafit.m A function to fit (x,y) data with a polynomial

tdatafit.m A function to fit (x,y) data with a trigonometric polynomial

pointinputter.m A script file for inputting points with the mouse (right-click when done)

parametriccurve.fig An example of approximating a parametric curve, using a 10^th degree trigonometric polynomial. The data is in parametricdata.mat