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.

integrate.solve_ivp(fun, t_span, y0,...)

fun(t, y), where t is a scalar and y is an ndarray with len(y) = len(y0)

Double pendulum: complete derivation

from IPython.display import Image
Image(url='http://upload.wikimedia.org/wikipedia/commons/c/c9/Double-compound-pendulum-dimensioned.svg')

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)==state.'''
    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, 100.0, 4000)

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: 76511
solution shape: (4, 4000)

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_joint = ell*np.sin(theta1)
y_joint = -ell*np.cos(theta1)
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)  # only every 4-th step is plotted

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")  # need to fix the plot size, so that it is stable
ax.grid()
line, = ax.plot([], [], "o-", linewidth=2)                                  # initial plot of empty rods, that we will be updating
time_text = ax.text(0.05, 0.92, "", transform=ax.transAxes)                 # initial text that we will be updating

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: 2.7209344761785854e-08

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.

Estimating the largest Lyapunov exponent

Double pendulum is known to be chaotic for large enough initial displacements. See for example the plot of Lyapunov exponents at wikipedia.

The Lyapunov exponent measures the exponential growth of a small displacement between two nearby trajectories. If their initial separation is $\delta_0$, we expect, while the trajectories remain close,

$$\delta(t) \simeq \delta_0 e^{\lambda t}.$$

Why not evolve the two trajectories only once?

At first, we could integrate one reference trajectory and one slightly perturbed trajectory over the entire time interval and estimate $\lambda$ from the slope of $\log[\delta(t)/\delta_0]$. This works only during the initial exponential-growth regime. Eventually, the two trajectories become separated by a distance comparable to the size of the accessible phase space. Their separation then saturates and can no longer grow exponentially. If we include this saturated region, we underestimate the Lyapunov exponent. Conversely, if the separation becomes extremely small, roundoff error can dominate it.

Periodic rescaling

We therefore divide the evolution into short intervals of duration $\Delta t$. At the beginning of every interval, the perturbed state is a small distance $\delta_0$ from the reference state. For each interval we:

1. Evolve both states through the same interval $\Delta t$.
2. Form their final displacement $\Delta\mathbf y_k=\mathbf y_k^{\mathrm{pert}}-\mathbf y_k^{\mathrm{ref}}$ and measure $\delta_k=\lVert\Delta\mathbf y_k\rVert$. Angular differences are wrapped into $[-\pi,\pi)$ before the norm is calculated.
3. Record the logarithmic growth $\log(\delta_k/\delta_0)$.
4. Keep the direction of $\Delta\mathbf y_k$, but reset its length to $\delta_0$:

$$\Delta\mathbf y_k^{\mathrm{new}}=\delta_0\frac{\Delta\mathbf y_k}{\delta_k}, \qquad \mathbf y_{k+1}^{\mathrm{pert}}=\mathbf y_k^{\mathrm{ref}}+\Delta\mathbf y_k^{\mathrm{new}}.$$

Only the length is reset; the direction is retained. Repetition naturally aligns the displacement with the most rapidly expanding direction in phase space, which is why the method approaches the largest Lyapunov exponent. The growth removed during every rescaling is not lost: it was recorded in step 3. After $N$ intervals, the estimate is

$$\lambda_{\max} \simeq \frac{1}{N\Delta t}\sum_{k=1}^{N}\log\left(\frac{\delta_k}{\delta_0}\right).$$

The function below returns this running finite-time estimate. A value that approaches zero indicates regular motion, whereas convergence toward a positive value indicates sensitive dependence on initial conditions and chaotic motion.

def largest_lyapunov(rhs, initial_state, rhs_args=(), delta0=1e-8, segment_time=0.2, nsegments=500):
    """Estimate the largest Lyapunov exponent by periodic rescaling."""
    reference = np.asarray(initial_state, dtype=float).copy() # initial state copy, so that we don't change initial state

    # Initially perturb theta1.
    displacement = np.zeros_like(reference)
    displacement[0] = delta0                # start with first angle be little different
    perturbed = reference + displacement    # this is the second initial condition

    log_growth = np.empty(nsegments)

    for i in range(nsegments):
        reference_solution = integrate.solve_ivp(rhs, (0.0, segment_time), reference, args=rhs_args, method="DOP853",rtol=1e-10, atol=1e-12)
        perturbed_solution = integrate.solve_ivp(rhs, (0.0, segment_time), perturbed, args=rhs_args, method="DOP853",rtol=1e-10, atol=1e-12)

        reference = reference_solution.y[:, -1] # solution after segment_time 
        perturbed = perturbed_solution.y[:, -1] # the second solution after segment time

        displacement = perturbed - reference # difference between the two solutions

        # Compare angular coordinates modulo 2*pi.
        displacement[:2] = (displacement[:2] + np.pi) % (2*np.pi) - np.pi
        # euclidian norm of both displacements
        separation = np.linalg.norm(displacement)
        # log_growth[i] = log(delta_k/delta0)
        log_growth[i] = np.log(separation/delta0)

        # Restore the separation to delta0 without changing its direction.
        displacement *= delta0/separation
        perturbed = reference + displacement

    times = segment_time*np.arange(1, nsegments + 1)
    running_exponent = np.cumsum(log_growth)/times

    return times, running_exponent

initial_state = [0.4*np.pi, 0.1*np.pi, 0, 0]
times, lyapunov = largest_lyapunov(double_pendulum_rhs, initial_state, rhs_args=(m, ell, g))

plt.plot(times, lyapunov)
plt.xlabel("time")
plt.ylabel(r"finite-time $\lambda_{\max}$")
plt.grid()
plt.show()

print("Estimated largest Lyapunov exponent:", lyapunov[-1])

Estimated largest Lyapunov exponent: 0.3414560992043457

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.

Basic idea of FFT

Define

$$W_N=e^{-2\pi i/N}.$$

The discrete Fourier transform is $$X_k=\sum_{j=0}^{N-1}x_jW_N^{jk}.$$

For even (N), separate the input into its even- and odd-indexed elements:

$$ X_k = \sum_{m=0}^{N/2-1}x_{2m}W_N^{2mk} + \sum_{m=0}^{N/2-1}x_{2m+1}W_N^{(2m+1)k}. $$

Using

$$W_N^{2mk}=W_{N/2}^{mk},$$

we obtain

$$ X_k = \sum_{m=0}^{N/2-1}x_{2m}W_{N/2}^{mk} + W_N^k \sum_{m=0}^{N/2-1}x_{2m+1}W_{N/2}^{m k}. $$ or compactly $$X_k=E_k+W_N^kO_k,$$

where

$$E_k=\sum_{m=0}^{N/2-1}x_{2m}W_{N/2}^{mk}$$

and

$$O_k=\sum_{m=0}^{N/2-1}x_{2m+1}W_{N/2}^{mk}$$

are two Fourier transforms of length $N/2$. Notice that this $X_k$ is valid for $k\in [0,N)$.

But we can simplify the calculation so that the same operations processing the first part $k\in [0,N/2)$ also gives the second half of the spectrum $k\in [N/2,N)$. Since

$$E_{k+N/2}=E_k,\qquad O_{k+N/2}=O_k$$

because

$$W^{m(k+N/2)}_{N/2}=W^{m k}_{N/2} (W^{N/2}_{N/2})^m=W^{m k}_{N/2}$$

and $$W_N^{k+N/2}=-W_N^k,$$

we obtain

\begin{eqnarray} X_k=E_k+W_N^kO_k \\ X_{k+N/2}=E_k-W_N^kO_k, \end{eqnarray} valid for $k\in [0,N/2)$.

This pair of operations is called an FFT butterfly.

For $N=2^p$, the even–odd decomposition can be repeated $\log_2N$ times. Each recursion level requires $O(N)$ operations, so the complete transform requires $O(N\log N)$ operations rather than $O(N^2)$.

Simple example

Consider simple analytic function: $$f(t) = a \cos(2\pi \alpha t) + b \sin(2\pi \beta t).$$

Since the function is real, we can use rfft and irfft for Fourier and inverse Fourier transform. For complex functions, we need to use fft and ifft.

The Fourier transform of a real functions has a property $X_{-k}^* = X_k$ hence we don't need negative frequencies and consequently rfftfreq omits negative frequencies, while fftfreq includes them.

We expect the result: $$X(\omega) = \frac{a}{2} \delta(\omega-2\pi\alpha) - i\frac{\beta}{2} \delta(\omega-2\pi\beta) + \frac{a}{2} \delta(\omega+2\pi\alpha) + i\frac{\beta}{2} \delta(\omega+2\pi\beta)$$

The last two terms don't show up because they are assumed property of rfft because $X(-\omega)=X^*(\omega)$

from scipy.fft import rfft, irfft, rfftfreq


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)

ramplitude=spectrum.real *2./n_samples
iamplitude=spectrum.imag *2./n_samples

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, ramplitude, '-bo')
axes[1].stem(frequency, iamplitude, '-ro')
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

If the signal contains frequency beyond Nyquist frequency, it is errorneously folded into spectrum below Nyquist frequency. Frequencies above Nyquist appear at incorrect frequencies below Nyquist because of aliasing.

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*110*t) + 0.7*np.sin(2*np.pi*120*t)
spectrum = rfft(signal)
frequency = rfftfreq(n_samples, d=1/sampling_rate)

reconstructed = irfft(spectrum, n=n_samples)

ramplitude=spectrum.real *2./n_samples
iamplitude=spectrum.imag *2./n_samples



fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(t, signal-reconstructed)
#axes[0].plot(t, reconstructed)
axes[0].set(xlim=(0, 0.5), xlabel="time (s)", ylabel="signal")
axes[1].stem(frequency, ramplitude, '-bo')
axes[1].stem(frequency, iamplitude, '-ro')
axes[1].set(xlim=(0, 100), xlabel="frequency (Hz)", ylabel="amplitude")
plt.show()

Why does the reconstructed spectrum look identical to the signal up to machine precision?

Explanation in this cell:

t_fine = np.linspace(0.0, 0.05, 5000)

original_continuous = 1.2*np.cos(2*np.pi*110*t_fine) + 0.7*np.sin(2*np.pi*120*t_fine)
aliased_continuous  = 1.2*np.cos(2*np.pi*90*t_fine) - 0.7*np.sin(2*np.pi*80*t_fine)

plt.figure(figsize=(10, 4))

plt.plot(t_fine,original_continuous, label="original: 110 Hz and 120 Hz")

plt.plot(t_fine, aliased_continuous, "--", label="aliased: 90 Hz and −80 Hz")

mask = t <= 0.05
plt.plot(t[mask], signal[mask], "ko", label="samples at 200 Hz")

plt.xlabel("time [s]")
plt.ylabel("signal")
plt.legend()
plt.show()

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 periodic 
signal_leaking = np.sin(2*np.pi*12.3*t)     # not periodic on the interval [0,duration]. FFT forces periodic function.
amp_periodic = 2*np.abs(rfft(signal_periodic))/n_samples
amp_leaking = 2*np.abs(rfft(signal_leaking))/n_samples

plt.stem(frequency, amp_periodic,'C0')
plt.stem(frequency, amp_leaking, 'C1')
plt.show()

window = hann(n_samples, sym=False)

plt.plot(t, signal_leaking, label='original signal')
plt.plot(t, window, label='Hann window')
plt.plot(t, window*signal_leaking, label='signal * Hann window') # makes it periodic.
plt.legend(loc='best')
plt.show()

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)

plt.stem(frequency, amp_leaking, 'C0', label='non-periodic function')
plt.stem(frequency, amp_windowed, 'C1', label='non-periodic*window')

plt.xlim([0, 35])
plt.xlabel("frequency (Hz)")
plt.ylabel("amplitude")
plt.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")
result = optimize.minimize(potential, x0=[0.8], bounds=[(-2.0, 2.0)], method="L-BFGS-B")

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 closest to the starting point")
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: [0.93351692]
V(x_min): -0.1338469407825757

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.

Both CubicSpline and PchipInterpolator join neighboring data points with cubic polynomials, and both pass exactly through the supplied values. They differ in how the derivatives at the joining points are chosen.

Cubic spline. CubicSpline requires the interpolating function, its first derivative, and its second derivative to be continuous at every data point. In other words, it is a $C^2$ interpolant. Two additional boundary conditions complete the construction; Continuous curvature produces a very smooth curve, which is desirable when the underlying function is known to be smooth. However, the smoothness conditions can force the polynomial to overshoot the data or create artificial maxima and minima.

PCHIP. PchipInterpolator requires the function and its first derivative to be continuous, so it is a $C^1$ interpolant, but its second derivative is generally discontinuous at the data points. Instead of imposing continuous curvature, PCHIP chooses the first derivative from the slopes of the neighboring intervals. If the two neighboring slopes have opposite signs, indicating a sampled maximum or minimum, the derivative is set to zero. If they have the same sign, a weighted harmonic mean is used. These conditions prevent overshoot and preserve the monotonicity of each monotonic portion of the data.

PCHIP does not require the complete data set to be monotonic. It can interpolate functions such as $\sin x$ or $\cos x$, treating their increasing and decreasing portions separately. For a very smooth function such as a sine wave, however, a cubic spline will usually give a smoother and more accurate approximation.

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.

x_data = np.linspace(0.0, 2*np.pi, 8)
y_data = np.cos(x_data)

x = np.linspace(0.0, 2*np.pi, 1000)
y_exact = np.cos(x)

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

y_cubic = cubic(x)
y_pchip = pchip(x)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)

# Interpolated functions
ax1.plot(x, y_exact, "k", linewidth=2, label="exact")
ax1.plot(x, y_cubic, "--", label="cubic spline")
ax1.plot(x, y_pchip, ":", linewidth=2, label="PCHIP")
ax1.plot(x_data, y_data, "ko", label="sampled points")
ax1.set(ylabel=r"$\cos x$")
ax1.legend()
ax1.grid()

# Absolute interpolation errors
ax2.semilogy(x,np.maximum(np.abs(y_cubic - y_exact), 1e-16),label="cubic-spline error")
ax2.semilogy(x,np.maximum(np.abs(y_pchip - y_exact), 1e-16),label="PCHIP error")
ax2.set(xlabel=r"$x$",ylabel="absolute error")
ax2.legend()
ax2.grid()

plt.show()

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.