Use NumPy arrays, plots, fitting, numerical methods, units and instrument data files to do real engineering and science calculations in Python.

Scientific Python

## Why this matters

Measurements become useful when they pin down the numbers in a model: the spring constant of a sample, the rate constant of a reaction, the time constant of a sensor, the Michaelis constant of an enzyme, the gravitational acceleration from a falling object. **Curve fitting** finds the parameter values that make the model match the data best, and, just as important, tells you **how uncertain** those values are. A fitted value without an uncertainty is only half an answer.

## Core idea: least squares

For a model `f(x; p)` with parameters `p` and data points `(x_i, y_i)`, least squares chooses `p` to minimise the sum of squared residuals:

```text
S(p) = Σ (y_i - f(x_i; p))^2
```

If the model is **linear in its parameters** (a straight line, a polynomial, `a sin t + b cos t`), there is a direct solution with a design matrix, as you saw with `np.linalg.lstsq`. If not (`A e^(-t/τ)`, `Vmax S / (Km + S)`), an iterative optimiser starts from a guess and walks downhill. SciPy's `curve_fit` handles both.

## scipy.optimize.curve_fit

```python
import numpy as np
from scipy.optimize import curve_fit

def decay(t, A, tau, C):
    return A * np.exp(-t / tau) + C

t = np.linspace(0, 10, 30)                      # seconds
rng = np.random.default_rng(0)
y = decay(t, 5.0, 2.5, 0.8) + rng.normal(0, 0.05, t.size)

p0 = [4.0, 2.0, 0.5]                            # starting guess
popt, pcov = curve_fit(decay, t, y, p0=p0)
perr = np.sqrt(np.diag(pcov))                   # 1-sigma uncertainties

for name, val, err in zip(["A", "tau", "C"], popt, perr):
    print(f"{name} = {val:.3f} ± {err:.3f}")
```

The pieces:

- The **model function** takes `x` first, then each parameter as a separate argument.
- `p0` is the starting guess. For non-linear models a sensible guess matters: estimate `A` from the first point, `C` from the last, `tau` from where the curve has fallen about two thirds of the way.
- `popt` holds the best-fit parameters.
- `pcov` is the **covariance matrix**. Its diagonal holds the variances, so `sqrt(diag(pcov))` gives standard errors. Off-diagonal terms show how parameters trade off against each other.

If you know each point's measurement uncertainty, pass it: `curve_fit(f, x, y, sigma=y_err, absolute_sigma=True)`. Without `absolute_sigma=True`, SciPy rescales the uncertainties to match the scatter of the residuals.

## Where the uncertainties come from

For a linear model with design matrix `X` (one column per parameter), the covariance of the fitted parameters is

```text
cov(p) = s^2 (X^T X)^(-1),      s^2 = S_min / (n - p)
```

where `s^2` is the residual variance and `n - p` the degrees of freedom. This is exactly what `curve_fit` reports for linear models and approximates for non-linear ones. You can compute it yourself with NumPy, which is what the exercise does.

### Worked example: g from a dropped ball

A video tracker gives the height of a tossed ball every 0.05 s. The model is `h(t) = h0 + v0 t - g t^2 / 2`: linear in `h0`, `v0` and `g`, so the design matrix has columns `[1, t, -t^2/2]`.

```python
X = np.column_stack([np.ones_like(t), t, -0.5 * t**2])
p, *_ = np.linalg.lstsq(X, h, rcond=None)
res = h - X @ p
s2 = res @ res / (len(t) - 3)
cov = s2 * np.linalg.inv(X.T @ X)
h0, v0, g = p
g_err = np.sqrt(cov[2, 2])
```

With 20 points and 1 cm of tracking noise, a typical result is `g = 9.78 ± 0.08 m/s²`: consistent with 9.81, and the ± tells you that a second toss could easily give 9.86.

## Judging the fit

The optimiser always returns numbers. Check them:

1. **Plot the residuals** against `x`. They should look like random scatter around zero. A curve or trend means the model is wrong or incomplete.
2. **Compare the residual size** with your measurement uncertainty. With known errors, the reduced chi-square `χ²/ν = Σ((y - f)/σ)² / (n - p)` should be near 1. Much larger: the model misses something or errors are underestimated. Much smaller: errors are overestimated.
3. **Look at the uncertainties.** A parameter with error larger than its value is not determined by your data.
4. **Check correlations.** Strongly correlated parameters (correlation near ±1) mean the data cannot separate them; consider fixing one or collecting data over a wider range.

