SciPy: scientific algorithms for Python

NumPy provides arrays and fast array operations. SciPy builds on NumPy and supplies numerical algorithms for integration, differential equations, special functions, Fourier transforms, optimization, root finding, interpolation, sparse matrices, statistics, and more.

This lecture concentrates on a small set of tools that recur throughout computational physics.

import numpy as np
import matplotlib.pyplot as plt

from scipy import integrate, optimize, special
from scipy.fft import rfft, irfft, rfftfreq
from scipy.interpolate import CubicSpline, PchipInterpolator
from scipy.signal.windows import hann

%matplotlib inline

Import SciPy by submodule

SciPy is organized into submodules such as scipy.integrate, scipy.optimize, and scipy.special. Explicit imports keep the origin of each algorithm visible and avoid overwriting NumPy or Python names.

The numerical arrays passed to and returned by SciPy are ordinarily NumPy arrays.

1. Special functions

Special functions arise as solutions of differential equations with particular geometries or boundary conditions. For example, Bessel functions appear in cylindrical problems such as vibrating circular membranes, cylindrical waveguides, and scattering from circular objects.

x = np.linspace(0.0, 15.0, 500)

fig, ax = plt.subplots(figsize=(8, 4))
for n in range(4):
    ax.plot(x, special.jv(n, x), label=rf"$J_{n}(x)$")

ax.axhline(0.0, color="0.7", linewidth=0.8)
ax.set(xlabel="$x$", ylabel="$J_n(x)$", title="Bessel functions of the first kind")
ax.legend()
plt.show()

Boundary conditions often select the zeros of a special function. For a circular membrane fixed at its edge, the allowed radial wave numbers are proportional to zeros of $J_n$.

for n in range(3):
    zeros = special.jn_zeros(n, 4)
    print(f"first four zeros of J_{n}: {zeros}")

first four zeros of J_0: [ 2.40482556  5.52007811  8.65372791 11.79153444]
first four zeros of J_1: [ 3.83170597  7.01558667 10.17346814 13.32369194]
first four zeros of J_2: [ 5.1356223   8.41724414 11.61984117 14.79595178]

Spherical harmonics

Spherical harmonics $Y_l^m(\theta,\phi)$ describe angular dependence in central-potential quantum mechanics, multipole expansions, radiation patterns, and crystal-field theory.

SciPy uses

• theta for the polar angle in $[0,\pi]$;
• phi for the azimuthal angle in $[0,2\pi)$.

The call is special.sph_harm_y(l, m, theta, phi).

theta = np.linspace(0.0, np.pi, 100)
phi = np.linspace(0.0, 2*np.pi, 200, endpoint=False)
theta_grid, phi_grid = np.meshgrid(theta, phi, indexing="ij")

l, m = 3, 2
Ylm = special.sph_harm_y(l, m, theta_grid, phi_grid)

# Use |Y_l^m| as the radius and Re(Y_l^m) as the color.
radius = np.abs(Ylm)
x = radius * np.sin(theta_grid) * np.cos(phi_grid)
y = radius * np.sin(theta_grid) * np.sin(phi_grid)
z = radius * np.cos(theta_grid)

colors = plt.cm.RdBu_r(
    plt.Normalize(vmin=-np.max(np.abs(Ylm.real)), vmax=np.max(np.abs(Ylm.real)))(Ylm.real)
)

fig = plt.figure(figsize=(7, 6))
ax = fig.add_subplot(projection="3d")
ax.plot_surface(x, y, z, facecolors=colors, linewidth=0, antialiased=True)
ax.set_box_aspect((1, 1, 1))
ax.set_title(rf"$Y_{{{l}}}^{{{m}}}$: radius $|Y|$, color $\mathrm{{Re}}(Y)$")
ax.set_axis_off()
plt.show()

2. Numerical integration

integrate.quad computes

$$ I=\int_a^b f(x)\,dx $$

using adaptive quadrature. It returns both the numerical result and an estimate of the absolute numerical error. The estimate is not the difference from the exact answer, which is normally unknown.

