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

Scientific and engineering data comes in bulk: 10,000 thermocouple readings, a 2,048-point spectrum, a strain gauge sampled at 1 kHz for an hour. Plain Python lists and `for` loops can handle that, but they are slow and wordy. **NumPy** gives you one object, the `ndarray`, that holds numbers of one type in a single block of memory and runs whole-array operations in compiled code. The result is code that is shorter, closer to the maths on paper, and often 10 to 100 times faster.

Almost every scientific Python library (SciPy, pandas, matplotlib, scikit-learn) is built on NumPy arrays, so this is the foundation for the rest of the course.

## Core idea: an array is a grid of one type

```python
import numpy as np

t = np.array([0.0, 0.5, 1.0, 1.5, 2.0])   # from a list
print(t.shape, t.dtype)                     # (5,) float64

x = np.linspace(0, 10, 101)   # 101 evenly spaced points from 0 to 10 inclusive
k = np.arange(0, 1, 0.25)     # [0.  0.25 0.5  0.75] (end excluded)
grid = np.zeros((3, 4))       # 3 rows, 4 columns of 0.0
```

Key attributes:

- `shape`: the size along each axis, such as `(3, 4)`.
- `ndim`: the number of axes.
- `dtype`: the element type, usually `float64` for measurements or `int64` for counts.

## Vectorisation: operate on the whole array

Arithmetic on arrays works **element by element**, with no loop in your code:

```python
celsius = np.array([20.0, 25.5, 30.1, 18.9])

# Loop version
kelvin_loop = []
for c in celsius:
    kelvin_loop.append(c + 273.15)

# Vectorised version
kelvin = celsius + 273.15
```

NumPy's "universal functions" (ufuncs) such as `np.sqrt`, `np.exp`, `np.sin` and `np.log` also work on every element. Reductions collapse an array to a summary: `x.sum()`, `x.mean()`, `x.std()`, `x.min()`, `x.max()`, `np.argmax(x)`.

### Worked example: RMS voltage

The root-mean-square of a signal is `sqrt(mean(v^2))`. For a sine wave of amplitude 10 V, it should be `10 / sqrt(2) ≈ 7.071` V.

```python
t = np.linspace(0, 1, 10_001)               # one second
v = 10 * np.sin(2 * np.pi * 50 * t)         # 50 Hz, 10 V peak
v_rms = np.sqrt(np.mean(v ** 2))
print(round(v_rms, 3))                       # 7.071
```

Three lines, no loop, and each line reads like the formula.

## Boolean masks: select by condition

Comparisons produce arrays of `True`/`False`. Use them to filter or count:

```python
readings = np.array([4.98, 5.02, 7.40, 5.01, 4.99, -1.00])
valid = (readings > 0) & (readings < 6)     # use & and |, not "and"/"or"
print(readings[valid])                       # [4.98 5.02 5.01 4.99]
print(valid.sum())                           # 4 valid readings
```

Brackets matter: `(readings > 0) & (readings < 6)`. Without them Python's operator precedence evaluates `0 & readings` first.

## Broadcasting: combine different shapes

When two arrays have different shapes, NumPy tries to **broadcast** them: it lines the shapes up from the right, and any axis of size 1 is stretched to match the other array. No data is copied.

```python
data = np.array([[1.0, 2.0, 3.0],
                 [4.0, 5.0, 6.0]])           # shape (2, 3)
col_means = data.mean(axis=0)                # shape (3,)
centred = data - col_means                   # (2, 3) - (3,) -> (2, 3)
```

`axis=0` means "collapse the rows", giving one value per column. To subtract **row** means you need the shape `(2, 1)`, which `keepdims=True` or `[:, None]` gives you:

```python
row_means = data.mean(axis=1, keepdims=True) # shape (2, 1)
centred_rows = data - row_means
```

### Worked example: all pairwise distances

Given `n` points in 2D (shape `(n, 2)`), the distance matrix `D[i, j]` can be built with no loops:

```python
pts = np.array([[0.0, 0.0], [3.0, 4.0], [6.0, 8.0]])
diff = pts[:, None, :] - pts[None, :, :]     # (3,1,2) - (1,3,2) -> (3,3,2)
D = np.sqrt((diff ** 2).sum(axis=-1))        # (3,3)
print(D[0, 1], D[0, 2])                      # 5.0 10.0
```

