Scipy Integration

Introduction

SciPy provides numerical integration functions for definite integrals, double integrals, and solving ordinary differential equations.

Numerical Integration

from scipy import integrate
import numpy as np

# Single integral
result, error = integrate.quad(lambda x: x**2, 0, 1)
print(f"Result: {result}, Error estimate: {error}")

# Integral with infinite bounds
result, error = integrate.quad(lambda x: np.exp(-x), 0, np.inf)
print(result)  # 1.0

# Integrand with parameters
def f(x, a, b):
    return a * x + b
result, error = integrate.quad(f, 0, 1, args=(2, 3))

Multiple Integration

from scipy import integrate

# Double integral
def f(x, y):
    return x + y

result, error = integrate.dblquad(f, 0, 2, lambda x: 0, lambda x: 1)
print(result)

# Triple integral
result = integrate.tplquad(lambda x, y, z: x*y*z, 0, 1, 0, 2, lambda x, y: 0, lambda x, y: 1)

Ordinary Differential Equations

from scipy.integrate import solve_ivp
import numpy as np

# Define ODE system
def dy_dt(t, y):
    return [-y[0] + y[1], -y[0] - y[1]]

# Solve
t_span = (0, 10)
y0 = [1, 0]
sol = solve_ivp(dy_dt, t_span, y0, t_eval=np.linspace(0, 10, 100))

print(sol.t)      # Time points
print(sol.y[0])   # First variable solution

ODE with Parameters

from scipy.integrate import solve_ivp

def lorenz(t, state, sigma, rho, beta):
    x, y, z = state
    return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]

params = (10, 28, 8/3)
sol = solve_ivp(lorenz, (0, 50), [1, 1, 1], args=params, dense_output=True)

# Interpolate at specific times
t_vals = [5, 10, 15]
y_vals = sol.sol(t_vals)

Practice Problems

1. Calculate definite integral of trigonometric function
2. Solve first-order ODE system
3. Implement double integration
4. Solve Lorenz system with proper parameters
5. Compare quad vs trapz accuracy