Summation and Integration¶

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.

Exercise A: 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 B: the Midpoint Rule, 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 C: 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), \quad h = \frac{b-a}{n}. $$

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

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

Testing and demonstrations within a module file¶

Create each of these three functions in a module integration.

Option A: using Spyder and 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 D: present your work so far with a notebook¶

Once the module is working for the above three methods, create a notebook which imports those functions and runs various examples with them. (That is, copy all the stuff from within the above if __name__ == "__main__": blocks to the notebook.)

The notebook should also describe the mathematical background; in particular, the formulas for all the methods.

Bonus Exercise E: 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.

Bonus Exercise F: 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$.

Final exercise: complete the notebook¶

Add further test cases; one interesting type of example is a periodic function.

Python for Scientific Computing