SciPy – Scientific Computing

Advanced scientific computing library built on NumPy, providing algorithms for optimization, integration, interpolation, and more.

When to Use

• Integrating functions (numerical integration, ODEs)
• Optimizing functions (minimization, root finding, curve fitting)
• Interpolating data (1D, 2D, splines)
• Advanced linear algebra (sparse matrices, decompositions)
• Signal processing (filtering, Fourier transforms, wavelets)
• Statistical analysis (distributions, hypothesis tests)
• Image processing (filters, morphology, measurements)
• Spatial algorithms (distance matrices, clustering, Voronoi)
• Special mathematical functions (Bessel, gamma, error functions)
• Solving differential equations (ODEs, PDEs)

Reference Documentation

Official docs: https://docs.scipy.org/
Search patterns: scipy.integrate.quad, scipy.optimize.minimize, scipy.interpolate, scipy.stats, scipy.signal

Core Principles

Use SciPy For

Do NOT Use For

• Basic array operations (use NumPy)
• Machine learning (use scikit-learn)
• Deep learning (use PyTorch, TensorFlow)
• Symbolic mathematics (use SymPy)
• Data manipulation (use pandas)

Quick Reference

Installation

# pip
pip install scipy

# conda
conda install scipy

# With NumPy
pip install numpy scipy

Standard Imports

import numpy as np
from scipy import integrate, optimize, interpolate
from scipy import linalg, signal, stats
from scipy.integrate import odeint, solve_ivp
from scipy.optimize import minimize, root
from scipy.interpolate import interp1d, UnivariateSpline

Basic Pattern – Integration

from scipy import integrate
import numpy as np

# Define function
def f(x):
    return x**2

# Integrate from 0 to 1
result, error = integrate.quad(f, 0, 1)
print(f"Integral: {result:.6f} ± {error:.2e}")

Basic Pattern – Optimization

from scipy import optimize
import numpy as np

# Function to minimize
def f(x):
    return (x - 2)**2 + 1

# Minimize
result = optimize.minimize(f, x0=0)
print(f"Minimum at x = {result.x[0]:.6f}")
print(f"Minimum value = {result.fun:.6f}")

Basic Pattern – Interpolation

from scipy import interpolate
import numpy as np

# Data points
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

# Create interpolator
f = interpolate.interp1d(x, y, kind='cubic')

# Interpolate at new points
x_new = np.linspace(0, 4, 100)
y_new = f(x_new)

Critical Rules

✅ DO

• Check convergence – Always verify optimization converged
• Specify tolerances – Set appropriate rtol and atol
• Use appropriate methods – Choose algorithm for problem type
• Validate inputs – Check array shapes and values
• Handle edge cases – Deal with singularities and discontinuities
• Set integration limits carefully – Watch for infinite limits
• Use vectorization – Functions should accept arrays
• Check statistical assumptions – Verify distribution assumptions
• Specify degrees of freedom – For interpolation and fitting
• Use sparse matrices – For large, sparse systems

❌ DON’T

• Ignore convergence warnings – They indicate problems
• Use inappropriate tolerances – Too loose or too tight
• Apply wrong distribution – Check data characteristics
• Forget initial guesses – Optimization needs good starting points
• Integrate discontinuous functions – Without special handling
• Extrapolate beyond data – Interpolation is not extrapolation
• Mix incompatible units – Keep consistent units
• Ignore error estimates – They provide confidence levels
• Use wrong coordinate system – Check Cartesian vs polar
• Overfit with high-degree polynomials – Causes oscillations

Anti-Patterns (NEVER)

from scipy import integrate, optimize
import numpy as np

# ❌ BAD: Ignoring convergence
result = optimize.minimize(f, x0=0)
optimal_x = result.x  # Didn't check if converged!

# ✅ GOOD: Check convergence
result = optimize.minimize(f, x0=0)
if result.success:
    optimal_x = result.x
else:
    print(f"Optimization failed: {result.message}")

# ❌ BAD: Non-vectorized function for integration
def bad_func(x):
    if x < 0.5:
        return x
    else:
        return 1 - x

# ✅ GOOD: Vectorized function
def good_func(x):
    return np.where(x < 0.5, x, 1 - x)

# ❌ BAD: Poor initial guess
result = optimize.minimize(complex_func, x0=[1000, 1000])
# May converge to local minimum or fail!

