16. Summation and Integration

16.1. Introduction

This unit starts with the methods for approximating definite integrals seen in a calculus course like the (Composite) Midpoint Rule and (Composite) Trapezoidal Rule.

As bonus exercises, we will then work through some more advanced methods as seen in a numerical methods course; the algorithms will be stated below.

We will also learn more about working with modules: how to run tests while developing a module, and how to incorporate demonstrations into a notebook once a module is working.

16.2. Exercises

Exercise 16.1 (The Midpoint Rule, using a for loop)

Write a Python function to be used as

M_n = midpoint0(f, a, b, n)

which returns the approximation of \(I = \int_a^b f(x) \, dx\) given by the [Composite] Midpoint Rule with \(n\) intervals of equal width.

That is,

\[ M_n = h \left( f(a+h/2) + f(a+3h/2) + f(a+5h/2)+ \cdots + f(b-h/2)\right), \quad h = \frac{b-a}{n}. \]

In this first version, accumulate the sum using a for loop.

Test this with \(\displaystyle I_1 := \int_1^e \frac{dx}{x}\) and several choices of \(n\), such as 10 and 100.

Work out the exact value, and use this to display the size of the errors in each approximation.

Exercise 16.2 (The Midpoint Rule again, using function sum)

Express the above Midpoint Rule in summation notation “\(\Sigma\)” and reimplement as

M_n = midpoint1(f, a, b, n)

using the Python function sum to avoid any loops. (Doing it this way give a code that is more concise and closer to mathematical form, so hopefully more readable; it will also probably run faster.)

Exercise 16.3 (The Trapezoidal Rule)

Write a Python function to be used as

T_n = trapezoidal(f, a, b, n)

which returns the approximation of \(I = \int_a^b f(x) \, dx\) given by the [Composite] Trapezoidal Rule with \(n\) intervals of equal width.

That is,

\[ T_n = h \left( \frac{f(a)}{2} + f(a+h) + \cdots + f(b-h) + \frac{f(b)}{2} \right), \]

Use the “sum” method as in Exercise 16.2 above; no loops allowed!

Again test this with \(I_1 = \int_1^e \frac{dx}{x}\).

Exercise 16.4 (Testing and demonstrating within a module file)

Create each of the above three three functions in a module integration. Then

Option A: using Spyder and the IPython command line

Place the test cases below the function definitions within one or more if blocks starting

if __name__ == "__main__":

(Each of those long dashes is a pair of underscores.)

The contents of such an if block are executed when you run the file directly (as if it were a normal Python program file) but are ignored when the module is used with import from another file. Thus, these blocks can be used while developing the module, and later to provide a demonstration of the module’s capabilities.

Option B: using just notebooks and module files

As you may have seen, importing a module a second or subsequent time from another Python file or notebook does not get the updated version of the module’s contents unless you restart the Python kernel. (Python does this to avoid duplicated effort when the same module is mentioned in several import statements in a file.) Thus while revising a module, it is more convenient to treat it like a normal Python file, testing from within by this method rather than by importing elsewhere.

Exercise 16.5 (Bonus Exercise: a more accurate result)

Define another Python function which uses the above functions to compute the “weighted average” \(\displaystyle S_n = \frac{2 M_n + T_n}{3}.\)

Yes, this is the Composite Simpson Rule, so make it

S_n = simpson(f, a, b, n)

Compare its accuracy to that of \(M_n\) and \(T_n\).

Still do testing within the module file, but also add demonstrations to the above notebook once you are sure that this function is working.

Exercise 16.6 (Bonus Exercise: The Romberg Method)

The Romberg method generates a triangular array of approximations \(R_{i,j}, j \leq i\) of an integral \(I = \int_a^b f(x)\; dx\) with the end-of-row values \(R_{0,0}, R_{1,1}, \dots R_{n,n}\) being the main, successively better (usually) approximations.

It starts with the trapezoid rule \(R_{0,0} := T_1 = \frac{f(a) + f(b)}{2}(b-a)\); the basic trapezoid rule.

Then using \(T_{2n} = \frac{T_n + M_n}{2}\), one defines

\[ R_{i, 0} := T_{2_i} = \frac{T_{2^{i}/2} + M_{2^{i}/2}}{2} = \frac{T_{2^{i-1}} + M_{2^{i-1}}}{2} = \frac{R_{i-1, 0} + M_{2^{i-1}}}{2}, \quad i \geq 1 \]

Finally, Richardson extrapolation leads to

\[ R_{i,1} := S_{2^i} = \frac{4 T_{2^i} - T_{2^i/2}}{4 - 1} = \frac{4 R_{i,0} - R_{i-1,0}}{4 - 1} = R_{i,0} + \frac{R_{i,0} - R_{i-1,0}}{4 - 1}, \; i \geq 1 \]

and with further extrapolations to the more general formula

\[ R_{i,j} := \frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1} = R_{i,j-1} + \frac{R_{i,j-1} - R_{i-1,j-1}}{4^i - 1}, \; 1 \leq j \leq i \]

Implement this, using these three formulas plus the above function for the composite midpoint rule.

One natural data structure is a 2D array with unused entries above the main diagonal. However, you might consider how to store this triangular collection of data as a list of lists, succesively of lengths 1, 2 and so on up to \(n\).