MATH2071: LAB 1(b): Selected topics in Matlab

Introduction Exercise 1
Adaptive Quadrature Exercise 2
Roundoff errors Exercise 3
  Exercise 4
  Exercise 5
  Exercise 6
  Exercise 7
  Exercise 8
  Exercise 9


There are two versions of the first lab. This version discusses some special topics and is intended for students who took Math 2070. If you have not already taken Math 2070, please see Lab 1(a). That version of the first lab introduces the Matlab environment and programming language, and presents the general format of the work you need to hand in. Students who take Lab 1(a) will not be in no way disadvantaged in this course because they ``missed'' Lab 1(b). If you prefer, you will find a version of this lab in Adobe pdf format here.

This version of the first lab is intended only for students who have already taken Math 2070.
This lab is concerned with several different topics. It covers material that is supplemental for students in Math 2071, but new students will not be shortchanged when they miss this material in favor of the introductory material presented in Lab 1(a).

The first topic discussed in this lab is a simple approach in a two-dimensional adaptive integration routine. The difference between the work in this lab and in last semester's homework assignment is that this lab will be using square mesh elements instead of triangular ones and an elementary technique for determining which mesh elements need to be refined in order to meet the error requirements.

The second topic is a demonstration of how roundoff error arises in a matrix calculation.

The third topic is a brief introduction to ordinary differential equations. This topic is also covered in Lab 1(a). The methods presented here will be used to motivate some of the linear algebra algorithms in the following labs.

Adaptive Quadrature

In this section you will construct a Matlab function to compute the integral of a given mathematical function over a square region in the plane. One way to do such a task would be to regard the square to be the Cartesian product of two one-dimensional lines and integrate using a one-dimensional adaptive quadrature routine such as adaptquad from last semester. Instead, in this lab you will be looking at the square as a region in the plane and you will be dividing it up into many (square) subregions, computing the integral of the given function over each subregion, and adding them up to get the integral over the given square.

The basic outline of the method used in this lab is the following:

  1. Start with a list containing a single item: the square region of integration.
  2. Use a Gaußian integration rule to integrate the function over each subregion in the list and estimate the resulting error of integration. The integral over the whole region is the sum of the integrals over the subregions, as is the estimated error.
  3. If the total estimated error of the integral is too large, find the subregion with largest error, divide it into four subregions, and replace the single subregion with four smaller subregions, and return to the previous step.
The way the notion of a ``list'' is implemented will introduce a data structure (discussed in detail below) that is more versatile than arrays or matrices.

Two-dimensional Gauß quadrature

One simple way of deriving a two-dimensional integration formula over a square is to use iterated integration. Recall that a one-dimensional Gauß integration rule can be written as

$\displaystyle \int_x^{x+h} f(x) dx \approx \sum_{n=1}^N w_nf(x_n).$ (1)