# ✅ GOOD: Reasonable initial guess
x0 = np.array([0.0, 0.0])  # Near expected minimum
result = optimize.minimize(complex_func, x0=x0)

# ❌ BAD: Extrapolation with interpolation
f = interpolate.interp1d(x_data, y_data)
y_new = f(100)  # x_data max is 10, this will crash!

# ✅ GOOD: Check bounds or use extrapolation
f = interpolate.interp1d(x_data, y_data, fill_value='extrapolate')
y_new = f(100)  # Now works (but be cautious!)

# ❌ BAD: Wrong statistical test
# Using t-test for non-normal data
stats.ttest_ind(non_normal_data1, non_normal_data2)

# ✅ GOOD: Use appropriate test
# Use Mann-Whitney U for non-normal data
stats.mannwhitneyu(non_normal_data1, non_normal_data2)

Integration (scipy.integrate)

Numerical Integration (Quadrature)

from scipy import integrate
import numpy as np

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

result, error = integrate.quad(f, 0, np.inf)
print(f"∫exp(-x²)dx from 0 to ∞ = {result:.6f}")
print(f"Error estimate: {error:.2e}")

# Integral with parameters
def g(x, a, b):
    return a * x**2 + b

result, error = integrate.quad(g, 0, 1, args=(2, 3))
print(f"Result: {result:.6f}")

# Integral with singularity
def h(x):
    return 1 / np.sqrt(x)

# Specify singularity points
result, error = integrate.quad(h, 0, 1, points=[0])

Double and Triple Integrals

from scipy import integrate
import numpy as np

# Double integral: ∫∫ x*y dx dy over [0,1]×[0,2]
def f(y, x):  # Note: y first, x second
    return x * y

result, error = integrate.dblquad(f, 0, 1, 0, 2)
print(f"Double integral: {result:.6f}")

# Triple integral
def g(z, y, x):
    return x * y * z

result, error = integrate.tplquad(g, 0, 1, 0, 1, 0, 1)
print(f"Triple integral: {result:.6f}")

# Variable limits
def lower(x):
    return 0

def upper(x):
    return x

result, error = integrate.dblquad(f, 0, 1, lower, upper)

Solving ODEs

from scipy.integrate import odeint, solve_ivp
import numpy as np

# Solve dy/dt = -k*y (exponential decay)
def exponential_decay(y, t, k):
    return -k * y

# Initial condition and time points
y0 = 100
t = np.linspace(0, 10, 100)
k = 0.5

# Solve with odeint (older interface)
solution = odeint(exponential_decay, y0, t, args=(k,))

# Solve with solve_ivp (modern interface)
def decay_ivp(t, y, k):
    return -k * y

sol = solve_ivp(decay_ivp, [0, 10], [y0], args=(k,), t_eval=t)

print(f"Final value (odeint): {solution[-1, 0]:.6f}")
print(f"Final value (solve_ivp): {sol.y[0, -1]:.6f}")

System of ODEs

from scipy.integrate import solve_ivp
import numpy as np

# Lotka-Volterra equations (predator-prey)
def lotka_volterra(t, z, a, b, c, d):
    x, y = z
    dxdt = a*x - b*x*y
    dydt = -c*y + d*x*y
    return [dxdt, dydt]

# Parameters
a, b, c, d = 1.5, 1.0, 3.0, 1.0

# Initial conditions
z0 = [10, 5]  # [prey, predator]

# Time span
t_span = (0, 15)
t_eval = np.linspace(0, 15, 1000)

# Solve
sol = solve_ivp(lotka_volterra, t_span, z0, args=(a, b, c, d), 
                t_eval=t_eval, method='RK45')

prey = sol.y[0]
predator = sol.y[1]

print(f"Prey population at t=15: {prey[-1]:.2f}")
print(f"Predator population at t=15: {predator[-1]:.2f}")

Optimization (scipy.optimize)

Function Minimization

from scipy import optimize
import numpy as np

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

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

# Minimize with Nelder-Mead
result = optimize.minimize(rosenbrock, x0, method='Nelder-Mead')
print(f"Nelder-Mead: x = {result.x}, f(x) = {result.fun:.6f}")

