Code
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from mechanicskit import fplot
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.dpi'] = 1207.6 Numerical Integration
Introduction
When we evaluate definite integrals, we naturally prefer analytical solutions: finding an antiderivative and evaluating it at the boundaries. However, many integrals that arise in engineering practice cannot be evaluated analytically, requiring numerical approximation instead. In this chapter, we examine several numerical integration schemes that balance accuracy against computational cost, allowing us to approximate integrals with controlled error.
We begin by exploring a straightforward approach: constructing integration rules that are easy to evaluate. This investigation reveals both the power and limitations of simple schemes, motivating more sophisticated methods.
Consider the definite integral of a function \(f(x)\) over an arbitrary interval:
\[ \int_a^b f(x) d x \]
To develop numerical integration rules that apply to any interval \([a,b]\), we face a fundamental challenge: designing rules that work independently of the specific limits. Rather than creating separate integration formulas for each possible interval, we use coordinate transformation to transform any interval \([a,b]\) to a standard reference interval. This approach separates the geometry of the integration domain from the integration rule itself, enabling us to develop general-purpose methods once and reuse them for arbitrary intervals.
We map any range \(x \in[a, b]\) to the standard range \(\xi \in[0,1]\) using the linear transformation \(x(\xi):=(1-\xi) a+\xi b\). The derivative of this mapping gives \(\frac{d x}{d \xi}=b-a\), which appears as the Jacobian in the change of variables:
\[ \boxed{\int_a^b f(x) d x=\int_0^1 f(x(\xi)) \frac{d x}{d \xi} d \xi=(b-a) \int_0^1 f(x(\xi)) d \xi} \]
This transformation reveals that integration over \([a,b]\) can be computed by scaling the result from the standard interval \([0,1]\) by the interval length \(b-a\).
Example:
Let’s say we want to integrate some polynomial of order \(n\)
\[ \int_a^b x^n d x=(b-a) \int_0^1 x(\xi)^n d \xi=(b-a) \int_0^1((1-\xi) a+\xi b)^n d \xi=-\frac{a^{n+1}-b^{n+1}}{n+1} \]
This gives a rule for integrating polynomials of order \(n\) over any range \([a,b]\), which can be implemented in code as a lookup table. For example,
\[ \int_2^5 2 x-4 x^2 d x=-2 \frac{a^{1+1}-b^{1+1}}{1+1}--4 \frac{a^{2+1}-b^{2+1}}{2+1}=-135 \]
But in general it is awkward to use this rule to integrate e.g., \(\int_{x_1}^{x_2}(1-x)(1-x) d x\) or \(\int_{x_1}^{x_2} \varphi^T \varphi d x\) without additional preparation which makes it cumbersome to implement in a computer program and run efficiently.
The same goes for the integration rulse such as
\[ \int_{-1}^1 x^n d x= \begin{cases}\frac{2}{n+1} & \text { if } n \text { is even } \\ 0 & \text { if } n \text { is odd }\end{cases} \]
and
\[ \int_0^1 x^n d x=\frac{1}{n+1} \]
We want a more general approach to numerical integration
sympy.integrate
The sympy library in Python provides a powerful tool for symbolic mathematics, including integration. The integrate function can be used to compute definite and indefinite integrals of functions. This is great for analytical solutions!
Let’s consider a polynomial function \(f(x)=1+3 x^2-x^3\) for \(x \in[0,3]\)
x = sp.Symbol('x')
f_expr = 1 + 3*x**2 - x**3
x1 = 0
x2 = 3Code
fig, ax = plt.subplots(figsize=(6, 4))
fplot(f_expr, range=(x1, x2), color='blue', linewidth=2,
label=f'$f(x) = {sp.latex(f_expr)}$')
f_numeric = sp.lambdify(x, f_expr, 'numpy')
x_fill = np.linspace(x1, x2, 500)
y_fill = f_numeric(x_fill)
ax.fill_between(x_fill, y_fill, color='lightblue', alpha=0.5, label=f'Area from {x1} to {x2}')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title('Plot of $f(x)$ and Shaded Integral Area')
ax.axhline(0, color='black', linewidth=0.5) # Add x-axis
ax.axvline(0, color='black', linewidth=0.5) # Add y-axis
ax.grid(True, linestyle='--', alpha=0.7)
ax.legend(loc='upper left')
ax.set_xlim(x1, x2)
ax.set_ylim(0, max(y_fill)*1.1)
plt.show()
We can integrate this analytically
integral_value = sp.integrate(f_expr, (x, x1, x2))
integral_value\(\displaystyle \frac{39}{4}\)
integral_value.evalf()\(\displaystyle 9.75\)
Non-trivial integrals
This works very well even for functions that do not have simple symbolic integrals.
f_expr = 1/sp.sqrt(sp.sin(x))
x1 = 0
x2 = 1Code
fig, ax = plt.subplots(figsize=(8, 5))
x1 = 0.01
fplot(f_expr, range=(x1, x2), color='blue', linewidth=2,
label=f'$f(x) = {sp.latex(f_expr)}$')
f_numeric = sp.lambdify(x, f_expr, 'numpy')
x_fill = np.linspace(x1, x2, 500)
y_fill = f_numeric(x_fill)
ax.fill_between(x_fill, y_fill, color='lightblue', alpha=0.5, label=f'Area from {x1} to {x2}')
ax.set_xlim(x1, x2)
ax.set_ylim(0, max(y_fill)*1.1)
ax.legend(loc='upper right')
ax.grid(False)
plt.show()
Note that we avoid the singularity by shifting the lower limit away from zero.
sp.integrate(f_expr, (x, x1, x2))\(\displaystyle \int\limits_{0.01}^{1} \frac{1}{\sqrt{\sin{\left(x \right)}}}\, dx\)
x1 = 0
sp.integrate(f_expr, (x, x1, x2)).evalf()\(\displaystyle 2.03480531920757\)
Another example is the function \(f(x)=e^{\cos(x)}\) for \(x \in[0,3]\)
f_expr = sp.exp(sp.cos(x))
x1 = 0; x2 = 3Code
fig, ax = plt.subplots(figsize=(8, 5))
fplot(f_expr, range=(x1, x2), color='blue', linewidth=2,
label=f'$f(x) = {sp.latex(f_expr)}$')
f_numeric = sp.lambdify(x, f_expr, 'numpy')
x_fill = np.linspace(x1, x2, 500)
y_fill = f_numeric(x_fill)
ax.fill_between(x_fill, y_fill, color='lightblue', alpha=0.5, label=f'Area from {x1} to {x2}')
ax.set_xlim(x1, x2)
ax.set_ylim(0, max(y_fill)*1.1)
ax.legend(loc='upper right')
ax.grid(False)
plt.show()
sp.integrate(f_expr, (x, x1, x2))\(\displaystyle \int\limits_{0}^{3} e^{\cos{\left(x \right)}}\, dx\)
Notice how integrate has failed to find a symbolic integral. It took a long time, and actually fell back on sympy.Integral.
sympy.Integral
sympy.integrate can fail when the integral is too complex. In such cases, we can use the sympy.Integral class to represent the integral symbolically without evaluating it immediately. This allows us to manipulate the integral expression further or evaluate it later when more information is available.
sp.Integral(f_expr, (x, 0, 3))\(\displaystyle \int\limits_{0}^{3} \left(2 x \sin{\left(12 x \right)} + 3\right)\, dx\)
We see here that it simply returns an unevaluated integral object.
We can however evaluate it numerically using the evalf method.
sp.Integral(f_expr, (x, 0, 3)).evalf()\(\displaystyle 9.05020713851588\)
There are two options here: use sp.integrate(f, (x, x1, x2)).evalf() or use sympy.Integral(f, (x, x1, x2)).evalf().
sp.integrate(f_expr, (x, 0, 3)).evalf()\(\displaystyle 9.05020713851588\)
sp.Integral(f_expr, (x, 0,3 )).evalf()\(\displaystyle 9.05020713851588\)
sympy.integrate vs sympy.Integral
The outcome is the same, but the second option is often faster since it does not try to find a symbolic expression for the integral. If sp.integrate fails it will fall back to sp.Integral. Relying on this fallback is not good practice.
For reliably computing the numerical value of an integral, always use sympy.Integral(...).evalf(). It is the direct, explicit, and robust method.
Use sympy.integrate() only when you need the symbolic antiderivative itself.
For engineering applications requiring only occasional numerical integration, symbolic methods like sympy.integrate or adaptive quadrature suffice. However, finite element analysis and similar computational methods require evaluating millions of integrals: element stiffness matrices, force vectors, and residuals at every iteration of every time step. In these contexts, computational efficiency becomes paramount. We therefore examine numerical integration schemes designed for efficient implementation and repeated evaluation, focusing on methods that balance accuracy against cost in a predictable, controllable manner.
Riemann sum
We begin our examination of numerical integration with the most fundamental approach: the Riemann sum. This method approximates the integral by dividing the interval \([a,b]\) into \(N\) subintervals of equal width \(\Delta x = (b-a)/N\) and evaluating the function at one point within each subinterval. Using the left endpoint of each subinterval gives
\[ \int_a^b f(x) \approx \frac{b-a}{N} \sum_{i=1}^N f\left(x_i\right) \]
where \(x_i\) denotes the evaluation point in the \(i\)-th subinterval. This formula treats the integral as the limit of a sum of rectangular areas, each with width \(\Delta x\) and height \(f(x_i)\). While conceptually simple, we shall see that this approach converges slowly, motivating more sophisticated methods.
Code
x = sp.Symbol('x')
f_expr = 1 + 3*x**2 - x**3
x1 = 0; x2 = 3
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 50 # Increase if needed
from matplotlib.colors import to_rgba
# Define N values
N_values = [1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 150, 200]
N_values = list(range(2, 20, 1)) + list(range(20, 201, 5))
def animate(frame_idx):
"""Animation function for each frame"""
ax.clear()
N = N_values[frame_idx]
# Compute Riemann sum data
f_numeric = sp.lambdify(x, f_expr, 'numpy')
dx = (x2 - x1) / N
x_points = np.linspace(x1 + dx/2, x2 - dx/2, N) # Center points for bars
y_points = f_numeric(x_points - dx/2) # Left Riemann sum (use left edge)
riemann_sum = dx * np.sum(y_points)
# Compute exact integral
exact_integral = sp.Integral(f_expr, (x, x1, x2)).evalf()
error = abs(float(exact_integral) - riemann_sum)
# Determine edge opacity based on N
if N <= 10:
edge_alpha = 0.5
edge_color = to_rgba('black', alpha=edge_alpha)
linewidth = 1
elif N < 75:
edge_alpha = 0.5 - ((N - 10) / (75 - 10))*0.5
edge_color = to_rgba('black', alpha=edge_alpha)
linewidth = 1
else:
edge_color = 'lightcoral'
# edge_color = to_rgba('black', alpha=0)
linewidth = 1
# Plot Riemann sum as bars
ax.bar(x_points, y_points, width=dx, align='center',
facecolor='lightcoral',
# alpha=1,
edgecolor=edge_color, linewidth=linewidth,
label='Riemann Sum')
# Plot the function on top
x_plot = np.linspace(x1, x2, 500)
y_plot = f_numeric(x_plot)
ax.plot(x_plot, y_plot, color='blue', linewidth=2,
label=f'$f(x) = {sp.latex(f_expr)}$', zorder=5)
# Set labels and limits
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$f(x)$', fontsize=12)
ax.set_xlim(x1, x2)
y_max = max(f_numeric(np.linspace(x1, x2, 500)))
ax.set_ylim(0, y_max * 1.2)
ax.grid(True, linestyle='--', alpha=0.7)
ax.axhline(0, color='black', linewidth=0.5)
ax.axvline(0, color='black', linewidth=0.5)
# Add text box with results
textstr = f'$N = {N}$\n'
textstr += f'Exact Integral: ${exact_integral:.6f}$\n'
textstr += f'Riemann Sum: ${riemann_sum:.6f}$\n'
textstr += f'Error: ${error:.6f}$'
props = dict(boxstyle='round', facecolor='wheat', alpha=0.8)
ax.text(0.02, 0.98, textstr, transform=ax.transAxes, fontsize=11,
verticalalignment='top', bbox=props)
ax.legend(loc='upper right')
# Create figure and axis (close to prevent static display)
fig, ax = plt.subplots(figsize=(6, 4))
plt.close(fig) # Prevent static image from showing
# Create animation
anim = FuncAnimation(fig, animate, frames=len(N_values),
interval=800, repeat=True)
# Display as interactive HTML
# HTML(anim.to_jshtml())
# Get the HTML and wrap it with responsive CSS
anim_html = anim.to_jshtml()
# Convert to responsive HTML
import mechanicskit as mk
mk.to_responsive_html(anim, container_id='riemann_sum_animation')To assess the Riemann sum’s accuracy quantitatively, we measure the absolute integration error as the difference between the Riemann approximation and the exact integral value:
\[ \epsilon=\left|I_{\text {Riemann }}-I_{\text {exact }}\right| \]
For our polynomial example with \(N=200\) function evaluations, we obtain:
N = 200
f_numeric = sp.lambdify(x, f_expr, 'numpy')
dx = (x2 - x1) / N
x_points = np.linspace(x1, x2, N)
y_points = f_numeric(x_points)
I_riemann = dx * np.sum(y_points)
I_exact = sp.Integral(f_expr, (x, x1, x2)).evalf()error_riemann = abs( I_exact - I_riemann)
error_riemann\(\displaystyle 0.03391959798995\)
Note that the absolute integration error is still significantly large after 200 function evaluations!
This slow convergence reveals the central challenge in numerical integration: developing schemes that reduce error more rapidly as we increase the number of function evaluations. Since computational cost scales directly with the number of evaluations \(N\), methods that achieve a given accuracy with fewer evaluations provide substantial practical advantage. The trapezoidal rule, which we examine next, demonstrates how a simple refinement can dramatically improve convergence.
Trapezoidal rule
The Riemann sum approximates the function as piecewise constant over each subinterval. A natural improvement uses piecewise linear approximation instead, connecting function values at adjacent points with straight line segments. This geometric insight leads to the trapezoidal rule.
Over a single interval \([a,b]\), we approximate the function linearly using shape functions that interpolate between the endpoint values: \(f(x) \approx \varphi_1(x) f(a)+\varphi_2(x) f(b)\) where \(\varphi_1(x)=(b-x)/(b-a)\) and \(\varphi_2(x)=(x-a)/(b-a)\). Integrating this linear approximation over the interval gives
\[ I \approx \frac{f(a)}{b-a} \int_a^b b-x d x+\frac{f(b)}{b-a} \int_a^b x-a d x=\frac{f(a)}{b-a}\left[b x-\frac{1}{2} x^2\right]_a^b+\frac{f(b)}{b-a}\left[\frac{1}{2} x^2-a x\right]_a^b \]
\[ =\frac{f(a)}{b-a}\left(b^2-\frac{1}{2} b^2-a b-\frac{1}{2} a^2\right)+\frac{f(b)}{b-a}\left(\frac{1}{2} b^2-a b-\frac{1}{2} a^2-a^2\right) \\ \]
\[ =\frac{f(a)}{b-a} \frac{1}{2}(b-a)^2+\frac{f(b)}{b-a} \frac{1}{2}(b-a)^2 \\ \]
\[ =\frac{1}{2}(b-a)(f(a)+f(b)) \]
This is called the trapezoidal rule. And for a chained approximation with \(N+1\) evenly spaced points the approximation becomes:
\[ \boxed{\int_a^b f(x) d x \approx \frac{b-a}{2 N} \sum_{i=1}^N\left(f\left(x_i\right)+f\left(x_{i+1}\right)\right)} \]
For example for the function above using \(N=200\) we have:
f_expr = 1 + 3*x**2 - x**3
x1 = 0; x2 = 3
I_exact = sp.Integral(f_expr, (x, x1, x2)).evalf()
N = 200
f_numeric = sp.lambdify(x, f_expr, 'numpy')
x_points = np.linspace(x1, x2, N+1)
x_i = x_points[0:-1]
x_i1 = x_points[1:]
I_trapz = np.sum(f_numeric(x_i) + f_numeric(x_i1)) * (x2-x1)/(2*N)
I_trapznp.float64(9.74983125)error_trapz = abs(I_exact - I_trapz)
error_trapz\(\displaystyle 0.000168750000000273\)
Trapezoidal rule gives much better accuracy than the Riemann sum for the same number of function evaluations! There is a built-in function in numpy for the trapezoidal rule: numpy.trapezoid
np.trapezoid(f_numeric(x_points), x_points)np.float64(9.74983125)error_riemann/error_trapz\(\displaystyle 201.005025125305\)
We see a 200x improvement in accuracy for the same number of function evaluations!
The dramatic improvement in accuracy reveals different error scaling laws. Riemann sum error decreases as \(O(h^2)\) where \(h=(b-a)/N\) is the interval width, meaning that halving the interval reduces error by a factor of 4. Trapezoidal rule error decreases as \(O(h^3)\) for smooth functions, giving eightfold error reduction when halving the interval. For our problem with \(N=200\) subintervals, the ratio of error scaling terms \((1/N)^2 / (1/N)^3 = N\) explains the observed 200-fold improvement. This \(O(h^3)\) convergence makes the trapezoidal rule practical for many applications, though Gaussian quadrature achieves even faster convergence through optimal node placement.
Code
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 50 # Increase if needed
from matplotlib.colors import to_rgba
f_expr = 1 + 3*x**2 - x**3
x1 = 0; x2 = 3
f_numeric = sp.lambdify(x, f_expr, 'numpy')
I_exact = sp.Integral(f_expr, (x, x1, x2)).evalf()
N_values = list(range(4, 20, 1)) + list(range(20, 201, 5))
# N_values = [4,8,15]
def makePlot(N):
x_points = np.linspace(x1, x2, N+1)
dx = (x2 - x1) / (N+1)
y_points = f_numeric(x_points)
I_riemann = dx * np.sum(y_points)
error_riemann = abs( I_exact - I_riemann)
x_i = x_points[0:-1]
x_i1 = x_points[1:]
I_trapz = np.sum(f_numeric(x_i) + f_numeric(x_i1)) * (x2-x1)/(2*N)
error_trapz = abs(I_exact - I_trapz)
if N <= 10:
edge_alpha = 0.5
edge_color = to_rgba('black', alpha=edge_alpha)
elif N < 50:
edge_alpha = 0.5 - ((N - 10) / (50 - 10))*0.5
edge_color = to_rgba('black', alpha=edge_alpha)
else:
edge_color = to_rgba('cyan', alpha=0)
# fig, ax = plt.subplots(figsize=(6, 4))
ax.grid(False)
for iel in range(N):
xs = [x_i[iel], x_i[iel], x_i1[iel], x_i1[iel]]
ys = [0, f_numeric(x_i[iel]), f_numeric(x_i1[iel]), 0]
ax.fill(xs, ys, facecolor = to_rgba('cyan', alpha=0.7),
edgecolor = edge_color)
# Plot the function on top
x_plot = np.linspace(x1, x2, 500)
y_plot = f_numeric(x_plot)
ax.plot(x_plot, y_plot, color='blue', linewidth=2,
label=f'$f(x) = {sp.latex(f_expr)}$', zorder=5)
# Set labels and limits
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$f(x)$', fontsize=12)
ax.set_xlim(x1, x2)
y_max = max(f_numeric(np.linspace(x1, x2, 500)))
ax.set_ylim(0, y_max * 1.2)
ax.grid(True, linestyle='--', alpha=0.7)
ax.axhline(0, color='black', linewidth=0.5)
ax.axvline(0, color='black', linewidth=0.5)
# Add text box with results
textstr = f'$N = {N}$\n'
textstr += f'I_e: ${I_exact:.6f}$\n'
textstr += f'I_R: ${I_riemann:.6f}$, error:${error_riemann:.4f}$\n'
textstr += f'I_T: ${I_trapz:.6f}$, error:${error_trapz:.4f}$'
props = dict(boxstyle='round', facecolor='w', alpha=0.8)
ax.text(0.02, 0.98, textstr, transform=ax.transAxes, fontsize=10,
verticalalignment='top', bbox=props)
ax.legend(loc='upper right')
def animate(frame_idx):
"""Animation function for each frame"""
ax.clear()
N = N_values[frame_idx]
makePlot(N)
# Create figure and axis (close to prevent static display)
fig, ax = plt.subplots(figsize=(6, 4))
plt.close(fig) # Prevent static image from showing
# Create animation
anim = FuncAnimation(fig, animate, frames=len(N_values),
interval=800, repeat=True)
mk.to_responsive_html(anim, container_id='trapz_animation')