Riemanns Integral¶

The simplest method for approximating integrals is by summing the area of rectangles that are defined for each subinterval. The width of the rectangle is \(x_{i+1} - x_i = h\), and the height is defined by a function value \(f(x)\) for some \(x\) in the subinterval. An obvious choice for the height is the function value at the left endpoint, \(x_i\), or the right endpoint, \(x_{i+1}\), because these values can be used even if the function itself is not known. This method gives the Riemann Integral approximation, which is

\[\int_a^b f(x) dx \approx \sum_{i = 0}^{n-1} hf(x_i),\]

or

\[\int_a^b f(x) dx \approx \sum_{i = 1}^{n} hf(x_i),\]

depending on whether the left or right endpoint is chosen.

As with numerical differentiation, we want to characterize how the accuracy improves as \(h\) gets small. To determine this characterizing, we first rewrite the integral of \(f(x)\) over an arbitrary subinterval in terms of the Taylor series. The Taylor series of \(f(x)\) around \(a = x_i\) is

\[f(x) = f(x_i) + f^{\prime}(x_i)(x-x_i) + \cdots\]

Thus

\[\int_{x_i}^{x_{i+1}} f(x) dx = \int_{x_i}^{x_{i+1}} (f(x_i) + f^{\prime}(x_i)(x-x_i) + \cdots)\ dx\]

by substitution of the Taylor series for the function. Since the integral distributes, we can rearrange the right side into the following form:

\[\int_{x_i}^{x_{i+1}} f(x_i) dx + \int_{x_i}^{x_{i+1}} f^{\prime}(x_i)(x-x_i)dx + \cdots.\\]

Solving each integral separately results in the approximation

\[\int_{x_i}^{x_{i+1}} f(x) dx = hf(x_i) + \frac{h^2}{2}f^{\prime}(x_i) + O(h^3),\]

which is just

\[\int_{x_i}^{x_{i+1}} f(x) dx = hf(x_i) + O(h^2).\]

Since the \(hf(x_i)\) term is our Riemann integral approximation for a single subinterval, the Riemann integral approximation over a single interval is \(O(h^2)\).

If we sum the \(O(h^2)\) error over the entire Riemann sum, we get \(nO(h^2)\). The relationship between \(n\) and \(h\) is

\[h = \frac{b - a}{n},\]

and so our total error becomes \(\frac{b - a}{h}O(h^2) = O(h)\) over the whole interval. Thus the overall accuracy is \(O(h)\).

The Midpoint Rule takes the rectangle height of the rectangle at each subinterval to be the function value at the midpoint between \(x_i\) and \(x_{i+1}\), which for compactness we denote by \(y_i = \frac{x_{i+1} + x_i}{2}\). The Midpoint Rule says

\[\int_a^b f(x)dx \approx \sum_{i = 0}^{n-1} hf(y_i).\]

Similarly to the Riemann integral, we take the Taylor series of \(f(x)\) around \(y_i\), which is

\[f(x) = f(y_i) + f^{\prime}(y_i)(x - y_i) + \frac{f''(y_i)(x - y_i)^2}{2!} + \cdots\]

Then the integral over a subinterval is

\[\int_{x_i}^{x_{i+1}} f(x) dx= \int_{x_i}^{x_{i+1}} \left(f(y_i) + f^{\prime}(y_i)(x - y_i) + \frac{f''(y_i)(x - y_i)^2}{2!} + \cdots\right) dx,\]

which distributes to

\[\int_{x_i}^{x_{i+1}} f(x) dx= \int_{x_i}^{x_{i+1}} f(y_i)dx + \int_{x_i}^{x_{i+1}} f^{\prime}(y_i)(x - y_i)dx + \int_{x_i}^{x_{i+1}} \frac{f''(y_i)(x - y_i)^2}{2!}dx + \cdots.\]

Recognizing that since \(x_i\) and \(x_{i+1}\) are symmetric around \(y_i\), then \(\int_{x_i}^{x_{i+1}} f^{\prime}(y_i)(x - y_i)dx = 0\). This is true for the integral of \((x - y_i)^p\) for any odd \(p\). For the integral of \((x - y_i)^p\) and with \(p\) even, it suffices to say that \(\int_{x_i}^{x_{i+1}} (x - y_i)^p dx = \int_{-\frac{h}{2}}^{\frac{h}{2}} x^p dx\), which will result in some multiple of \(h^{p+1}\) with no lower order powers of \(h\).

Utilizing these facts reduces the expression for the integral of \(f(x)\) to

\[\int_{x_i}^{x_{i+1}} f(x) dx= hf(y_i) + O(h^3).\]

Since \(hf(y_i)\) is the approximation of the integral over the subinterval, the Midpoint Rule is \(O(h^3)\) for one subinterval, and using similar arguments as for the Riemann Integral, is \(O(h^2)\) over the whole interval. Since the Midpoint Rule requires the same number of calculations as the Riemann Integral, we essentially get an extra order of accuracy for free! However, if \(f(x_i)\) is given in the form of data points, then we will not be able to compute \(f(y_i)\) for this integration scheme.

TRY IT! Use the left Riemann Integral, right Riemann Integral, and Midpoint Rule to approximate \(\int_{0}^{\pi} \text{sin}(x) dx\) wtih 11 evenly spaced grid ponts over the whole interval. Compare this value to the exact value of 2.

import numpy as np

a = 0
b = np.pi
n = 11
h = (b - a) / (n - 1)
x = np.linspace(a, b, n)
f = np.sin(x)

I_riemannL = h * sum(f[:n-1])
err_riemannL = 2 - I_riemannL

I_riemannR = h * sum(f[1::])
err_riemannR = 2 - I_riemannR

I_mid = h * sum(np.sin((x[:n-1] \
        + x[1:])/2))
err_mid = 2 - I_mid

print(I_riemannL)
print(err_riemannL)

print(I_riemannR)
print(err_riemannR)

print(I_mid)
print(err_mid)

9835235375094546
01647646249054535
9835235375094546
01647646249054535
0082484079079745
-0.008248407907974542

< 21.1 Numerical Integration Problem Statement | Contents | 21.3 Trapezoid Rule >

Python Numerical Methods

Riemanns Integral¶