# Minimize with BFGS (uses gradient)
result = optimize.minimize(rosenbrock, x0, method='BFGS')
print(f"BFGS: x = {result.x}, f(x) = {result.fun:.6f}")

# Minimize with bounds
bounds = [(0, 2), (0, 2)]
result = optimize.minimize(rosenbrock, x0, method='L-BFGS-B', bounds=bounds)
print(f"L-BFGS-B with bounds: x = {result.x}")

Root Finding

from scipy import optimize
import numpy as np

# Find roots of f(x) = 0
def f(x):
    return x**3 - 2*x - 5

# Root finding with scalar function
root = optimize.brentq(f, 0, 3)  # Search in [0, 3]
print(f"Root: {root:.6f}")

# Root finding with newton method (needs derivative)
def f_prime(x):
    return 3*x**2 - 2

root = optimize.newton(f, x0=2, fprime=f_prime)
print(f"Root (Newton): {root:.6f}")

# System of equations
def system(x):
    return [x[0]**2 + x[1]**2 - 1,
            x[0] - x[1]]

result = optimize.root(system, [1, 1])
print(f"Solution: {result.x}")

Curve Fitting

from scipy import optimize
import numpy as np

# Generate data with noise
x_data = np.linspace(0, 10, 50)
y_true = 2.5 * np.exp(-0.5 * x_data)
y_data = y_true + 0.2 * np.random.randn(len(x_data))

# Model function
def exponential_model(x, a, b):
    return a * np.exp(b * x)

# Fit model to data
params, covariance = optimize.curve_fit(exponential_model, x_data, y_data)
a_fit, b_fit = params

print(f"Fitted parameters: a = {a_fit:.4f}, b = {b_fit:.4f}")
print(f"True parameters: a = 2.5, b = -0.5")

# Standard errors
perr = np.sqrt(np.diag(covariance))
print(f"Parameter errors: ±{perr}")

Constrained Optimization

from scipy import optimize
import numpy as np

# Minimize f(x) = x[0]^2 + x[1]^2
# Subject to: x[0] + x[1] >= 1
def objective(x):
    return x[0]**2 + x[1]**2

# Constraint: g(x) >= 0
def constraint(x):
    return x[0] + x[1] - 1

con = {'type': 'ineq', 'fun': constraint}

# Minimize
x0 = np.array([2, 2])
result = optimize.minimize(objective, x0, constraints=con)

print(f"Optimal point: {result.x}")
print(f"Objective value: {result.fun:.6f}")

Interpolation (scipy.interpolate)

1D Interpolation

from scipy import interpolate
import numpy as np

# Data points
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([0, 0.8, 0.9, 0.1, -0.8, -1])

# Linear interpolation
f_linear = interpolate.interp1d(x, y, kind='linear')

# Cubic interpolation
f_cubic = interpolate.interp1d(x, y, kind='cubic')

# Evaluate at new points
x_new = np.linspace(0, 5, 100)
y_linear = f_linear(x_new)
y_cubic = f_cubic(x_new)

print(f"Value at x=2.5 (linear): {f_linear(2.5):.4f}")
print(f"Value at x=2.5 (cubic): {f_cubic(2.5):.4f}")

Spline Interpolation

from scipy import interpolate
import numpy as np

# Data
x = np.linspace(0, 10, 11)
y = np.sin(x)

# B-spline interpolation
tck = interpolate.splrep(x, y, s=0)  # s=0 for exact interpolation

# Evaluate
x_new = np.linspace(0, 10, 100)
y_new = interpolate.splev(x_new, tck)

# Or use UnivariateSpline
spl = interpolate.UnivariateSpline(x, y, s=0)
y_spl = spl(x_new)

# Smoothing spline (s > 0)
spl_smooth = interpolate.UnivariateSpline(x, y, s=1)
y_smooth = spl_smooth(x_new)

2D Interpolation

from scipy import interpolate
import numpy as np

# Create 2D data
x = np.linspace(0, 4, 5)
y = np.linspace(0, 4, 5)
X, Y = np.meshgrid(x, y)
Z = np.sin(X) * np.cos(Y)

# Flatten for irregular grid
x_flat = X.flatten()
y_flat = Y.flatten()
z_flat = Z.flatten()

# 2D interpolation
f = interpolate.interp2d(x_flat, y_flat, z_flat, kind='cubic')

