M50 · SciPy for Science & Engineering

The SciPy map

SciPy is a collection of subpackages, each wrapping well-studied algorithms (many in Fortran and C libraries such as LAPACK, QUADPACK, and HiGHS). You import the subpackage you need: from scipy import optimize. You already used scipy.linalg and scipy.sparse in M49.

SciPy's main subpackages. Each takes and returns NumPy arrays.

A habit for the whole module: every numerical answer is an approximation, so check it. Plug the root back into the function, compare a fitted curve with the data, integrate a case you know by hand, and read the solver's status flags.

scipy.optimize: minimize, find roots, fit, and plan

Minimization

"Find the inputs that make this function as small as possible" covers a huge range of problems: the cheapest design, the best-fitting model (minimum error), the lowest-energy molecule. optimize.minimize(f, x0) starts at a guess x0 and walks downhill. Think of standing on a foggy hillside: you can feel the slope under your feet (the gradient) and step downhill until the ground is flat. Methods like BFGS also learn the curvature as they go, so they take smarter steps. Try it with the visualizer: the step size (learning rate) decides whether you crawl, converge, or overshoot.

import numpy as np
from scipy import optimize

# A box with a square base and no lid must hold 32 litres. Minimize the material used.
def surface(side):
    height = 32 / side ** 2
    return side ** 2 + 4 * side * height

res = optimize.minimize_scalar(surface, bounds=(0.5, 10), method="bounded")
print(f"best side {res.x:.4f}, height {32 / res.x**2:.4f}, area {res.fun:.4f}")

# Rosenbrock's banana valley: the minimum is at (1, 1)
def rosen(p):
    x, y = p
    return (1 - x) ** 2 + 100 * (y - x ** 2) ** 2

def rosen_grad(p):                         # the slope in each direction (calculus; optional)
    x, y = p
    return np.array([-2 * (1 - x) - 400 * x * (y - x ** 2), 200 * (y - x ** 2)])

res = optimize.minimize(rosen, x0=[-1.5, 2.0], method="BFGS", jac=rosen_grad)
print("BFGS:", res.x.round(5), " f =", f"{res.fun:.2e}", " iterations:", res.nit, " success:", res.success)
res = optimize.minimize(rosen, x0=[-1.5, 2.0], method="BFGS")   # gradient estimated numerically
print("BFGS, no jac:", res.x.round(5), " success:", res.success, "-", res.message)
res = optimize.minimize(rosen, x0=[-1.5, 2.0], method="Nelder-Mead")
print("Nelder-Mead (no gradients):", res.x.round(4), " evaluations:", res.nfev)
res = optimize.minimize(rosen, x0=[-1.5, 2.0], bounds=[(-2, 0.5), (-1, 3)], method="L-BFGS-B")
print("with bounds x<=0.5:", res.x.round(4))

Notice the second BFGS run: without an exact gradient it estimates slopes numerically, lands almost on (1, 1), and still reports success: False with "precision loss". The answer is fine here, but the lesson is general: read res.success and res.message, and judge the result yourself.

Root finding with brentq

A root of f is an x with f(x) = 0: the break-even price, the time a ball hits the ground, the interest rate that makes a loan's payments add up. If f is continuous, f(a) < 0, and f(b) > 0, the graph must cross zero somewhere between (the intermediate value theorem). Bisection halves the interval repeatedly, keeping the half where the sign changes, just like binary search in M30. brentq combines bisection's guarantee with faster interpolation steps. You must give it a bracket [a, b] with opposite signs.

Fitting models with curve_fit

M49's least squares fits models that are linear in their parameters. Many real models aren't: exponential decay A·e^(−kt), logistic growth, a damped spring. optimize.curve_fit(model, x, y, p0) adjusts the parameters to minimize the sum of squared errors, starting from your guess p0 (a sensible guess matters for non-linear fits). It returns the best parameters and their covariance matrix; the square roots of its diagonal are the parameters' standard errors, a measure of how precisely the data pins each one down.

