Scientific Computing (Python)

Question 1 (20)

(a) (6) By Vieta: $x_1+x_2 = -b/a = 10^8$, $x_1 x_2 = c/a = 1$. So one root is large $\approx 10^8$, the other small $\approx 10^{-8}$. Exactly: $\sqrt{b^2-4ac}=\sqrt{10^{16}-4}\approx 10^8(1-2\times10^{-16})$. (2)

• Large root: $x_+ = \frac{10^8 + \sqrt{\cdots}}{2}\approx 10^8$ — computed by adding two nearly-equal positive quantities → no cancellation. (2)
• Small root using $-b-\sqrt{b^2-4ac} = 10^8 - 10^8(1-2\times10^{-16})$: subtraction of two nearly-equal numbers → catastrophic cancellation; only ~0 significant digits of the difference survive because the true difference is at the $10^{-16}$ relative scale. Relative error blows up to $O(1)$. (2)

(b) (5) Define $q=-\tfrac12\big(b+\operatorname{sgn}(b)\sqrt{b^2-4ac}\big)$. Here the two terms $b$ and $\operatorname{sgn}(b)\sqrt{\cdots}$ have the same sign, so the sum is an addition of like-signed numbers → no cancellation. (2) Then $x_1=q/a$ is the well-conditioned large root. The small root follows from the product $x_1 x_2=c/a$, i.e. $x_2=c/q$ — obtained by division, which is backward-stable and introduces no cancellation. (2) Since here $b<0$, $\operatorname{sgn}(b)=-1$, $q=-\tfrac12(b-\sqrt{b^2-4ac})=\tfrac12(|b|+\sqrt{\cdots})\approx10^8$; $x_2=1/q\approx10^{-8}$. (1)

(c) (6)

import numpy as np
def stable_roots(a, b, c):
    disc = np.sqrt(b*b - 4*a*c)
    q = -0.5*(b + np.sign(b)*disc)
    return q/a, c/q                      # (large, small)
 
def newton(f, df, x0, tol=1e-14, maxit=100):
    x = x0
    for _ in range(maxit):
        fx = f(x)
        dx = fx/df(x)
        x -= dx
        if abs(dx) < tol:
            break
    return x
 
a,b,c = 1.0, -1e8, 1.0
big, small = stable_roots(a,b,c)
f  = lambda x: a*x*x + b*x + c
df = lambda x: 2*a*x + b
root = newton(f, df, small)              # refine the small root

(3 for stable_roots, 2 for correct newton, 1 for refining the smaller-magnitude root.) Order of convergence: quadratic ($|e_{n+1}|\le C|e_n|^2$) near a simple root. At a double root $f'=0$ there, convergence degrades to linear (rate $1/2$). (1 embedded)

(d) (3) Stable formula: relative error $\approx$ a few $\varepsilon \sim 10^{-16}$ (only rounding in division/sqrt). (1) Naive formula for the small root: cancellation amplifies relative error by $\sim |b|/\text{true difference}\approx 10^8/10^{-16}$... concretely the subtraction $10^8-10^8$ loses ~16 digits, so the naive small root has relative error $O(1)$ (often completely wrong / even $0$). (2)

Question 2 (22)

(a) (5) $\lambda_j = 2-2\cos(j\pi/4)$:

• $j=1$: $\cos(\pi/4)=\tfrac{\sqrt2}{2}$, $\lambda_1 = 2-\sqrt2 \approx 0.5858$
• $j=2$: $\cos(\pi/2)=0$, $\lambda_2 = 2$
• $j=3$: $\cos(3\pi/4)=-\tfrac{\sqrt2}{2}$, $\lambda_3 = 2+\sqrt2 \approx 3.4142$ (2)

$\omega_j=\sqrt{\lambda_j}$: $\omega_1\approx0.7654$, $\omega_2=\sqrt2\approx1.4142$, $\omega_3\approx1.8478$. (1)

Numerically: w2, V = np.linalg.eigh(K). eigh exploits that $K$ is real symmetric → guarantees real eigenvalues, orthonormal eigenvectors, and uses a more stable/efficient symmetric algorithm. eig (general) may return complex values with round-off imaginary parts and non-orthogonal vectors. (2)

(b) (6)