# Evaluate at new points
x_new = np.linspace(0, 4, 50)
y_new = np.linspace(0, 4, 50)
Z_new = f(x_new, y_new)

print(f"Interpolated shape: {Z_new.shape}")

Linear Algebra (scipy.linalg)

Advanced Matrix Operations

from scipy import linalg
import numpy as np

# Matrix
A = np.array([[1, 2], [3, 4]])

# Matrix exponential
exp_A = linalg.expm(A)

# Matrix logarithm
log_A = linalg.logm(A)

# Matrix square root
sqrt_A = linalg.sqrtm(A)

# Matrix power
A_power_3 = linalg.fractional_matrix_power(A, 3)

print(f"exp(A):\n{exp_A}")
print(f"A^3:\n{A_power_3}")

Matrix Decompositions

from scipy import linalg
import numpy as np

A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 10]])

# LU decomposition
P, L, U = linalg.lu(A)
print(f"A = P @ L @ U: {np.allclose(A, P @ L @ U)}")

# QR decomposition
Q, R = linalg.qr(A)
print(f"A = Q @ R: {np.allclose(A, Q @ R)}")

# Cholesky decomposition (for positive definite)
A_pos_def = np.array([[4, 2], [2, 3]])
L = linalg.cholesky(A_pos_def, lower=True)
print(f"A = L @ L.T: {np.allclose(A_pos_def, L @ L.T)}")

# Schur decomposition
T, Z = linalg.schur(A)
print(f"A = Z @ T @ Z.H: {np.allclose(A, Z @ T @ Z.T.conj())}")

Sparse Matrices

from scipy import sparse
import numpy as np

# Create sparse matrix
row = np.array([0, 0, 1, 2, 2, 2])
col = np.array([0, 2, 2, 0, 1, 2])
data = np.array([1, 2, 3, 4, 5, 6])

# COO format (coordinate)
A_coo = sparse.coo_matrix((data, (row, col)), shape=(3, 3))

# Convert to CSR (efficient row operations)
A_csr = A_coo.tocsr()

# Convert to dense
A_dense = A_csr.toarray()
print(f"Dense matrix:\n{A_dense}")

# Sparse matrix operations
B = sparse.eye(3)  # Sparse identity
C = A_csr + B
D = A_csr @ B  # Matrix multiplication

print(f"Number of non-zeros: {A_csr.nnz}")
print(f"Sparsity: {1 - A_csr.nnz / (3*3):.2%}")

Signal Processing (scipy.signal)

Filter Design

from scipy import signal
import numpy as np

# Butterworth filter design
b, a = signal.butter(4, 0.5, btype='low')

# Apply filter to signal
t = np.linspace(0, 1, 1000)
sig = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*50*t)
filtered = signal.filtfilt(b, a, sig)

# Chebyshev filter
b_cheby, a_cheby = signal.cheby1(4, 0.5, 0.5)

# Frequency response
w, h = signal.freqz(b, a)
print(f"Filter order: {len(b)-1}")

Convolution and Correlation

from scipy import signal
import numpy as np

# Signals
sig1 = np.array([1, 2, 3, 4, 5])
sig2 = np.array([0, 1, 0.5])

# Convolution
conv = signal.convolve(sig1, sig2, mode='same')
print(f"Convolution: {conv}")

# Correlation
corr = signal.correlate(sig1, sig2, mode='same')
print(f"Correlation: {corr}")

# 2D convolution (for images)
image = np.random.rand(100, 100)
kernel = np.array([[1, 0, -1],
                   [2, 0, -2],
                   [1, 0, -1]]) / 4
filtered_img = signal.convolve2d(image, kernel, mode='same')

Peak Finding

from scipy import signal
import numpy as np

# Signal with peaks
t = np.linspace(0, 10, 1000)
sig = np.sin(2*np.pi*t) + 0.5*np.sin(2*np.pi*5*t)

# Find peaks
peaks, properties = signal.find_peaks(sig, height=0.5, distance=50)

print(f"Found {len(peaks)} peaks")
print(f"Peak positions: {peaks[:5]}")
print(f"Peak heights: {properties['peak_heights'][:5]}")

# Find local maxima with width
peaks_width, properties_width = signal.find_peaks(sig, width=20)
widths = properties_width['widths']