value, estimated_error = integrate.quad(lambda x: x**2, 0.0, 1.0)

print(f"integral       = {value:.16f}")
print(f"exact value    = {1/3:.16f}")
print(f"estimated error= {estimated_error:.2e}")
print(f"actual error   = {abs(value - 1/3):.2e}")

integral       = 0.3333333333333333
exact value    = 0.3333333333333333
estimated error= 3.70e-15
actual error   = 0.00e+00

Parameters can be passed to the integrand with args. Here we verify the Gaussian integral

$$ \int_{-\infty}^{\infty}e^{-a x^2}\,dx=\sqrt{\frac{\pi}{a}}. $$

def gaussian(x, a):
    return np.exp(-a*x**2)

a = 2.5
value, estimated_error = integrate.quad(gaussian, -np.inf, np.inf, args=(a,))
exact = np.sqrt(np.pi/a)

print(f"numerical = {value:.12f}")
print(f"exact     = {exact:.12f}")
print(f"estimated absolute error = {estimated_error:.2e}")

numerical = 1.120998243280
exact     = 1.120998243280
estimated absolute error = 5.02e-11

For a two-dimensional integral, dblquad expects the integrand as f(y, x). The limits of the inner integral may depend on the outer variable. For the triangular region $0\le x\le1$ and $0\le y\le1-x$,

$$ \int_0^1 dx\int_0^{1-x}dy\,(x+y)=\frac13. $$

value, estimated_error = integrate.dblquad(
    lambda y, x: x + y,
    0.0, 1.0,
    lambda x: 0.0,
    lambda x: 1.0 - x,
)

print(f"integral = {value:.12f}, estimated error = {estimated_error:.2e}")

integral = 0.333333333333, estimated error = 5.55e-15

Adaptive quadrature is powerful, but not magic. Singularities, discontinuities, rapid oscillations, and slowly decaying tails should be communicated to the algorithm or treated analytically. Always inspect warnings and test convergence by changing tolerances or reformulating the integral.

3. Ordinary differential equations

integrate.solve_ivp solves an initial-value problem

$$ \frac{d\mathbf y}{dt}=\mathbf f(t,\mathbf y),\qquad \mathbf y(t_0)=\mathbf y_0. $$

A higher-order equation is first written as a system of first-order equations. The solver chooses internal time steps adaptively; t_eval only requests the times at which the solution is returned.

Double compound pendulum: complete derivation

Consider two identical uniform rods, each with mass $m$ and length $\ell$. The first is pivoted at one end, and the second is hinged to the end of the first. Their angles from the downward vertical are $\theta_1$ and $\theta_2$.

The centers of mass are

$$ \mathbf r_1= \begin{bmatrix} \frac{\ell}{2}\sin\theta_1\\ -\frac{\ell}{2}\cos\theta_1 \end{bmatrix}, \qquad \mathbf r_2= \begin{bmatrix} \ell\sin\theta_1+\frac{\ell}{2}\sin\theta_2\\ -\ell\cos\theta_1-\frac{\ell}{2}\cos\theta_2 \end{bmatrix}. $$

For a uniform rod about its center of mass, $I=m\ell^2/12$.

The translational and rotational kinetic energies give

$$ T=\frac12m\ell^2\left[ \frac43\dot\theta_1^2+ \frac13\dot\theta_2^2+ \cos(\theta_1-\theta_2)\dot\theta_1\dot\theta_2 \right]. $$

Taking zero gravitational potential at the pivots,

$$ V=-mg\ell\left(\frac32\cos\theta_1+\frac12\cos\theta_2\right), $$

and therefore

$$ L=T-V= \frac12m\ell^2\left[ \frac43\dot\theta_1^2+ \frac13\dot\theta_2^2+ \cos(\theta_1-\theta_2)\dot\theta_1\dot\theta_2 \right] +mg\ell\left(\frac32\cos\theta_1+\frac12\cos\theta_2\right). $$

The generalized momenta are

