Riemann Sums

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


A Riemann sum of a function $f(x)$ over a partition

$$ x_0 = a < x_1 < \cdots < x_{N-1} < x_N = b $$

is a sum of the form

$$ \sum_{i=1}^N f(x_i^ * ) (x_i - x_{i-1}) \ \ , \ x_i^* \in [x_{i-1},x_i] $$

where each value $x_i^* \in [x_{i-1},x_i]$ in each subinterval is arbitrary.

Riemann sums are important because they provide an easy way to approximate a definite integral

$$ \int_a^b f(x) \, dx \approx \sum_{i=1}^N f(x_i^ * ) (x_i - x_{i-1}) \ \ , \ x_i^* \in [x_{i-1},x_i] $$

Notice that the product $f(x_i^ * ) (x_i - x_{i-1})$ for each $i$ is the area of a rectangle of height $f(x_i^ * )$ and width $x_i - x_{i-1}$. We can think of a Riemann sum as the area of $N$ rectangles with heights determined by the graph of $y=f(x)$.

The value $x_i^*$ chosen in each subinterval is arbitrary however there are certain obvious choices:

  • A left Riemann sum is when each $x_i^* = x_{i-1}$ is the left endpoint of the subinterval $[x_{i-1},x_i]$
  • A right Riemann sum is when each $x_i^* = x_i$ is the right endpoint of the subinterval $[x_{i-1},x_i]$
  • A midpoint Riemann sum is when each $x_i^* = (x_{i-1} + x_i)/2$ is the midpoint of the subinterval $[x_{i-1},x_i]$

Let's visualize rectangles in the left, right and midpoint Riemann sums for the function

$$ f(x) = \frac{1}{1 + x^2} $$

over the interval $[0,5]$ with a partition of size $N=10$.

f = lambda x : 1/(1+x**2)
a = 0; b = 5; N = 10
n = 10 # Use n*N+1 points to plot the function smoothly

x = np.linspace(a,b,N+1)
y = f(x)

X = np.linspace(a,b,n*N+1)
Y = f(X)


x_left = x[:-1] # Left endpoints
y_left = y[:-1]
plt.title('Left Riemann Sum, N = {}'.format(N))

x_mid = (x[:-1] + x[1:])/2 # Midpoints
y_mid = f(x_mid)
plt.title('Midpoint Riemann Sum, N = {}'.format(N))

x_right = x[1:] # Left endpoints
y_right = y[1:]
plt.title('Right Riemann Sum, N = {}'.format(N))



Notice that when the function $f(x)$ is decreasing on $[a,b]$ the left endpoints give an overestimate of the integral $\int_a^b f(x) dx$ and right endpoints give an underestimate. The opposite is true is when the function is increasing.

Let's compute the value of each of the Riemann sums:

dx = (b-a)/N
x_left = np.linspace(a,b-dx,N)
x_midpoint = np.linspace(dx/2,b - dx/2,N)
x_right = np.linspace(dx,b,N)

print("Partition with",N,"subintervals.")
left_riemann_sum = np.sum(f(x_left) * dx)
print("Left Riemann Sum:",left_riemann_sum)

midpoint_riemann_sum = np.sum(f(x_midpoint) * dx)
print("Midpoint Riemann Sum:",midpoint_riemann_sum)

right_riemann_sum = np.sum(f(x_right) * dx)
print("Right Riemann Sum:",right_riemann_sum)
Partition with 10 subintervals.
Left Riemann Sum: 1.613488696614725
Midpoint Riemann Sum: 1.373543428316664
Right Riemann Sum: 1.1327194658454942

We know the exact value

$$ \int_0^5 \frac{1}{1 + x^2} dx = \arctan(5) $$

and we can compare the Riemann sums to the value