Spectral Analysis

from scipy import signal
import numpy as np

# Generate signal
fs = 1000  # Sampling frequency
t = np.linspace(0, 1, fs)
sig = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)

# Periodogram
f, Pxx = signal.periodogram(sig, fs)

# Welch's method (smoother estimate)
f_welch, Pxx_welch = signal.welch(sig, fs, nperseg=256)

# Spectrogram
f_spec, t_spec, Sxx = signal.spectrogram(sig, fs)

print(f"Frequency resolution: {f[1] - f[0]:.2f} Hz")
print(f"Number of frequency bins: {len(f)}")

Statistics (scipy.stats)

Probability Distributions

from scipy import stats
import numpy as np

# Normal distribution
mu, sigma = 0, 1
norm = stats.norm(loc=mu, scale=sigma)

# PDF, CDF, PPF
x = np.linspace(-3, 3, 100)
pdf = norm.pdf(x)
cdf = norm.cdf(x)
ppf = norm.ppf(0.975)  # 97.5th percentile

print(f"P(X ≤ 1.96) = {norm.cdf(1.96):.4f}")
print(f"97.5th percentile: {ppf:.4f}")

# Generate random samples
samples = norm.rvs(size=1000)

# Other distributions
exponential = stats.expon(scale=2)
poisson = stats.poisson(mu=5)
binomial = stats.binom(n=10, p=0.5)

Hypothesis Testing

from scipy import stats
import numpy as np

# Generate two samples
np.random.seed(42)
sample1 = stats.norm.rvs(loc=0, scale=1, size=100)
sample2 = stats.norm.rvs(loc=0.5, scale=1, size=100)

# t-test (independent samples)
t_stat, p_value = stats.ttest_ind(sample1, sample2)
print(f"t-test: t = {t_stat:.4f}, p = {p_value:.4f}")

# Mann-Whitney U test (non-parametric)
u_stat, p_value_mw = stats.mannwhitneyu(sample1, sample2)
print(f"Mann-Whitney: U = {u_stat:.4f}, p = {p_value_mw:.4f}")

# Chi-square test
observed = np.array([10, 20, 30, 40])
expected = np.array([25, 25, 25, 25])
chi2, p_chi = stats.chisquare(observed, expected)
print(f"Chi-square: χ² = {chi2:.4f}, p = {p_chi:.4f}")

# Kolmogorov-Smirnov test (normality)
ks_stat, p_ks = stats.kstest(sample1, 'norm')
print(f"KS test: D = {ks_stat:.4f}, p = {p_ks:.4f}")

Correlation and Regression

from scipy import stats
import numpy as np

# Data
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = 2*x + 1 + np.random.randn(10)*0.5

# Pearson correlation
r, p_value = stats.pearsonr(x, y)
print(f"Pearson r = {r:.4f}, p = {p_value:.4f}")

# Spearman correlation (rank-based)
rho, p_spear = stats.spearmanr(x, y)
print(f"Spearman ρ = {rho:.4f}, p = {p_spear:.4f}")

# Linear regression
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
print(f"y = {slope:.4f}x + {intercept:.4f}")
print(f"R² = {r_value**2:.4f}")

Fast Fourier Transform (scipy.fft)

Basic FFT

from scipy import fft
import numpy as np

# Time-domain signal
fs = 1000  # Sampling rate
T = 1/fs
N = 1000   # Number of samples
t = np.linspace(0, N*T, N)

# Signal: 50 Hz + 120 Hz
signal = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t)

# Compute FFT
yf = fft.fft(signal)
xf = fft.fftfreq(N, T)

# Magnitude spectrum
magnitude = np.abs(yf)

# Get positive frequencies only
positive_freq_idx = xf > 0
xf_positive = xf[positive_freq_idx]
magnitude_positive = magnitude[positive_freq_idx]

print(f"Peak frequencies: {xf_positive[magnitude_positive > 100]}")

Inverse FFT

from scipy import fft
import numpy as np

# Original signal
t = np.linspace(0, 1, 1000)
signal = np.sin(2*np.pi*10*t)

# FFT
signal_fft = fft.fft(signal)

# Modify in frequency domain (e.g., remove high frequencies)
signal_fft[100:] = 0

# Inverse FFT
signal_filtered = fft.ifft(signal_fft).real