$$ p_1=m\ell^2\left[\frac43\dot\theta_1+ \frac12\cos(\theta_1-\theta_2)\dot\theta_2\right], $$

$$ p_2=m\ell^2\left[\frac13\dot\theta_2+ \frac12\cos(\theta_1-\theta_2)\dot\theta_1\right]. $$

Inverting these two equations gives

$$ \dot\theta_1=\frac{6}{m\ell^2} \frac{2p_1-3\cos(\theta_1-\theta_2)p_2} {16-9\cos^2(\theta_1-\theta_2)}, $$

$$ \dot\theta_2=\frac{6}{m\ell^2} \frac{8p_2-3\cos(\theta_1-\theta_2)p_1} {16-9\cos^2(\theta_1-\theta_2)}. $$

The Euler–Lagrange equations are

$$ \frac{d}{dt}\frac{\partial L}{\partial\dot\theta_i} =\frac{\partial L}{\partial\theta_i}. $$

For the angular derivatives of the Lagrangian,

$$ \frac{\partial L}{\partial\theta_1} =-\frac12m\ell^2\sin(\theta_1-\theta_2)\dot\theta_1\dot\theta_2 -\frac32mg\ell\sin\theta_1, $$

$$ \frac{\partial L}{\partial\theta_2} =\frac12m\ell^2\sin(\theta_1-\theta_2)\dot\theta_1\dot\theta_2 -\frac12mg\ell\sin\theta_2. $$

Since $p_i=\partial L/\partial\dot\theta_i$, the remaining two first-order equations are

$$ \dot p_1=-\frac12m\ell^2\left[ \dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2) +3\frac{g}{\ell}\sin\theta_1\right], $$

$$ \dot p_2=-\frac12m\ell^2\left[ -\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2) +\frac{g}{\ell}\sin\theta_2\right]. $$

Thus the state vector $\mathbf y=(\theta_1,\theta_2,p_1,p_2)$ obeys four coupled first-order equations, exactly the form required by solve_ivp.

def double_pendulum_rhs(t, state, m, ell, g):
    '''Time derivative of (theta1, theta2, p1, p2).'''
    theta1, theta2, p1, p2 = state
    delta = theta1 - theta2
    c = np.cos(delta)
    denominator = 16.0 - 9.0*c**2
    prefactor = 6.0/(m*ell**2*denominator)

    theta1_dot = prefactor*(2.0*p1 - 3.0*c*p2)
    theta2_dot = prefactor*(8.0*p2 - 3.0*c*p1)
    coupling = theta1_dot*theta2_dot*np.sin(delta)

    p1_dot = -0.5*m*ell**2*(coupling + 3.0*g/ell*np.sin(theta1))
    p2_dot = -0.5*m*ell**2*(-coupling + g/ell*np.sin(theta2))

    return np.array([theta1_dot, theta2_dot, p1_dot, p2_dot])

m, ell, g = 1.0, 1.0, 9.81
initial_state = np.array([np.pi - 0.1, -np.pi/2, 0.0, 0.0])
t_eval = np.linspace(0.0, 20.0, 801)

sol = integrate.solve_ivp(
    double_pendulum_rhs,
    t_span=(t_eval[0], t_eval[-1]),
    y0=initial_state,
    t_eval=t_eval,
    args=(m, ell, g),
    method="DOP853",
    rtol=1e-9,
    atol=1e-11,
)

print("success:", sol.success)
print("message:", sol.message)
print("right-hand-side evaluations:", sol.nfev)
print("solution shape:", sol.y.shape)

success: True
message: The solver successfully reached the end of the integration interval.
right-hand-side evaluations: 14498
solution shape: (4, 801)

sol.y has shape (number of variables, number of requested times). Solver success should always be checked before using the trajectory.

theta1, theta2, p1, p2 = sol.y