Linear programming with linprog

A linear program optimizes a linear objective subject to linear inequality constraints. A small workshop makes tables (profit $70, 3 hours of labour, 4 m of wood) and chairs (profit $50, 2 hours, 1.5 m). This week it has 120 hours and 140 m of wood, and can sell at most 30 chairs. How many of each maximizes profit?

Each constraint cuts the plane with a straight line; the allowed plans form a polygon. A linear objective always reaches its best value at a corner, which is why LP solvers can be so fast.

linprog minimizes c @ x subject to A_ub @ x ≤ b_ub (and optional equalities and bounds). To maximize profit, minimize its negative.

Remember: minimize for smooth objectives (watch for local minima), brentq for roots inside a sign-changing bracket, curve_fit for non-linear models (check residuals and standard errors), linprog for linear objectives with linear constraints.

scipy.integrate: areas and differential equations

quad: the area under a curve

An integral ∫ₐᵇ f(x) dx is the area under the graph of f between a and b: total distance from a speed curve, total energy from a power curve, a probability from a density. You can approximate it by slicing the area into thin rectangles and adding them up. integrate.quad does this adaptively (more slices where the curve bends) and returns the value and an estimate of its error.

solve_ivp: simulate how things change

Many laws of nature describe rates of change. A differential equation says how fast a quantity changes given its current state. A falling object with air drag: its velocity v changes at rate dv/dt = g − (c/m)·v² (gravity pulls, drag pushes back harder the faster it goes). Population growth with a limited food supply: dP/dt = r·P·(1 − P/K). Given the starting state (the initial value), a solver takes many small time steps, each time using the rate to predict the next state. solve_ivp chooses the step sizes automatically to stay within your tolerances (rtol, atol).

The state can be a vector. For the falling object, the state is [height, velocity], and the rate function returns [dh/dt, dv/dt]. An event function lets the solver stop exactly when something happens, like height reaching zero.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

g, m, c = 9.81, 80.0, 0.27                       # gravity, mass (kg), drag coefficient (kg/m)

def fall(t, state):
    h, v = state                                 # v is downward speed
    return [-v, g - (c / m) * v * abs(v)]

def hit_ground(t, state):
    return state[0]
hit_ground.terminal = True
hit_ground.direction = -1

sol = solve_ivp(fall, (0, 120), [1000.0, 0.0], events=hit_ground, max_step=0.5, rtol=1e-8)
print(f"lands after {sol.t_events[0][0]:.1f} s at {sol.y_events[0][0][1]:.1f} m/s;"
      f" terminal velocity √(mg/c) = {np.sqrt(m * g / c):.1f} m/s")

r, K = 0.4, 1000
logistic = solve_ivp(lambda t, P: r * P * (1 - P / K), (0, 30), [10], t_eval=np.linspace(0, 30, 200))
exact = K / (1 + (K / 10 - 1) * np.exp(-r * logistic.t))
print("max error vs exact logistic solution:", np.abs(logistic.y[0] - exact).max().round(4))

fig, axes = plt.subplots(1, 3, figsize=(12, 3.4))
axes[0].plot(sol.t, sol.y[0]); axes[0].set_title("height (m)")
axes[1].plot(sol.t, sol.y[1]); axes[1].axhline(np.sqrt(m * g / c), ls="--", color="gray"); axes[1].set_title("speed (m/s)")
axes[2].plot(logistic.t, logistic.y[0]); axes[2].axhline(K, ls="--", color="gray"); axes[2].set_title("logistic population")
for ax in axes:
    ax.set_xlabel("t (s)" if ax is not axes[2] else "t (years)"); ax.grid(alpha=0.3)
fig.tight_layout()

Remember: quad returns (value, error_estimate). solve_ivp(fun, (t0, t1), y0) with fun(t, y) returning rates; use t_eval for output times and events to stop on a condition.