print(f"Signal reconstructed: {np.allclose(signal[:100], signal_filtered[:100])}")

2D FFT (for images)

from scipy import fft
import numpy as np

# Create 2D signal (image)
x = np.linspace(0, 2*np.pi, 128)
y = np.linspace(0, 2*np.pi, 128)
X, Y = np.meshgrid(x, y)
image = np.sin(5*X) * np.cos(5*Y)

# 2D FFT
image_fft = fft.fft2(image)
image_fft_shifted = fft.fftshift(image_fft)

# Magnitude spectrum
magnitude = np.abs(image_fft_shifted)

# Inverse 2D FFT
image_reconstructed = fft.ifft2(image_fft).real

print(f"Image shape: {image.shape}")
print(f"FFT shape: {image_fft.shape}")

Spatial Algorithms (scipy.spatial)

Distance Computations

from scipy.spatial import distance
import numpy as np

# Two points
p1 = np.array([0, 0])
p2 = np.array([3, 4])

# Euclidean distance
eucl = distance.euclidean(p1, p2)
print(f"Euclidean distance: {eucl:.4f}")

# Manhattan distance
manh = distance.cityblock(p1, p2)
print(f"Manhattan distance: {manh:.4f}")

# Cosine distance
cosine = distance.cosine(p1, p2)

# Pairwise distances
points = np.random.rand(10, 2)
dist_matrix = distance.cdist(points, points, 'euclidean')
print(f"Distance matrix shape: {dist_matrix.shape}")

Convex Hull

from scipy.spatial import ConvexHull
import numpy as np

# Random points
points = np.random.rand(30, 2)

# Compute convex hull
hull = ConvexHull(points)

print(f"Number of vertices: {len(hull.vertices)}")
print(f"Hull area: {hull.area:.4f}")
print(f"Hull volume (perimeter): {hull.volume:.4f}")

# Get hull points
hull_points = points[hull.vertices]

Delaunay Triangulation

from scipy.spatial import Delaunay
import numpy as np

# Points
points = np.random.rand(30, 2)

# Triangulation
tri = Delaunay(points)

print(f"Number of triangles: {len(tri.simplices)}")

# Check if point is in triangulation
test_point = np.array([0.5, 0.5])
simplex_index = tri.find_simplex(test_point)
print(f"Point inside: {simplex_index >= 0}")

KDTree for Nearest Neighbors

from scipy.spatial import KDTree
import numpy as np

# Build tree
points = np.random.rand(100, 3)
tree = KDTree(points)

# Query nearest neighbors
query_point = np.array([0.5, 0.5, 0.5])
distances, indices = tree.query(query_point, k=5)

print(f"5 nearest neighbors:")
print(f"Distances: {distances}")
print(f"Indices: {indices}")

# Query within radius
indices_radius = tree.query_ball_point(query_point, r=0.3)
print(f"Points within r=0.3: {len(indices_radius)}")

Practical Workflows

Numerical Integration of Physical System

from scipy.integrate import odeint
import numpy as np

# Damped harmonic oscillator: m*x'' + c*x' + k*x = 0
def damped_oscillator(y, t, m, c, k):
    x, v = y
    dxdt = v
    dvdt = -(c/m)*v - (k/m)*x
    return [dxdt, dvdt]

# Parameters
m = 1.0  # mass
c = 0.5  # damping
k = 10.0  # spring constant

# Initial conditions
y0 = [1.0, 0.0]  # [position, velocity]

# Time points
t = np.linspace(0, 10, 1000)

# Solve
solution = odeint(damped_oscillator, y0, t, args=(m, c, k))

position = solution[:, 0]
velocity = solution[:, 1]

print(f"Final position: {position[-1]:.6f}")
print(f"Final velocity: {velocity[-1]:.6f}")

Parameter Estimation from Data

from scipy import optimize
import numpy as np

# Generate synthetic data
x_true = np.linspace(0, 10, 50)
params_true = [2.5, 1.3, 0.8]
y_true = params_true[0] * np.exp(-params_true[1] * x_true) + params_true[2]
y_data = y_true + 0.2 * np.random.randn(len(x_true))

# Model
def model(x, a, b, c):
    return a * np.exp(-b * x) + c

