Past exam of the short options paper shared by Parts A and B of the physics course of the University of Oxford / 2019 / S31

Round-off error arises from representing numbers and carrying out arithmetic at finite precision: storing $1/3$ or subtracting nearly equal numbers can lose information even when the mathematical algorithm is exact. Truncation error arises from replacing a continuous or infinite mathematical operation by a finite approximation, such as discarding higher Taylor series terms in a finite difference. Refining a grid usually reduces truncation error but may increase accumulated round-off error; the two have different origins.

Write $h=\Delta x$, $k=\Delta t$ and ${\delta }_{xx}{u}_i=(u_{i+1}-2u_i+u_{i-1})/h^2$. Insert a smooth exact heat equation solution into the Crank-Nicolson method. Expanding about the midpoint $(x_i,t_n+k/2)$ givesThus the residual, defined as left side minus right side, isSince $u_t=D{u}_{xx}$ implies $u_{ttt}=D{u}_{xxtt}$,This proves second-order numerical consistency in both space and time. About the old time level, uncombined Taylor remainders can include $O(k^3)$; the midpoint expression displays the symmetric leading errors most cleanly. Multiplying by $k$ gives the one-step defect $O(k^3+k{h}^2)$. Stable evolution with suitable initial and boundary data then gives global error $O(k^2+h^2)$; consistency alone does not prove global convergence.

For balanced global error $\epsilon$, a method second order in both variables can take $h\sim {\epsilon }^{1/2}$ and $k\sim {\epsilon }^{1/2}$. With $O(1/h)$ work per time step in one spatial dimension, its cost scales as ${\epsilon }^{-1}$. A method second order in space but first order in time instead needs $k\sim \epsilon$, so its cost is $\boxed{O({\epsilon }^{-3/2})}$ under the same per-step cost assumption. Equivalently at fixed $h$, using $k\sim h^2$ rather than $k\sim h$ needs about $1/h$ more steps. A restrictive stability condition or different linear-solver cost can dominate these accuracy-based estimates.

Use a smooth exact solution such as $u(x,t)=e^{-D{q}^{2}t}\sin (qx)$ on a domain with matching boundary conditions, and compare numerical and exact values at a fixed physical time. To measure spatial order, halve $h$ while making temporal error negligible; the error norm should fall by a factor of four. To measure temporal order, halve $k$ on a sufficiently fine fixed spatial grid and again expect a factor of four. Report observed order $p={\log }_2(E_h/E_{h/2})$ and use several refinements to identify the asymptotic regime.

Check more than one smooth mode, boundary implementation and expected conservation or decay properties. A manufactured solution with an added known source can test nontrivial boundaries and forcing. Repeat in maximum and root-mean-square norms, stop before finite-precision or solver tolerances dominate, and separate stability failures from errors in convergence order. A single small error or a single visually smooth solution does not establish accuracy.

Let the step be $h$. Forward Euler method gives $x_1=x_0+h{v}_0$, $v_1=v_0+hF(x_0)$. Applying the same map with $-h$ gives $x_2=x_1-h{v}_1=x_0-h^{2}F(x_0)$, which generally differs from $x_0$. It is therefore not its own inverse under time reversal.

A kick–drift–kick leapfrog integrator instead hasThe step $-h$ first recovers $v_{1/2}$, then $x_0$, then $v_0$, so $\boxed{{\Phi }_{-h}={\Phi }_h^{-1}}$. This assumes a position-dependent force and the same fixed step. Newtonian mechanics is reversible; respecting that structure helps avoid artificial directional bias and often improves long-term qualitative behaviour. The usual leapfrog also has a structure preserving symplectic geometry, which supports bounded oscillatory energy error, but reversibility by itself does not imply exact conservation of energy. Adaptive stepping requires special care to preserve these properties.

For a second-order scalar partial differential equation, classification refers to its principal second-derivative form. An elliptic partial differential equation has a definite principal form, exemplified by Poisson's equation; it typically determines a spatial field from boundary data. A parabolic partial differential equation has a distinguished first-order evolution direction and diffusive second spatial derivatives, exemplified by $u_t=D{u}_{xx}$. Its initial data relax or spread with time. Thus the pseudo-time equation is parabolic, while its steady-state spatial equation is elliptic.

At a steady state in pseudo-time, ${\partial }_{\tau }\Phi =0$ and the evolution becomes ${\Phi }_{xx}=4\pi G\rho$, the desired Poisson equation. Deviations from a steady solution obey a diffusion equation, so nonconstant error modes decay. This turns an elliptic boundary-value solve into repeated relaxation steps using familiar time-evolution methods. Boundary data and source compatibility still matter; pseudo-time does not remove a nonexistent steady solution.

