- Introduction
- The MATLAB Sum Method
- The Midpoint Method
- Precision
- The Trapezoid Method
- Newton-Cotes Rules
- Gauss Quadrature
- Composite Newton-Cotes Rules
- Composite Gauss Rules
- Error Estimation
- ASSIGNMENT

In numerical analysis, **Quadrature** is the estimation of an integral.
We may want to integrate some function **f(x)** or a set of tabulated
data. The domain might be a finite or infinite interval, it might be a
rectangle or irregular shape, it might be a three dimensional cell.

We consider some simple tests of the precision of a rule, and then describe (simple versions of) the midpoint and trapezoidal rules.

We consider the Newton-Cotes and Gauss-Legendre families of rules, and discuss how to get more accurate approximations by increasing the order of the rule or subdividing the interval. Finally, we take a look at error estimation.

The "MATLAB Sum" method (I made this name up) for estimating an
integral is to evaluate the function at **N** evenly spaced
points in **[A,B]** and multiply by **h**, the spacing
between the points.

Write a routine called *matsum.m*, which looks like

The MATLAB Sum method might look like something like this:function quad = matsum ( f, a, b, n )

h = ( b - a ) / ( n - 1 ) xvec = ? fvec = ? quad = h * sum ( fvec )

**Exercise** - Use your MATLAB sum routine to estimate the integral
of **5*x ^{4}** over

N h "MATLAB Sum" Result Error 2 1 ---------------------- ---------------------- 11 0.1 ---------------------- ---------------------- 101 0.01 ---------------------- ---------------------- 1001 0.001 ---------------------- ----------------------How fast is the error decreasing with

The MATLAB Sum rule is easy to use. It's not very accurate -
it doesn't even integrate **f(x)=1** exactly!
We will keep this method in mind as a basis for comparison.

The *midpoint rule* on a single interval approximates an integral
as the product of the interval width **(b-a)** times the
function value at the midpoint:

h = ( b - a )

x = ( a + b ) / 2

quad = h * f(x)

It's a little tricky to do a composite midpoint rule with MATLAB, because
we need to compute the locations of the midpoints. The easiest
way to compute the midpoints is to compute the **N**
*endpoints* and then average pairs of them, giving a new vector of
**N-1** midpoints.

Write a routine called *midsum.m*, which is declared as:

The coding might look like:function quad = midsum ( f, a, b, n )

h = ( b - a ) / ( n - 1 ) xvec = ? xvec = 0.5 * ( xvec(1:n-1) + xvec(2:n) )replaces ends by midpointsfvec = ? quad = h * sum ( fvec )

**Exercise** - Use your midpoint method to estimate the integral
of **5*x ^{4}** over

N h Midpoint Result Error 2 1 ---------------------- ---------------------- 11 0.1 ---------------------- ---------------------- 101 0.01 ---------------------- ---------------------- 1001 0.001 ---------------------- ----------------------How fast is the error decreasing with

If a quadrature rule can compute exactly the integral of any polynomial
up to some specific degree, we will call this its
*degree of precision*.
Thus a rule that can correctly integrate any cubic, but not quartics,
has precision 3.

To determine the degree of precision of a rule, we might look at the approximations of the integrals from 0 to 1 of the functions

since the exact value of every one of these integrals is 1.p_{0}(x) = 1

p_{1}(x) = 2 * x

p_{2}(x) = 3 * x^{2}

...

p_{k}(x) = ( k + 1 ) * x^{k}

**Exercise** - To study the precision of the midpoint method, use a single
interval ( **N** = 2 ), and estimate the integrals of the test functions
over **[0,1]**. The correct answer should be 1 each time.

f Midpoint Result Error 1 ---------------------- ---------------------- 2 * x ---------------------- ---------------------- 3 * x^2 ---------------------- ---------------------- 4 * x^3 ---------------------- ----------------------When do the results stop being essentially exact? So what is the precision of the midpoint rule?

The *trapezoidal rule* over one interval multiplies the interval
width times the average of the function values at the endpoints:

h = ( b - a )

Q = h * ( f(a) + f(b) )/2

To apply the composite trapezoid rule, we need to generate
**N** points and evaluate the function there. Then we *almost*
repeat the "MATLAB Sum" formula, except that the first and last points
have a factor of (1/2) instead of 1. This amounts to a small correction
to the "MATLAB Sum" formula.

