Trapezoidal Rule for Numerical Integration in Python

1. From Rectangles to Trapezoids: A Better Approximation

Recall from our last lesson that the concept of an integral is built on the idea of summing the areas of infinitely many, infinitesimally thin rectangles under a curve (a Riemann sum). Numerical integration takes this idea but uses a finite number of shapes to approximate the area.

While rectangles are simple, they can lead to significant errors, either over- or under-shooting the curve. We can get a much better approximation by using trapezoids instead.

As the video explains, by fitting the top of our shape to the curve at both ends of the interval, we hug the function more closely, leading to a more accurate estimate of the area.

2. The Formula for the Trapezoidal Rule

Let's formalize this idea. To approximate the integral ${\int }_a^{b}f(x)dx$, we divide the interval $[a,b]$ into $N$ smaller subintervals, each of width $\Delta x=(b-a)/N$. Each of these subintervals forms the base of a trapezoid.

The area of a single trapezoid over a subinterval $[x_{i-1},x_i]$ is given by the average of the heights (the function values at the endpoints) multiplied by the width:

To find the total approximate area, we sum the areas of all $N$ trapezoids.

After summing all the trapezoid areas and collecting terms, we arrive at the general formula for the trapezoid rule:

Notice that the first and last function values are used once, while all the interior values are used twice (as they are part of two adjacent trapezoids).

3. Implementing the Trapezoidal Rule in Python

Now let's translate this formula into code. Your background in Python and NumPy from the first module will be very helpful here. We'll write a function that takes a function f, the interval [a, b], and the number of trapezoids N as arguments.

Let's break down how to build this function step-by-step.

import numpy as np

def trapezoid_rule(f, a, b, N):
    """
    Approximates the definite integral of f(x) from a to b
    using the trapezoidal rule with N subintervals.
    """
    # 1. Calculate the width of the trapezoids
    dx = (b - a) / N
    
    # 2. Create an array of N+1 points from a to b
    x = np.linspace(a, b, N + 1)
    
    # 3. Evaluate the function at these points
    y = f(x)
    
    # 4. Apply the trapezoidal rule formula
    # T = (dx/2) * (y_0 + 2y_1 + ... + 2y_{N-1} + y_N)
    # This can be implemented efficiently by summing all y values,
    # and then adding the interior points once more.
    
    integral_approx = (dx / 2) * (y[0] + y[-1] + 2 * np.sum(y[1:-1]))
    
    return integral_approx

Let's test it on an integral we can solve analytically. From our last lesson, we know ${\int }_0^{1}3x^{2}dx=[x^3{]}_0^1=1^3-0^3=1$.

# Define the function to integrate
def my_func(x):
    return 3 * x**2

# Approximate the integral
approx_val = trapezoid_rule(my_func, 0, 1, 100)
print(f"Approximation: {approx_val}")
print(f"Exact value: 1")

Running this code, you'll see the approximation is very close to 1!

# First, define the function using numpy's sin
def sin_func(x):
    return np.sin(x)

# With N=10
approx_10 = trapezoid_rule(sin_func, 0, np.pi/2, 10)
print(f"Approximation with N=10: {approx_10}")

# With N=1000
approx_1000 = trapezoid_rule(sin_func, 0, np.pi/2, 1000)
print(f"Approximation with N=1000: {approx_1000}")

4. Error Analysis and Using Library Functions

A key part of numerical methods is understanding the error. As we saw, increasing N reduces the error. The following video provides a complete walkthrough of implementing the trapezoid rule and, importantly, visualizing how the error decreases as N increases.

While implementing algorithms from scratch is fantastic for learning, in practice, you'll often use highly optimized functions from libraries like NumPy or SciPy. Both provide a trapz function.

Here's how you would use NumPy's trapz function to do the same calculation. It's even simpler because the function deduces the interval width from the x values you provide.

import numpy as np

# Define the function, limits, and N
f = lambda x: 3 * x**2
a, b, N = 0, 1, 100

# Create the x and y data points
x_vals = np.linspace(a, b, N + 1)
y_vals = f(x_vals)

# Use numpy's built-in function
# Note: The order is y, then x.
numpy_approx = np.trapz(y_vals, x_vals)

print(f"NumPy's trapz approximation: {numpy_approx}")

As you can see, libraries make the process very straightforward, but now you understand exactly what's happening under the hood. The resource "Implementing The Trapezium Rule With Python" (ID: LINK) even includes a performance benchmark showing how much faster these library functions are compared to a pure Python implementation, especially for large N.

Conclusion

Congratulations on completing this lesson and the module on Calculus! You now have a solid foundation in both the analytical and numerical sides of integration.

Key Takeaways:

• Numerical integration is used to approximate definite integrals when an analytical solution is difficult or impossible, or when you only have data points.
• The trapezoidal rule improves on simple rectangle-based approximations by fitting the curve more closely, reducing error.
• The formula for the trapezoidal rule is a weighted sum of the function's values at evenly spaced points:
  $T_N(f)=\frac{\Delta x}{2}\left[y_0+2y_1+\cdots +2y_{N-1}+y_N\right]$
• Increasing the number of trapezoids (N) increases the accuracy of the approximation.
• We can implement the rule from scratch in Python, but for practical applications, optimized functions like numpy.trapz are preferred for their simplicity and performance.

Preview of the Next Module:

We have now covered the essentials of Python programming, Linear Algebra, and Calculus—three major pillars of data science and AI. You're well-equipped to move on to the next one: Probability and Statistics. In our next lesson, we will start with the fundamentals of counting by learning to apply fundamental counting principles, including permutations and combinations. These concepts form the bedrock for calculating probabilities.