Skip to main content

Numerical Integration & Optimization: NumPy

Numerical integration and optimization are core techniques for solving real-world problems: computing areas under curves, finding roots of equations, fitting models to data, and minimizing functions. While NumPy provides foundational array operations, scipy.integrate and scipy.optimize extend it with specialized algorithms. This article teaches you to implement and use numerical methods efficiently, when to apply each technique, and how to validate results.

Numerical Integration: Computing Areas​

Integration computes the area under a curve. Numerical methods approximate the integral when analytical solutions don't exist.

Trapezoidal Rule (Simple Approximation)​

The trapezoidal rule breaks the interval into trapezoids and sums their areas—simple but effective for smooth functions.

import numpy as np
from scipy import integrate

# Define a function to integrate
def f(x):
    return np.sin(x)

# Trapezoidal rule: NumPy built-in
x = np.linspace(0, np.pi, 100)  # 100 points from 0 to π
y = f(x)
integral_trap = np.trapz(y, x)  # Trapezoidal approximation
print(f"Trapezoidal integral: {integral_trap:.6f}")
print(f"Exact value (2): {2.0:.6f}")
print(f"Error: {abs(integral_trap - 2.0):.6f}")

# Finer grid → better approximation
x_fine = np.linspace(0, np.pi, 1000)
y_fine = f(x_fine)
integral_fine = np.trapz(y_fine, x_fine)
print(f"Finer trapezoidal integral: {integral_fine:.6f}")

Scipy's quad(): Higher-Accuracy Integration​

scipy.integrate.quad() uses adaptive quadrature (Gauss-Kronrod) for high accuracy with fewer function evaluations.

import numpy as np
from scipy import integrate

# Function to integrate
def f(x):
    return np.exp(-x**2)  # Gaussian

# quad() returns (result, error_estimate)
result, error = integrate.quad(f, 0, 1)
print(f"Integral of exp(-x^2) from 0 to 1: {result:.10f}")
print(f"Estimated error: {error:.2e}")

# Practical: integrate a complex function
def integrand(x):
    return np.sin(x) / (1 + x**2)

integral, error = integrate.quad(integrand, 0, np.pi)
print(f"Integral of sin(x)/(1+x^2): {integral:.6f}")

Multiple Integration​

For double and triple integrals, use dblquad or tplquad:

from scipy import integrate

# Double integral: ∫∫ x*y dxdy over unit square
def integrand_2d(y, x):
    return x * y

result, error = integrate.dblquad(integrand_2d, 0, 1, 0, 1)
print(f"Double integral: {result:.6f}")
# Analytical: (∫x dx)(∫y dy) = (1/2)(1/2) = 1/4

Root Finding: Solving f(x) = 0​

Finding roots is essential for solving equations and parameter fitting.

Newton-Raphson Method (Derivative-Based)​

Fast but requires the derivative:

from scipy import optimize
import numpy as np

# Find root of f(x) = x^3 - 2
def f(x):
    return x**3 - 2

def f_prime(x):
    return 3 * x**2

# Newton-Raphson (requires derivative)
root = optimize.newton(f, x0=1.5, fprime=f_prime)
print(f"Root of x^3 - 2: {root:.10f}")
print(f"Verification: f(root) = {f(root):.2e}")

# Secant method (no derivative)
root_secant = optimize.newton(f, x0=1.5)  # fprime omitted
print(f"Root (secant): {root_secant:.10f}")

Bracketed Root Finding (Robust)​

When you know the root is between a and b, use bracketing methods:

from scipy import optimize
import numpy as np

def f(x):
    return np.sin(x) - x/2

# Brent's method (fast, robust)
root = optimize.brentq(f, 0, 2)  # Requires f(a) and f(b) to have opposite signs
print(f"Root of sin(x) - x/2: {root:.10f}")
print(f"Verification: f(root) = {f(root):.2e}")

# Multiple roots: use fsolve for systems
def system(x):
    return [x[0]**2 + x[1]**2 - 1,  # Circle
            x[0] - x[1]]             # Line

solution = optimize.fsolve(system, x0=[0, 0])
print(f"Intersection of circle and line: {solution}")

Optimization: Minimizing Functions​

Univariate Minimization​

For single-variable optimization:

from scipy import optimize
import numpy as np

# Function with a minimum at x = sqrt(2)
def objective(x):
    return (x - np.sqrt(2))**2

# Minimize
result = optimize.minimize_scalar(objective, method='brent')
print(f"Minimum at x = {result.x:.10f}")
print(f"Minimum value: {result.fun:.2e}")

# Bounded minimization
result_bounded = optimize.minimize_scalar(objective, bounds=(1, 2), method='bounded')
print(f"Bounded minimum at x = {result_bounded.x:.10f}")

Multivariate Optimization​

