Lab 7—Numerical Integration

Purpose

This lesson introduces you to two classic algorithms for evaluating definite integrals. In the process, it also develops your ability to program with for loops, and reinforces your ability to work with vectors and plotting.

Background

Numerical integration is the process of evaluating a definite integral for its numeric value. A number of algorithms for numerical integration have been developed, starting in the days when people had to evaluate integrals by hand, and more recently intended for implementation on computers. The Wikipedia article on numerical integration provides a good sense of the variety of methods that are employed, although it covers too much material in too little space to give a very good understanding of any one method.

You will explore two classical simple numerical integration algorithms in this lesson: Riemann sums and the trapezoid rule. For an introduction to Riemann sums see the Khan Academy videos on “Simple Riemann Approximation Using Rectangles” and “Generalizing a Left Riemann Sum with Equally Spaced Rectangles.” For the trapezoid rule, see “Trapezoidal Approximation of Area under Curve.”

This lesson also introduces a new Matlab feature, for loops. Our textbook discusses for loops in section 13.4, and we will talk about them in class on March 30.

Riemann sums and the trapezoid rule are both algorithms for evaluating definite integrals. Any time you have multiple algorithms for solving the same problem the question of which algorithm is “best” arises. The answer usually depends on context, i.e., sometimes you need a fast solution above all else, sometimes you need one that is very accurate, sometimes you need one that fits in a small amount of computer memory, and very often you need one that strikes a balance between many criteria. In this lesson you will compare Riemann sums and the trapezoid rule to see which one most quickly converges on a close approximation to the actual value of an integral. The comparison will be qualitative rather than quantitative, and, as you will see, not entirely clear cut.

Activity

Do the following.

Step 1

Write two Matlab functions, one of which uses Riemann sums to compute a definite integral and the other of which uses the trapezoid rule.

Both functions should take a function handle for the function being integrated, the bounds of the interval over which to integrate, and the number of sub-intervals to divide that interval into as arguments. Both functions should return the estimated integral as their result. Finally, both functions should draw a graph showing the function they are integrating and the rectangles or trapezoids that they use to approximate the area under it.

Use your functions to visually explore how Riemann sums and the trapezoid rule approximate the area under the function y = x³ between x = 0 and x = 1. Try this for a number of sub-intervals (usually called n in discussions of numerical integration) ranging from 4 or 5 up to 100 or more. Then repeat this exploration for the function y = sinx between x = 0 and x = π. Based on what you see, which integration algorithm do you think will produce the most accurate results with the fewest sub-intervals?

Step 2

To test your insights from Step 1, write a Matlab script that visualizes how close the estimated integral computed by each of the functions from Step 1 is to the actual integral. This script should step through a number of n values, call the Riemann and trapezoid functions for each n, and plot the returned values against n. Your script should also plot a straight line whose Y value is the actual value of the integral being estimated by the Riemann and trapezoid functions.

Even though the instructions above say to call your Riemann and trapezoid functions for various n values, you may actually call modified versions of those functions that don’t do the plotting—you may well call the functions hundreds or thousands of times, and you probably don’t want to see that many plots.

Run your script for the function y = x³ with x ranging from 0 to 1. Also run your script for y = sinx with x ranging from 0 to π. Are the visualizations consistent with what you expected from inspecting plots in Step 1? In general, would you favor Riemann sums over the trapezoid rule for computing definite integrals, or vice versa, or can you not say? What caveats come with your answer?

Extra Credit

The error in integrating f(x) from a to b by the trapezoid rule (i.e., the difference between the actual integral and the value found using the trapezoid rule) is bounded by the inequality

|E| ≤ M(b-a)³ / (12n²)

where E is the error and M is an upper bound on the absolute value of the second derivative of f in the interval [a, b]. As usual, n is the number of intervals used by the trapezoid rule.

For up to 2 points extra credit, extend the plots you draw in Step 2 to also plot the error envelope for the trapezoid rule, i.e., curves showing the actual value of the integral plus and minus the maximum error as calculated by the above relation. How large is this error envelope? Are the values produced by the trapezoid rule generally well inside it, or are they close to one edge of the envelope?

Comment

I will demonstrate the kinds of plots and functions I want you to write for this exercise in class.

Follow-Up

I will grade this exercise in a face-to-face meeting with you. During this meeting I will look at your solution, ask you any questions I have about it, answer questions you have, etc. Please bring your solution on your computer, as that will speed the process along.

Sign up for a meeting via Google calendar. Please make the meeting 15 minutes long, and schedule it to finish before the end of the “Grade By” date above.