Lab 3: Programming in Matlab and basic techniques in quadrature (solving integrals numerically)

Discussion

A program is basically a list of MATLAB commands saved in a "m-file" (a file with extension .m)

There are  scripts (a list of commands that executes in order) and functions, which allow arguments to be passed and values to be returned. It is the latter kind of program that we will concentrate on.

Preliminaries:

Start by creating a directory on drive C; with your_name as its name (it will be automatically deleted tonight so don't leave anything in it you want later).

Start MATLAB and open a diary file. For this lab you will pass in your diary as well as some of your programs.

It will be useful to add your directory to MATLABs path. Type:
 

>> path(path,'c:\your_name')

You can see the whole path

>> path

A very simple example:

Use word pad to create a file called "c2f.m" with the following two lines and save it in your directory

function f = c2f(c)
f=(9/6)*c+32;

As an alternative to word pad, you can type "edit" at the MATLAB prompt to use the MATLAB editor. It has the nice feature that it "colorizes" your program.

Remember, the name of the function must be the same as the name of the m-file!

After saving, go back to MATLAB and  type something like

>> c2f(22)
>> c2f(0)
>> c2f(100)

etc.

A function can have several arguments and can return several values.

Example 2:

The next function takes two 1x2 vectors as arguments and returns two numbers. Can you figure out what it does?

function [m,b]=tpline(p1,p2)
m=(p2(2)-p1(2))/(p2(1)-p1(1));
b=p2(2)-m*p2(1);

Try

>> a=[0 1]
>> b=[1 0]
>> line(a,b)

Note that you only get one of the two outputs. To get them both type

>> [m,b]=line(a,b)

This is what we will  use to pass the integrand to our quadrature (i.e., numerical integration) programs.

We can pass built-in functions or functions of our own design. For example, suppose we want to integrate sin^2(x).  The the function we will pass (sin2.m) can be defined:

function y = sin2(x)
y=sin(x).^2;

Please define this function now.

Next, we need a  function that will take sin2 (or any function) as its input. We start with something really simple (not really a quadrature).

function s=summ(fun,n)
s=0;
for i=1:n
    s=s+feval(fun,i);
end
 

Before running this, several comments:

1. Note the use of the "for" loop. MATLAB also has an "if...else..." command and many other usual and familiar  control structures. Look in your MATLAB book or one of the web tutorials for details.

As an alternative, if you want to know how "while" works, type

>> help while

2. Having introduced control structures, let me comment that in MATLAB it's not always best to use loops.

For example, to compute the sum of the squares of the entries of a vector x, it is probably faster (as well as simpler to program) to compute

y=x*x';

instead of the more conventional

y=0;
for i=1:length(x)
    y=y+x(i)^2;
end

This takes some getting used to if you've done a lot of functional programming.

In this sense MATLAB is psychologically similar to parallel computing (although it's probably not really executed in parallel)

3. Important: the "feval" command allows the passed function name "fun" to be evaluated in the program. See below when we execute "summ" for one more crucial feature of this "name passing" process.

Here is the evaluation:

>> summ('sin2',10)

What is the second crucial feature I'm talking about?

Assignment

part 1.

We want to estimate

     / b
     |
     |    f(x)dx
     |
    /  a
 

using a quadrature with n steps.

Write m-files that take arguments f, a, b and n and perform the following rules:

* left endpoint
* right endpoint
* trapezoid
* Simpson's rule

Hint: Try to use vector operations rather than loops whenever possible. Can you come up with a snappy way to compute a vector of the weights [1 2 2 2 2...2 1] or [1 2 4 2 4 .. 4 2 1] for the trapezoid rule or Simpson's rule.

Now compute the following integrals for various rules and values of n.

Compare to the actual values (well known numbers)

     / 1
     |
a.  |    (x^2+3*x+5)dx
     |
    /  0
 

     / 4
     |
b.  |    sqrt(16-x^2)dx
     |
    /  0

     / 1
     |
c.  |    1/x dx
     |
    /  0.1
 

Hint: Try long format.

Plot or print out a table showing how the error depends on n for different rules.

You might also want to measure how many flops the computation uses.

Can you find n large enough to see the effects round-off error? (Note: because MATLAB uses double precision this is exceedingly difficult, but if you choose a bad enough function you might be able to see it).

We will be able to see round off error using the variable precision available in  Maple).

---

 Part 2.

Now we let's compare our homegrown routines with MATLAB's built-in ones.

MATLAB uses adaptive quadrature and always aims for a stepwise that results in a specific very small error (how would you modify your programs to aim for a specific error?)

Try

>> help quad
>> help quad8

Try each of these out on your functions. Compare the number of flops necessary in each to get 8 place accuracy.

Discussion (Continued)

Now let's consider double integrals. Make an m-file

function z=saddle(x,y)
z=1-(x.^2-y.^2);

>> dblquad('saddle',0,1,0,1)

Assignment

Part 3. (Option 1)

     / 1   / 1
     |      |
     |      |    log[ (t^2+1)^2-2t(1-t^2){cos(2 pi x)+cos(2 pi y)}]dxdy
     |      |
    /  0  / 0

It is Onsager's 1944 famous analytic solution of the Ising model (see B. Cipra, An Introduction to the Ising Model, American Mathematical Monthly, Vol. 94, No. 10, Dec., 1987). To interpret this it is important to see how it behaves as a function of t.
 

In MATLAB we can make the parameter t a global variable. Look at the MATLAB book or type

>> help global

We can use this to pass t to a function file that is a function of x and y. Make t global for the session too, and try integrating the function for several values of t.

Now we can write a function that takes T as an argument, passes it to the global variable t and outputs the result, say as z. Call this function "phase.m"

It should then be possible to plot z (see the plots at the end of the last lab).

When you get a graph you like, you can choose "export" from the file menu, save your picture, and post it to your web page. A good format for the web is "jpeg". You will need to use "ftp" to get the picture to your web page. Be sure to send it in "binary" rather than text format.

By the way, if you were passing in a paper report, say in TeX, you might want to export your graph in encapsulated postscript, another option on the export menu.

Comment: This will be challenging but please try. Also note, the program will be SLOW!!!

part 3: (Option 2)

Let p(u₁,...,u_n) be a nonzero polynomial with integer coefficients in the d variables (u_i)^(+-1).

The integral
 

                             / 2pi  / 2pi
                             |       |
 (2pi)^(-n) * exp [   | ...   |    log |p( exp(i t₁),...exp(i t_n)| dt₁...dt_n ]
                             |       |
                            /  0   / 0

This is called the Mahler measure of p and denoted M(p).

When d=1 it follows from "Jensen's equality" that if p=c₀+c₁ u+...+c_m u^m
then M(p)=|c_m| * (the product of all roots bigger than 1).

Try this with p(u)=u²-u-1.

Try this with p(u)=1+u-u³-u⁴-u⁵-u⁶-u⁷+u⁹+u¹⁰.

This should give M(p) approximately 1.176... as computed in 1933 by D. H. Lehmer. This is called the Lehmer number.

No polynomial has ever been found with a smaller Mahler measure (bigger than 1). Lehmer's question is, does such a polynomial exist? This is still an unsolved problem.

Can you write a function that will compute Mahler measures? (i.e. that takes the polynomial as an input and outputs the Mahler measure). Note: the polynomial should be entered as a function.

How about for p a polynomial in several variables? The challenge is to make the number of variables an input parameter.

Try this with p(u,v)=1+u-v+2uv+u²+v²

Try with a polynomial in 3 variables.

How close can you get to the Lehmer number? Can you automate your search?

You can find many references on this problem if you serach the web under "Mahler Measure".

What to pass in (on your secret web page):
 

1. A diary of the entire lab. You can patch together several sessions. Please edit it for readability and include comments as appropriate.
2. Your m-files for the trapezoid rule and Simpson's rule.
3. Your m-file phase.m (or your mahler.m if you do option 2)
4. The graph in jpeg format (if you do option 1). In option 2 some polynomials with good Lehmer numbers.