Here, $ N$ is the index of the rule. For the case $ N=2$, the points $ x_n$ are $ x+h/2\pm h/(2\sqrt{3})$ and the weights are $ w_1=w_2=h/2$. The degree of precision is 3, and the error is proportional to $ h^5 max\vert f''''\vert$. (If you look up the error in a reference somewhere, you will notice that the error is usually given as proportional to $ h^4$, not $ h^5$. The extra power of $ h$ appearing in (1) comes from the fact that the region of integration is $ [x,x+h]$.) Applying (1) twice in the iterated integral

$\displaystyle \int_x^{x+h}\int_y^{y+h} f(x,y) dxdy \approx \sum_{n=1}^N \sum_{m=1}^Mw_nw_mf(x_n,y_m).$ (2)

For the case $ N=M=2$, (2) becomes

$\displaystyle \int_x^{x+h}\int_y^{y+h} f(x,y) dxdy \approx \sum_{n=1}^4(h^2/4)f(x_n,y_n),$ (3)

where the points $ (x_n,y_n)=(x+h/2\pm h/(2\sqrt{3}),y+h/2\pm h/(2\sqrt{3}))$. The error is $ O(h^6)$, and (3) is exact for monomials $ x^ny^m$ with $ n\leq3$ and $ m\leq3$, and for sums of such monomials. In the following exercise, you will implement this method in Matlab.

Exercise 1:
  1. Write a Matlab function to compute the integral of a function over a single square element using (3) with $ (x_n,y_n)=(x+h/2\pm h/(2\sqrt{3}),y+h/2\pm h/(2\sqrt{3}))$. Name the function m-file q_elt.m and have it begin
    function q=q_elt(f,x,y,h)
    % q=q_elt(f,x,y,h)
    % INPUT
    % f=???
    % x=???
    % y=???
    % h=???
    % OUTPUT
    % q=???
    % your name and the date
  2. Test qelt on the functions $ 1$, $ xy$, $ x^2y$, $ x^2y^2$, and $ x^3y^3$ over the square $ [0,1]\times[0,1]$ and show that the result is exact, up to roundoff.
  3. Test qelt on the function $ x^4y^4$ to see that it is not exact, thus showing the degree of precision is 3.

Error estimation

In order to do any sort of adaptive quadrature, you need to be able to estimate the error in one element. Remember, this is only an estimate because without the true value of the quadrature, you cannot get the true error.

Suppose you have a square element with side of length $ h$. If you divide it into four sub-squares with sides of length $ h/2$, then you can compute the quadrature twice: once on the single square with side of length $ h$ and once by adding up the four quadratures over the four squares with sides of length $ h/2$. Consider the following figure.

% big square
\put(0,0){\lin... 1-4
Denote the true integral over this square as $ q$ and its approximation over the square with side of length $ h$ as $ q_h$. Denote the four approximate integrals over the four squares with sides of length $ h/2$ as $ q^1_{h/2}$, $ q^2_{h/2}$, $ q^3_{h/2}$, and $ q^4_{h/2}$. Assuming that the fourth derivative of $ f$ is roughly constant over the squares, the following expressions can be written.
$\displaystyle q_{h}$ $\displaystyle =$ $\displaystyle q+Ch^6$  
$\displaystyle q^1_{h/2}+q^2_{h/2}+q^3_{h/2}+q^4_{h/2}$ $\displaystyle =$ $\displaystyle q+4C(h/2)^6=q+(4/64)Ch^6.$  

This system of equations can be solved for the error in $ q_h$ as

error$\displaystyle =Ch^6=\frac{64}{60}\left(q_{h}- (q^1_{h/2}+q^2_{h/2}+q^3_{h/2}+q^4_{h/2})\right)$ (4)

In the following exercise you will write a Matlab function to estimate both the integral and the error over a single element.

Exercise 2:
  1. Write an m-file named qerr_elt.m to estimate the error according to (4). Use q_elt.m to evaluate the integrals. qerr_elt.m should begin
    function [q,errest]=qerr_elt(f,x,y,h)
    % q=qerr_elt(f,x,y,h)
    % more comments
    % your name and the date
  2. Use qerr_elt to estimate the integral and error for the function $ 16x^3y^3$ over the square $ [0,1]\times[0,1]$. Since the exact integral is 1, and the method as degree of precision equal to 3, both the error estimate and the true error should be zero or roundoff.
  3. Use qerr_elt to estimate the integral and error for the function $ 25x^4y^4$ over the square $ [0,1]\times[0,1]$. You should observe that the estimated error is very close to the true error because the fourth derivative of the function is constant.

A versatile data storage method

Up to now all the Matlab programs you have used involved fairly simple ways of storing data, involving variables, vectors, and matrices. Another valuable programming tool for storing data is the so-called ``structure.'' In many programming languages, such as Java, C and C++, it is called a ``struct,'' in Pascal it is called a ``record,'' and in Fortran, it is called a ``defined type.'' In Matlab, the term ``structure'' is used, although everyone will understand if you call it a ``struct.''

You can find detailed information about structures at this Matlab page

A structure consists of several named sub-fields, each containing a value, and separated from the variable name with a dot. Rather than going into full detail here, consider just the simple concept of a square region in space, with sides parallel to the coordinate axes. Such a square can be specified with three numerical quantities: the x and y coordinates of the lower left corner point, and the length of a side. These three quantities will be called x, y, and h for the purposes of this lab. Thus, if a Matlab (structure) variable named elt were to refer to the square $ [-1,1]\times[-1,1]$, it could be given as

It is important to realize that the value of elt.x is unrelated to the value of x that might appear elsewhere in a program. For the purpose of this lab, two other quantities will be included in this structure: the approximate integral of the function over this element, called q, and the estimated error, called errest.

In the following exercises, you will be using a subscripted array (vector) of structures to implement the notion of a ``list of elements.'' Structures can be indexed, and the resulting syntax for the $ k^$th entry of the array of structures named elt would be elt(k). The sub-fields of elt(k) are denoted


Exercise 3: In this exercise you will build up a function that estimates the integral of a function and its error over a square by dividing the square into many smaller squares, all the same size.
  1. Begin a function named q_total.m with the following code.
    function [q,errest]=q_total(f,x,y,H,n)
    % [q,errest]=q_total(f,x,y,H,n)
    % more comments
    % n=number of intervals along one side
    % your name and the date
    h=( ??? )/n;
    for k=1:n
      for j=1:n
        elt(elt_count).x= ???
        elt(elt_count).y= ???
        elt(elt_count).h= ???
        % the following two values will be corrected later
        elt(elt_count).q= h^2;
    for k=1:length(elt);
  2. Test the partially-written function q_total by choosing any function f (since it is unused so far, it does not matter) and using it to estimate the integral over the square $ [0,1]\times[0,1]$ using $ n=10$. Since it actually is computing the area of the square, you should get 1.0. If you do not, you have either computed the value of h incorrectly or you have somehow generated the wrong number of elements. The length of the vector elt should be precisely $ n^2$.
  3. As a second test, apply it to the square $ [-1,1]\times[-1,1]$ using $ n=13$. Again, you should get the area of the square.
  4. Use the function qerr_elt to estimate the values of q and errest based on the elemental values of x, y, and h and place them into elt(elt_count).q and elt(elt_count).errest.
  5. Estimate the integral and error of the function $ 16x^3y^3$ over the square $ [0,1]\times[0,1]$ for the value $ n=1$. If you do not get 1.0, you have computed elt(elt_count).x or elt(elt_count).y or elt(elt_count.h incorrectly or used qerr_elt incorrectly. Note: length(elt) is precisely 1 in this case.
  6. Estimate the integral and error of the function $ 16x^3y^3$ over the square $ [0,1]\times[0,1]$ for the value $ n=2$. If you do not get 1, review your changes to q_total.m carefully.
  7. Fill in the following table for the integral of the function $ 25x^4y^4$ over the square $ [0,1]\times[0,1]$.
      n   integral    estimated error   true error
      2   ________     ____________     __________
      4   ________     ____________     __________
      8   ________     ____________     __________
     16   ________     ____________     __________
  8. Are your results consistent with the global order of accuracy of $ O(h^4)$?

In order to further test the integration and error estimation, a more complicated function is needed. One function that is neither too easy nor too hard to integrate is the following function.

f(x,y)= \frac{1}{\sqrt{1+100 x^2+100(y-0.5)^2}}+\\

This function has three peaks on the square $ [-1,1]\times[-1,1]$: two at $ (\pm 0.5,-0.5)$ and a much sharper one at $ (0,0.5)$. Its integral over the square $ [-1,1]\times[-1,1]$ is 1.35231319384. Matlab code to effect this function is:
function z=three_peaks(x,y)
% z=three_peaks(x,y)
% three peaks at (-.5,-.5), (+.5,-.5), (0,.5)
% the integral of this function over 
%   [-1,1]X[-1,1] is 1.75522375591726

z=1./sqrt(1+ 100*(x+0.5).^2+ 100*(y+0.5).^2)+ ...
  1./sqrt(1+ 100*(x-0.5).^2+ 100*(y+0.5).^2)+ ...
  1./sqrt(1+ 100*(x    ).^2+ 100*(y-0.5).^2);
A perspective plot of the function is:

Exercise 4:
  1. Use cut-and-paste to copy the above code to a function m-file named three_peaks.m.
  2. Integrate three_peaks over the square $ [-1,1]\times[-1,1]$ and fill in the following table. The true value of the integral is 1.75522375591726. Warning: the larger values of $ n$ will take some time-be patient.
      n   integral    estimated error   true error
     10   ________     ____________     __________
     20   ________     ____________     __________
     40   ________     ____________     __________
     80   ________     ____________     __________
    160   ________     ____________     __________
  3. Are the true error values consistent with the convergence rate of $ O(h^4)$?
  4. Notice that the estimated errors are much larger than the true errors, especially for larger values of $ n$. This is because the elemental errors are sometimes positive and sometimes negative and should cancel each other, but we take absolute values in the code for q_total. Make a copy of q_total.m called q_total_noabs.m and remove the absolute value from the summation of the elemental error estimates. Using q_total_noabs, compute the integral of three_peaks over $ [-1,1]\times[-1,1]$ using $ n=80$. You should observe that the true and estimated errors agree much more closely. Nonetheless, use of absolute value is adequate for smaller values of $ n$ and is more conservative in all cases.

An adaptive strategy

The objective of this discussion of quadrature is to present an adaptive strategy for quadrature. You have seen all the pieces and now it is time to put them together. The adaptive strategy used here is the following.

  1. Start with a vector named elt of structures as used in q_total above. This vector will have only one subscript:
    and the values represent the given square region over which the integral is to be taken, with elt(1).q and elt(1).errest computed using qerr_elt.
  2. Add up all the elemental values of q and absolute values of errest to get the total q and errest. If errest is smaller than the tolerance, stop and return the values of q and errest.
  3. If the total estimated error of the integral is too large, find the value of k for which abs(elt(k).errest) is largest and divide it into four subregions.
  4. Replace elt(k) with values from the upper right of the four smaller square subelements. You can use code similar to the following.
    % new values for this element
    [elt(k).q,elt(k).errest]=qerr_elt( ??? )
  5. Add three more elements to the vector of elements using code similar to the following for each one.
    elt(K+1).x= ???
    elt(K+1).y= ???
    [elt(K+1).q,elt(K+1).errest]=qerr_elt( ??? )
  6. Return to the second step above.

Exercise 5:
  1. Write a Matlab function m-file named q_adaptive.m that implements the preceeding algorithm. Your function could use the following outline.
    function [q,errest,elt]=q_adaptive(f,x,y,H,tolerance)
    % [q,errest,elt]=q_adaptive(f,x,y,H,tolerance)
    % more comments
    % your name and the date
    % initialize elt
    ??? more code ???
    for passes=1:MAX_PASSES
      % compute q by adding up elemental values
      % and compute errest by adding up absolute elemental values
      % use a loop for this because the "sum" function doesn't
      % work for structures.
      % if error meets tolerance, return
      % find the element with largest abs(errest)
      % replace that element with a quarter-sized element
      % add three more quarter-sized elements
  2. Test q_adaptive by computing the integral of the function $ 16x^3y^3$ over the square $ [0,1]\times[0,1]$ to a tolerance of 1.e-2. The result should be exactly correct because the degree of precision is 3.
  3. Test q_adaptive by computing the integral of the function $ 25x^4y^4$ over the square $ [0,1]\times[0,1]$ to a tolerance of 1.e-2. You should see that length(elt) is precisely 4 because only a single refinement pass was required.

    If you do not get the correct length, you can debug by temporarily setting MAX_PASSES=2 in the code and look at qelt. Is length(qelt equal to 4? If not, look at the coordinates of each of the elements in elt. There should be no duplicates or omissions. When you have corrected your error, do not forget to reset MAX_PASSES=500

In the following exercise you will see how the adaptive strategy worked.

Exercise 6:
  1. Download a plotting function plotelt.m that displays the elements. Elements colored green have small estimated error, elements colored amber have mid-sized error estimates and elements colored red have the largest error estimates. Red elements are candidates for the next mesh refinement.
  2. Use q_adaptive to estimate the integral of the function $ 25x^4y^4$ over the square $ [0,1]\times[0,1]$ to an accuracy of 1.e-5. What are the integral, the estimated error, and the true error? How many elements were used? You should observe that the exact and estimated errors are close in size.
  3. Use plotelt to plot the final mesh used. Please include a copy of this plot when you send me your work. You should observe that the error is largest near the upper and right edges and that several of the larger elements near the upper right corner have already been refines.
  4. Again estimate the integral of $ 25x^4y^4$ over the square $ [0,1]\times[0,1]$, but to an accuracy of 9.e-6, a little bit smaller than before. Plot the resulting mesh, and include a copy with your summary. You can see that a few more of the red elements have been refined.

Exercise 7: Use q_adaptive to estimate the integral and error of the function three_peaks over the square $ [-1,1]\times[-1,1]$ to a tolerance of 1.e-6. Include the result, the estimated error, the true error, and the number of mesh elements (length(elt)) with your summary. Plot the resulting mesh and include plot with your summary.

You should be able to see from the plot that elements near the peaks themselves have been refined, but also elements in the areas between the peaks.

Roundoff errors

Last term you saw some effects of roundoff errors. Later in this term you will look at roundoff errors again. Right now, though, is a good time to look at how some roundoff errors come about.

In the exercise below we will have occasion to use a special matrix called the Frank matrix. Row $ k$ of the $ n\times n$ Frank matrix has the formula:

$\displaystyle {\bf A}_{k,j}=\left\{ \begin{array}{cl}
0 & \mbox{ for }j<k-2\\
n+1-k & \mbox{ for }j=k-1\\
n+1-j & \mbox{ for }j\ge k

The Frank matrix for $ n=5$ looks like:

$\displaystyle \left(\begin{array}{ccccc}

The determinant of the Frank matrix is 1, but is difficult to compute numerically. This matrix has a special form called Hessenberg form wherein all elements below the first subdiagonal are zero. Matlab provides the Frank matrix in its ``gallery'' of matrices, gallery('frank',n), but we will use an m-file frank.m. The inverse of the Frank matrix also consists of integer entries and an m-file for it can be downloaded as frank_inv.m. You can find more information about the Frank matrix from the Matrix Market, and the references therein.

Exercise 8: Let's look carefully at the Frank matrix and its inverse. For convenience, define A to be the Frank matrix of order 6, and Ainv its inverse, computed using frank and frank_inv. Similarly, let B and Binv be the Frank matrix of order 24 and its inverse. Do not use the Matlab inv function for this exercise! You know that both A*Ainv and B*Binv should equal the identity matrices of order 6 and 24 respectively. Let's just look at the top left entry. Compute
  A(1,:)*Ainv(:,1)= __________
  B(1,:)*Binv(:,1)= __________
Both of these answers should equal 1. The first does and the second does not. To see what goes right, compute the terms:
 A(1,6)*Ainv(6,1)= __________
 A(1,5)*Ainv(5,1)= __________
 A(1,4)*Ainv(4,1)= __________
 A(1,3)*Ainv(3,1)= __________
 A(1,2)*Ainv(2,1)= __________
 A(1,1)*Ainv(1,1)= __________
           sum   = __________
Note that the signs alternate, so that when you add them up, each term tends to cancel part of the preceeding term. Now, to see what goes wrong, compute the terms:
 B(1,24)*Binv(24,1)= __________
 B(1,23)*Binv(23,1)= __________
 B(1,22)*Binv(22,1)= __________
 B(1,21)*Binv(21,1)= __________
 B(1,20)*Binv(20,1)= __________
 B(1,16)*Binv(16,1)= __________
 B(1,11)*Binv(11,1)= __________
 B(1,6) *Binv(6,1) = __________
 B(1,1) *Binv(1,1) = __________
You can see what happens to the sum. The first few terms are huge compared with the final sum of 1. Matlab uses 64-bit floating point numbers, so you can only rely on the first thirteen or fourteen significant digits in numbers like B(1,24)*Binv(24,1). Further, they are of opposing signs so that there is extensive cancellation. There simply are not enough bits in the calculation to get anything like the correct answer. (Remark: It would not have been productive to compute each of the products B(1,k)*Binv(k,1) for each k, so I had you do the five largest and then sampled the rest. I chose to sample the terms with an odd-sized interval between adjacent terms. Had I chosen an even interval-say every other term-the alternating sign pattern would have been obscured. When you are sampling errors or residuals for any reason, never take every other term!)

Ordinary differential equations

In this section you will see a brief introduction to solving differential equations.

In general, a first-order ordinary differential equation can be written in the form

$\displaystyle y'=f(x,y)$ (5)

where $ y'=\frac{dy}{dx}$. Such an equation needs an initial condition $ y(x_0)=x_0$. Perhaps the simplest method for numerically finding a solution of (5) is to use the ``explicit Euler'' method wherein a discrete selection of points $ x_k=x_0+h*(k-1)$ where $ h$ is some fixed step size and $ k=1,2,3,\dots$. Writing the approximate value of $ y(x_k)$ as $ y_k$, then Euler's explicit method can be derived by approximating the derivative $ (y')_k\approx(y_{k+1}-y_k)/h$ and writing

$\displaystyle y_{k+1}=y_k+hf(x_k,y_k)$ (6)

The differential equation

$\displaystyle y'=-y+\sin x$

with initial condition $ y(0)=0$ has an exact solution $ y(x)=.5*(e^{-x}+\sin x-\cos x)$. It also has an approximate numerical solution defined by Euler's formula as

$\displaystyle y_{k+1}= y_k+h*(-y_k+\sin x_k).$ (7)

In some sense, $ y(x_{k+1})\approx y_{k+1}$ We are going to look at how this expression evolves for $ x>0$.

Exercise 9: Copy the following text into a file named exer9.m and then answer the questions about the code.
function error=exer9(nsteps)
% error=exer9(nsteps)
% compute the solution of the differential equation
% y'+y=sin(x)
% starting at y=0 at x=0 using Euler's method
% and ending at x=25
% nsteps=number of steps taken

% Your name and the date


clear x y exactSolution
for k=1:nsteps
plot(x,y);  % default line color is blue
hold on
plot(x,exactSolution,'g'); % g for green line
legend('Euler solution','Exact solution')
hold off

  1. Add your name and the date to the comments at the beginning of the file.
  2. Run exer9 with nsteps=160. You can see that the approximate and exact solutions are quite close. Please include this plot with your summary.
  3. Run exer9 with nsteps=10. You can see that the approximate and exact solutions are not close and are diverging from each other rapidly. Please include this plot with your summary.
  4. Fill in the following table using exer9.
      nsteps      error       ratio
        10      _________    _______
        20      _________    _______
        40      _________    _______
        80      _________    _______
       160      _________    _______
    The convergence rate should be $ O((\mathtt{nsteps})^{-p})$ for some integer $ p$. What is your estimate of $ p$?

A slightly more complicated method for (5) is Euler's implicit method, which can be written in the following way.

$\displaystyle y_{k+1}=y_k+hf(x_{k+1},y_{k+1})$ (8)

This method is ``implicit'' because (8) is not a formula for $ y_{k+1}$ but instead must be solved to find $ y_{k+1}$. In the case under consideration, where $ f(x,y)=-y+\sin x$, it is easy to solve (8) for $ y_{k+1}$ because (8) is a linear equation. You have learned (or will soon learn) in class that the payoff for the greater complication of solving a system is that a larger stepsize can be taken while maintaining stability.

Exercise 10:
  1. Copy exer9.m to a new m-file named exer10.m. Modify it so that Euler's implict method is used instead of Euler's explicit method. Be sure to modify your comments as necessary.
  2. Run exer10 with nsteps=10. You should observe that, while the solution is not very accurate, at least it is not diverging. Please include a copy of this plot with your summary.
  3. Fill in the following table using exer10.
      nsteps      error       ratio
        10      _________    _______
        20      _________    _______
        40      _________    _______
        80      _________    _______
       160      _________    _______
    The convergence rate should be $ O((\mathtt{nsteps})^{-p})$ for some integer $ p$. What is your estimate of $ p$?
  4. You should be able to see that Euler's implict and explicit methods are about the same accuracy for large nsteps but that the implict method is much better when nsteps is small.

Back to MATH2071 page.

Mike Sussman 2009-01-06