## Common mistakes

- **No starting guess** for a non-linear model: `curve_fit` defaults every parameter to 1, which may converge to nonsense or fail with "Optimal parameters not found".
- **Reporting `pcov` diagonals as uncertainties**: they are variances; take the square root.
- **Fitting too many parameters** for the data, which gives huge, correlated uncertainties.
- **Extrapolating** the fitted model far outside the measured range.
- **Ignoring units**: the units of a parameter follow from the model (if `t` is in seconds, `tau` is in seconds).

## Quick recap

- `popt, pcov = curve_fit(model, x, y, p0=guess)`; uncertainties are `np.sqrt(np.diag(pcov))`.
- For linear-in-parameter models, `cov = s² (XᵀX)⁻¹` with `s² = SSR / (n - p)`.
- Always look at residuals, the uncertainties and parameter correlations before trusting a fit.

## Try it

In the **Do** tab you will fit the dropped-ball model with NumPy, compute the uncertainty in `g` from the covariance matrix, and test it on seeded noisy data.


Fit a model to data with scipy.optimize.curve_fit, read parameter uncertainties from the covariance matrix, and judge the fit with residuals.

Curve fitting

Curve Fitting in Python (2022)

Curve fitting in Python with curve_fit

## Why this matters

Most physical models are written as rates of change. A capacitor's voltage changes at a rate set by the current; a drug's concentration falls at a rate proportional to how much is left; a mass on a spring accelerates according to the spring force; a reactor's temperature rises with the heat released. These are **ordinary differential equations (ODEs)**. Only the simplest have neat formula solutions. For everything else, including the realistic versions with friction, non-linear springs and coupled reactions, you solve them **numerically** by stepping forward in time.

## Core idea: state, rate, step

Write the system in first-order form:

```text
dy/dt = f(t, y),    y(t0) = y0
```

`y` is the **state** (a number or a NumPy array), and `f` returns its rate of change. A second-order equation becomes two first-order ones by adding the velocity to the state. For a mass on a spring, `m x'' = -k x - c x'`:

```python
import numpy as np

def spring(t, y, k=4.0, c=0.0, m=1.0):
    x, v = y
    return np.array([v, (-k * x - c * v) / m])     # [dx/dt, dv/dt]
```

## Euler's method

The simplest step: follow the current slope for a time `h`.

```text
y_{n+1} = y_n + h f(t_n, y_n)
```

```python
def euler(f, y0, t):
    y = np.zeros((len(t), np.size(y0)))
    y[0] = y0
    for n in range(len(t) - 1):
        h = t[n + 1] - t[n]
        y[n + 1] = y[n] + h * f(t[n], y[n])
    return y
```

The loop over time steps is unavoidable (each step needs the previous one), but each step is vectorised over the state. Euler is easy to understand and badly inaccurate: its error is proportional to `h`, so halving the step only halves the error. On the undamped spring it slowly **gains energy**, and the oscillation grows when it should not.

## Fourth-order Runge-Kutta (RK4)

RK4 samples the slope four times per step and averages them:

```text
k1 = f(t,       y)
k2 = f(t + h/2, y + h k1 / 2)
k3 = f(t + h/2, y + h k2 / 2)
k4 = f(t + h,   y + h k3)
y_next = y + h (k1 + 2 k2 + 2 k3 + k4) / 6
```

Its error scales like `h^4`: halving the step cuts the error by 16. For the same accuracy it is usually far cheaper than Euler, despite four function calls per step.

### Worked example: comparing the two

Exponential decay `dy/dt = -0.5 y`, `y(0) = 1`, has the exact solution `e^(-0.5 t)`. Over 10 s with `h = 0.1`:

```text
method   y(10)        error
exact    0.0067379
Euler    0.0059205    8.2e-4   (12% relative error)
RK4      0.0067379    1.8e-9
```

RK4 is accurate to about eight significant digits with the same 100 steps.

## The real tool: solve_ivp

In practice you use SciPy, which picks step sizes automatically to meet an error tolerance:

```python
from scipy.integrate import solve_ivp

sol = solve_ivp(spring, t_span=(0, 10), y0=[1.0, 0.0],
                method="RK45", rtol=1e-8, atol=1e-10,
                t_eval=np.linspace(0, 10, 501), args=(4.0, 0.3, 1.0))
x = sol.y[0]          # sol.y has shape (n_states, n_times)
print(sol.success, sol.message)
```

Points to know:

- The function signature is `f(t, y)`: **time first**. (The older `odeint` uses `f(y, t)`, a common source of bugs.)
- `rtol` and `atol` control accuracy. The defaults (`1e-3`, `1e-6`) are loose; tighten them and check the answer does not change.
- `t_eval` chooses where you get output; the solver's internal steps are separate.
- `events=` stops the integration when something happens, such as a projectile hitting the ground.
- **Stiff** problems (chemical kinetics with very fast and very slow reactions, circuits with tiny time constants) make explicit methods take absurdly small steps. Switch to `method="LSODA"`, `"Radau"` or `"BDF"`.

## Checking a numerical solution

1. **Convergence**: halve the step (or tighten `rtol`) and confirm the answer stops changing.
2. **Conserved quantities**: for an undamped spring, energy `0.5 k x² + 0.5 m v²` should stay constant. Drift exposes a poor method.
3. **Known limits**: set damping to zero, or compare with an analytic case.
4. **Units and signs**: a decay that grows usually means a sign error in `f`.

## Common mistakes

- **Wrong argument order** (`f(y, t)` versus `f(t, y)`).
- **Forgetting to convert** a second-order equation to first-order form.
- **Trusting Euler** with a step that "looks small". Always do a convergence check.
- **Returning a list** from `f` and then adding it to arrays: return a NumPy array.
- **Loose default tolerances** in `solve_ivp` for long simulations.

## Quick recap

- Write `dy/dt = f(t, y)` with a state vector; second-order becomes first-order by adding velocities.
- Euler: `y + h f`, error ~ `h`. RK4: four slopes, error ~ `h^4`.
- Use `solve_ivp(f, (t0, t1), y0, rtol=..., t_eval=...)` in real work; use a stiff method when steps collapse.
- Check convergence and conserved quantities.

## Try it

In the **Do** tab you will implement Euler and RK4 steppers in NumPy and test them on exponential decay and a spring.


Step ordinary differential equations forward with Euler and Runge-Kutta in NumPy, then solve them properly with scipy.integrate.solve_ivp.

Solving ODEs

How to Solve Differential Equations in PYTHON

Coding a Fourth-Order Runge-Kutta Integrator in Python and Matlab

## Why this matters

Two questions come up constantly in engineering and science work:

- **Where is it zero?** The water depth in a channel that carries a given flow, the friction factor in the Colebrook equation, the temperature at which a reaction's equilibrium constant equals 1, the operating point of a diode circuit. These are **root-finding** problems: solve `f(x) = 0` when no formula gives `x` directly.
- **What is the total?** Energy from a power-versus-time log, work from a force-displacement curve, the area under a chromatography peak, the distance from sampled velocity. These are **integrals**, often of data rather than of a formula.

Both have simple, reliable numerical methods you should understand, plus SciPy functions that do them robustly.

## Root finding: bisection

If `f` is continuous and `f(a)` and `f(b)` have **opposite signs**, there is a root between them. Bisection halves the bracket each step:

```python
def bisect(f, a, b, tol=1e-10):
    fa = f(a)
    if fa * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    while b - a > tol:
        m = 0.5 * (a + b)
        fm = f(m)
        if fa * fm <= 0:
            b = m                 # root is in [a, m]
        else:
            a, fa = m, fm         # root is in [m, b]
    return 0.5 * (a + b)
```

It is slow (one binary digit per step, about 34 steps for `1e-10` on a unit bracket) but **guaranteed** to converge once you have a bracket. That reliability is why it is the fallback in every serious solver.

## Root finding: Newton's method

Newton follows the tangent line to where it crosses zero:

```text
x_{n+1} = x_n - f(x_n) / f'(x_n)
```

```python
def newton(f, df, x0, tol=1e-12, max_iter=50):
    x = x0
    for _ in range(max_iter):
        step = f(x) / df(x)
        x -= step
        if abs(step) < tol:
            return x
    raise RuntimeError("Newton did not converge")
```

