The Ultimate SciPy Cheatsheet

A comprehensive guide to SciPy's core modules for scientific computing, including optimization, signal processing, and numerical methods.

Equation Solving

SciPy provides powerful tools for solving systems of linear and nonlinear equations through scipy.linalg and scipy.optimize modules.

Linear Equation Systems

Solve systems of linear equations using LU factorization:

import numpy as np
from scipy import linalg as la

# Define the coefficient matrix and right-hand side vector
A = np.array([[2, 3], [5, 4]])
b = np.array([4, 3])

# Perform LU factorization
P, L, U = la.lu(A)

# Display the L and U matrices
print("L matrix:")
print(L)
# Output: array([[ 1. ,  0. ],
#                 [ 0.4,  1. ]])

print("U matrix:")
print(U)
# Output: array([[ 5. ,  4. ],
#                 [ 0. ,  1.4]])

# Solve the system Ax = b
x = la.solve(A, b)
print("Solution:", x)
# Output: array([-1.,  2.])

Linear Least Squares (Overdetermined Systems)

Fit a quadratic model to noisy data using linear least-squares:

# Define true model parameters
x = np.linspace(-1, 1, 100)
a, b, c = 1, 2, 3
y_exact = a + b * x + c * x**2

# Simulate noisy data
m = 100
X = 1 - 2 * np.random.rand(m)
Y = a + b * X + c * X**2 + np.random.randn(m)

# Fit the data to the model using linear least squares
A = np.vstack([X**0, X**1, X**2])
sol, r, rank, sv = la.lstsq(A.T, Y)

# Use the solution to plot the fitted curve
y_fit = sol[0] + sol[1] * x + sol[2] * x**2

Nonlinear Univariate Equations (Root Finding)

Find the root of the equation $e^x - 2 = 0$ using several methods:

from scipy import optimize

# Define the function
f = lambda x: np.exp(x) - 2

# Using the bisection method
result = optimize.bisect(f, -2, 2)
print("Bisection result:", result)
# Output: 0.6931471805592082

# Using Newton's method (secant variant)
x_root_guess = 2
result = optimize.newton(f, x_root_guess)
print("Newton result:", result)
# Output: 0.69314718056

# Newton's method with derivative
fprime = lambda x: np.exp(x)
result = optimize.newton(f, x_root_guess, fprime=fprime)
print("Newton with derivative:", result)
# Output: 0.69314718056

# Using Brent's methods (generally preferred)
result1 = optimize.brentq(f, -2, 2)
result2 = optimize.brenth(f, -2, 2)
print("Brent methods:", result1, result2)

Systems of Nonlinear Equations

Solve the system of nonlinear equations $y = x^3 + 2x^2 - 1$ and $y = -x^2 + 1$:

def f(x):
    return [x[1] - x[0]**3 - 2 * x[0]**2 + 1,
            x[1] + x[0]**2 - 1]

# Solve using fsolve
solution = optimize.fsolve(f, [1, 1])
print("Solution:", solution)
# Output: array([ 0.73205081,  0.46410162])

# Providing the Jacobian can improve performance
def f_jacobian(x):
    return [[-3*x[0]**2 - 4*x[0], 1],
            [2*x[0], 1]]

solution = optimize.fsolve(f, [1, 1], fprime=f_jacobian)
print("Solution with Jacobian:", solution)

Interpolation

Construct interpolation functions from discrete data points using scipy.interpolate.

Spline Interpolation

Use a cubic spline to interpolate Runge's function:

from scipy import interpolate

# Define Runge's function
def runge(x):
    return 1/(1 + 25 * x**2)

# Generate sample points
x = np.linspace(-1, 1, 11)
y = runge(x)

# Create a cubic spline interpolation function
f_i = interpolate.interp1d(x, y, kind=3)

# Evaluate at finer points
xx = np.linspace(-1, 1, 100)
plt.figure(figsize=(8, 4))
plt.plot(xx, runge(xx), 'k', lw=1, label="Runge's function")
plt.plot(x, y, 'ro', label='sample points')
plt.plot(xx, f_i(xx), 'r--', lw=2, label='spline order 3')
plt.legend()
plt.show()

Multivariate Interpolation

Interpolate 2D data on both regular and irregular grids:

# 2D interpolation on a regular grid
x = y = np.linspace(-2, 2, 10)

def f(x, y):
    return (_exp(-(x + .5)**2 - 2*(y + .5)**2) -
            _exp(-(x - .5)**2 - 2*(y - .5)**2))

X, Y = np.meshgrid(x, y)
Z = f(X, Y) + 0.05 * np.random.randn(*X.shape)

# Create interpolation function
f_i = interpolate.interp2d(x, y, Z, kind='cubic')

# 2D interpolation for unstructured data
def f(x, y):
    return _exp(-x**2 - y**2) * _cos(4*x) * _sin(6*y)

N = 500
xdata = np.random.uniform(-1, 1, N)
ydata = np.random.uniform(-1, 1, N)
zdata = f(xdata, ydata)

# Regular grid for interpolation
x_grid = y_grid = np.linspace(-1, 1, 100)
X_grid, Y_grid = np.meshgrid(x_grid, y_grid)

# Interpolate unstructured data
Zi = interpolate.griddata((xdata, ydata), zdata,
                          (X_grid, Y_grid), method='cubic')

Integration

Numerical integration (quadrature) for functions and tabulated data using scipy.integrate.

Numerical Integration of Functions

Evaluate the integral $\int_{-1}^{1} e^{-x^2} dx$:

from scipy import integrate

def f(x):
    return np.exp(-x**2)

# Compute the integral
val, err = integrate.quad(f, -1, 1)

print("Integral value:", val)
# Output: 1.493648265624854

print("Estimated error:", err)
# Output: 1.6582826951881447e-14

Integration of Tabulated Data

Integrate sampled data points of the function $f(x) = \sqrt{x}$:

# Define the function and integration limits
f = lambda x: np.sqrt(x)
a, b = 0, 2

# Generate sample points
x = np.linspace(a, b, 25)
y = f(x)

# Use trapezoid rule
val_trapz = integrate.trapz(y, x)
print("Trapezoid rule:", val_trapz)
# Output: 1.88082171605

# Use Simpson's rule
val_simps = integrate.simps(y, x)
print("Simpson's rule:", val_simps)
# Output: 1.88366510245

Multiple Integration

Evaluate the double integral $\int_{0}^{1} \int_{0}^{1} e^{-x^2 - y^2} dx dy$:

def f(x, y):
    return np.exp(-x**2 - y**2)

# Define integration limits
a, b = 0, 1
g = lambda x: 0  # Lower limit for y
h = lambda x: 1  # Upper limit for y

# Compute double integral
result, error = integrate.dblquad(f, a, b, g, h)
print("Double integral result:", result)
# Output: (0.5577462853510337, 6.1922276789587025e-15)

Signal Processing

Spectral analysis with Fast Fourier Transform (FFT) and digital filters using scipy.signal and scipy.fftpack.

Spectral Analysis with FFT

Compute and plot the frequency spectrum of a noisy signal:

from scipy import fftpack, signal as sig
import matplotlib.pyplot as plt

# Define signal parameters
B = 30.0      # bandwidth
f_s = 2 * B   # sampling frequency
delta_f = 0.01 # frequency resolution
N = int(f_s / delta_f)
T = N / f_s
t = np.linspace(0, T, N)

# Signal generation function
def signal_samples(t):
    return (2 * np.sin(2 * np.pi * t) +
            3 * np.sin(22 * 2 * np.pi * t) +
            2 * np.random.randn(*np.shape(t)))

# Create signal
f_t = signal_samples(t)

# Compute FFT
F = fftpack.fft(f_t)
f = fftpack.fftfreq(N, 1.0/f_s)

# Plot spectrum for positive frequencies
mask = np.where(f >= 0)

fig, axes = plt.subplots(3, 1, figsize=(8, 6))
axes[0].plot(f[mask], np.log(abs(F[mask])))
axes[1].plot(f[mask], abs(F[mask])/N)
axes[1].set_xlim(0, 2)
axes[2].plot(f[mask], abs(F[mask])/N)
axes[2].set_xlim(21, 23)
plt.show()

Digital Filters

Design and apply a low-pass FIR filter to temperature data:

# Load temperature data
df = pd.read_csv('temperature_outdoor_2014.tsv', delimiter="\t",
                 names=["time", "temperature"])
df.time = (pd.to_datetime(df.time.values, unit="s")
           .tz_localize('UTC').tz_convert('Europe/Stockholm'))
df = df.set_index("time").resample("H").ffill()
temperature = df.temperature.values

# Design FIR filter
n = 101       # Filter length
f_s = 1 / 3600  # Sampling frequency (1 sample per hour)
nyq = f_s / 2   # Nyquist frequency

# Create low-pass filter
b = sig.firwin(n, cutoff=nyq/12, nyq=nyq, window="hamming")

# Apply filter to signal
temperature_filt = sig.lfilter(b, 1, temperature)

# Plot results
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(df.index, temperature, label="original", alpha=0.5)
ax.plot(df.index, temperature_filt, color="red", lw=2, label="FIR filtered")
ax.legend(loc=0)
plt.show()

Essential SciPy Imports

from scipy import linalg as la
from scipy import optimize
from scipy import interpolate
from scipy import integrate
from scipy import fftpack, signal as sig
from scipy.optimize import root, minimize

Module	Purpose	Common Functions
`scipy.linalg`	Linear algebra	`solve`, `inv`, `det`, `eig`
`scipy.optimize`	Optimization	`minimize`, `fsolve`, `bisect`, `brentq`
`scipy.interpolate`	Interpolation	`interp1d`, `interp2d`, `griddata`
`scipy.integrate`	Integration	`quad`, `dblquad`, `trapz`, `simps`
`scipy.signal`	Signal processing	`firwin`, `lfilter`, `fftconvolve`
`scipy.stats`	Statistics	`norm`, `t`, `chi2`, `f_oneway`

SciPy Best Practices

Performance Tips:

Provide derivatives when available to improve optimization convergence
Use appropriate solvers based on problem characteristics
Handle numerical precision appropriately for your domain
Choose optimal interpolation methods based on data regularity

Error Handling:

try:
    result = optimize.brentq(f, a, b)
except ValueError as e:
    print(f"Root finding failed: {e}")
    # Try alternative method
    result = optimize.newton(f, (a + b) / 2)

Validation:

# Always validate numerical results
def validate_solution(A, b, x):
    residual = np.dot(A, x) - b
    return np.max(np.abs(residual)) < 1e-10

SciPy provides robust scientific computing tools that complement NumPy's array operations. This cheatsheet covers the most important capabilities for mathematical computing, signal processing, and statistical analysis.

Updated: January 15, 2025
Author: Danial Pahlavan
Category: Scientific Computing & Data Science