# Objective function (residual sum of squares)
def objective(params):
    y_pred = model(x_true, *params)
    return np.sum((y_data - y_pred)**2)

# Optimize
params_init = [1.0, 1.0, 1.0]
result = optimize.minimize(objective, params_init)

print(f"True parameters: {params_true}")
print(f"Estimated parameters: {result.x}")
print(f"Relative errors: {np.abs(result.x - params_true) / params_true * 100}%")

Signal Filtering Pipeline

from scipy import signal
import numpy as np

def filter_pipeline(noisy_signal, fs):
    """Complete signal processing pipeline."""
    # 1. Design Butterworth lowpass filter
    fc = 10  # Cutoff frequency
    w = fc / (fs / 2)  # Normalized frequency
    b, a = signal.butter(4, w, 'low')
    
    # 2. Apply filter (zero-phase)
    filtered = signal.filtfilt(b, a, noisy_signal)
    
    # 3. Remove baseline drift
    baseline = signal.medfilt(filtered, kernel_size=51)
    detrended = filtered - baseline
    
    # 4. Normalize
    normalized = (detrended - np.mean(detrended)) / np.std(detrended)
    
    return normalized

# Example usage
fs = 1000
t = np.linspace(0, 1, fs)
clean_signal = np.sin(2*np.pi*5*t)
noise = 0.5 * np.random.randn(len(t))
drift = 0.1 * t
noisy_signal = clean_signal + noise + drift

processed = filter_pipeline(noisy_signal, fs)
print(f"Original SNR: {10*np.log10(np.var(clean_signal)/np.var(noise)):.2f} dB")

Interpolation and Smoothing

from scipy import interpolate
import numpy as np

# Noisy data
x = np.linspace(0, 10, 20)
y_true = np.sin(x)
y_noisy = y_true + 0.3 * np.random.randn(len(x))

# Smoothing spline
spl = interpolate.UnivariateSpline(x, y_noisy, s=2)

# Evaluate on fine grid
x_fine = np.linspace(0, 10, 200)
y_smooth = spl(x_fine)

# Compare with true function
y_true_fine = np.sin(x_fine)
rmse = np.sqrt(np.mean((y_smooth - y_true_fine)**2))

print(f"RMSE: {rmse:.6f}")

# Can also get derivatives
y_prime = spl.derivative()(x_fine)
y_double_prime = spl.derivative(n=2)(x_fine)

Statistical Analysis Workflow

from scipy import stats
import numpy as np

def analyze_experiment(control, treatment):
    """Complete statistical analysis of experiment data."""
    results = {}
    
    # 1. Descriptive statistics
    results['control_mean'] = np.mean(control)
    results['treatment_mean'] = np.mean(treatment)
    results['control_std'] = np.std(control, ddof=1)
    results['treatment_std'] = np.std(treatment, ddof=1)
    
    # 2. Normality tests
    _, results['control_normal_p'] = stats.shapiro(control)
    _, results['treatment_normal_p'] = stats.shapiro(treatment)
    
    # 3. Choose appropriate test
    if results['control_normal_p'] > 0.05 and results['treatment_normal_p'] > 0.05:
        # Both normal: use t-test
        t_stat, p_value = stats.ttest_ind(control, treatment)
        results['test'] = 't-test'
        results['statistic'] = t_stat
        results['p_value'] = p_value
    else:
        # Non-normal: use Mann-Whitney U
        u_stat, p_value = stats.mannwhitneyu(control, treatment)
        results['test'] = 'Mann-Whitney U'
        results['statistic'] = u_stat
        results['p_value'] = p_value
    
    # 4. Effect size (Cohen's d)
    pooled_std = np.sqrt((results['control_std']**2 + 
                          results['treatment_std']**2) / 2)
    results['cohens_d'] = ((results['treatment_mean'] - 
                            results['control_mean']) / pooled_std)
    
    return results

# Example
control = stats.norm.rvs(loc=10, scale=2, size=30)
treatment = stats.norm.rvs(loc=12, scale=2, size=30)

results = analyze_experiment(control, treatment)
print(f"Test: {results['test']}")
print(f"p-value: {results['p_value']:.4f}")
print(f"Effect size: {results['cohens_d']:.4f}")

Performance Optimization

Choosing the Right Method

from scipy import integrate, optimize
import numpy as np
import time

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