Inserting a new axis with `None` (also called `np.newaxis`) is the most common broadcasting trick.

## Common mistakes

- **Using lists where you meant arrays.** `[1, 2] * 2` is `[1, 2, 1, 2]`; `np.array([1, 2]) * 2` is `[2, 4]`.
- **Integer arrays.** `np.array([1, 2, 3])` is `int64`. Assigning `0.5` into it silently truncates to `0`. Create float arrays (`np.array([1.0, 2, 3])` or `dtype=float`) for measurements.
- **Views, not copies.** Slicing (`a[2:5]`) returns a view; changing it changes the original. Call `.copy()` when you need an independent array.
- **Shape surprises.** A `(3,)` array and a `(3, 1)` array broadcast to `(3, 3)`. Print `.shape` whenever a result looks too big.
- **Writing a loop and wrapping it in `np.vectorize`.** It is convenient but runs a Python loop underneath; it is not real vectorisation.

## Quick recap

- An `ndarray` holds one dtype in a fixed shape; arithmetic and ufuncs act on every element.
- Reductions (`sum`, `mean`, `std`, `max`) take an `axis` argument.
- Boolean masks filter and count without loops.
- Broadcasting aligns shapes from the right and stretches size-1 axes; `[:, None]` adds an axis.

## Try it

Open the **Do** tab and replace three loop-style functions with vectorised NumPy versions: a unit conversion, an RMS calculation and a pairwise distance matrix.


Replace slow Python loops with NumPy array math, and use broadcasting to combine arrays of different shapes.

Arrays and vectorisation

Python NumPy Tutorial for Beginners

Numpy Array Broadcasting In Python Explained

## Why this matters

A surprising amount of engineering and science ends up as a set of linear equations. Kirchhoff's laws for a circuit, force balance at the joints of a truss, mass balances in a chemical process, calibrating a spectrometer against several standards: each gives you **n equations in n unknowns**. By hand, three unknowns is tedious and ten is hopeless. With NumPy it is one line, and the harder part becomes writing the system down correctly and checking the answer.

## Core idea: write it as A x = b

Take a two-loop circuit. A 12 V source drives loop 1, which shares a 4 Ω resistor with loop 2. Writing Kirchhoff's voltage law around each loop with mesh currents `i1` and `i2`:

```text
loop 1:  (2 + 4) i1 -       4 i2 = 12
loop 2:      -4 i1 + (4 + 6) i2 =  0
```

Collect the coefficients into a matrix `A`, the unknowns into a vector `x`, and the right-hand sides into `b`:

```text
A = [[ 6, -4],      x = [i1,      b = [12,
     [-4, 10]]           i2]           0]
```

Each **row** of `A` is one equation, each **column** belongs to one unknown.

```python
import numpy as np

A = np.array([[6.0, -4.0],
              [-4.0, 10.0]])
b = np.array([12.0, 0.0])

x = np.linalg.solve(A, b)
print(x)          # [2.7273 1.0909]  amps
```

## Matrix products and the @ operator

`A @ x` is the matrix-vector product (`np.dot` does the same). Plain `A * x` is element-wise and broadcasts, which is almost never what you want here.

```python
print(A @ x)      # [12.  0.]  reproduces b
```

Other tools you will use:

- `A.T`: transpose.
- `np.linalg.det(A)`: determinant (zero means no unique solution).
- `np.linalg.inv(A)`: the inverse. Avoid it for solving; `solve` is faster and more accurate.
- `np.linalg.norm(v)`: the length of a vector.
- `np.linalg.eig(A)`: eigenvalues and eigenvectors (vibration modes, principal stresses).

## Always check: the residual

A solver always returns numbers. Whether they are trustworthy is your job. The **residual** `r = A x - b` should be tiny compared with `b`:

```python
r = A @ x - b
print(np.linalg.norm(r))                 # about 1e-15
print(np.allclose(A @ x, b))             # True
```

If you typed a coefficient into the wrong row, the residual against your **original** equations will still be small (the solver solved what you gave it). So also sanity-check physics: currents should satisfy Kirchhoff's current law, forces should sum to zero, concentrations should be positive.

## Conditioning: when small errors blow up