A backward Euler method with spatial second-order central differences givesThe source is fixed in pseudo-time; if it were evolving, its evaluation would also need specifying. Each pseudo-time step therefore requires a linear system solve rather than a purely local explicit update.

Let $r=\Delta \tau /h^2$ and impose periodic indices modulo $J$. With ${\mathbf{\Phi }}^n=({\Phi }_0^n,\dots ,{\Phi }_{J-1}^n{)}^T$, the linear system iswhere $S$ shifts one grid cell cyclically. For $J\ge 3$, $A$ has diagonal $1+2r$, adjacent and corner entries $-r$, and all other entries zero. It is a cyclic tridiagonal matrix, rather than a nonperiodic tridiagonal one.

The periodic Laplacian annihilates the constant vector, so its steady equation requires $\boxed{\sum _j{\rho }_j=0}$. A strictly positive gravitational density violates this condition: summing the update gives ${\bar{\Phi }}^{n+1}={\bar{\Phi }}^n-4\pi G\Delta \tau \bar{\rho }$, an indefinitely drifting mean. Periodic gravitational calculations normally solve for $\rho -\bar{\rho }$ and fix the arbitrary mean gravitational potential. The per-step matrix $A$ itself remains invertible; the issue is existence of a periodic steady solution of Poisson's equation.

For an error Fourier mode $e_j^n={\hat{e}}^n{e}^{i\theta j}$, the periodic discrete Laplacian has eigenvalue $-4{\sin }^2(\theta /2)/h^2$. The Backward Euler method therefore has amplification factorvon Neumann stability analysis thus imposes no upper limit on positive $\Delta \tau$. The constant mode has $g=1$ and is not damped; compatible forcing and a gauge condition remain necessary as discussed in part (e). Large steps can accelerate relaxation of nonconstant modes, but require a suitably accurate linear solver and can worsen its conditioning. Unconditional linear stability is not a statement about time accuracy or cost.

Put $y^{(0)}=y$ and $y^{(1)}=\dot{y}$. ThenThe second row is the original damped harmonic oscillator equation and the first row defines the auxiliary velocity. This converts a second-order ordinary differential equation into two first-order equations without changing the required two initial values.

A numerical evolution has numerical stability if small perturbations in initial data, forcing or arithmetic are not amplified without control over the time interval of interest. For this constant-coefficient decaying linear system, an absolute-stability test requires the homogeneous numerical amplification factors to have magnitude at most one, with no growing Jordan-block contribution on the unit circle. Stability is distinct from accuracy: a strongly overdamped numerical solution can be stable while representing the physical oscillation poorly.

The characteristic polynomial is ${\lambda }^2+2\gamma \lambda +k^2=0$. Since $0<\gamma <|k|$,Forward Euler method has amplification eigenvalues $g_{\pm }=1+h{\lambda }_{\pm }$. Their squared modulus isThusStrict inequality gives numerical damping; at the upper boundary the amplitude is constant despite physical decay, so that boundary is a poor accuracy choice. In the undamped limit $\gamma =0$, $|g{|}^2=1+k^2{h}^2>1$ for every $h>0$: forward Euler method is unstable for any nonzero step on an undamped harmonic oscillator.

Backward Euler method has update ${\mathbf{y}}^{n+1}=(I-hM{)}^{-1}{\mathbf{y}}^n$ and amplification eigenvalues $g_{\pm }=(1-h{\lambda }_{\pm }{)}^{-1}$. Hencefor every $h\ge 0$. It is therefore unconditionally stable for this linear system. Even with $\gamma =0$, it damps every nonzero numerical mode: stability is obtained at the cost of artificial damping and phase error. Large $h$ does not accurately resolve the oscillation period.

An explicit numerical method normally has a simple, inexpensive step and needs no simultaneous solve, but stability can force very small steps, especially for stiff differential equations. An implicit numerical method requires a linear or nonlinear solve, storage and suitable solver tolerances; its larger stability region can permit larger steps when stability is the main restriction. Accuracy, phase fidelity and preservation of physical structure still constrain the step. For this oscillator, backward Euler method is stable but dissipative, while forward Euler method is cheap but unstable without damping. Neither statement makes all implicit methods superior to all explicit ones; method order, problem structure and total cost determine the useful choice.