x1 = 0.5*ell*np.sin(theta1)
y1 = -0.5*ell*np.cos(theta1)
x_joint = ell*np.sin(theta1)
y_joint = -ell*np.cos(theta1)
x2 = x_joint + 0.5*ell*np.sin(theta2)
y2 = y_joint - 0.5*ell*np.cos(theta2)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(sol.t, theta1, label=r"$\theta_1$")
axes[0].plot(sol.t, theta2, label=r"$\theta_2$")
axes[0].set(xlabel="time", ylabel="angle (rad)")
axes[0].legend()

axes[1].plot(x1, y1, label="center of rod 1")
axes[1].plot(x2, y2, label="center of rod 2")
axes[1].set(xlabel="$x$", ylabel="$y$", aspect="equal")
axes[1].legend()
plt.show()

Animation

The numerical solution supplies the angles at each requested time. The animation converts those angles into the positions of the pivot, hinge, and free end. We display every fourth calculated frame to keep the notebook output manageable.

from matplotlib.animation import FuncAnimation
from IPython.display import HTML

x_end = x_joint + ell*np.sin(theta2)
y_end = y_joint - ell*np.cos(theta2)
frame_indices = np.arange(0, sol.t.size, 4)

fig, ax = plt.subplots(figsize=(5, 5))
ax.set(xlim=(-2.1*ell, 2.1*ell), ylim=(-2.1*ell, 2.1*ell), aspect="equal")
ax.grid()
line, = ax.plot([], [], "o-", linewidth=2)
time_text = ax.text(0.05, 0.92, "", transform=ax.transAxes)

def update(frame):
    i = frame_indices[frame]
    line.set_data([0.0, x_joint[i], x_end[i]], [0.0, y_joint[i], y_end[i]])
    time_text.set_text(f"t = {sol.t[i]:.2f} s")
    return line, time_text

animation = FuncAnimation(
    fig, update, frames=len(frame_indices), interval=75, blit=True
)
plt.close(fig)
HTML(animation.to_jshtml())

A physical diagnostic: energy conservation

The exact system is conservative. Energy drift therefore provides a useful numerical diagnostic. This is more informative than checking only whether the solver returned without an exception.

def double_pendulum_energy(state, m, ell, g):
    theta1, theta2, p1, p2 = state
    derivatives = np.column_stack([
        double_pendulum_rhs(0.0, state[:, i], m, ell, g)
        for i in range(state.shape[1])
    ])
    theta1_dot, theta2_dot = derivatives[0], derivatives[1]
    delta = theta1 - theta2

    kinetic = 0.5*m*ell**2*(
        (4/3)*theta1_dot**2
        + (1/3)*theta2_dot**2
        + np.cos(delta)*theta1_dot*theta2_dot
    )
    potential = -m*g*ell*(1.5*np.cos(theta1) + 0.5*np.cos(theta2))
    return kinetic + potential

energy = double_pendulum_energy(sol.y, m, ell, g)
relative_drift = (energy - energy[0])/max(abs(energy[0]), 1.0)

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.plot(sol.t, relative_drift)
ax.set(xlabel="time", ylabel="relative energy change")
ax.ticklabel_format(axis="y", style="sci", scilimits=(0, 0))
plt.show()

print("maximum relative energy drift:", np.max(np.abs(relative_drift)))

maximum relative energy drift: 9.647012183291229e-09

Smaller tolerances normally reduce the error but require more function evaluations. Adaptive high-order solvers are excellent general-purpose tools, although they do not preserve Hamiltonian structure exactly over arbitrarily long times.

4. Fast Fourier transform

For samples $x_j=x(t_j)$ with $t_j=j\Delta t$, the discrete Fourier transform is

$$ X_k=\sum_{j=0}^{N-1}x_j e^{-2\pi i jk/N}. $$

Direct evaluation takes $O(N^2)$ operations. The FFT exploits symmetries of the exponential factors to evaluate the same transform in $O(N\log N)$ operations.

For a real signal, rfft stores only nonnegative frequencies. rfftfreq returns frequencies in cycles per unit time; multiply by $2\pi$ when angular frequency is required.

sampling_rate = 200.0                 # samples per second
duration = 2.0                        # seconds
n_samples = int(sampling_rate*duration)
t = np.arange(n_samples)/sampling_rate