scipy.interpolate: filling in between measurements

You measured temperature every 3 hours but need it at 10:30. Interpolation builds a function that passes exactly through your data points and estimates values between them. The simplest is straight lines between neighbours (np.interp). Splines join smooth cubic pieces so the curve and its slope have no kinks, which looks and behaves more like physical data.

np.interp(x_new, x, y): fast piecewise-linear, 1-D.
make_interp_spline(x, y, k=3): an interpolating B-spline of degree k, callable, with derivatives available. CubicSpline is similar, with boundary-condition options.
interp1d still exists but SciPy now documents it as legacy; prefer the two above in new code.

Remember: interpolants pass through every point. Linear for speed, splines for smoothness; don't extrapolate.

scipy.signal and scipy.fft: cleaning and reading signals

Frequencies and the FFT

Any signal sampled over time (sound, vibration, heartbeat, electricity use) can be described as a sum of sine waves of different frequencies (cycles per second, Hz). The Fourier transform measures how much of each frequency is present; the FFT (fast Fourier transform) computes it in O(n log n). The result, the spectrum, shows peaks at the frequencies that make up the signal. Use scipy.fft.rfft for real-valued signals and rfftfreq(n, d=1/fs) for the matching frequencies, where fs is the sampling rate.

Filtering

A low-pass filter keeps slow changes and removes fast jitter (noise); a high-pass filter does the opposite. signal.butter(order, cutoff, fs=fs, output="sos") designs a smooth Butterworth filter, and signal.sosfiltfilt applies it forwards and backwards so the output isn't shifted in time. signal.find_peaks locates local maxima, with options like height, distance (minimum spacing), and prominence (how much a peak stands out) to ignore noise bumps.

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import rfft, rfftfreq

fs = 200                                         # samples per second
t = np.arange(0, 5, 1 / fs)
rng = np.random.default_rng(0)
clean = np.sin(2 * np.pi * 1.2 * t) + 0.4 * np.sin(2 * np.pi * 3 * t)
noisy = clean + 0.5 * np.sin(2 * np.pi * 50 * t) + rng.normal(0, 0.3, t.size)   # 50 Hz hum + noise

freqs = rfftfreq(t.size, d=1 / fs)
spectrum = np.abs(rfft(noisy)) / t.size * 2
top = freqs[np.argsort(spectrum)[-3:]]
print("strongest frequencies (Hz):", np.sort(top).round(2))

sos = signal.butter(4, 8, btype="low", fs=fs, output="sos")     # keep < 8 Hz
smooth = signal.sosfiltfilt(sos, noisy)
print("RMS error vs clean: noisy", np.sqrt(np.mean((noisy - clean) ** 2)).round(3),
      " filtered", np.sqrt(np.mean((smooth - clean) ** 2)).round(3))

