Scientific Computing (Python)

Question 1 (12)

(a) RK4 stages: (4, 1 each) $k_1 = f(t_n, y_n)$ $k_2 = f\!\left(t_n+\tfrac{h}{2},\, y_n+\tfrac{h}{2}k_1\right)$ $k_3 = f\!\left(t_n+\tfrac{h}{2},\, y_n+\tfrac{h}{2}k_2\right)$ $k_4 = f\!\left(t_n+h,\, y_n+h\,k_3\right)$ $y_{n+1} = y_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$ Why: weighted average of slopes sampled at endpoints (weight 1) and midpoint (weight 2) matches Taylor expansion to $O(h^4)$.

(b) (5)

import numpy as np
def rk4(f, y0, t0, t1, n):
    t = np.linspace(t0, t1, n+1)
    h = (t1 - t0) / n
    y0 = np.atleast_1d(np.asarray(y0, dtype=float))
    y = np.empty((n+1, y0.size))
    y[0] = y0
    for i in range(n):
        k1 = f(t[i], y[i])
        k2 = f(t[i] + h/2, y[i] + h/2*k1)
        k3 = f(t[i] + h/2, y[i] + h/2*k2)
        k4 = f(t[i] + h,   y[i] + h*k3)
        y[i+1] = y[i] + h/6*(k1 + 2*k2 + 2*k3 + k4)
    return t, y

Marks: linspace/h (1), init/atleast_1d (1), 4 k-stages (2), update+return (1).

(c) (3) Local truncation error per step is $O(h^5)$; accumulated over $\sim 1/h$ steps gives global error $O(h^4)$ — hence "4th order." Halving $h$ reduces error ~16×.

Question 2 (10)

(a) (3) $\int_a^b f\,dx \approx \frac{h}{3}\Big[f(x_0) + 4\!\!\sum_{i\,odd}\!\!f(x_i) + 2\!\!\sum_{i\,even}\!\!f(x_i) + f(x_n)\Big]$ (odd indices weight 4, interior even weight 2, endpoints weight 1).

(b) (5)

import numpy as np
def simpson(f, a, b, n):
    if n % 2 != 0:
        raise ValueError("n must be even")
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    w = np.ones(n+1)
    w[1:-1:2] = 4
    w[2:-1:2] = 2
    return h/3 * np.dot(w, y)

Marks: even check (1), linspace/h (1), weight array (2), dot/return (1).

(c) (2) $n=2,h=1$: nodes $0,1,2$, $f=0,1,8$. Estimate $=\frac{1}{3}(0+4\cdot1+8)=\frac{12}{3}=4$. Exact $\int_0^2 x^3 = 4$. Simpson is exact for cubics.

Question 3 (10)

(a) (3) $x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}$. Fails when $f'(x_n)=0$ (division by zero) / near-flat derivative, or poor initial guess causing divergence/cycling.

(b) (4)

def newton(f, df, x0, tol=1e-10, maxit=100):
    x = x0
    for _ in range(maxit):
        fx = f(x)
        if abs(fx) < tol:
            return x
        x = x - fx / df(x)
    raise RuntimeError("did not converge")

Marks: loop (1), convergence check (1), update (1), raise (1).

(c) (3) $f=x^2-2,\ f'=2x$. $x_1 = 1 - \frac{1-2}{2} = 1 + \frac12 = \frac{3}{2}$. $x_2 = \frac32 - \frac{(9/4)-2}{3} = \frac32 - \frac{1/4}{3} = \frac32 - \frac{1}{12} = \frac{17}{12}$. ($\frac{17}{12}\approx1.4167$, close to $\sqrt2$.)

Question 4 (10)

(a) (3) (1) Align shapes right-to-left; prepend 1's to the shorter shape. (2) In each dimension sizes must be equal or one of them must be 1 (the size-1 dim is stretched/broadcast).

(b) (2) Result shape (3,4). Element $[i,j] = a[i,0] + b[0,j]$.

(c) (5)

D2 = ((X[:, None, :] - X[None, :, :])**2).sum(-1)

Explain: X[:,None,:] is (m,1,d), X[None,:,:] is (1,m,d); broadcasting produces all $(i,j)$ differences as an (m,m,d) array in a single C-level operation, then summing over the last axis gives the (m,m) squared-distance matrix — no Python-level loop, so it runs in vectorized compiled code.

Question 5 (10)

(a) (4) Catastrophic cancellation: subtracting two nearly-equal floating-point numbers loses leading significant digits, hugely amplifying relative error. (2) Stable smaller root: compute $q = -\tfrac12\big(b + \operatorname{sign}(b)\sqrt{b^2-4ac}\big)$, then roots are $x_1=q/a$ and $x_2=c/q$; choosing the $+\operatorname{sign}(b)$ avoids subtracting near-equal quantities. (2)

(b) (3) (1) Forming inv(A) is more expensive and less accurate (extra rounding); solve uses LU factorization directly. (2) solve is numerically more stable / better conditioned and preserves structure. (Either two valid reasons.)

(c) (3) Returns the solution vector $x$. Solve: $2x_1+x_2=3,\ x_1+3x_2=5$. $\det=5$. $x_1=(3\cdot3-1\cdot5)/5=4/5$, $x_2=(2\cdot5-1\cdot3)/5=7/5$. So $x=(0.8,\ 1.4)$.

Question 6 (8)

(a) (3) np.fft.fft computes the Discrete Fourier Transform: $X_k=\sum_{n}x_n e^{-2\pi i kn/N}$, returning complex amplitudes for each frequency bin. fftfreq(n, d) returns the sample frequencies (cycles per unit) in the FFT bin order, with spacing $1/(n\cdot d)$, including negative frequencies in the second half.

(b) (5, 1 each) (i) scipy.integrate — (ii) scipy.optimize — (iii) scipy.linalg — (iv) scipy.stats — (v) scipy.sparse.

[
  {"claim": "Simpson n=2 estimate of int_0^2 x^3 = 4", "code": "h=1; xs=[0,1,2]; f=lambda x:x**3; est=h/3*(f(xs[0])+4*f(xs[1])+f(xs[2])); result = (est==4)"},
  {"claim": "Newton x1=3/2, x2=17/2... =17/12 for x^2-2 from x0=1", "code": "from sympy import Rational; f=lambda x:x**2-2; df=lambda x:2*x; x0=Rational(1); x1=x0-f(x0)/df(x0); x2=x1-f(x1)/df(x1); result = (x1==Rational(3,2)) and (x2==Rational(17,12))"},
  {"claim": "Solution of the 2x2 system is (4/5, 7/5)", "code": "A=Matrix([[2,1],[1,3]]); b=Matrix([3,5]); x=A.solve(b); result = (x[0]==Rational(4,5)) and (x[1]==Rational(7,5))"},
  {"claim": "Broadcasting (3,1)+(1,4) gives shape (3,4)", "code": "import numpy as np; a=np.zeros((3,1)); b=np.zeros((1,4)); result = ((a+b).shape==(3,4))"}
]