Python & Scientific Computing

Question 1

(a) [5] — R[:, None, :] has shape (N,1,2); R[None,:,:] has shape (1,N,2). (1) Broadcasting rule: compare trailing dimensions; a dimension of size 1 is stretched to match. Align: (N,1,2) vs (1,N,2) → dims (2 vs 2) match, (1 vs N)→N, (N vs 1)→N, giving (N,N,2). (2)

D = R[:, None, :] - R[None, :, :]   # D[i,j] = R[i]-R[j]

Rule stated correctly + code = (2).

(b) [8]

import numpy as np
def total_energy(R, q):
    D = R[:, None, :] - R[None, :, :]      # (N,N,2)
    dist = np.sqrt((D**2).sum(-1))         # (N,N)
    qq = q[:, None] * q[None, :]           # (N,N)
    iu = np.triu_indices(len(q), k=1)      # upper triangle, i<j
    return np.sum(qq[iu] / dist[iu])

• Displacement + distance vectorized (3)
• Charge outer product q[:,None]*q[None,:] (2)
• Diagonal (self-distance = 0) avoided by selecting triu_indices(k=1), which also enforces $i<j$ so no double counting (3). Alternative: mask diagonal with np.fill_diagonal(dist, np.inf) then halve full sum — award full if correct.

(c) [6]

def forces(R, q):
    D = R[:, None, :] - R[None, :, :]      # (N,N,2)
    dist = np.sqrt((D**2).sum(-1))         # (N,N)
    np.fill_diagonal(dist, np.inf)         # zero out self term
    inv3 = dist**(-3)                       # (N,N), diagonal -> 0
    qq = q[:, None] * q[None, :]           # (N,N)
    F = (qq * inv3)[..., None] * D          # (N,N,2)
    return F.sum(axis=1)                     # sum over j

• Correct inv3 with self-term neutralized via inf (3)
• Broadcasting scalar coefficient (N,N,1) against D and summing over j (3).

(d) [5] Proof: Let $\mathbf{f}_{ij} = q_iq_j (\mathbf{r}_i-\mathbf{r}_j)/\|\mathbf{r}_i-\mathbf{r}_j\|^3$. Then $\mathbf{f}_{ji} = q_jq_i(\mathbf{r}_j-\mathbf{r}_i)/\|\cdot\|^3 = -\mathbf{f}_{ij}$ (antisymmetric). (2) $\sum_i \mathbf{F}_i = \sum_i\sum_{j\ne i}\mathbf{f}_{ij} = \sum_{i<j}(\mathbf{f}_{ij}+\mathbf{f}_{ji}) = \sum_{i<j}\mathbf{0} = \mathbf{0}$. (2) Test:

assert np.allclose(forces(R,q).sum(axis=0), 0, atol=1e-9)

(1)

Question 2

(a) [5]

def partial_sums():
    s, n = 0.0, 1
    while True:
        s += 1.0 / n**2
        yield s
        n += 1

Generator holds only s, n (O(1) memory) and produces values lazily; a list of all partials would require O(k) storage and forces a fixed truncation. (2) for explanation, (3) for correct generator.

(b) [5]

import itertools, math
target = math.pi**2 / 6
ks = {10, 100, 1000, 10000}
kmax = max(ks)
errs = {k: abs(s - target)
        for k, s in enumerate(itertools.islice(partial_sums(), kmax), start=1)
        if k in ks}

Dict comprehension (3); single pass through generator using islice (efficient) (2).

(c) [7] Integral comparison: $1/n^2$ is decreasing, so $\int_{k+1}^{\infty}\frac{dx}{x^2} \le \sum_{n=k+1}^\infty \frac1{n^2} \le \int_{k}^{\infty}\frac{dx}{x^2}.$ (3) (setting up upper/lower rectangle bounds against the integral). Evaluate: $\int_a^\infty x^{-2}dx = 1/a$, giving $\frac{1}{k+1}\le \varepsilon_k \le \frac{1}{k}$. (2) For $k=1000$: $\varepsilon_{1000}\approx 1/1000 -$ small correction $\approx 9.995\times10^{-4}$, and $1/1001 \approx 9.990\times10^{-4} \le \varepsilon \le 10^{-3}=1/1000$. ✓ (2)

(d) [3] Summing small terms first (large $n$ → small $1/n^2$) accumulates them into a small running total before adding the dominant early terms; this preserves low-order bits that would be lost to rounding if a large partial sum were already present (the small terms would fall below the ULP of the big accumulator). Hence ascending-magnitude summation reduces catastrophic loss of significance. (3)

Question 3

(a) [6]

import pandas as pd
df = pd.read_csv("measurements.csv", parse_dates=["timestamp"])
stats = df.groupby("sensor_id")["value"].agg(["mean", "std"])

parse_dates (2), groupby.agg giving mean+std ignoring NaN automatically (4).

(b) [6]

g = df.groupby("sensor_id")["value"]
mu = g.transform("mean")
sd = g.transform("std")
z = (df["value"] - mu) / sd
flagged = df[z.abs() > 3]

transform broadcasts group stats back to row shape (3); vectorized z-score + boolean mask (3).

(c) [4] Steps:

• python -m venv .venv && source .venv/bin/activate (isolated interpreter) (1)
• pip freeze > requirements.txt to pin exact versions (or conda env export > environment.yml) (1)
• Commit: source .py/.ipynb, requirements.txt, README. (1)
• .gitignore: .venv/, __pycache__/, large data files / *.csv, .ipynb_checkpoints/ (1)

[
  {"claim": "total_energy matches pairwise loop for a random config", "code": "import numpy as np\nrng=np.random.default_rng(0)\nR=rng.standard_normal((5,2)); q=rng.standard_normal(5)\nD=R[:,None,:]-R[None,:,:]; dist=np.sqrt((D**2).sum(-1))\nqq=q[:,None]*q[None,:]; iu=np.triu_indices(5,1)\nU=np.sum(qq[iu]/dist[iu])\nUloop=sum(q[i]*q[j]/np.linalg.norm(R[i]-R[j]) for i in range(5) for j in range(i+1,5))\nresult=bool(abs(U-Uloop)<1e-9)"},
  {"claim": "net force is zero", "code": "import numpy as np\nrng=np.random.default_rng(1)\nR=rng.standard_normal((6,2)); q=rng.standard_normal(6)\nD=R[:,None,:]-R[None,:,:]; dist=np.sqrt((D**2).sum(-1))\nnp.fill_diagonal(dist,np.inf); inv3=dist**-3\nqq=q[:,None]*q[None,:]; F=((qq*inv3)[...,None]*D).sum(1)\nresult=bool(np.allclose(F.sum(0),0,atol=1e-9))"},
  {"claim": "Basel error at k=1000 lies in [1/1001, 1/1000]", "code": "import math\nS=sum(1.0/n**2 for n in range(1,1001))\neps=math.pi**2/6 - S\nresult=bool(1/1001 <= eps <= 1/1000)"},
  {"claim": "integral bound evaluates to 1/a", "code": "from sympy import integrate, oo, symbols\nx,a=symbols('x a',positive=True)\nresult=bool(integrate(x**-2,(x,a,oo))==1/a)"}
]