Write a routine called *trapsum.m*, which is declared as:

The coding might look like something like this:function quad = trapsum ( f, a, b, n )

h = ( b - a ) / ( n - 1 ) xvec = ? fvec = ? quad = h * ( sum ( fvec ) - 0.5 * ( fvec(1) + fvec(n) ) )

**Exercise** - Use the trapezoid method to estimate the integral
of **5*x ^{4}** over

N h Trapezoid Result Error 2 1 ---------------------- ---------------------- 11 0.1 ---------------------- ---------------------- 101 0.01 ---------------------- ---------------------- 1001 0.001 ---------------------- ----------------------How fast is the error decreasing with

A Newton-Cotes formula uses the idea that if it's reasonable to approximate a function by its interpolant, then it's reasonable to approximate the integral of that function with the integral of its interpolant. For instance, the trapezoidal rule is the Newton-Cotes formula of order 2, derived by integrating the linear interpolant to a function.

If our Newton-Cotes formula (over a single interval) samples the
function at **N** points, then the interpolant will be exact
for any polynomial of degree **N-1** or less, hence the
integral will be approximated exactly. Therefore, a Newton-Cotes
rule of order **N** has precision at least **N-1**.

We didn't actually prove it can't be better. And in fact,
the rules of order **N** have precision **N** when N is odd.
Hence, you'd expect the midpoint rule to have precision....???

Write a routine called *ncsum.m*, which is declared as:

The coding might look like something like this:function quad = ncsum ( f, a, b, n )

h = ( b - a )The routineWatch out! This is not divided by (n-1)!x = linspace ( a, b, n ) wvec = nc_weight ( n ) fvec = ? quad = h * wvec' * fvec'Transpose both.

**Exercise** - Use Newton-Cotes methods to estimate the integral
of **5*x ^{4}** over

N h Newton-Cotes Result Error 2 1.0 ---------------------- ---------------------- 3 0.5 ---------------------- ---------------------- 4 0.33 ---------------------- ---------------------- 5 0.25 ---------------------- ----------------------

Unfortunately, "more precise" does not guarantee "more accurate"!
If one rule is "more precise" than another, then it can integrate
exactly polynomials of higher degree. So you might think that,
for any function **f(x)**, if we want to approximate the integral
really well, we could just use Newton-Cotes methods of higher and
higher precision.

But the Newton-Cotes rules do **not** have the very desirable
property that the sequence of approximations to the integral is
guaranteed to converge to the integral. It is time to examine
one of our favorite examples.

**Exercise**: consider the Runge function **1/(1+x^2)** on the
interval [-4,4]. The exact integral is **2*atan(4)**. Use rules of
increasing order, each time recording the quadrature error.

N h Newton-Cotes Result Error 3 4.0 ---------------------- ---------------------- 5 2.0 ---------------------- ---------------------- 9 1.0 ---------------------- ---------------------- 17 0.5 ---------------------- ----------------------

Using more precise Newton-Cotes rules means using increasing degree polynomial interpolation with equally spaced nodes. When did we last try this idea, and what was the result then?

Back in the old days, arithmetic was done with a device called a "pencil", and people did everything they could to find the most effective way of reducing the pain of using this device. Gauss, in particular, had to compute many approximate integrals. He knew about the Newton-Cotes formulas, but discovered that he could approximate an integral using far fewer points (saving many, many pencils), if he was willing to replace the equally spaced sampling points by points whose location was determined in another way.

*(Gauss's points are actually the places where the Legendre polynomial
is zero. If you get a chance, plot the Legendre polynomial of degree
5 and then look at the values of the 5 Gauss-Legendre abscissas in
the routine gl_space.m which I am about to tell you about.)*

Write a routine called *glsum.m*, which is declared as:

The coding might look like something like this:function quad = glsum ( f, a, b, n )

h = ( b - a ) / 2The weights are computed by the routineNotice this h is different!wvec = gl_weight ( n ) xvec = gl_space ( a, b, n ) fvec = ? quad = h * wvec' * fvec'Transpose both

**Exercise** - Use Gauss-Legendre methods to estimate the integral
of **5*x ^{4}** over