I = np.arctan(5)
print("Left Riemann Sum Error:",np.abs(left_riemann_sum - I))
print("Midpoint Riemann Sum:",np.abs(midpoint_riemann_sum - I))
print("Right Riemann Sum:",np.abs(right_riemann_sum - I))
Left Riemann Sum Error: 0.24008792966970915
Midpoint Riemann Sum: 0.00014266137164820059
Right Riemann Sum: 0.24068130109952168

Error Formulas

A Riemann sum is an approximation of a definite integral. A natural question arises: how good of an approximation is a Riemann sum?

Theorem. Let $L_N(f)$ denote the left Riemann sum

$$ L_N(f) = \sum_{i=1}^N f(x_{i-1} ) \Delta x $$

where $\Delta x = (b-a)/N$ and $x_i = a + i \Delta x$. The error bound is

$$ E_N^{L}(f) = \left| \ \int_a^b f(x) \ dx - L_N(f) \ \right| \leq \frac{(b-a)^2}{2 N} K_1 $$

where $\left| \, f'(x) \, \right| \leq K_1$ for all $x \in [a,b]$.

Theorem. Let $R_N(f)$ denote the right Riemann sum

$$ R_N(f) = \sum_{i=1}^N f(x_{i} ) \Delta x $$

where $\Delta x = (b-a)/N$ and $x_i = a + i \Delta x$. The error bound is

$$ E_N^{R}(f) = \left| \ \int_a^b f(x) \ dx - R_N(f) \ \right| \leq \frac{(b-a)^2}{2 N} K_1 $$

where $\left| \, f'(x) \, \right| \leq K_1$ for all $x \in [a,b]$.

Theorem. Let $M_N(f)$ denote the midpoint Riemann sum

$$ M_N(f) = \sum_{i=1}^N f(x_i^* ) \Delta x $$

where $\Delta x = (b-a)/N$ and $x_i^* = (x_{i-1} + x_i)/2$ for $x_i = a + i \Delta x$. The error bound is

$$ E_N^{M}(f) = \left| \ \int_a^b f(x) \ dx - M_N(f) \ \right| \leq \frac{(b-a)^3}{24 N^2} K_2 $$

where $\left| \, f''(x) \, \right| \leq K_2$ for all $x \in [a,b]$.

There are several points to notice:

  • Left and right Riemann sums have the same error bound which depends on the first derivative $f'(x)$.
  • Midpoint Riemann sum error bound depends on the second derivative $f''(x)$.
  • We expect the midpoint Riemann sum to give a better approximation as $N \to \infty$ since its error bound is inversely proportional to $N^2$ but left/right Riemann sum error bound is inversely proportional only to $N$.


Let's write a function called riemann_sum which takes 5 input parameters f, a, b, N and method and returns the Riemann sum

$$ \sum_{i=1}^N f(x_i^*) \Delta x $$

where $\Delta x = (b-a)/N$ and $x_i = a + i\Delta x$ defines a partition with $N$ subintervals of equal length, and method determines whether we use left endpoints, right endpoints or midpoints (with midpoints as the default method).

