Double Integrals with SciPy

Welcome to the lesson on double integrals, an essential concept in calculus. In previous lessons, you have learned how to perform definite integration using the SciPy library in Python. Today, we're going to expand on that knowledge by diving into double integrals.

Double integrals allow us to integrate functions of two variables over a specific region. This is incredibly useful in fields like physics and engineering, where you often need to calculate areas, volumes, or other properties across two-dimensional spaces.

Now, let's define the mathematical function we want to integrate.

Here's how you can define a simple function that involves two variables, x and y:

import numpy as np

def f(x, y):
    return x * y + np.sin(x)

Our task will be to compute the following integral:

$\int_{x=0}^{2} \int_{y=0}^{1} \left( x \cdot y + \sin(x) \right) \, dy \, dx$

Visualizing the function and area of integration can provide significant insight into the nature of the task. We'll use Matplotlib to create a 3D plot for this.

First, we create grid data for x, y, and the function values z:

x_lower = 0
x_upper = 2
y_lower = 0
y_upper = 1

x = np.linspace(x_lower, x_upper, 100)
y = np.linspace(y_lower, y_upper, 100)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

x_lower, x_upper, y_lower, y_upper are the integration limits. np.meshgrid helps in generating a grid of x and y values, while Z contains the corresponding function values.

Now, plot the surface using Matplotlib:

fig = plt.figure()  # Create a new figure
ax = fig.add_subplot(111, projection='3d')  # Add a 3D subplot to the figure
ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)  # Plot the surface with a colormap and transparency

# Highlight the area of integration
x_fill = np.linspace(x_lower, x_upper, 20)  # Create a finer grid for x
y_fill = np.linspace(y_lower, y_upper, 20)  # Create a finer grid for y
X_fill, Y_fill = np.meshgrid(x_fill, y_fill)  # Generate a meshgrid for the fill area
Z_fill = f(X_fill, Y_fill)  # Compute the function values for the fill area
ax.plot_surface(X_fill, Y_fill, np.zeros_like(Z_fill), color='r', alpha=0.3)  # Plot the integration area

ax.set_xlabel('X')  # Set the x-axis label
ax.set_ylabel('Y')  # Set the y-axis label
ax.set_zlabel('Integrand Value')  # Set the z-axis label
ax.set_title('Visualization of the Integrand Function with Integration Area')  # Set the plot title
plt.show()  # Display the plot

Here is the output plot:

• The surface plot (ax.plot_surface) displays the function across the integration domain.
• The highlighted area (in red) represents the integration limits.

We will use SciPy's dblquad function to compute the double integral of our defined function. Let's break this down step by step.

First, you need to import dblquad from SciPy's integration module:

from scipy.integrate import dblquad

Now, use dblquad to calculate the double integral:

double_integral, error = dblquad(f, x_lower, x_upper, lambda x: y_lower, lambda x: y_upper)

• dblquad takes the function f, the x limits, and two lambda functions defining the y limits dependent on x.
• double_integral stores the value of the integral, while error estimates the numerical approximation error.

Here's the method call with respective arguments for better demonstration:

double_integral, error = dblquad(
    f,  # The function to integrate
    x_lower,  # Lower limit for x
    x_upper,  # Upper limit for x
    lambda x: y_lower,  # Lower limit for y as a function of x
    lambda x: y_upper  # Upper limit for y as a function of x
)

Finally, print the result:

print(f"Double integral: {double_integral}") # 1.92

In the dblquad function, the y_bounds are defined with lambda functions because the limits for the inner integral (with respect to y) can depend on the variable of the outer integral (x). Using lambda functions allows for flexibility in specifying these y limits as functions of x. In simple cases where the integration limits are constant, the lambda function simply returns those constants.

In certain cases, the integration bounds for y can depend on the variable x. This is particularly useful when dealing with regions like circles or other complex shapes. Let's illustrate how to define these dependent bounds for a circular region centered at (0, 0) with a radius of 1. The equation for this circle is $x^2 + y^2 \leq 1$, the bounds for x will range from -1 to 1. The bounds for y need to be defined as functions of x to represent the circle:

# Define x bounds
x_lower = -1
x_upper = 1

# Define y bounds as functions of x
y_lower = lambda x: -np.sqrt(1 - x**2)
y_upper = lambda x: np.sqrt(1 - x**2)

Here, y_lower(x) is $-\sqrt{1 - x^2}$ and y_upper(x) is $\sqrt{1 - x^2}$, capturing the upper and lower semicircles for any given x value.

To demonstrate how to call the method with these bounds:

double_integral_circle, error_circle = dblquad(
    f,  # The function to integrate
    x_lower,  # Lower limit for x
    x_upper,  # Upper limit for x
    y_lower,  # Lower limit for y as a function of x
    y_upper  # Upper limit for y as a function of x
)

In this lesson, you have learned how to define a function of two variables, compute its double integral using SciPy, and visualize the results with Matplotlib. These steps are crucial in understanding and applying double integrals in various scientific and engineering contexts.

As you move forward to the hands-on practice, experiment with different functions and integration bounds to reinforce your understanding. You're progressing well in this course, and with these foundations, you're well-equipped to tackle more complex problems. Keep up the great work!