signal = 1.2*np.cos(2*np.pi*5.0*t) + 0.7*np.sin(2*np.pi*12.5*t)
spectrum = rfft(signal)
frequency = rfftfreq(n_samples, d=1/sampling_rate)
amplitude = 2.0*np.abs(spectrum)/n_samples
amplitude[0] /= 2.0                   # DC is not paired with a negative frequency

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(t, signal)
axes[0].set(xlim=(0, 0.5), xlabel="time (s)", ylabel="signal")
axes[1].stem(frequency, amplitude, basefmt=" ")
axes[1].set(xlim=(0, 30), xlabel="frequency (Hz)", ylabel="amplitude")
plt.show()

For a record of duration $T=N\Delta t$:

• the frequency spacing is $\Delta f=1/T$;
• the highest representable frequency is the Nyquist frequency $f_{\rm Nyquist}=1/(2\Delta t)$;
• frequencies above Nyquist are aliased to lower frequencies.

The inverse transform reconstructs the sampled data, up to floating-point roundoff.

reconstructed = irfft(spectrum, n=n_samples)

print(f"frequency resolution = {frequency[1] - frequency[0]:.3f} Hz")
print(f"Nyquist frequency    = {sampling_rate/2:.1f} Hz")
print(f"maximum reconstruction error = {np.max(np.abs(reconstructed-signal)):.2e}")

frequency resolution = 0.500 Hz
Nyquist frequency    = 100.0 Hz
maximum reconstruction error = 8.88e-16

Spectral leakage

The DFT treats the finite data record as one period of a periodic signal. If the record does not contain an integer number of oscillations, joining repeated copies creates a discontinuity and spreads power over many frequency bins. A window reduces this leakage at the cost of broadening spectral peaks.

sampling_rate = 200.0
duration = 1.0
n_samples = int(sampling_rate*duration)
t = np.arange(n_samples)/sampling_rate
frequency = rfftfreq(n_samples, 1/sampling_rate)

signal_periodic = np.sin(2*np.pi*12.0*t)     # exactly 12 cycles
signal_leaking = np.sin(2*np.pi*12.3*t)     # not an integer number of cycles
window = hann(n_samples, sym=False)

amp_periodic = 2*np.abs(rfft(signal_periodic))/n_samples
amp_leaking = 2*np.abs(rfft(signal_leaking))/n_samples
amp_windowed = 2*np.abs(rfft(window*signal_leaking))/np.sum(window)

fig, ax = plt.subplots(figsize=(9, 4))
ax.semilogy(frequency, np.maximum(amp_periodic, 1e-8), label="12.0 Hz: integer cycles")
ax.semilogy(frequency, np.maximum(amp_leaking, 1e-8), label="12.3 Hz: leakage")
ax.semilogy(frequency, np.maximum(amp_windowed, 1e-8), label="12.3 Hz: Hann window")
ax.set(xlim=(0, 35), ylim=(1e-6, 2), xlabel="frequency (Hz)", ylabel="amplitude")
ax.legend()
plt.show()

Zero-padding can make a plotted spectrum look smoother, but it does not add information or improve the fundamental resolution $1/T$. A longer observation time is needed for finer frequency resolution.

5. Optimization

An optimizer returns a result object containing more than the proposed minimum. Inspect success, message, the location x, and the function value fun.

def potential(x):
    return 0.25*x**4 - 0.5*x**2 + 0.12*x

result = optimize.minimize_scalar(potential, bounds=(-2.0, 0.0), method="bounded")

x_plot = np.linspace(-2.0, 2.0, 500)
plt.plot(x_plot, potential(x_plot))
plt.plot(result.x, result.fun, "ro", label="minimum in the chosen interval")
plt.xlabel("$x$")
plt.ylabel("$V(x)$")
plt.legend()
plt.show()

print("success:", result.success)
print("x_min:", result.x)
print("V(x_min):", result.fun)