Near a root it converges **quadratically**: the number of correct digits roughly doubles each step. For `√2` (root of `x² - 2`) from `x0 = 1`: 1.5, 1.41667, 1.4142157, 1.41421356237. But it needs the derivative, and from a bad start it can jump away or cycle. A good practice is to bracket first, then polish with Newton.

### Worked example: open-channel flow

A rectangular channel 3 m wide must carry `Q = 6 m³/s`. With Manning's equation, the flow depth `y` satisfies `f(y) = Q_manning(y) - 6 = 0`. You cannot solve that for `y` by algebra, but `f(0.1) < 0` and `f(5) > 0`, so bisection on `[0.1, 5]` finds the depth to millimetre precision in about 13 halvings. The same pattern solves the Colebrook friction factor, bubble-point temperatures and diode operating points.

## Integration: the trapezoid rule

For samples `y_i` at points `x_i`, join neighbours with straight lines and add up the trapezoids:

```text
∫ y dx ≈ Σ (x_{i+1} - x_i) (y_i + y_{i+1}) / 2
```

Vectorised, with no loop:

```python
import numpy as np

def trapezoid(y, x):
    y = np.asarray(y, float); x = np.asarray(x, float)
    return float(np.sum(np.diff(x) * (y[1:] + y[:-1]) / 2))
```

It works for **unevenly spaced** data, which is what loggers often produce. Its error falls like `h²`. NumPy has it built in: `np.trapezoid(y, x)` in NumPy 2.0 and later (older versions call it `np.trapz`).

### Worked example: energy from a power log

A motor's power is logged every 10 s for a minute: `[0, 1.2, 2.5, 2.5, 2.4, 1.0, 0]` kW. The energy is `trapezoid(P, t)` = 96 kJ (with `t` in seconds and `P` in kW, the result is in kJ). Units of an integral are the units of `y` times the units of `x`.

## Integration: Simpson's rule

Simpson fits parabolas through each pair of intervals. With an **even** number `n` of equal intervals of width `h`:

```text
∫ f dx ≈ h/3 [ f0 + 4 f1 + 2 f2 + 4 f3 + ... + 4 f_{n-1} + f_n ]
```

Its error falls like `h⁴`. For `∫₀^π sin x dx = 2` with 10 intervals, the trapezoid rule gives 1.9835 and Simpson gives 2.00011.

## The SciPy toolbox

```python
from scipy.optimize import brentq, root_scalar, fsolve
from scipy.integrate import quad, simpson, trapezoid

root = brentq(lambda y: f(y), 0.1, 5)          # bracketing, fast and safe
value, err = quad(np.sin, 0, np.pi)              # adaptive, for functions you can call
area = simpson(y, x=x)                           # for sampled data
```

- `brentq` combines bisection's safety with faster interpolation; it is the default choice for a bracketed scalar root.
- `fsolve` and `root` handle systems of non-linear equations.
- `quad` adapts its sampling to the function and returns an error estimate. Use it when you can evaluate the function; use `trapezoid` or `simpson` when you only have samples.

## Common mistakes

- **No bracket check**: bisection on an interval without a sign change silently returns garbage. Raise an error instead.
- **Newton from a bad start** or with a derivative near zero: it diverges. Bracket first.
- **Odd number of intervals** for Simpson's rule.
- **Forgetting units**: the integral of kW over seconds is kJ, not kWh.
- **Integrating noisy data with high-order rules**: Simpson amplifies noise less than differentiation does, but smoothing or a fitted model may be better.

## Quick recap

- Bisection: needs a sign change, always converges, slowly. Newton: fast near a root, needs `f'` and a good start.
- Trapezoid: `sum(diff(x) * (y[1:] + y[:-1]) / 2)`, any spacing, error ~ `h²`. Simpson: even intervals, error ~ `h⁴`.
- In real work: `brentq` for roots, `quad` for functions, `trapezoid` and `simpson` for data.

## Try it

The **Do** tab asks you to write bisection, a vectorised trapezoid rule and Simpson's rule, and test them on known answers.


Find roots with bisection and Newton's method, integrate sampled and analytic functions with the trapezoid and Simpson rules, and know the SciPy equivalents.

Numerical integration and roots

Finding Zeros of Functions In Python (Bisection Method and Scipy)

Numerical Integration: Discrete Riemann Integrals and Trapezoid Rule

Fit models to data and read parameter uncertainties, step ODEs forward in time, and find roots and integrals numerically, with the SciPy functions that do it for you.