N Gauss-Legendre Result Error 1 ---------------------- ---------------------- 2 ---------------------- ---------------------- 3 ---------------------- ---------------------- 4 ---------------------- ----------------------How fast is the error decreasing with

Unlike Newton-Cotes, Gauss-Legendre rules **do** have the important
property that, in general, using the sequence of increasing-order
approximations to the integral will converge to the integral.

**Exercise**: Vanquish the beast! Consider the Runge function
**1/(1+x^2)** on the interval **[-4,4]**. The exact integral is
**2*atan(4)**. Use rules of increasing order, each time recording
the quadrature error.

N Gauss-Legendre Result Error 1 ---------------------- ---------------------- 2 ---------------------- ---------------------- 4 ---------------------- ---------------------- 8 ---------------------- ----------------------

As long as we keep the order low, there's nothing wrong with a
Newton-Cotes rule. But if we want a more accurate approximation
to an integral, we can't increase the order of the rule, because
the results go bad. The solution is to use a composite rule.
We break the interval **[A,B]** up into **M** subintervals
and use the rule repeatedly over each subinterval. If the rule
uses **N** points, then we want a total of **M*(N-1)+1**
evenly spaced points. The easiest way to organize this computation
is simply to compute the endpoints of the subintervals, and
call the routine we already wrote.

Write a routine, *comp_ncsum.m*, which looks like

which approximates the integral offunction quad = comp_ncsum ( f, a, b, n, m )

ends = linspace ( a, b, m+1) quad = 0 for i = 1 to m quad = quad + ncsum ( f, ends(i), ends(i+1), n )Note that we are computing

**Exercise**: consider the Runge function **1/(1+x^2)** on the
interval **[-4,4]**. The exact integral is **2*atan(4)**. Try the
composite rules of order **N**, with **M** intervals,
each time recording the error.

N M Composite Newton-Cotes Error 2 4 ---------------------- ---------------------- 2 8 ---------------------- ---------------------- 2 16 ---------------------- ---------------------- 3 4 ---------------------- ---------------------- 3 8 ---------------------- ---------------------- 3 16 ---------------------- ----------------------

A composite Gauss rule can be put together in much the same way as a Newton-Cotes method, except that the spacing of the points is not even.

Write a routine, *comp_glsum.m*, which looks like

which approximates the integral offunction quad = comp_glsum ( f, a, b, n, m )

ends = linspace ( a, b, m+1) quad = 0 for i = 1 to m quad = quad + glsum ( f, ends(i), ends(i+1), n )

**Exercise**: consider the Runge function **1/(1+x^2)** on the
interval [-4,4]. The exact integral is **2*atan(4)**. Try the
composite Gauss-Legendre rules of order **N**, with **M** intervals,
each time recording the quadrature error.

N M Composite Gauss-Legendre Error 2 4 ---------------------- ---------------------- 2 8 ---------------------- ---------------------- 2 16 ---------------------- ---------------------- 3 4 ---------------------- ---------------------- 3 8 ---------------------- ---------------------- 3 16 ---------------------- ----------------------

In quadrature, we know we're only approximating an integral. It is often useful to have an estimate of the error made in approximation. While this estimate can't be taken completely seriously, it is often used as a guide to deciding when an approximation is probably accurate enough, or when it is likely that using more points would produce a better answer.

To approximate the error, we simply make two estimates of the integral, in such a way that we expect the second estimate to be significantly more accurate. We could, for instance compare two rules of different order, or two composite rules, one of which uses twice as many subintervals. In both cases, the difference between the two approximations to the integral is an estimate of the error.

**Exercise**: consider the Runge function **1/(1+x^2)** on the
interval **[-4,4]**. The exact integral is **2*atan(4)**. Compute
the MATLAB Sum estimate for the integral, using the following number
of points. Use the difference of the first and second sums
as an error estimate for the first sum. Use the difference of
the second and third sums to estimate the error in the second sum,
and so on.

N MATLAB Sum Estimated Error True Error 2 ----------------- -------------- ----------- 5 ----------------- -------------- ----------- 17 ----------------- -------------- ----------- 65 ----------------- -------------- ----------- 257 -----------------

Using any method discussed in this lab, estimate the integral
over **[-1,1]** of

Mail me the name of the quadrature method you used, the order0.92 * cosh ( x ) - cos ( x )