Introduction | Exercise 1 |

Adaptive Quadrature | Exercise 2 |

Roundoff errors | Exercise 3 |

Exercise 4 | |

Exercise 5 | |

Exercise 6 | |

Exercise 7 | |

Exercise 8 | |

Exercise 9 |

Introduction

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.

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:

- Start with a list containing a single item: the square region of integration.
- 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.
- 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.

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

Here, is the index of the rule. For the case , the points are and the weights are . The degree of precision is 3, and the error is proportional to . (If you look up the error in a reference somewhere, you will notice that the error is usually given as proportional to , not . The extra power of appearing in (1) comes from the fact that the region of integration is .) Applying (1) twice in the iterated integral

For the case , (2) becomes

where the points . The error is , and (3) is exact for monomials with and , and for sums of such monomials. In the following exercise, you will implement this method in Matlab.

**Exercise 1**:- Write a Matlab function to compute the integral of a function
over a single square element using (3) with
. Name the
function m-file
`q_elt.m`and have it beginfunction 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

- Test
`qelt`on the functions , , , , and over the square and show that the result is exact, up to roundoff. - Test
`qelt`on the function to see that it is not exact, thus showing the degree of precision is 3.

- Write a Matlab function to compute the integral of a function
over a single square element using (3) with
. Name the
function m-file

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

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

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

**Exercise 2**:- 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 beginfunction [q,errest]=qerr_elt(f,x,y,h) % q=qerr_elt(f,x,y,h) % more comments % your name and the date

- Use
`qerr_elt`to estimate the integral and error for the function over the square . 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. - Use
`qerr_elt`to estimate the integral and error for the function over the square . You should observe that the estimated error is very close to the true error because the fourth derivative of the function is constant.

- Write an m-file named

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
, it could be given as

elt.x=-1; elt.y=-1; elt.h=2;It is important to realize that the value of

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
th entry of the
array of structures named `elt` would be `elt(k)`. The
sub-fields of `elt(k)` are denoted

elt(k).x elt(k).y elt(k).h elt(k).q elt(k).errest

**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.- 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; elt_count=0; for k=1:n for j=1:n elt_count=elt_count+1; 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; elt(elt_count).errest=0; end end q=0; errest=0; for k=1:length(elt); q=q+elt(k).q; errest=errest+abs(elt(k).errest); end

- 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 using . 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 . - As a second test, apply it to the square using . Again, you should get the area of the square.
- 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`. - Estimate the integral and error of the function
over the square
for the value .
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. - Estimate the integral and error of the function
over the square
for the value . If you do not
get 1, review your changes to
`q_total.m`carefully. - Fill in the following table for the integral of the function
over the square
.
n integral estimated error true error 2 ________ ____________ __________ 4 ________ ____________ __________ 8 ________ ____________ __________ 16 ________ ____________ __________

- Are your results consistent with the global order of accuracy of ?

- Begin a function named

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.

This function has three peaks on the square : two at and a much sharper one at . Its integral over the square 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**:- Use cut-and-paste to copy the above code to a function m-file
named
`three_peaks.m`. - Integrate
`three_peaks`over the square and fill in the following table. The true value of the integral is 1.75522375591726. Warning: the larger values of will take some time-be patient.n integral estimated error true error 10 ________ ____________ __________ 20 ________ ____________ __________ 40 ________ ____________ __________ 80 ________ ____________ __________ 160 ________ ____________ __________

- Are the true error values consistent with the convergence rate of ?
- Notice that the estimated errors are much larger than the true
errors, especially for larger values of . 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 using . You should observe that the true and estimated errors agree much more closely. Nonetheless, use of absolute value is adequate for smaller values of and is more conservative in all cases.

- Use cut-and-paste to copy the above code to a function m-file
named

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.

- Start with a vector named
`elt`of structures as used in`q_total`above. This vector will have only one subscript:elt(1).x elt(1).y elt(1).h elt(1).q elt(1).errest

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`. - 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`. - 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. - Replace
`elt(k)`with values from the upper right of the four smaller square subelements. You can use code similar to the following.x=elt(k).x; y=elt(k).y; h=elt(k).h; % new values for this element elt(k).x=x+h/2; elt(k).y=y+h/2; elt(k).h=h/2; [elt(k).q,elt(k).errest]=qerr_elt( ??? )

- Add three more elements to the vector of elements using code
similar to the following for each one.
K=length(elt); elt(K+1).x= ??? elt(K+1).y= ??? elt(K+1).h=h/2; [elt(K+1).q,elt(K+1).errest]=qerr_elt( ??? )

- Return to the second step above.

**Exercise 5**:- 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 MAX_PASSES=500; % initialize elt elt(1).x=??? ??? 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 end

- Test
`q_adaptive`by computing the integral of the function over the square to a tolerance of`1.e-2`. The result should be exactly correct because the degree of precision is 3. - Test
`q_adaptive`by computing the integral of the function over the square 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`

- Write a Matlab function m-file named

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

**Exercise 6**:- 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. - Use
`q_adaptive`to estimate the integral of the function over the square 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. - 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. - Again estimate the integral of over the square , 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.

- Download a plotting function

**Exercise 7**: Use`q_adaptive`to estimate the integral and error of the function`three_peaks`over the square 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 of the Frank matrix has the formula:

**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*You know that both`inv`function for this exercise!`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. ComputeA(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

where . Such an equation needs an initial condition . Perhaps the simplest method for numerically finding a solution of (5) is to use the ``explicit Euler'' method wherein a discrete selection of points where is some fixed step size and . Writing the approximate value of as , then Euler's explicit method can be derived by approximating the derivative and writing

The differential equation

In some sense, We are going to look at how this expression evolves for .

**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 FINAL_TIME=25.0; stepsize=FINAL_TIME/nsteps; clear x y exactSolution y(1)=0; x(1)=0; exactSolution(1)=0; for k=1:nsteps x(k+1)=x(k)+stepsize; y(k+1)=y(k)+stepsize*(-y(k)+sin(x(k))); exactSolution(k+1)=.5*(exp(-x(k+1))+sin(x(k+1))-cos(x(k+1))); end 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 error=norm(y-exactSolution)/norm(exactSolution);

- Add your name and the date to the comments at the beginning of the file.
- Run
`exer9`with`nsteps=160`. You can see that the approximate and exact solutions are quite close. Please include this plot with your summary. - 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. - Fill in the following table using
`exer9`.nsteps error ratio 10 _________ _______ 20 _________ _______ 40 _________ _______ 80 _________ _______ 160 _________ _______

The convergence rate should be for some integer . What is your estimate of ?

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

This method is ``implicit'' because (8) is not a formula for but instead must be solved to find . In the case under consideration, where , it is easy to solve (8) for 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**:- 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. - 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. - Fill in the following table using
`exer10`.nsteps error ratio 10 _________ _______ 20 _________ _______ 40 _________ _______ 80 _________ _______ 160 _________ _______

The convergence rate should be for some integer . What is your estimate of ? - 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.

- Copy

Back to MATH2071 page.

Mike Sussman 2009-01-06