Matematik och Datavetenskap, Chalmers Tekniska Högskola och Göteborgs Universitet

ALA-B, 2003, studio 2.1

Do not print this file!! Read it from the computer screen.

How small is the error?

The goal of this exercise is to study how fast the error in your approximate integral decreases when the stepsize h decreases. The theory is found in AMBS chapter 30.

1. The algorithm in The Fundamental Theorem is called the RECTANGLE RULE.

Compute the integral of cos(x) over the interval [0, 1] by hand. Then compute the same integral with your program my_int. Note that the value of the integral is the last component of the vector U. How do you get this component after you have computed the vector U?

Hint: length(U) is the number of components. Put the following commands in a script file studio21.m and execute them by typing studio21 in the matlab command window.

Iexact=??
h=0.1
[x,U]=my_int('cos', [0, 1], 0, h)
I=U(length(U))
Q=abs( I-Iexact )

Does the computation stop exactly at x=1? It is very important that the computation stops exactly at x=1. In order to guarantee that the program does not take one step too much we must change the stopping criterion in my_int as follows:

while x(i) < (b-eps)

Do this now and check that it works.

Create a row vector h of stepsizes: h(1)=2^(-1), ..., h(10)=2^(-10). Do this by putting a for loop in the file studio21.m:

for i=1:10
h(i)=2^(-i);
end

Test the program by typing

>> studio21
>> h

on the command line. Then for each step in the loop compute the integral I(i) with steplength h(i) and the absolute value of the quadrature error Q(i). Make a table by putting I, h and Q as columns in a matrix:

A=[I' h' Q'].

Then change to long format:

format long

and type

A.

Can you see that the numbers in the first column of the table form a Cauchy sequence (=decimal expansion=real number)?
How much does the error decrease each time the step is decreased by a factor 0.5?

Save the table for the documentation.

Plot the error versus the stepsize: plot(h,Q,'o-'). It is difficult to see the details. Therefore it is better to plot the logarithms of h and Q:

plot(log(h),log(Q),'o-')
axis equal

The argument 'o-' puts a ring at each point and draws a straight line between the points.

Note: The natural logarithm ln(x) is called log(x) in English and in matlab.

2. THE MIDPOINT RULE. Now change the program my_int so that it uses the midpoint values of f, f(x(i-1)+h/2), instead of the left end values f(x(i-1)). (This algorithm is called the midpoint rule.) Then repeat all the computations in 1. Use the command hold on so that the new plot comes in the same window as the old one. Use another marker than 'o-' so that you can distinguish between the two curves. For example: 'd-'. Type help plot to see a list of markers.

Save or print the figure for the documentation.

Can you see that the numbers in the first column of table form a Cauchy sequence (=decimal expansion=real number)?
How much does the error decrease each time the step is decreased by a factor 0.5?

3. The theory in Chapter 30 says that Q < Kh^s for some positive number s. For example, according to formula (30.3) we have s=1 for the rectangle rule. If we take the logarithm of this inequality we get

log (Q) < log (K) + s log (h)

so the slope in your diagram is s. Read off s from your diagrams. What do you get for the rectangle rule and for the midpoint rule?

This is an important way to to test a numerical computer program. If you do not see the right rate of convergence (slope s) in the test example, then there is a bug in your program.

4. A function that changes rapidly. (If you have time.) Do all the computations with the function cos(20 x) instead of cos(x).

/stig