⬡ Hub
Skip to content

SciPy: Optimization and Integration

SciPy builds on NumPy to provide a collection of routines for common scientific and engineering tasks. Among its most powerful modules are scipy.optimize for finding minima, maxima, and roots of functions, and scipy.integrate for numerical integration.

1. Optimization (scipy.optimize)

The scipy.optimize module provides various optimization algorithms, including minimization (finding the lowest point of a function), curve fitting, and root finding.

a. Minimization of Scalar Functions (minimize_scalar)

For finding the minimum of a single-variable function.

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

# Define the function to minimize: f(x) = x^2 + 10*sin(x)
def f(x):
    return x**2 + 10 * np.sin(x)

# Plot the function to visualize minima
x_vals = np.linspace(-10, 10, 500)
y_vals = f(x_vals)

plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals, label='f(x) = x^2 + 10*sin(x)')
plt.title('Function for Minimization')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)

# Find the minimum using bounded optimization
result = minimize_scalar(f, bounds=(-10, 10), method='bounded')

print("Minimization Result (Bounded):")
print(f"Success: {result.success}")
print(f"Status: {result.status}")
print(f"Message: {result.message}")
print(f"Optimal x: {result.x:.4f}")
print(f"Minimum f(x): {result.fun:.4f}")

plt.scatter(result.x, result.fun, color='red', marker='o', s=100, label='Global Minimum Found')
plt.legend()
plt.show()

# You can also use other methods like 'golden' or 'brent' for unbounded search (within an interval)
# result_golden = minimize_scalar(f, bracket=(-5, 0, 5), method='golden')
# print(f"\nOptimal x (Golden): {result_golden.x:.4f}, f(x): {result_golden.fun:.4f}")

b. Minimization of Multi-variable Functions (minimize)

For finding the minimum of functions with multiple variables. This is crucial for training machine learning models.

import numpy as np
from scipy.optimize import minimize

# Define the function to minimize: f(x, y) = (x-1)^2 + (y-2)^2
def objective_function(params):
    x, y = params
    return (x - 1)**2 + (y - 2)**2

# Initial guess
x0 = np.array([0, 0])

# Perform minimization using BFGS method
result = minimize(objective_function, x0, method='BFGS')

print("\nMinimization Result (Multi-variable, BFGS):")
print(f"Success: {result.success}")
print(f"Message: {result.message}")
print(f"Optimal x, y: {result.x}")
print(f"Minimum f(x, y): {result.fun:.4f}")

# For constrained optimization, you can pass constraints and bounds
# from scipy.optimize import Bounds, LinearConstraint
# bounds = Bounds([-10, -10], [10, 10])
# linear_constraint = LinearConstraint([1, 1], -np.inf, 3) # x + y <= 3
# result_constrained = minimize(objective_function, x0, method='SLSQP', bounds=bounds, constraints=linear_constraint)
# print(f"\nOptimal x, y (Constrained): {result_constrained.x}, f(x, y): {result_constrained.fun}")

c. Curve Fitting (curve_fit)

curve_fit uses non-linear least squares to fit a function, f, to data.

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

# Define the function we want to fit (e.g., a power law)
def power_law(x, a, b):
    return a * x**b

# Generate synthetic data with noise
x_data = np.linspace(1, 10, 50)
y_true = 2 * x_data**1.5 # True parameters: a=2, b=1.5
y_data = y_true + 1.5 * np.random.normal(size=len(x_data)) # Add some noise

# Fit the function to the data
# popt: optimal values for the parameters
# pcov: covariance matrix of the parameters
popt, pcov = curve_fit(power_law, x_data, y_data)

print("\nCurve Fitting Results:")
print(f"Optimal parameters (a, b): {popt}")

# Plot the original data and the fitted curve
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='Noisy Data')
plt.plot(x_data, power_law(x_data, *popt), color='red', label='Fitted Curve')
plt.title('Curve Fitting with scipy.optimize.curve_fit')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.grid(True)
plt.show()

2. Integration (scipy.integrate)

The scipy.integrate module provides several techniques to compute definite integrals, both for functions and for samples.

a. Integrating Functions (quad, dblquad, tplquad)

For numerical integration of single, double, or triple integrals.

import numpy as np
from scipy.integrate import quad, dblquad
import matplotlib.pyplot as plt

# Define a function to integrate: f(x) = e^(-x^2)
def f_single(x):
    return np.exp(-x**2)

# Integrate f(x) from -infinity to +infinity (Gaussian integral)
result_quad, error_quad = quad(f_single, -np.inf, np.inf)
print(f"Single Integral (e^(-x^2) from -inf to inf): {result_quad:.4f} (Error: {error_quad:.2e})")
print(f"Expected: {np.sqrt(np.pi):.4f}")

# Define a function for double integral: f(x, y) = x*y
def f_double(y, x): # Note: order of arguments is (y, x) for dblquad
    return x * y

# Integrate f(x, y) = x*y over a region
# x from 0 to 1, y from 0 to 2
result_dblquad, error_dblquad = dblquad(f_double, 0, 1, lambda x: 0, lambda x: 2)
print(f"\nDouble Integral (x*y, x=[0,1], y=[0,2]): {result_dblquad:.4f} (Error: {error_dblquad:.2e})")
# Expected: (1/2) * (1/2) * (2^2 - 0^2) = 1/4 * 4 = 1.0

b. Integrating Sampled Data (simpson, trapz)

For integrating data points provided as arrays, rather than a function.

import numpy as np
from scipy.integrate import simpson, trapz

# Sample data
x_data = np.linspace(0, np.pi, 100)
y_data = np.sin(x_data) # Function is sin(x)

# Integrate using Simpson's rule
integral_simpson = simpson(y_data, x_data)
print(f"Integral of sin(x) from 0 to pi (Simpson): {integral_simpson:.4f}")

# Integrate using the trapezoidal rule
integral_trapz = trapz(y_data, x_data)
print(f"Integral of sin(x) from 0 to pi (Trapezoidal): {integral_trapz:.4f}")

# The true integral of sin(x) from 0 to pi is 2.0

Further Topics:

  • Root finding (fsolve, root)
  • Linear programming (linprog)
  • Global optimization (differential_evolution, basinhopping)
  • Ordinary Differential Equation (ODE) solvers (solve_ivp)
  • Numerical differentiation.

SciPy's optimization and integration tools are foundational for many scientific and engineering computations, enabling you to solve complex mathematical problems numerically where analytical solutions might be difficult or impossible to find.