Level 4 — ApplicationScientific Computing (Python)

Scientific Computing (Python)

⏱ 60 minutes60 marksprintable — key stays hidden on paper

Level 4 — Application (novel problems, no hints) Time limit: 60 minutes Total marks: 60

Answer all questions. Show working. Code should be idiomatic NumPy/SciPy. Numerical answers to 4 significant figures unless stated.

Question 1 — Broadcasting & Vectorization (12 marks)

You are given points, a NumPy array of shape (N, 3) holding 3D coordinates.

(a) Write a single vectorized expression (no Python loops) that produces the (N, N) matrix D of squared Euclidean distances between all pairs of points. State the shapes of the intermediate arrays your expression creates. (5)

(b) For N = 2000, estimate the peak temporary memory (in MB, float64) allocated by the broadcasting approach in (a), and explain why this can be a "gotcha". (3)

(c) Consider this attempted normalization:

X = np.random.rand(4, 3)
col_means = X.mean(axis=0)      # shape (3,)
Xc = X - col_means
row_norms = np.linalg.norm(Xc, axis=1)   # shape (4,)
Y = Xc / row_norms

The last line raises no error but gives the wrong result. Explain precisely why, and give the corrected line. (4)

Question 2 — Linear Algebra & Floating Point (12 marks)

(a) For the function $f(x) = \dfrac{1 - \cos x}{x^2}$ , evaluate it naively in float64 at $x = 10^{-8}$ . State the computed value, the true limiting value as $x\to 0$ , and name the phenomenon responsible. (4)

(b) Rewrite $f(x)$ into an algebraically equivalent form that is numerically stable for small $x$ , and justify why it avoids the problem. (4)

(c) You must solve $A\mathbf{x} = \mathbf{b}$ where $A$ is a $500\times500$ symmetric positive-definite matrix, repeatedly for 1000 different right-hand sides $\mathbf{b}$ with the same $A$ . Describe the most efficient standard approach (name the factorization and the SciPy calls) and give the reason it beats calling np.linalg.solve 1000 times. (4)

Question 3 — Numerical Methods From Scratch (14 marks)

(a) Implement newton(f, df, x0, tol=1e-12, maxit=50) returning the root and the iteration count. Include a fallback: if a step lands where df is near zero (|df| < 1e-14), raise ValueError. (6)

(b) Apply Newton's method by hand to $f(x) = x^3 - 2x - 5$ starting at $x_0 = 2$ . Compute $x_1$ and $x_2$ exactly (as fractions or 6-dp decimals). (4)

(c) Composite Simpson's rule requires an even number of intervals $n$ . Explain what your implementation should do if the caller passes odd $n$ , and state the order of the error term (in $h$ ) for Simpson vs the trapezoidal rule. (4)

Question 4 — SciPy Integration & Optimization (12 marks)

Consider the ODE system (damped pendulum): $\ddot\theta + 0.3\,\dot\theta + \sin\theta = 0,\qquad \theta(0)=1,\ \dot\theta(0)=0.$

(a) Write the first-order state form and a Python rhs(t, y) function suitable for scipy.integrate.solve_ivp. State exactly how you would request the solution sampled at t = np.linspace(0, 20, 200) using the DOP853 method. (6)

(b) You have noisy data (t_data, y_data) believed to follow $y = A e^{-kt}\cos(\omega t + \phi)$ . Give the scipy.optimize.curve_fit call to estimate $(A, k, \omega, \phi)$ , and explain why supplying a p0 initial guess matters for this particular model. (4)

(c) curve_fit returns popt, pcov. Explain what pcov contains and how to extract the standard error (1σ) of parameter $k$ from it. (2)

Question 5 — Signal / FFT / Sparse (10 marks)

(a) A signal is sampled at fs = 1000 Hz for 1 second (1000 samples). After X = np.fft.rfft(x), what is the length of X, and what physical frequency (Hz) does index k = 50 correspond to? Give the general formula. (4)

(b) You build a $10^6 \times 10^6$ matrix where each row has at most 5 nonzero entries. Explain which sparse format (CSR vs CSC) you'd choose for repeated matrix–vector products A @ v, and roughly how many float64 + index values are stored versus the dense count. (4)

(c) State one reason scipy.linalg.solve may be preferred over numpy.linalg.solve. (2)

Answer keyMark scheme & solutions

Question 1 (12)

(a) (5)

diff = points[:, None, :] - points[None, :, :]   # shape (N, N, 3)
D = (diff**2).sum(axis=-1)                        # shape (N, N)

points[:,None,:] → shape (N,1,3); points[None,:,:] → (1,N,3). (1)
Broadcasting difference → (N,N,3). (2)
Sum over last axis → (N,N). (1) Correct final expression (1).

(b) (3) The intermediate diff has N*N*3 = 2000·2000·3 = 1.2×10^7 float64 = 1.2e7 × 8 bytes = 9.6×10^7 B ≈ 91.55 MB (or ~96 MB using 10^6). (2) Gotcha: the temporary is much larger than input/output; broadcasting silently materializes a big 3D array → possible MemoryError / cache thrash. (1)

(c) (4) Xc is (4,3), row_norms is (4,). Broadcasting (4,3)/(4,) aligns row_norms against the last axis (length 3 ≠ 4) → error unless by coincidence... actually it raises? No: shapes (4,3) and (4,) → trailing dims 3 vs 4 mismatch → broadcasting error. The intended per-row division needs a column vector. Corrected: Y = Xc / row_norms[:, None] (shape (4,1) broadcasts across columns). (2) for diagnosis, (2) for fix. (Accept: student notes the raise; key point is the reshape fix.)