success: True
x_min: -1.0553231462371546
V(x_min): -0.3734064603277733

The bounds are part of the mathematical problem, not merely numerical decoration. This tilted double-well potential has more than one local minimum; changing the interval can select a different one.

6. Root finding

Bracketed methods require an interval whose endpoints have opposite signs. They are robust for a continuous function because the intermediate-value theorem guarantees a root inside the bracket. Open methods use one or more initial guesses and can be faster, but convergence depends more strongly on those guesses.

def root_function(x):
    return np.cos(x) - x

bracketed = optimize.root_scalar(root_function, bracket=(0.0, 1.0), method="brentq")
newton = optimize.root_scalar(
    root_function,
    x0=0.7,
    fprime=lambda x: -np.sin(x) - 1.0,
    method="newton",
)

print(f"Brent:  root={bracketed.root:.14f}, converged={bracketed.converged}")
print(f"Newton: root={newton.root:.14f}, converged={newton.converged}")
print(f"residual={root_function(bracketed.root):.2e}")

Brent:  root=0.73908513321516, converged=True
Newton: root=0.73908513321516, converged=True
residual=7.88e-15

A sign change can also indicate a pole

Bracketing assumes continuity. The tangent function changes sign across $\pi/2$, but it has a pole there, not a zero. A solver can shrink the interval around the discontinuity and report apparent convergence, so always evaluate the residual and understand the function being solved.

false_root = optimize.root_scalar(np.tan, bracket=(1.4, 1.7), method="brentq")

print("reported location:", false_root.root)
print("distance from pi/2:", false_root.root - np.pi/2)
print("function value there:", np.tan(false_root.root))

reported location: 1.570796326794342
distance from pi/2: -5.546674231027282e-13
function value there: 1802683150020.3972

7. Interpolation

Interpolation constructs a function that passes through supplied data. It does not create new physical information between measurements, and a visually smooth curve need not be physically reasonable.

CubicSpline produces a twice-continuously differentiable piecewise cubic. PchipInterpolator preserves monotonicity within each interval and is often safer for data with sharp transitions.

x_data = np.arange(6.0)
y_data = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0])
x_dense = np.linspace(x_data[0], x_data[-1], 500)

cubic = CubicSpline(x_data, y_data)
pchip = PchipInterpolator(x_data, y_data)

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(x_data, y_data, "ko", label="data")
ax.plot(x_dense, cubic(x_dense), label="cubic spline")
ax.plot(x_dense, pchip(x_dense), label="PCHIP")
ax.axhspan(0.0, 1.0, color="0.9", zorder=-1)
ax.set(xlabel="$x$", ylabel="interpolated value")
ax.legend()
plt.show()

print(f"cubic range: [{cubic(x_dense).min():.3f}, {cubic(x_dense).max():.3f}]")
print(f"PCHIP range: [{pchip(x_dense).min():.3f}, {pchip(x_dense).max():.3f}]")

cubic range: [-0.128, 1.128]
PCHIP range: [0.000, 1.000]

The cubic spline overshoots the range of these data, whereas PCHIP remains between 0 and 1. Neither is universally superior:

• use a cubic spline when smooth derivatives are important and overshoot is acceptable;
• use PCHIP when monotonicity and shape preservation matter;
• use a fitted or smoothing model, rather than an interpolant, when measurements contain noise.

Avoid extrapolation unless the underlying model supplies a physical reason to trust it.

Summary

Problem	SciPy tool	bul also check...
Special functions	scipy.special	conventions and domains
Quadrature	integrate.quad, dblquad	error estimate and convergence
Initial-value ODE	integrate.solve_ivp	success, tolerances, conserved quantities
Fourier analysis	scipy.fft	sampling, Nyquist, leakage, normalization
Minimization	optimize.minimize_scalar	bounds, success, local versus global minima
Root finding	optimize.root_scalar	bracket continuity and residual
Interpolation	CubicSpline, PchipInterpolator	overshoot and extrapolation

The central habit is to inspect the numerical result rather than trusting a returned number merely because no exception occurred.