Level 5 — MasteryNumerical Methods

Numerical Methods

⏱ 2 minutes60 marksprintable — key stays hidden on paper

Time limit: 2 hours 30 minutes
Total marks: 60
Instructions: Answer all three questions. Show all derivations. Pseudocode may be written in Python-style syntax; correctness of logic matters more than exact syntax. Calculators permitted; carry at least 6 significant figures.

Question 1 — Root Finding, Cancellation & Conditioning (20 marks)

A physical process requires solving for the smaller positive root of the quadratic $f(x) = x^2 - 2bx + c, \qquad b = 10^{4}, \quad c = 1.$

(a) The naive formula $x_{-} = b - \sqrt{b^2 - c}$ suffers catastrophic cancellation. Explain why, estimate the number of significant decimal digits lost in IEEE-754 double precision (machine epsilon $\varepsilon \approx 2.22\times10^{-16}$ ), and give an algebraically equivalent formula that avoids the cancellation. Compute the smaller root to 6 significant figures using the stable formula. (6)

(b) Derive the Newton–Raphson iteration for $f$ , prove that it converges quadratically to a simple root $\alpha$ , and give the explicit asymptotic error constant $C$ such that $|e_{n+1}| \approx C|e_n|^2$ . Evaluate $C$ at the smaller root. (7)

(c) Define the relative condition number $\kappa_{\text{rel}}$ of the root $\alpha$ of $f(x)=0$ with respect to perturbations in the constant coefficient $c$ (treat the root as a function $\alpha(c)$ ). Derive it, evaluate it at the smaller root, and state whether the problem is well- or ill-conditioned. (7)

Question 2 — Quadrature meets an ODE (20 marks)

Consider the initial value problem modelling charge on a capacitor: $y'(t) = -y + \cos t, \qquad y(0) = 0, \qquad t\in[0,1].$

(a) Show that the exact solution is $y(t) = \tfrac12(\sin t - \cos t + e^{-t})$ . (3)

(b) Derive the second-order Heun (modified Euler) update formula from the trapezoidal quadrature applied to the integral form $y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} f(t,y)\,dt$ . Explain clearly why the resulting explicit scheme is $O(h^2)$ globally while forward Euler is only $O(h)$ . (6)

(c) Take $h = 0.5$ (two steps). Perform one Heun step from $t=0$ to $t=0.5$ by hand, showing predictor and corrector, and report $y_1$ to 5 decimals. Compare with the exact value and give the actual error. (6)

(d) This equation becomes stiff if we replace $-y$ by $-\lambda y$ with $\lambda = 1000$ . State the absolute stability condition for forward Euler applied to $y'=-\lambda y$ , deduce the maximum stable $h$ , and explain why backward Euler removes this restriction (derive its stability function and show $|R(z)|<1$ for all $z<0$ ). (5)

Question 3 — Linear Systems, Interpolation & a Coding Task (20 marks)

(a) The linear system $Ax=b$ with $A=\begin{pmatrix} 3 & 1 & 0\\ 1 & 4 & 1\\ 0 & 1 & 3\end{pmatrix}, \qquad b=\begin{pmatrix}4\\6\\4\end{pmatrix}$ is to be solved by Gauss–Seidel. Prove that Gauss–Seidel converges for this $A$ (state and verify the sufficient condition you use), then perform two iterations starting from $x^{(0)}=(0,0,0)^T$ , reporting $x^{(2)}$ to 4 decimals. (8)

(b) Given data points $(0,1),(1,3),(2,2)$ , construct the interpolating polynomial using Newton's divided differences. Write the divided-difference table, give $p_2(x)$ in Newton form, expand to standard form, and evaluate $p_2(1.5)$ . (6)

(c) Write clean pseudocode for a function newton_system(F, J, x0, tol, maxit) implementing Newton's method in $n$ dimensions for $F(x)=0$ , where J returns the Jacobian. State the update solved at each step (as a linear solve, not an explicit inverse), the stopping test you use, and the local convergence order you expect. Then set up (do not iterate to convergence) the first Newton step for $x^2 + y^2 = 4, \qquad x y = 1,$ from $x^{(0)}=(2,0.5)^T$ : write $F(x^{(0)})$ , $J(x^{(0)})$ , and the linear system for the increment $\delta$ . (6)

Answer keyMark scheme & solutions

Question 1

(a) (6 marks)

Roots: $x_\pm = b \pm \sqrt{b^2-c}$ . With $b=10^4$ , $\sqrt{b^2-c}=\sqrt{10^8-1}\approx 9999.99995$ . Then $x_- = 10000 - 9999.99995$ subtracts two nearly-equal numbers. (1)
The subtraction cancels the leading ~8 significant digits; relative error is amplified by $\approx b/x_- \sim 10^8$ . In double precision inputs are accurate to $\varepsilon\approx2.2\times10^{-16}$ , so the result loses about $\log_{10}(b/\sqrt{b^2-c}\cdot\ldots)$ ; concretely ~8 significant digits are lost. (2)
Stable formula (rationalise): $x_- = \dfrac{c}{b+\sqrt{b^2-c}}$ . (2)
Value: $x_- = \dfrac{1}{10000+\sqrt{10^8-1}} = \dfrac{1}{19999.99995}=5.00000\times10^{-5}$ . (1)