# Compare integration methods
methods = ['quad', 'romberg', 'simpson']
times = []

for method in methods:
    start = time.time()
    if method == 'quad':
        result, _ = integrate.quad(f, 0, 10)
    elif method == 'romberg':
        x = np.linspace(0, 10, 1000)
        result = integrate.romberg(f, 0, 10)
    elif method == 'simpson':
        x = np.linspace(0, 10, 1000)
        y = f(x)
        result = integrate.simpson(y, x)
    elapsed = time.time() - start
    times.append(elapsed)
    print(f"{method}: {result:.8f} ({elapsed*1000:.2f} ms)")

Vectorization

from scipy import interpolate
import numpy as np

# Bad: Loop over points
def interpolate_loop(x, y, x_new):
    results = []
    for xi in x_new:
        # Create interpolator for each point (very slow!)
        f = interpolate.interp1d(x, y)
        results.append(f(xi))
    return np.array(results)

# Good: Vectorized
def interpolate_vectorized(x, y, x_new):
    f = interpolate.interp1d(x, y)
    return f(x_new)  # Pass entire array

# Benchmark
x = np.linspace(0, 10, 100)
y = np.sin(x)
x_new = np.linspace(0, 10, 1000)

# Only use vectorized version (loop version is too slow)
result = interpolate_vectorized(x, y, x_new)

Common Pitfalls and Solutions

Integration Convergence

from scipy import integrate
import numpy as np

# Problem: Oscillatory function
def oscillatory(x):
    return np.sin(100*x) / x if x != 0 else 100

# Bad: Default parameters may not converge well
result, error = integrate.quad(oscillatory, 0, 1)

# Good: Increase limit and tolerance
result, error = integrate.quad(oscillatory, 0, 1, 
                               limit=100, epsabs=1e-10, epsrel=1e-10)

# Good: For highly oscillatory, use special methods
from scipy.integrate import quad_vec
result, error = quad_vec(oscillatory, 0, 1)

Optimization Local Minima

from scipy import optimize
import numpy as np

# Function with multiple minima
def multi_minima(x):
    return np.sin(x) + np.sin(10*x/3)

# Bad: Single starting point may find local minimum
result = optimize.minimize(multi_minima, x0=0)
print(f"Local minimum: {result.x[0]:.4f}")

# Good: Try multiple starting points
x0_list = np.linspace(0, 10, 20)
results = [optimize.minimize(multi_minima, x0=x0) for x0 in x0_list]
global_min = min(results, key=lambda r: r.fun)
print(f"Global minimum: {global_min.x[0]:.4f}")

# Best: Use global optimization
result_global = optimize.differential_evolution(multi_minima, bounds=[(0, 10)])
print(f"Global (DE): {result_global.x[0]:.4f}")

Statistical Test Assumptions

from scipy import stats
import numpy as np

# Generate non-normal data
non_normal = stats.expon.rvs(size=30)

# Bad: Using t-test without checking normality
t_stat, p_value = stats.ttest_1samp(non_normal, 1.0)

# Good: Check normality first
_, p_normal = stats.shapiro(non_normal)
if p_normal < 0.05:
    print("Data is not normal, using Wilcoxon test")
    stat, p_value = stats.wilcoxon(non_normal - 1.0)
else:
    t_stat, p_value = stats.ttest_1samp(non_normal, 1.0)

Sparse Matrix Efficiency

from scipy import sparse
import numpy as np
import time

# Large sparse matrix
n = 10000
density = 0.01
data = np.random.rand(int(n*n*density))
row = np.random.randint(0, n, len(data))
col = np.random.randint(0, n, len(data))

# Bad: Convert to dense
A_coo = sparse.coo_matrix((data, (row, col)), shape=(n, n))
# A_dense = A_coo.toarray()  # Uses huge memory!

# Good: Keep sparse and use appropriate format
A_csr = A_coo.tocsr()  # Efficient for row operations
A_csc = A_coo.tocsc()  # Efficient for column operations

# Sparse operations
x = np.random.rand(n)
y = A_csr @ x  # Fast sparse matrix-vector product

print(f"Sparse matrix: {A_csr.nnz / (n*n) * 100:.2f}% non-zero")

This comprehensive SciPy guide covers 50+ examples across all major scientific computing workflows!