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:

Python
Copy to clipboard
1import numpy as np
2
3def f(x, y):
4    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:

Python
Copy to clipboard
1x_lower = 0
2x_upper = 2
3y_lower = 0
4y_upper = 1
5
6x = np.linspace(x_lower, x_upper, 100)
7y = np.linspace(y_lower, y_upper, 100)
8X, Y = np.meshgrid(x, y)
9Z = 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:

Python
Copy to clipboard
1fig = plt.figure()  # Create a new figure
2ax = fig.add_subplot(111, projection='3d')  # Add a 3D subplot to the figure
3ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)  # Plot the surface with a colormap and transparency
4
5# Highlight the area of integration
6x_fill = np.linspace(x_lower, x_upper, 20)  # Create a finer grid for x
7y_fill = np.linspace(y_lower, y_upper, 20)  # Create a finer grid for y
8X_fill, Y_fill = np.meshgrid(x_fill, y_fill)  # Generate a meshgrid for the fill area
9Z_fill = f(X_fill, Y_fill)  # Compute the function values for the fill area
10ax.plot_surface(X_fill, Y_fill, np.zeros_like(Z_fill), color='r', alpha=0.3)  # Plot the integration area
11
12ax.set_xlabel('X')  # Set the x-axis label
13ax.set_ylabel('Y')  # Set the y-axis label
14ax.set_zlabel('Integrand Value')  # Set the z-axis label
15ax.set_title('Visualization of the Integrand Function with Integration Area')  # Set the plot title
16plt.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:

Python
Copy to clipboard
1from scipy.integrate import dblquad

Now, use dblquad to calculate the double integral:

Python
Copy to clipboard
1double_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:

Python
Copy to clipboard
1double_integral, error = dblquad(
2    f,  # The function to integrate
3    x_lower,  # Lower limit for x
4    x_upper,  # Upper limit for x
5    lambda x: y_lower,  # Lower limit for y as a function of x
6    lambda x: y_upper  # Upper limit for y as a function of x
7)

Finally, print the result:

Python
Copy to clipboard
1print(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:

Python
Copy to clipboard
1# Define x bounds
2x_lower = -1
3x_upper = 1
4
5# Define y bounds as functions of x
6y_lower = lambda x: -np.sqrt(1 - x**2)
7y_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:

Python
Copy to clipboard
1double_integral_circle, error_circle = dblquad(
2    f,  # The function to integrate
3    x_lower,  # Lower limit for x
4    x_upper,  # Upper limit for x
5    y_lower,  # Lower limit for y as a function of x
6    y_upper  # Upper limit for y as a function of x
7)

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!