from scipy import optimize
import numpy as np

# Rosenbrock function: a classic optimization benchmark
def rosenbrock(x):
    return 100*(x[1] - x[0]**2)**2 + (1 - x[0])**2

# Start from initial guess
x0 = np.array([-1, -1])
result = optimize.minimize(rosenbrock, x0, method='BFGS')
print(f"Optimum at: {result.x}")
print(f"Optimal value: {result.fun:.2e}")

# With gradient (faster)
def rosenbrock_grad(x):
    return np.array([
        -400*x[0]*(x[1] - x[0]**2) - 2*(1 - x[0]),
        200*(x[1] - x[0]**2)
    ])

result_grad = optimize.minimize(rosenbrock, x0, method='BFGS', jac=rosenbrock_grad)
print(f"Optimum with gradient: {result_grad.x}")

Curve Fitting​

Fit a model to data using least-squares:

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

# Generate synthetic data: y = 2*x + 3 + noise
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 2*x + 3 + np.random.randn(50) * 2

# Define model: y = a*x + b
def model(x, a, b):
    return a*x + b

# Fit using curve_fit
popt, pcov = optimize.curve_fit(model, x, y)
a_fit, b_fit = popt
print(f"Fitted parameters: a = {a_fit:.4f}, b = {b_fit:.4f}")
print(f"True parameters: a = 2.0, b = 3.0")

# Predictions
y_fit = model(x, *popt)
residuals = y - y_fit
mse = np.mean(residuals**2)
print(f"Mean squared error: {mse:.6f}")

# Covariance matrix: diagonal elements are parameter variances
perr = np.sqrt(np.diag(pcov))
print(f"Standard errors: a ± {perr[0]:.4f}, b ± {perr[1]:.4f}")

# Exponential fit
def exp_model(x, a, k):
    return a * np.exp(-k * x)

# Data: exponential decay
x_exp = np.linspace(0, 5, 50)
y_exp = 10 * np.exp(-0.5 * x_exp) + np.random.randn(50) * 0.1

popt_exp, _ = optimize.curve_fit(exp_model, x_exp, y_exp, p0=[10, 0.5])
print(f"Exponential fit: a = {popt_exp[0]:.4f}, k = {popt_exp[1]:.4f}")

Practical Example: Fitting a Sigmoid​

Sigmoidal curves appear in biology, neural networks, and learning curves. Fit one using curve_fit:

import numpy as np
from scipy import optimize

# Sigmoid function
def sigmoid(x, x0, k, L):
    return L / (1 + np.exp(-k * (x - x0)))

# Generate synthetic S-curve data
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = sigmoid(x, x0=5, k=1, L=1) + np.random.randn(50) * 0.05

# Fit
popt, pcov = optimize.curve_fit(sigmoid, x, y, p0=[5, 1, 1])
x0_fit, k_fit, L_fit = popt
print(f"Fitted sigmoid: x0={x0_fit:.2f}, k={k_fit:.2f}, L={L_fit:.2f}")

# Smooth prediction
x_smooth = np.linspace(0, 10, 1000)
y_smooth = sigmoid(x_smooth, *popt)
print(f"Mid-point (x0): {x0_fit:.2f}")
print(f"Steepness (k): {k_fit:.2f}")
print(f"Maximum value (L): {L_fit:.2f}")

Key Takeaways​

  • Use np.trapz() for simple trapezoidal integration; use scipy.integrate.quad() for high accuracy.
  • For multiple integrals, use dblquad or tplquad from scipy.integrate.
  • Root finding: Newton-Raphson is fast with derivatives; bracketing methods are robust without derivatives.
  • scipy.optimize.minimize() handles multivariate optimization with optional gradients (faster).
  • curve_fit() fits models to data; returns parameters and covariance (uncertainty estimates).
  • Always validate: check residuals, condition numbers, and sensitivities.

Frequently Asked Questions​

Should I always provide a derivative to scipy.optimize?​

Providing derivatives (via jac or fprime) usually improves convergence speed and stability. If computing the derivative is expensive or error-prone, numerical differentiation is acceptable.

How do I know if my optimization converged?​

Check result.success and result.message. Verify by plotting or comparing objective function values. Re-run from different initial guesses to check consistency.

What if my function has multiple roots or minima?​

Root finding and optimization find local solutions. For global solutions, use scipy.optimize.basinhopping() or differential_evolution(), or try multiple initial guesses.

How do I constrain optimization (bounds, inequality constraints)?​

Use bounds parameter in minimize() for box constraints. For general constraints, use constraints parameter with dictionaries specifying type and functions.

Why does curve_fit sometimes fail?​

Common causes: poor initial guess (provide reasonable p0), ill-conditioned problem, or numerically unstable model. Rescale data to [-1, 1] range for stability.

Further Reading​