Riemann integrals

import numpy as np
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Left and Right Riemann Integrals¶

Riemann’'s integral is the simplest of the Newton-Cotes formulae with a polynomial of degree 0, i.e. a constant over each subdivision. Variants exist depending on which point in the subdomain $[x_i, x_{i+1}]$ is chosen to evaluate the constant.

If we take the function value at the left of the subdomain $[x_i, x_{i+1}]$ , we obtain the left Riemann integral:

\int_a^b f(x) \,\mathrm{d}x \approx \sum_{i=0}^{n-1} f(x_i) h

(1)

Conversely, if we take the function value at the right limit, we obtain the right Riemann integral:

\int_a^b f(x) \,\mathrm{d}x \approx \sum_{i=1}^{n} f(x_i) h

(2)

numerical_integration_module.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def left_riemann(f, a, b, n):
    h = (b - a) / n
    x_left = np.linspace(a, b - h, n)
    y_left = f(x_left)
    integral_approx = np.sum(y_left * h)

    return {
        "x_left": x_left,
        "y_left": y_left,
        "n": n,
        "integral_approx": integral_approx,
        "a": a,
        "b": b,
        "h": h,
    }

exe.display_source(exe.left_riemann)

def f(x):
    return x**2
result = exe.left_riemann(f, a=0, b=1, n=8)
fig = exe.visualize_left_riemann(f, result)
fig.show()

numerical_integration_module.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def right_riemann(f, a, b, n):
    h = (b - a) / n
    x_right = np.linspace(a + h, b, n)
    y_right = f(x_right)

    integral_approx = np.sum(y_right * h)

    return {
        "x_right": x_right,
        "y_right": y_right,
        "n": n,
        "integral_approx": integral_approx,
        "a": a,
        "b": b,
        "h": h,
    }

exe.display_source(exe.right_riemann)

def f(x):
    return x**2
result = exe.right_riemann(f, a=0, b=1, n=8)
fig = exe.visualize_right_riemann(f, result)
fig.show()

Error Analysis of One-Sided Riemann Sums¶

To analyze the error, consider integrating the Taylor expansion of the left integral (the right being trivially similar):

\begin{aligned} f(x) &= f(x_i) + f'(x_i) h + O(h^2) \\ \int_{x_i}^{x_{i+1}} f(x) \,\mathrm{d}x &= \int_{x_i}^{x_{i+1}} [f(x_i) + f'(x_i) h + O(h^2)] \,\mathrm{d}x \\ &= f(x_i) h + f'(x_i) \frac{h^2}{2} + O(h^3) \\ &= f(x_i) h + O(h^2) \end{aligned}

(3)

The truncation error per element is $O(h^2)$ . Since there are $n = (b-a)/h$ elements, the global composite error scales as:

E_{composite} = \sum_{i=0}^{n-1} O(h^2) = n \cdot O(h^2) = O(h)

(4)

making left and right Riemann sums first-order accurate.

The Midpoint Riemann Integral¶

Rather than favor the left or right subdomain limit, let’'s take the midpoint $y_i = \frac{x_{i+1} + x_i}{2}$ . The Midpoint Rule states:

\int_a^b f(x) \,\mathrm{d}x \approx \sum_{i=0}^{n-1} f(y_i) h

(5)

numerical_integration_module.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def midpoint_riemann(f, a, b, n):
    h = (b - a) / n
    x_mid = np.linspace(a + h / 2, b - h / 2, n)
    y_mid = f(x_mid)

    integral_approx = np.sum(y_mid * h)

    return {
        "x_mid": x_mid,
        "y_mid": y_mid,
        "n": n,
        "integral_approx": integral_approx,
        "a": a,
        "b": b,
        "h": h,
    }

exe.display_source(exe.midpoint_riemann)

def f(x):
    return x**2
result = exe.midpoint_riemann(f, a=0, b=1, n=8)
fig = exe.visualize_midpoint_riemann(f, result)
fig.show()

Example: Midpoint Approximation of $\sin(x)$ ¶

Approximate $\int_0^\pi \sin(x) \,\mathrm{d}x = 2$ using the Midpoint rule across 8 subdivisions:

result = exe.midpoint_riemann(np.sin, a=0, b=np.pi, n=8)
fig = exe.visualize_midpoint_riemann(np.sin, result)
fig.show()

Error Analysis: Second-Order Midpoint Accuracy¶

The error can be deduced considering the Taylor series of $f(x)$ around the midpoint $y_i$ :

\begin{aligned} f(x) &= f(y_i) + f'(y_i)(x - y_i) + \frac{f''(y_i)}{2}(x - y_i)^2 + O(h^3) \\ \int_{x_i}^{x_{i+1}} f(x) \,\mathrm{d}x &= \int_{y_i - h/2}^{y_i + h/2} \left[f(y_i) + f'(y_i)(x - y_i) + \frac{f''(y_i)}{2}(x - y_i)^2 + O(h^3)\right] \,\mathrm{d}x \end{aligned}

(6)

Because the integration limits are symmetric around $y_i$ , the odd term vanishes: $\int_{-h/2}^{h/2} (x - y_i) \,\mathrm{d}x = 0$ .

Integrating the quadratic term yields:

\int_{x_i}^{x_{i+1}} f(x) \,\mathrm{d}x = f(y_i) h + \frac{f''(y_i)}{24} h^3 + O(h^5)

(7)

Summing across all $n$ subdivisions gives a composite error of $O(h^2)$ , a full order of accuracy higher than one-sided Riemann sums!

Note: If you only have discrete data measured at points $x_i$ , you cannot evaluate $f(y_i)$ without interpolating between samples.

Left and Right Riemann Integrals¶

Error Analysis of One-Sided Riemann Sums¶

The Midpoint Riemann Integral¶

Example: Midpoint Approximation of sin⁡(x)\sin(x)sin(x)¶

Error Analysis: Second-Order Midpoint Accuracy¶

Example: Midpoint Approximation of $\sin(x)$ ¶