Question 2 (12)

(a) (4) At $x=10^{-8}$ , $\cos x \approx 1 - 5\times10^{-17}$ , which rounds to 1.0 in float64, so 1 - cos x = 0.0, giving $f = 0/10^{-16} = 0.0$ . (2) True limit $\lim_{x\to0} f = 1/2 = 0.5$ . (1) Phenomenon: catastrophic cancellation (subtracting nearly-equal numbers). (1)

(b) (4) Use half-angle: $1-\cos x = 2\sin^2(x/2)$ , so $f(x)=\frac{2\sin^2(x/2)}{x^2}=\frac12\left(\frac{\sin(x/2)}{x/2}\right)^2.$ (2) Stable because sin(x/2) for small argument is computed accurately (no subtraction of near-equal quantities); the ratio sinc-like term → 1. (2)

(c) (4) Compute a Cholesky factorization once: c, low = scipy.linalg.cho_factor(A), then for each rhs x = scipy.linalg.cho_solve((c,low), b). (2) (Accept lu_factor/lu_solve for general A.) Factorization is $O(n^3)$ done once; each solve is only $O(n^2)$ . np.linalg.solve re-factorizes every call → 1000× the $O(n^3)$ cost. (2)

Question 3 (14)

(a) (6)

def newton(f, df, x0, tol=1e-12, maxit=50):
    x = x0
    for i in range(1, maxit+1):
        d = df(x)
        if abs(d) < 1e-14:
            raise ValueError("derivative too small")
        x_new = x - f(x)/d
        if abs(x_new - x) < tol:
            return x_new, i
        x = x_new
    return x, maxit

Marks: loop/update (2), derivative guard raising ValueError (2), convergence test on step size (1), returns root+count (1).

(b) (4) $f(x)=x^3-2x-5,\ f'(x)=3x^2-2$ .

$f(2)=8-4-5=-1$ , $f'(2)=10$ . $x_1 = 2 - (-1)/10 = 2.1$ . (2)
$f(2.1)=9.261-4.2-5=0.061$ , $f'(2.1)=3(4.41)-2=11.23$ . $x_2 = 2.1 - 0.061/11.23 = 2.094568...$ ≈ 2.094568. (2)

(c) (4) For odd $n$ : either raise an error / require even $n$ , or handle the last interval separately (e.g. add one interval treated by trapezoid or a 3/8 rule), or increment $n$ by 1. (2) Simpson error is $O(h^4)$ (globally); trapezoidal is $O(h^2)$ . (2)

Question 4 (12)

(a) (6) State $y=[\theta,\ \omega]$ with $\omega=\dot\theta$ : $\dot\theta=\omega,\qquad \dot\omega = -0.3\,\omega - \sin\theta.$ (2)

def rhs(t, y):
    th, w = y
    return [w, -0.3*w - np.sin(th)]
 
t_eval = np.linspace(0, 20, 200)
sol = solve_ivp(rhs, (0, 20), [1.0, 0.0],
                method='DOP853', t_eval=t_eval)

Marks: rhs function (2), correct t_span=(0,20), y0=[1,0], method='DOP853', t_eval (2).

(b) (4)

def model(t, A, k, w, phi):
    return A*np.exp(-k*t)*np.cos(w*t + phi)
popt, pcov = curve_fit(model, t_data, y_data, p0=[1,0.1,1,0])

(2). p0 matters because the model is nonlinear and oscillatory: the cost surface has many local minima (esp. in $\omega,\phi$ ); a poor start converges to a wrong minimum or fails. (2)

(c) (2) pcov is the estimated covariance matrix of popt. Standard error of $k$ (the 2nd parameter) = np.sqrt(pcov[1,1]). (2)

Question 5 (10)

(a) (4) rfft of a length-N=1000 real signal has length N//2 + 1 = 501. (2) Frequency of bin k: $f_k = k\,f_s/N = k\cdot 1000/1000 = k$ Hz, so k=50 → 50 Hz. (2)

(b) (4) Choose CSR for fast A @ v (row-wise access). (2) Nonzeros ≤ 5·10^6 → store ~5×10^6 float64 values + 5×10^6 column indices + 10^6+1 row pointers, i.e. ~10^7 numbers, versus dense (10^6)^2 = 10^{12} entries — a 10^5× saving. (2)

(c) (2) scipy.linalg is built on LAPACK with more/optimized routines and can be more numerically stable / offers structured solvers (assume_a, cho, etc.); any valid reason. (2)

[
 {"claim":"Q1b temporary memory 2000x2000x3 float64 = ~91.55 MB", "code":"result = abs(2000*2000*3*8/1024/1024 - 91.55) < 0.1"},
 {"claim":"Q2a naive f(1e-8) evaluates to 0.0 not 0.5", "code":"import math; x=1e-8; result = ((1-math.cos(x))/x**2 == 0.0)"},
 {"claim":"Q3b Newton x1=2.1 and x2 approx 2.094568", "code":"f=lambda x:x**3-2*x-5; df=lambda x:3*x**2-2; x1=2-f(2)/df(2); x2=x1-f(x1)/df(x1); result = abs(x1-2.1)<1e-12 and abs(x2-2.094568)<1e-5"},
 {"claim":"Q5 rfft length 501 and bin50 = 50Hz", "code":"N=1000; fs=1000; L=N//2+1; f50=50*fs/N; result = (L==501) and (abs(f50-50)<1e-9)"}
]