peaks, props = signal.find_peaks(smooth, height=0.5, distance=fs // 2)
print(f"{peaks.size} peaks at t =", t[peaks].round(2))

fig, axes = plt.subplots(1, 2, figsize=(12, 3.5))
axes[0].plot(t, noisy, lw=0.6, alpha=0.6, label="noisy")
axes[0].plot(t, smooth, lw=2, label="low-pass filtered")
axes[0].plot(t[peaks], smooth[peaks], "x", ms=10, mew=2, label="peaks")
axes[0].set_xlim(0, 3); axes[0].legend(); axes[0].set_xlabel("t (s)")
axes[1].plot(freqs, spectrum); axes[1].set_xlim(0, 60); axes[1].set_xlabel("frequency (Hz)"); axes[1].set_title("spectrum")
for ax in axes:
    ax.grid(alpha=0.3)
fig.tight_layout()

Remember: the FFT turns time into frequency; filter with butter + sosfiltfilt; tune find_peaks with distance and prominence.

scipy.spatial: nearest neighbours, distances, hulls

"Which store is closest to each customer?" Checking every store for every customer costs O(n·m). A KD-tree splits space into boxes, alternating between x and y (and further dimensions), like a binary search tree for points (M24). A query only visits the few boxes near the target, about O(log n) per query in low dimensions. Build it once with KDTree(points), then call tree.query(targets, k) for the k nearest, or tree.query_ball_point(x, r) for everything within radius r.

For all pairwise distances between two smallish sets, use scipy.spatial.distance.cdist(A, B, metric=...) (many metrics: "euclidean", "cityblock" (L1), "cosine", …). The convex hull is the smallest convex polygon containing all the points: imagine a rubber band snapped around nails. ConvexHull(points) gives its vertices, area, and perimeter.

import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import KDTree, ConvexHull
from scipy.spatial.distance import cdist

rng = np.random.default_rng(1)
stores = rng.uniform(0, 100, (500, 2))
customers = rng.uniform(0, 100, (5000, 2))

tree = KDTree(stores)
t0 = time.perf_counter(); dist, idx = tree.query(customers, k=1); t_tree = time.perf_counter() - t0
t0 = time.perf_counter(); brute = cdist(customers, stores).argmin(axis=1); t_brute = time.perf_counter() - t0
print(f"KDTree {t_tree * 1000:.1f} ms, brute force {t_brute * 1000:.1f} ms, same answer: {np.array_equal(idx, brute)}")
print("mean distance to nearest store:", dist.mean().round(2))
print("stores within 5 of the centre:", tree.query_ball_point([50, 50], r=5))

hull = ConvexHull(stores)
print(f"hull has {len(hull.vertices)} corner stores; area {hull.volume:.0f} (in 2-D, .volume is area and .area is perimeter)")

fig, ax = plt.subplots(figsize=(5.5, 5.5))
ax.scatter(*customers[:600].T, c=idx[:600], cmap="tab20", s=6, alpha=0.6)
ax.scatter(*stores.T, c="k", marker="^", s=12, label="stores")
loop = np.append(hull.vertices, hull.vertices[0])
ax.plot(*stores[loop].T, "r-", lw=1.5, label="convex hull")
ax.set_aspect("equal"); ax.legend(loc="upper right"); ax.set_title("customers coloured by nearest store")
fig.tight_layout()

Remember: build a KDTree once, query many times; cdist for full distance matrices; ConvexHull for outer boundaries (2-D volume = area).

scipy.stats: a preview

scipy.stats has about a hundred probability distributions, each with the same methods: pdf/pmf (density or probability), cdf (probability of being ≤ x), ppf (the inverse: which x has this cumulative probability), rvs (random samples), and fit. It also has descriptive statistics and hypothesis tests. M54 covers them properly; here is a taste.

Numerical pitfalls: floating point, conditioning, tolerance

Every result in this module is computed with 64-bit floating-point numbers, which have about 16 significant decimal digits (the IEEE-754 format; see the bits visualizer below). Three ideas explain almost every numerical surprise:

Rounding error: most decimals, like 0.1, can't be stored exactly. Each operation rounds, and errors accumulate. Worst of all is catastrophic cancellation: subtracting two nearly equal numbers wipes out the leading digits and leaves mostly noise.
Conditioning: some problems amplify small input errors no matter how good the algorithm (M49's condition number). An ill-conditioned problem needs reformulating, not a better solver.
Tolerance: iterative algorithms (roots, optimizers, ODE solvers) stop when a change is below a tolerance. Tolerances too loose give wrong answers; too tight waste time or fail to converge. Compare floats with np.isclose(a, b, rtol=..., atol=...), never ==, and remember that atol matters near zero, where relative tolerance is meaningless.

import math
import numpy as np
from scipy.integrate import solve_ivp

print("0.1 + 0.2 =", repr(0.1 + 0.2))
total = 0.0
for _ in range(10):
    total += 0.1                       # ten separate roundings
print("adding 0.1 ten times:", repr(total), " == 1.0?", total == 1.0)
print("sum() and math.fsum compensate for rounding:", sum([0.1] * 10) == 1.0, math.fsum([0.1] * 10) == 1.0)

# Catastrophic cancellation: (1 - cos x) / x^2 should approach 0.5 as x -> 0
for x in (1e-2, 1e-5, 1e-8):
    naive = (1 - math.cos(x)) / x ** 2
    stable = 2 * math.sin(x / 2) ** 2 / x ** 2         # same formula, rewritten
    print(f"x={x:.0e}  naive {naive:.10f}  stable {stable:.10f}")

print("expm1/log1p exist for this reason:", math.exp(1e-12) - 1, "vs", math.expm1(1e-12))

# Tolerance changes the answer of an ODE solver
exact = math.exp(-5)
for rtol in (1e-2, 1e-5, 1e-9):
    sol = solve_ivp(lambda t, y: -y, (0, 5), [1.0], rtol=rtol, atol=rtol * 1e-3)
    print(f"rtol={rtol:.0e}  y(5) error {abs(sol.y[0, -1] - exact):.1e}  steps {sol.t.size}")

print("isclose near zero:", np.isclose(1e-12, 0.0), " with atol=0:", np.isclose(1e-12, 0.0, atol=0.0))

Remember: expect rounding, avoid subtracting near-equal numbers, respect conditioning, and treat tolerances as part of the answer.

Practice

Use SciPy. Tests compare with tolerances.

A product's profit at price p is profit(p) = (p - unit_cost) * demand(p) - fixed_cost with demand(p) = 1000 * exp(-p / 40). Write break_even_prices(unit_cost, fixed_cost) returning a sorted tuple of the two prices where profit is exactly 0. First locate the price of maximum profit with optimize.minimize_scalar on -profit over bounds (unit_cost, 500); then use optimize.brentq on each side ([unit_cost, p_best] and [p_best, 500]). Raise ValueError if even the best price makes a loss.

Hint

res = optimize.minimize_scalar(lambda p: -profit(p), bounds=(unit_cost, 500), method="bounded"). If profit(res.x) <= 0, raise. Otherwise the signs flip on both sides of res.x.

Show solution

import numpy as np
from scipy import optimize

def break_even_prices(unit_cost, fixed_cost):
    def profit(p):
        return (p - unit_cost) * 1000 * np.exp(-p / 40) - fixed_cost

    res = optimize.minimize_scalar(lambda p: -profit(p), bounds=(unit_cost, 500), method="bounded")
    best = res.x
    if profit(best) <= 0:
        raise ValueError("no price is profitable")
    lo = optimize.brentq(profit, unit_cost, best)
    hi = optimize.brentq(profit, best, 500)
    return (lo, hi)

Write logistic(t, K, r, t0) returning K / (1 + exp(-r * (t - t0))), and fit_logistic(t, y) returning (K, r, t0) fitted with optimize.curve_fit. Use a data-driven initial guess: K ≈ max(y), r ≈ 1, t0 ≈ the time where y is closest to half its maximum.

Hint

t0_guess = t[np.argmin(np.abs(y - y.max() / 2))], then popt, _ = optimize.curve_fit(logistic, t, y, p0=[y.max(), 1.0, t0_guess]).

Show solution

import numpy as np
from scipy import optimize

def logistic(t, K, r, t0):
    return K / (1 + np.exp(-r * (t - t0)))

def fit_logistic(t, y):
    t, y = np.asarray(t, dtype=float), np.asarray(y, dtype=float)
    t0_guess = t[np.argmin(np.abs(y - y.max() / 2))]
    popt, _ = optimize.curve_fit(logistic, t, y, p0=[y.max(), 1.0, t0_guess])
    K, r, t0 = popt
    return float(K), float(r), float(t0)

Write cheapest_diet(cost, nutrients, minimums). cost[j] is the price per serving of food j; nutrients[i, j] is the amount of nutrient i per serving of food j; minimums[i] is the daily minimum of nutrient i. Return (servings, total_cost) minimizing cost with all servings ≥ 0. Raise ValueError if the problem is infeasible. (Remember linprog wants ≤ constraints: multiply ≥ rows by −1.)

Hint

linprog(c=cost, A_ub=-nutrients, b_ub=-minimums, bounds=[(0, None)] * len(cost)). Check res.success.

Show solution

import numpy as np
from scipy.optimize import linprog

def cheapest_diet(cost, nutrients, minimums):
    cost = np.asarray(cost, dtype=float)
    res = linprog(c=cost, A_ub=-np.asarray(nutrients, dtype=float), b_ub=-np.asarray(minimums, dtype=float),
                  bounds=[(0, None)] * cost.size, method="highs")
    if not res.success:
        raise ValueError(f"no feasible diet: {res.message}")
    return res.x, float(res.fun)

Write fall_time(height, mass, drag, g=9.81) that simulates a dropped object with dv/dt = g - (drag / mass) * v * |v| and dh/dt = -v, starting at rest, and returns the time when h reaches 0, using a terminal event. With drag=0 the answer must match sqrt(2 h / g).

Hint

Define def ground(t, s): return s[0], set ground.terminal = True and ground.direction = -1, integrate over a long span such as (0, 10_000) with rtol=1e-9, and return sol.t_events[0][0].

Show solution

import numpy as np
from scipy.integrate import solve_ivp

def fall_time(height, mass, drag, g=9.81):
    def rates(t, s):
        h, v = s
        return [-v, g - (drag / mass) * v * abs(v)]

    def ground(t, s):
        return s[0]
    ground.terminal = True
    ground.direction = -1

    sol = solve_ivp(rates, (0, 10_000), [height, 0.0], events=ground, rtol=1e-9, atol=1e-9)
    return float(sol.t_events[0][0])

Write assign_customers(stores, customers, max_distance). Return an integer array with the index of the nearest store for each customer, or -1 if the nearest store is farther than max_distance. Use scipy.spatial.KDTree. Then write store_loads(assignment, n_stores) returning how many customers each store serves (ignore -1).

Hint

dist, idx = KDTree(stores).query(customers, k=1), then np.where(dist <= max_distance, idx, -1). For loads, np.bincount(a[a >= 0], minlength=n_stores).

Show solution

import numpy as np
from scipy.spatial import KDTree

def assign_customers(stores, customers, max_distance):
    dist, idx = KDTree(stores).query(customers, k=1)
    return np.where(dist <= max_distance, idx, -1).astype(int)

def store_loads(assignment, n_stores):
    assignment = np.asarray(assignment)
    return np.bincount(assignment[assignment >= 0], minlength=n_stores)

Mini-project: model a cooling cup of coffee

Newton's law of cooling says an object's temperature T approaches the room temperature T_room at a rate proportional to the difference: dT/dt = −k (T − T_room). Its exact solution is T(t) = T_room + (T₀ − T_room)·e^(−kt). Here is the whole workflow in miniature: fit k from noisy thermometer readings with curve_fit, then simulate with solve_ivp and answer a practical question ("when is it drinkable at 60 °C?") with an event.

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy.integrate import solve_ivp

rng = np.random.default_rng(42)
minutes = np.arange(0, 61, 4.0)
true_k, room, start = 0.045, 21.0, 88.0
readings = room + (start - room) * np.exp(-true_k * minutes) + rng.normal(0, 0.8, minutes.size)

def newton(t, T0, k, T_room):
    return T_room + (T0 - T_room) * np.exp(-k * t)

popt, pcov = optimize.curve_fit(newton, minutes, readings, p0=[readings[0], 0.05, 20])
T0, k, T_room = popt
print(f"fitted T0={T0:.1f}°C  k={k:.4f}/min ± {np.sqrt(pcov[1, 1]):.4f}  room={T_room:.1f}°C")

def drinkable(t, T):
    return T[0] - 60
drinkable.terminal = True

sol = solve_ivp(lambda t, T: -k * (T - T_room), (0, 120), [T0], events=drinkable, dense_output=True)
t60 = sol.t_events[0][0]
print(f"reaches 60°C after {t60:.1f} min;  closed form: {np.log((T0 - T_room) / (60 - T_room)) / k:.1f} min")

tt = np.linspace(0, 90, 300)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(minutes, readings, "o", label="thermometer")
ax.plot(tt, newton(tt, *popt), label="fitted model")
ax.plot(tt[tt <= t60], sol.sol(tt[tt <= t60])[0], "--", lw=3, alpha=0.6, label="solve_ivp (until 60°C)")
ax.axhline(60, color="gray", ls=":"); ax.axvline(t60, color="gray", ls=":")
ax.set_xlabel("minutes"); ax.set_ylabel("°C"); ax.legend(); ax.grid(alpha=0.3)
fig.tight_layout()

Measure (or download) a real cooling curve: a mug of tea with a kitchen thermometer every 2 minutes, a laptop's CPU temperature after a benchmark, or a public dataset. Build a small package that fits Newton's law, simulates what-if scenarios, and reports results.

pixi init coolfit && cd coolfit
pixi add python=3.14 numpy scipy matplotlib pytest
pixi run python -m coolfit readings.csv --target 60      # prints a report, writes cooling.png
pixi run pytest -q

Spec

load_readings(path) reads a CSV of minutes,celsius into two arrays (use csv or np.loadtxt).
fit_newton(t, T) returns a frozen dataclass with T0, k, T_room and their standard errors from curve_fit.
simulate(params, minutes, k_scale=1.0) uses solve_ivp; k_scale models a what-if (an insulated mug might halve k).
time_to(params, target) uses a terminal event, and a test confirms it matches the closed-form ln((T0 − T_room)/(target − T_room)) / k.
The CLI prints the fitted parameters with uncertainties, the time to the target, and the residual RMS, and saves a figure with data, fit, residuals, and the insulated-mug scenario.

fit_newton recovers known parameters from synthetic data within 5% (tested)
simulate and the closed-form solution agree within 0.01 °C (tested)
time_to raises ValueError when the target is below room temperature (tested)
Residuals are plotted, and you wrote one sentence on whether they look like random noise
The report includes a ± uncertainty for k
pixi run pytest passes

Stretch goals

Add a second process: the mug itself absorbs heat (a two-variable ODE for coffee and mug temperatures) and compare fits.
Pick the best of several models with a residual-based criterion, then read about AIC in M54.
Smooth noisy readings first with a Savitzky–Golay filter (scipy.signal.savgol_filter) and compare the fitted k.

Quiz

Recap

SciPy is a toolbox of subpackages on top of NumPy: optimize, integrate, interpolate, linalg, sparse, signal, fft, spatial, stats, and more.
minimize walks downhill (local); brentq needs a sign-change bracket; curve_fit fits non-linear models with uncertainties; linprog solves linear programs at corners.
quad integrates with an error estimate; solve_ivp simulates ODEs, with events to stop on conditions.
Splines interpolate smoothly between points; don't extrapolate or interpolate noise.
FFT reveals frequencies; Butterworth + sosfiltfilt removes noise; find_peaks finds events.
KDTree for nearest neighbours, cdist for distance matrices, ConvexHull for boundaries.
Floating point rounds, problems can be ill-conditioned, and tolerances are part of every answer: verify.

Going further

SciPy user guide (one chapter per subpackage)
Optimization tutorial and solve_ivp reference
Interpolation tutorial (including the legacy status of interp1d)
Python tutorial: floating-point arithmetic, issues and limitations
Next: M51 · pandas: Series & DataFrames; for statistics, M54