Some systems are **ill-conditioned**: tiny changes in `b` (measurement noise) cause huge changes in `x`. The condition number measures this:

```python
print(np.linalg.cond(A))        # about 3.5: well conditioned

B = np.array([[1.0, 1.0],
              [1.0, 1.0001]])
print(np.linalg.cond(B))        # about 40,000: nearly singular
```

Rule of thumb: you can lose about `log10(cond)` significant digits. A condition number of `1e4` with data good to 3 digits means the answer may have none. Ill-conditioning usually signals two nearly identical equations, for example two calibration standards with almost the same concentration.

## Overdetermined systems: least squares

With more equations than unknowns (more measurements than parameters) there is usually no exact solution. **Least squares** finds the `x` that minimises `||A x - b||`.

### Worked example: a calibration line

A sensor's output `y` (mV) is measured at known temperatures `T` (°C). Fit `y = m T + c`:

```python
T = np.array([0.0, 10.0, 20.0, 30.0, 40.0])
y = np.array([0.02, 0.41, 0.79, 1.22, 1.60])

A = np.column_stack([T, np.ones_like(T)])   # columns: T and 1
coef, residuals, rank, sv = np.linalg.lstsq(A, y, rcond=None)
m, c = coef
print(round(m, 4), round(c, 4))             # 0.0397 0.014
```

The **design matrix** `A` has one column per parameter. The same trick fits any model that is linear in its parameters (polynomials, sums of known functions) and is the basis of curve fitting in Module 3.

SciPy adds more specialised tools (`scipy.linalg.lu_factor` to reuse a factorisation, sparse solvers in `scipy.sparse.linalg` for systems with millions of unknowns), but `numpy.linalg` covers most day-to-day work.

## Common mistakes

- **`*` instead of `@`.** Element-wise multiplication silently gives the wrong numbers.
- **Integer matrices.** Build matrices with floats (`6.0`, not `6`) so later edits do not truncate.
- **Using `inv(A) @ b`.** It works but is slower and less accurate than `solve(A, b)`.
- **Mixing up rows and columns.** One row per equation, one column per unknown, in the same order every row.
- **Ignoring `LinAlgError: Singular matrix`.** It means the equations are not independent: a missing equation or a duplicated one.

## Quick recap

- Write the problem as `A x = b`, then `x = np.linalg.solve(A, b)`.
- Use `@` for matrix products.
- Check `np.allclose(A @ x, b)` and the physics; check `np.linalg.cond(A)` when inputs are noisy.
- For more equations than unknowns, build a design matrix and use `np.linalg.lstsq`.

## Try it

In the **Do** tab you will solve a three-loop circuit, verify it with the residual, and fit a calibration line with least squares.


Set up a real engineering system as A x = b, solve it with numpy.linalg, and check the answer with residuals and the condition number.

Linear algebra with NumPy

Numpy Arrays - Fast Linear Algebra in Python

Inverse matrices, column space and null space | Chapter 7, Essence of linear algebra

## Why this matters

Every measurement has noise. A balance reads 12.503 g one time and 12.498 g the next; a voltmeter's last digit flickers; a thermocouple drifts. Before you trust a result you need to know how that noise flows through your calculation. Sometimes a formula answers that. Often it is easier, and safer, to **simulate**: generate thousands of plausible measurements, push each one through the calculation, and look at the spread of results. This is the **Monte Carlo** method, and NumPy makes it fast.

## Core idea: a seeded generator

Modern NumPy uses a `Generator` object. Create one with a **seed** so your results are reproducible:

```python
import numpy as np

rng = np.random.default_rng(42)      # same seed -> same numbers every run

rng.random(3)                # 3 uniform numbers in [0, 1)
rng.normal(10.0, 0.2, 5)     # 5 draws, mean 10.0, standard deviation 0.2
rng.uniform(-1, 1, (2, 3))   # a 2x3 array of uniform numbers in [-1, 1)
rng.integers(1, 7, 10)       # 10 dice rolls (upper bound excluded)
rng.poisson(4.0, 1000)       # counts, such as photons per detector window
```

Why seed? A colleague (or you in six months) should get **exactly** the same figures from the same code. Without a seed, every run differs a little and nobody can tell a bug from bad luck.