def riemann_sum(f,a,b,N,method='midpoint'):
    '''Compute the Riemann sum of f(x) over the interval [a,b].

    f : function
        Vectorized function of one variable
    a , b : numbers
        Endpoints of the interval [a,b]
    N : integer
        Number of subintervals of equal length in the partition of [a,b]
    method : string
        Determines the kind of Riemann sum:
        right : Riemann sum using right endpoints
        left : Riemann sum using left endpoints
        midpoint (default) : Riemann sum using midpoints

        Approximation of the integral given by the Riemann sum.
    dx = (b - a)/N
    x = np.linspace(a,b,N+1)

    if method == 'left':
        x_left = x[:-1]
        return np.sum(f(x_left)*dx)
    elif method == 'right':
        x_right = x[1:]
        return np.sum(f(x_right)*dx)
    elif method == 'midpoint':
        x_mid = (x[:-1] + x[1:])/2
        return np.sum(f(x_mid)*dx)
        raise ValueError("Method must be 'left', 'right' or 'midpoint'.")

Let's test our function with inputs where we know exactly what the output should be. For example, we know

$$ \int_0^{\pi/2} \sin(x) \, dx = 1 $$

and, since $\sin(x)$ is increasing on $[0,\pi/2]$, we know that left endpoints will give an under-estimate, and right endpoints will give an over-estimate.


We also know that $\int_0^1 x \, dx = 1/2$ and midpoint should give the result exactly for any $N$:

riemann_sum(lambda x : x,0,1,1)


Approximate Pi

Find a value $N$ which guarantees the right Riemann sum of $f(x)=\frac{4}{1 + x^2}$ over $[0,1]$ is within $10^{-5}$ of the exact value

$$ \int_0^1 \frac{4}{1 + x^2} dx = \pi $$


$$ f'(x) = -\frac{8x}{(1+x^2)^2} $$

Use brute force optimization to find a bound on $\left| f'(x) \right|$ on $[0,1]$:

x = np.linspace(0,1,1000)
y = np.abs(-8*x/(1 + x**2)**2)

Therefore, $\left| f'(x) \right| \leq 2.6$ for $x \in [0,1]$. Use the error bound

$$ \frac{(b-a)^2}{2 N} K_1 \leq 10^{-5} \ \Rightarrow \ \frac{1.3}{N} \leq 10^{-5} \ \Rightarrow \ 130000 \leq N $$

Let's compute the right Riemann sum for $N=130000$:

approximation = riemann_sum(lambda x : 4/(1 + x**2),0,1,130000,method='right')

Verify the accuracy of the approximation

np.abs(approximation - np.pi) < 10**(-5)

Approximate ln(2)

Find a value $N$ which guarantees the midpoint Riemann sum of $f(x)=\frac{1}{x}$ over $[1,2]$ is within $10^{-8}$ of the exact value

$$ \int_1^2 \frac{1}{x} dx = \ln(2) $$


$$ f''(x) = \frac{2}{x^3} $$

Since $f''(x)$ is decreasing for all $x>0$ we have $\left| \, f''(x) \, \right| \leq 2$ for all $x \in [1,2]$. Use the error bound:

$$ \frac{(b-a)^3}{24 N^2} K_2 \leq 10^{-8} \ \Rightarrow \ \frac{1}{12 N^2} \leq 10^{-8} \ \Rightarrow \frac{10^4}{\sqrt{12}} \leq N $$

10**4 / np.sqrt(12)

Therefore a partition of size $N=2887$ guarantees the desired accuracy:

approximation = riemann_sum(lambda x : 1/x,1,2,2887,method='midpoint')

Verify the accuracy of the approximation:

np.abs(approximation - np.log(2)) < 10**(-8)


  1. Consider the integral

    $$ \int_1^2 \frac{dx}{1+x^3} $$

    Without plotting the functions $f(x)$, $f'(x)$ or $f''(x)$, find a value $N$ such that $E_N^R(f) \leq 10^{-5}$ given

    $$ f(x) = \frac{1}{1 + x^3} \ , \ f'(x) = -\frac{3 x^{2}}{\left(x^{3} + 1\right)^{2}} \ , \ f''(x) = \frac{6 x \left(2 x^{3} - 1\right)}{\left(x + 1\right)^{3} \left(x^{2} - x + 1\right)^{3}} \\ $$

  2. Plot the function $f''(x)$ from the previous question on the interval $[1,2]$ and find a value $N$ such that $E_N^M(f) \leq 10^{-5}$ for the integral in the previous question.

  3. Let $f(x) = x^x$ and note that

    $$ f'(x) = x^{x} \left(\log{\left(x \right)} + 1\right) \ , \ f''(x) = x^{x} \left(\log{\left(x \right)} + 1\right)^{2} + x^{x-1} $$

    Plot the function $f''(x)$ and use that information to compute $T_N(f)$ for the integral

    $$ \int_1^2 x^x \, dx $$

    such that $E_N^T(f) \leq 10^{-3}$.