Numerical Methods

Question 1 (10)

(a) Taylor about $x_n$: $f(r) = f(x_n) + f'(x_n)(r - x_n) + O((r-x_n)^2)$. Set $f(r)=0$, drop quadratic term, solve for $r \approx x_{n+1}$: $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.$ (1 for expansion, 1 for setting $f(r)=0$ & truncating, 1 for formula.)

(b) Write $x_{n+1} - r = e_{n+1}$. Expand $f(r)$ and $f'$ around $x_n$: $0 = f(x_n) - e_n f'(x_n) + \tfrac{1}{2}e_n^2 f''(x_n) + O(e_n^3).$ Since $x_{n+1} = x_n - f(x_n)/f'(x_n)$, $e_{n+1} = e_n - \frac{f(x_n)}{f'(x_n)}.$ Substitute $f(x_n) = e_n f'(x_n) - \tfrac12 e_n^2 f''(x_n) + \dots$ (from expanding $f(r)=0$ rearranged) to get $e_{n+1} = \frac{f''(x_n)}{2 f'(x_n)} e_n^2 + O(e_n^3) \;\Rightarrow\; e_{n+1} \approx \frac{f''(r)}{2f'(r)} e_n^2.$ Asymptotic constant $C = \left|\dfrac{f''(r)}{2 f'(r)}\right|$; order $p=2$. (2 for expansion setup, 2 for substitution/cancellation of $e_n$ term, 1 for constant & stating quadratic.)

(c) If $f'(r)=0$ (multiple root) quadratic convergence is lost; convergence becomes linear (order 1). (1+1)

Question 2 (10)

(a) Catastrophic cancellation: subtracting two nearly equal floating-point numbers. Each operand already carries relative rounding error $\sim \varepsilon_{mach}$ (absolute error $\sim \varepsilon |x|$). Their difference is small but the absolute errors persist, so the relative error of the difference is magnified by roughly $\frac{|x|}{|x-y|}$, destroying significant digits. (1 def, 2 explanation.)

(b) $x = b - \sqrt{b^2 - c}$: when $b^2 \gg c$, $\sqrt{b^2-c} \approx b$, so we subtract nearly-equal quantities → cancellation. Multiply by conjugate: $x = \frac{(b-\sqrt{b^2-c})(b+\sqrt{b^2-c})}{b+\sqrt{b^2-c}} = \frac{c}{\,b + \sqrt{b^2 - c}\,},$ an addition — no cancellation. (2 identify problem, 2 stable form.)

(c) $\sqrt{b^2-c} \approx b\sqrt{1 - c/b^2} \approx b(1 - c/(2b^2))$. Cancellation removes $\log_{10}(b^2/c)$ digits $\approx \log_{10}(10^{16}/1) = 16$. Essentially ~16 digits lost — the naive result is worthless (true small root $\approx c/(2b) = 5\times10^{-9}$). Acceptable answer: about 15–16 digits. (3.)

Question 3 (12)

(a) With Lagrange quadratic $P_2$ through the three points and substitution $x = x_0 + sh$: $\int_{x_0}^{x_2} P_2\,dx = \frac{h}{3}\left(f_0 + 4f_1 + f_2\right).$ Derivation: integrate the interpolating polynomial (or Newton forward form) over $s\in[0,2]$; the coefficients $\tfrac13,\tfrac43,\tfrac13$ emerge. (3 setup/integration, 2 result.)

(b) Composite ($n$ even, $h=(b-a)/n$): $\int_a^b f\,dx \approx \frac{h}{3}\Big[f_0 + f_n + 4\!\!\sum_{i\text{ odd}}\!f_i + 2\!\!\sum_{i\text{ even}\atop i\neq 0,n}\!f_i\Big].$ Requirement: $n$ must be even. (2 formula, 1 requirement.)

(c) $h=0.25$; nodes $0,0.25,0.5,0.75,1$. $f_0=1,\ f_1=e^{-0.0625}=0.939413,\ f_2=e^{-0.25}=0.778801,\ f_3=e^{-0.5625}=0.569783,\ f_4=e^{-1}=0.367879.$ $I=\frac{0.25}{3}\big[1 + 0.367879 + 4(0.939413+0.569783) + 2(0.778801)\big]$ $=\frac{0.25}{3}[1.367879 + 4(1.509196) + 1.557602] = \frac{0.25}{3}(8.962265) = 0.74685.$ (True value $0.746824$.) $\approx 0.74685$. (2 evaluations, 2 answer.)

Question 4 (10)

(a) Pseudocode:

GaussSeidel(A, b, x0, tol, maxit):
    x = x0
    for k = 1 to maxit:
        x_old = copy(x)
        for i = 1 to n:
            s = b[i]
            for j = 1 to n, j != i:
                s = s - A[i][j] * x[j]   # uses updated x for j<i
            x[i] = s / A[i][i]
        if norm(x - x_old) < tol: return x
    return x   # (or flag non-convergence)

(2 loop structure, 2 in-place update using new values, 1 stopping test.)

(b) Sufficient condition: $A$ strictly diagonally dominant, $|a_{ii}| > \sum_{j\neq i}|a_{ij}|$ (also converges for SPD $A$). Gauss–Seidel uses already-updated components within the same sweep, so it incorporates newer information → generally faster than Jacobi. (2 condition, 1 reason.)

(c) $x_1 = 3/4 = 0.75$. Then $x_2 = (6 + 0.75)/4 = 1.6875$. Result: $(0.75,\ 1.6875)$. (2.) (True solution: $x_1=1.2,\ x_2=1.8$.)

Question 5 (10)

(a) $k_1 = f(t_n,y_n),\quad k_2 = f(t_n+\tfrac h2, y_n+\tfrac h2 k_1),$ $k_3 = f(t_n+\tfrac h2, y_n+\tfrac h2 k_2),\quad k_4 = f(t_n+h, y_n+h k_3),$ $y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4).$ (2 stages, 2 update.)

(b) RK4 local truncation error $O(h^5)$; global error $O(h^4)$. Euler: local $O(h^2)$, global $O(h)$. RK4 is 4th-order accurate vs Euler 1st-order. (1 LTE, 1 global, 1 Euler contrast.)

(c) $f=y$, $h=0.5$, $y_0=1$: $k_1=1$, $k_2=1+0.25(1)=1.25$, $k_3=1+0.25(1.25)=1.3125$, $k_4=1+0.5(1.3125)=1.65625$. $y_1 = 1 + \frac{0.5}{6}(1 + 2(1.25)+2(1.3125)+1.65625) = 1 + \frac{0.5}{6}(7.78125)=1.648438.$ $e^{0.5}=1.648721$; abs error $\approx 0.000283$. $y(0.5)\approx1.64844$, error $\approx2.8\times10^{-4}$. (2 stages, 1 answer+error.)

Question 6 (8)

(a) Relative condition number = (relative change in output)/(relative change in input): $\kappa = \lim_{\delta x\to 0}\left|\frac{(f(x+\delta x)-f(x))/f(x)}{\delta x / x}\right| = \left|\frac{f'(x)\,x}{f(x)}\right|.$ Using $f(x+\delta x)-f(x)\approx f'(x)\delta x$. (2 definition, 2 limit/derivative.)

(b) Ill-conditioned problem: inherent property — small input perturbations cause large output changes regardless of method ($\kappa \gg 1$). Unstable algorithm: the computational procedure introduces/amplifies rounding error beyond what the problem's conditioning warrants. A stable algorithm bounds added error, but the inherent sensitivity $\kappa$ still multiplies the unavoidable input rounding — so no algorithm can produce an accurate result for an ill-conditioned problem. (2 distinction, 2 why stability can't rescue.)

[
  {"claim":"Composite Simpson n=4 for integral e^{-x^2} on [0,1] ≈ 0.74685",
   "code":"h=Rational(1,4);f=lambda x: exp(-x**2);I=h/3*(f(0)+f(1)+4*(f(Rational(1,4))+f(Rational(3,4)))+2*f(Rational(1,2)));result=abs(float(I)-0.74685)<1e-4"},
  {"claim":"RK4 one step h=0.5 for y'=y,y0=1 gives 1.648438",
   "code":"h=0.5;y=1.0;f=lambda t,y:y;k1=f(0,y);k2=f(0,y+h/2*k1);k3=f(0,y+h/2*k2);k4=f(0,y+h*k3);y1=y+h/6*(k1+2*k2+2*k3+k4);result=abs(y1-1.648438)<1e-5"},
  {"claim":"RK4 abs error vs e^0.5 ≈ 2.8e-4",
   "code":"h=0.5;y=1.0;k1=y;k2=y+h/2*k1;k3=y+h/2*k2;k4=y+h*k3;y1=y+h/6*(k1+2*k2+2*k3+k4);result=abs(abs(y1-float(exp(Rational(1,2))))-0.000283)<5e-5"},
  {"claim":"Gauss-Seidel one sweep gives (0.75,1.6875)",
   "code":"x1=(3-0)/4;x2=(6+x1)/4;result=abs(x1-0.75)<1e-12 and abs(x2-1.6875)<1e-12"}
]