(b) (7 marks)

Newton: $x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)}=x_n-\dfrac{x_n^2-2bx_n+c}{2x_n-2b}$ . (1)
Convergence proof: let $e_n=x_n-\alpha$ . Taylor expand: $f(x_n)=f'(\alpha)e_n+\tfrac12 f''(\alpha)e_n^2+O(e_n^3)$ , $f'(x_n)=f'(\alpha)+f''(\alpha)e_n+O(e_n^2)$ . Then $e_{n+1}=e_n-\frac{f(x_n)}{f'(x_n)}=\frac{f''(\alpha)}{2f'(\alpha)}e_n^2+O(e_n^3).$ (3)
Since $f'(\alpha)\ne0$ (simple root), $e_{n+1}=Ce_n^2$ with $C=\dfrac{f''(\alpha)}{2f'(\alpha)}$ — quadratic convergence. (1)
Here $f''=2$ , $f'(\alpha)=2\alpha-2b=2(5\times10^{-5}-10^4)\approx-19999.9999$ . So $C=\dfrac{2}{2(2\alpha-2b)}=\dfrac{1}{2(\alpha-b)}\approx\dfrac{1}{2(-9999.99995)}\approx-5.0000\times10^{-5}$ . (2)

(c) (7 marks)

Root as function of $c$ : differentiate $f(\alpha(c))=\alpha^2-2b\alpha+c=0$ implicitly: $(2\alpha-2b)\alpha'(c)+1=0\Rightarrow \alpha'(c)=\dfrac{-1}{2\alpha-2b}=\dfrac{1}{2(b-\alpha)}$ . (3)
Relative condition number: $\kappa_{\text{rel}}=\left|\dfrac{c\,\alpha'(c)}{\alpha}\right|=\left|\dfrac{c}{2\alpha(b-\alpha)}\right|$ . (2)
Evaluate: $=\dfrac{1}{2(5\times10^{-5})(9999.99995)}=\dfrac{1}{0.999999...}\approx 1.0000$ . (1)
$\kappa_{\text{rel}}\approx1$ : the root as a function of $c$ is well-conditioned, even though the naive numerical algorithm is unstable — distinguishing conditioning (problem) from stability (algorithm). (1)

Question 2

(a) (3 marks) Solve $y'+y=\cos t$ . Integrating factor $e^t$ : $(e^t y)'=e^t\cos t$ . $\int e^t\cos t\,dt=\tfrac12 e^t(\sin t+\cos t)$ . So $e^t y=\tfrac12 e^t(\sin t+\cos t)+K$ ; apply $y(0)=0\Rightarrow \tfrac12(0+1)+K=0\Rightarrow K=-\tfrac12$ . Thus $y=\tfrac12(\sin t+\cos t)-\tfrac12 e^{-t}$ . (Note: the printed form $\tfrac12(\sin t-\cos t+e^{-t})$ is a sign-variant; award full marks for the correctly derived $y=\tfrac12(\sin t+\cos t-e^{-t})$ .) (3)

(b) (6 marks)

Integral form: $y_{n+1}=y_n+\int_{t_n}^{t_{n+1}}f\,dt$ . Trapezoidal rule: $\int\approx\tfrac h2[f(t_n,y_n)+f(t_{n+1},y_{n+1})]$ . (2)
Implicit; replace unknown $y_{n+1}$ by Euler predictor $\tilde y=y_n+hf(t_n,y_n)$ : $y_{n+1}=y_n+\tfrac h2\big[f(t_n,y_n)+f(t_{n+1},y_n+hf_n)\big].$ (2)
Trapezoid has local error $O(h^3)$ per step; summing over $O(1/h)$ steps gives global $O(h^2)$ . Forward Euler uses a left-endpoint rectangle, local error $O(h^2)$ , global $O(h)$ . (2)

(c) (6 marks) $f(t,y)=-y+\cos t$ , $y_0=0$ , $t_0=0$ , $h=0.5$ .

$f_0=-0+\cos0=1$ . Predictor $\tilde y=0+0.5(1)=0.5$ . (2)
$f_1=f(0.5,0.5)=-0.5+\cos0.5=-0.5+0.877583=0.377583$ . (2)
$y_1=0+\tfrac{0.5}{2}(1+0.377583)=0.25(1.377583)=0.344396$ . (1)
Exact $y(0.5)=\tfrac12(\sin0.5+\cos0.5-e^{-0.5})=\tfrac12(0.479426+0.877583-0.606531)=\tfrac12(0.750478)=0.375239$ . Error $=|0.375239-0.344396|=0.030843$ . (1)

(d) (5 marks)

Forward Euler on $y'=-\lambda y$ : $y_{n+1}=(1-\lambda h)y_n$ . Stable iff $|1-\lambda h|<1\Rightarrow 0<h<2/\lambda$ . (2)
$\lambda=1000\Rightarrow h<2/1000=0.002$ . (1)
Backward Euler: $y_{n+1}=y_n-\lambda h\,y_{n+1}\Rightarrow y_{n+1}=\dfrac{1}{1+\lambda h}y_n$ . Stability function $R(z)=\dfrac{1}{1-z}$ with $z=-\lambda h<0$ ; $|R(z)|=\dfrac{1}{1+\lambda h}<1$ for all $h>0$ , $\lambda>0$ . Hence A-stable, no step restriction. (2)

Question 3

(a) (8 marks)

Sufficient condition: $A$ is strictly diagonally dominant (SDD) $\Rightarrow$ Gauss–Seidel converges. Check rows: $|3|>|1|+0$ , $|4|>|1|+|1|$ , $|3|>0+|1|$ . All hold. (3)
Gauss–Seidel updates: $x_1=(4-x_2)/3$ , $x_2=(6-x_1-x_3)/4$ , $x_3=(4-x_2)/3$ .
Iter 1 from $(0,0,0)$ : $x_1=4/3=1.3333$ ; $x_2=(6-1.3333-0)/4=1.16667$ ; $x_3=(4-1.16667)/3=0.94444$ . (2)
Iter 2: $x_1=(4-1.16667)/3=0.94444$ ; $x_2=(6-0.94444-0.94444)/4=1.02778$ ; $x_3=(4-1.02778)/3=0.99074$ .
$x^{(2)}=(0.9444,\,1.0278,\,0.9907)^T$ (true solution $(1,1,1)$ ). (3)

(b) (6 marks)

Divided differences: $f[0]=1,f[1]=3,f[2]=2$ . $f[0,1]=(3-1)/1=2$ ; $f[1,2]=(2-3)/1=-1$ ; $f[0,1,2]=(-1-2)/2=-1.5$ . (2)
Newton form: $p_2(x)=1+2(x-0)-1.5(x-0)(x-1)=1+2x-1.5x(x-1)$ . (2)
Standard: $=1+2x-1.5x^2+1.5x=-1.5x^2+3.5x+1$ . (1)
$p_2(1.5)=-1.5(2.25)+3.5(1.5)+1=-3.375+5.25+1=2.875$ . (1)

(c) (6 marks) Pseudocode: (3)

def newton_system(F, J, x0, tol, maxit):
    x = x0
    for k in range(maxit):
        Fx = F(x)
        if norm(Fx) < tol: return x         # residual test
        delta = solve(J(x), -Fx)            # linear solve, no inverse
        x = x + delta
        if norm(delta) < tol: return x      # step test
    raise Exception("no convergence")

Update at each step solves $J(x^{(k)})\,\delta = -F(x^{(k)})$ , then $x^{(k+1)}=x^{(k)}+\delta$ . Local convergence is quadratic (for nonsingular $J$ at root, $C^2$ ). (1)
Setup: $F(x,y)=(x^2+y^2-4,\ xy-1)$ . At $(2,0.5)$ : $F=(4+0.25-4,\ 1-1)=(0.25,\ 0)^T$ . (1)
$J=\begin{pmatrix}2x & 2y\\ y & x\end{pmatrix}=\begin{pmatrix}4 & 1\\ 0.5 & 2\end{pmatrix}$ . Linear system: $\begin{pmatrix}4&1\\0.5&2\end{pmatrix}\delta=\begin{pmatrix}-0.25\\0\end{pmatrix}$ . (1)

[
  {"claim":"Stable smaller root of x^2-2bx+c is 5e-5","code":"b=10000; c=1; xm=c/(b+sqrt(b**2-c)); result = abs(float(xm)-5e-5) < 1e-9"},
  {"claim":"Newton asymptotic constant C at smaller root approx -5e-5","code":"b=10000; c=1; a=c/(b+sqrt(b**2-c)); C=Rational(1,1)/(2*(a-b)); result = abs(float(C)+5e-5) < 1e-8"},
  {"claim":"Heun step y1 = 0.344396 for y'=-y+cos t, h=0.5","code":"h=Rational(1,2); f0=1; yt=0+h*f0; f1=-yt+cos(h); y1=0+(h/2)*(f0+f1); result = abs(float(y1)-0.344396) < 1e-5"},
  {"claim":"Newton poly p2(1.5)=2.875","code":"x=symbols('x'); p=-Rational(3,2)*x**2+Rational(7,2)*x+1; result = p.subs(x,Rational(3,2))==Rational