Energy: $\dfrac{dE}{dt}= \dot{\mathbf u}^\top\ddot{\mathbf u} + \dot{\mathbf u}^\top K\mathbf u$ (using symmetry of $K$: $\tfrac{d}{dt}\tfrac12\mathbf u^\top K\mathbf u=\dot{\mathbf u}^\top K\mathbf u$). (2) Substitute $\ddot{\mathbf u}=-K\mathbf u$: $\dfrac{dE}{dt}= \dot{\mathbf u}^\top(-K\mathbf u)+\dot{\mathbf u}^\top K\mathbf u = 0.$ Energy conserved. (2)

(c) (6)

def rk4(f, y0, t0, tf, h):
    n = int(round((tf-t0)/h))
    t = t0; y = np.array(y0, float)
    T=[t]; Y=[y.copy()]
    for _ in range(n):
        k1=f(t,y); k2=f(t+h/2,y+h/2*k1)
        k3=f(t+h/2,y+h/2*k2); k4=f(t+h,y+h*k3)
        y = y + h/6*(k1+2*k2+2*k3+k4); t += h
        T.append(t); Y.append(y.copy())
    return np.array(T), np.array(Y)
 
K = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]], float)
A = np.block([[np.zeros((3,3)), np.eye(3)],[-K, np.zeros((3,3))]])
f = lambda t,y: A@y
y0 = np.array([1,0,0, 0,0,0], float)
T, Y = rk4(f, y0, 0, 20, 0.05)
 
from scipy.integrate import solve_ivp
sol = solve_ivp(f, (0,20), y0, method='DOP853', rtol=1e-10, atol=1e-12,
                dense_output=True)

(3 for correct RK4, 1 for solve_ivp call, 2 for diagnostic.) Diagnostic: plot $|E(t)-E(0)|$ (or relative drift) vs $t$. Fixed-step RK4 shows monotone/secular energy drift $O(h^4)$ per step accumulating; DOP853 with tight tol keeps drift near round-off. (2)

(d) (5) $u_1(t)$ is a superposition of the three modes → spectral peaks at $f_j=\omega_j/(2\pi)$: $f_1\approx0.1218$, $f_2\approx0.2251$, $f_3\approx0.2941$ cycles/unit time. (2) Record length $T=20\Rightarrow \Delta f = 1/T = 0.05$ Hz bin spacing. (1) Nyquist: np.fft/rfftfreq show max resolvable freq $=1/(2\Delta t)=1/0.1=10$ Hz. (2)

Question 3 (18)

(a) (4) The model is nonlinear in $\lambda$; the sum-of-squares surface has flat/curved valleys, so LM can stall or diverge from a poor start. Principled guesses: $C_0=\min y$ (or last-time plateau); $A_0=y(t_0)-C_0$; $\lambda_0$ from a log-linear fit: $\ln(y_i-C_0)\approx \ln A_0 -\lambda t_i$, take slope. (4)

(b) (4) pcov is the estimated covariance matrix of popt. Standard error on $\lambda$: sigma_lambda = np.sqrt(np.diag(pcov))[idx_lambda]. It relies on the assumption that residuals are independent, (approximately) Gaussian, and — if absolute_sigma=False — that the model is correct so the residual variance is used to scale $J^\top J)^{-1}$. (4)

(c) (6) Normal equations solve $A^\top A x = A^\top b$. Via SVD $A=U\Sigma V^\top$, $A^\top A = V\Sigma^2 V^\top$, so singular values get squared: $\kappa(A^\top A)=\kappa(A)^2$. (3) Squaring the condition number squares the error amplification; small singular values become negligible, losing accuracy (and $A^\top A$ may become numerically singular). QR/lstsq (SVD-based) work directly on $A$, keeping conditioning at $\kappa(A)$, hence more accurate/backward-stable. (3)

(d) (4) Two-sample Kolmogorov–Smirnov test: scipy.stats.ks_2samp. (Or Mann–Whitney/t-test — KS best for "same distribution".) $H_0$: the two samples are drawn from the same distribution. p $=0.002<\alpha=0.05\Rightarrow$ reject $H_0$: significant evidence the distributions differ. (4)

[
  {"claim":"Vieta small root of x^2-1e8 x+1 is ~1e-8 and stable formula c/q matches",
   "code":"a,b,c=1.0,-1e8,1.0; import math; disc=math.sqrt(b*b-4*a*c); q=-0.5*(b+ (1 if b>0 else -1)*disc); small=c/q; result = abs(small-1e-8) < 1e-12"},
  {"claim":"Eigenvalues of K are 2-