Older code uses `np.random.seed(0)` and `np.random.rand()`. It still works, but the `Generator` API is the recommended one: it has better statistical properties and keeps separate streams independent. Pass the `rng` object into functions instead of relying on a hidden global state.

## Simulating measurement noise

Suppose a length is truly 50.00 mm and the caliper has a random error with standard deviation 0.02 mm. Ten simulated readings:

```python
rng = np.random.default_rng(7)
true_length = 50.00
readings = true_length + rng.normal(0, 0.02, size=10)
print(readings.round(3))
print(readings.mean(), readings.std(ddof=1))
```

Use `ddof=1` for the **sample** standard deviation of measured data (it divides by `n - 1`). The standard error of the mean is `std / sqrt(n)`: averaging 10 readings shrinks the random error by a factor of about 3.

## Monte Carlo uncertainty propagation

### Worked example: density of a sample

You measure a cylinder's mass `m = 25.40 ± 0.05 g` and volume `V = 9.40 ± 0.10 cm³` (± is one standard deviation). The density is `ρ = m / V`. What is its uncertainty?

**By formula** (first-order propagation for a quotient):

```text
σρ / ρ = sqrt( (σm/m)^2 + (σV/V)^2 )
       = sqrt( (0.05/25.40)^2 + (0.10/9.40)^2 ) ≈ 0.0108
ρ = 2.702 g/cm³,   σρ ≈ 0.029 g/cm³
```

**By simulation:**

```python
rng = np.random.default_rng(2026)
n = 100_000
m = rng.normal(25.40, 0.05, n)
V = rng.normal(9.40, 0.10, n)
rho = m / V                       # 100,000 densities, vectorised
print(rho.mean(), rho.std())      # about 2.702 and 0.029
lo, hi = np.percentile(rho, [2.5, 97.5])
print(lo, hi)                     # a 95% interval
```

The two agree. The simulation wins when the formula gets messy: non-linear functions, correlated inputs, or non-normal errors (a uniform resolution error, a skewed drift). You change one line instead of re-deriving a derivative. `np.percentile` gives you an interval directly from the simulated results without assuming they are normal.

## How many samples?

Monte Carlo error shrinks like `1 / sqrt(n)`. With `n = 10,000` the estimated standard deviation is good to about 1%; with `n = 1,000,000` to about 0.1%. Vectorised NumPy handles a million samples in a fraction of a second, so start with 100,000 and increase if the answer still moves between seeds.

## Estimating a probability

Simulation also answers "how often" questions. If a part's diameter is `N(10.00, 0.03)` mm and the spec is 9.95 to 10.05 mm, what fraction is out of spec?

```python
d = np.random.default_rng(1).normal(10.00, 0.03, 200_000)
out = (d < 9.95) | (d > 10.05)
print(out.mean())          # about 0.096, the fraction of True values
```

The mean of a Boolean array is the fraction of `True` values, which is your estimated probability.

## Common mistakes

- **No seed.** Results change every run; figures cannot be reproduced. Seed once at the top and pass `rng` around.
- **Re-seeding inside a loop.** `default_rng(0)` inside a loop gives the same "random" numbers every iteration.
- **`std()` versus `std(ddof=1)`.** For a handful of measured values use `ddof=1`; for a large simulated population the difference is negligible.
- **Confusing standard deviation with standard error.** SD describes individual readings; SE (`SD / sqrt(n)`) describes the uncertainty of the mean.
- **Too few samples.** With 100 samples, a 95% interval can move noticeably between seeds.

## Quick recap

- `rng = np.random.default_rng(seed)` then `rng.normal`, `rng.uniform`, `rng.integers`, `rng.poisson`.
- Add simulated noise to a true value to test your analysis before real data arrives.
- Monte Carlo propagation: sample every input, compute the output array, report its mean, SD and percentiles.
- The mean of a Boolean array estimates a probability.

## Try it

The **Do** tab asks you to write a seeded noise simulator and a Monte Carlo density calculation, then compare it with the formula.


Use NumPy's seeded random generator to simulate measurement noise and propagate uncertainty through a calculation with Monte Carlo.

Random numbers and simulation

What is Monte Carlo Simulation?

6. Monte Carlo Simulation

Arrays and vectorised maths, broadcasting, linear systems with numpy.linalg, and seeded random numbers for simulating measurement noise.