Recently, I worked on a data science project that required calculating the area under a curve and solving differential equations. I faced the challenge of finding a reliable and efficient method for numerical integration in Python.
Fortunately, SciPy’s integrate module proved to be extremely helpful. In this article, I will demonstrate how to use the SciPy integrate module for various integration tasks, ranging from simple to complex. I’ll also share practical examples from my own Python journey.
Let’s get in!
What is SciPy Integrate?
SciPy integrate is a useful module within the SciPy library that provides functions for numerical integration. Whether you need to compute definite integrals, solve ordinary differential equations (ODEs), or work with multiple integrals, this module has you covered.
In my experience, SciPy integrates stands out because it combines ease of use with sophisticated numerical methods, making it suitable for both beginners and advanced users.
Read Python SciPy Pairwise Distance
Set Up Your Environment
Before we start using SciPy integrate, let’s make sure you have everything set up:
# Install SciPy if you haven't already
# pip install scipy
# Import the necessary modules
import numpy as np
from scipy import integrate
import matplotlib.pyplot as pltCheck out Python SciPy Smoothing
Method 1: Single Integration with quad()
Python quad() function is my go-to tool for calculating definite integrals of single-variable functions. Here’s how to use it:
def my_function(x):
return x**2
# Integrate x^2 from 0 to 1
result, error = integrate.quad(my_function, 0, 1)
print(f"Result: {result}") # Should be 1/3 (approximately 0.333...)
print(f"Estimated error: {error}")Output:
Result: 0.33333333333333337
Estimated error: 3.700743415417189e-15You can see the output in the screenshot below.
What I love about quad() is that it returns both the result and an error estimate, which helps assess the accuracy of your calculation.
Handle Infinite Limits
Sometimes you need to integrate over infinite intervals. SciPy makes this surprisingly easy:
# Integrate e^(-x^2) from -infinity to infinity
def gaussian(x):
return np.exp(-x**2)
result, error = integrate.quad(gaussian, -np.inf, np.inf)
print(f"Result: {result}") # Should be approximately sqrt(π) ≈ 1.77Method 2: Multiple Integration
When I work with multivariable functions, I use these specialized functions:
Double Integration with dblquad()
def integrand(y, x):
return x*y**2
# Integrate x*y^2 over the rectangle [0,1]×[0,1]
result, error = integrate.dblquad(integrand, 0, 1, # x limits
lambda x: 0, lambda x: 1) # y limits
print(f"Double integral result: {result}")Output:
Double integral result: 0.16666666666666669You can see the output in the screenshot below.
Notice that the order of arguments is important: the first argument after the function represents the outer integral variable.
Triple Integration with tplquad()
def integrand(z, y, x):
return x*y*z
# Integrate x*y*z over the cube [0,1]×[0,1]×[0,1]
result, error = integrate.tplquad(integrand, 0, 1, # x limits
lambda x: 0, lambda x: 1, # y limits
lambda x, y: 0, lambda x, y: 1) # z limits
print(f"Triple integral result: {result}") # Should be 1/8Read Python SciPy Ndimage Imread Tutorial
Method 3: Solve Ordinary Differential Equations (ODEs)
In many real-world applications, like predicting stock market trends or simulating physical systems, I’ve needed to solve differential equations. SciPy’s solve_ivp() function is perfect for this:
# Solve the differential equation dy/dt = -k*y
def model(t, y, k):
return -k * y
# Parameter k
k = 0.3
# Time points
t_span = (0, 10)
t_eval = np.linspace(0, 10, 100)
# Initial condition: y(0) = 5
y0 = [5]
# Solve the ODE
solution = integrate.solve_ivp(
lambda t, y: model(t, y, k),
t_span,
y0,
t_eval=t_eval,
method='RK45'
)
# Plot the solution
plt.figure(figsize=(10, 6))
plt.plot(solution.t, solution.y[0])
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.title('Solution to dy/dt = -0.3y with y(0) = 5')
plt.grid(True)Different ODE Solvers
SciPy offers various methods for solving ODEs. I’ve found these particularly useful:
RK45: A variable step-size Runge-Kutta method (my default choice)Radau: For stiff ODEsBDF: Another good choice for stiff problemsLSODA: Automatically selects between methods for stiff and non-stiff problems
# Using a different solver (BDF for stiff problems)
solution_bdf = integrate.solve_ivp(
lambda t, y: model(t, y, k),
t_span,
y0,
t_eval=t_eval,
method='BDF'
)Method 4: Numerical Integration Using Predefined Rules
When working with data points rather than functions, I use these methods:
# Simpson's rule for numerical integration
x = np.linspace(0, 1, 10)
y = x**2
result = integrate.simpson(y, x)
print(f"Simpson's rule result: {result}") # Should be close to 1/3
# Trapezoidal rule
result = integrate.trapezoid(y, x)
print(f"Trapezoidal rule result: {result}")Output:
Simpson's rule result: 0.33333333333333337
Trapezoidal rule result: 0.33539094650205764You can see the output in the screenshot below.
Check out Use Python SciPy Differential Evolution
Method 5: Practical Example – Compute Economic Indicators
Let’s look at a real-world example where I used SciPy integrate to calculate the Gini coefficient (a measure of income inequality) for hypothetical U.S. household income data:
# Generate sample income distribution (log-normal distribution is typical for incomes)
np.random.seed(42) # for reproducibility
incomes = np.random.lognormal(mean=10.5, sigma=1, size=1000) # in USD
# Sort incomes for Lorenz curve calculation
incomes_sorted = np.sort(incomes)
population = np.linspace(0, 1, len(incomes_sorted))
income_share = np.cumsum(incomes_sorted) / np.sum(incomes_sorted)
# Calculate Gini coefficient using the trapezoidal rule
# Gini = 1 - 2 * (area under Lorenz curve)
area_under_lorenz = integrate.trapezoid(income_share, population)
gini = 1 - 2 * area_under_lorenz
print(f"Gini coefficient: {gini:.4f}")
# Plot Lorenz curve
plt.figure(figsize=(10, 6))
plt.plot(population, income_share, label='Lorenz Curve')
plt.plot([0, 1], [0, 1], 'k--', label='Line of Equality')
plt.xlabel('Cumulative share of population')
plt.ylabel('Cumulative share of income')
plt.title(f'U.S. Income Distribution (Gini Coefficient: {gini:.4f})')
plt.legend()
plt.grid(True)Method 6: Solve Real-World Problems with odeint()
For systems of ODEs, I sometimes prefer the slightly older odeint() function due to its simpler interface:
# SIR Model for disease spread (e.g., COVID-19 in a U.S. city)
def sir_model(y, t, beta, gamma):
S, I, R = y
dSdt = -beta * S * I
dIdt = beta * S * I - gamma * I
dRdt = gamma * I
return [dSdt, dIdt, dRdt]
# Parameters
beta = 0.3 # infection rate
gamma = 0.1 # recovery rate
# Initial conditions: 1M population, 100 infected, 0 recovered
y0 = [0.999, 0.001, 0.0] # As proportions of total population
t = np.linspace(0, 100, 1000) # Time grid in days
# Solve the ODE system
solution = integrate.odeint(sir_model, y0, t, args=(beta, gamma))
S, I, R = solution.T
# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(t, S, 'b', label='Susceptible')
plt.plot(t, I, 'r', label='Infected')
plt.plot(t, R, 'g', label='Recovered')
plt.xlabel('Days')
plt.ylabel('Proportion of Population')
plt.title('SIR Model for Disease Spread')
plt.legend()
plt.grid(True)SciPy’s integrate module has been an essential tool in my Python toolkit for years. Whether I’m analyzing financial data, modeling physical systems, or solving complex mathematical problems, these integration techniques have helped me tackle challenging real-world scenarios.
The methods I’ve shared here range from simple definite integrals to sophisticated differential equation solvers, giving you a comprehensive toolbox for numerical integration in Python. If you’re working with any form of quantitative analysis, these techniques will prove invaluable.
You may read other SciPy-related articles:
I am Bijay Kumar, a Microsoft MVP in SharePoint. Apart from SharePoint, I started working on Python, Machine learning, and artificial intelligence for the last 5 years. During this time I got expertise in various Python libraries also like Tkinter, Pandas, NumPy, Turtle, Django, Matplotlib, Tensorflow, Scipy, Scikit-Learn, etc… for various clients in the United States, Canada, the United Kingdom, Australia, New Zealand, etc. Check out my profile.
