Past exam of the mathematics course of the University of Cambridge / 2014 / iii / Paper 69 (The Unified Method for Partial Differential Equations and Medical Imaging)

Use the spectral parameter for a linear boundary value problem $k\in \mathbb{C}$ and define the dispersion relation $\omega (k)=k^2-i\alpha k$. Direct differentiation gives the divergence formIndeed the coefficient of $u$ left over after expansion is $\omega -k^2+i\alpha k=0$, and the remaining factor is $u_t-u_{xx}-\alpha u_x$. Thus this is equivalent to the advection-diffusion equation for every $k$.

Introduce the Half-range Fourier transforms and finite-time spectral boundary transformswhere $g_1(t)=u_x(0,t)$ is the unknown normal derivative with the positive-$x$ convention. The outward normal at zero instead gives $-g_1$. The spatial transforms are analytic for $\mathrm{Im}k<0$ and continuous on the real axis under the stated decay assumptions. Integrating the divergence form on $(0,\infty )\times (0,t)$ gives the global relationThe sign follows from the lower spatial endpoint: the integrated spatial derivative is minus its value at zero. This sign will determine the boundary-forcing term in the solution.

Apply Fourier inversion to the global relation on the real axis, extending $u$ by zero to negative $x$. At an interior point $x>0$ this givesThe finite-time spectral boundary transforms are entire functions of $k$. Define the upper decay domainOrient $\partial D^{+}$ from its left infinite end through $i\alpha$ to its right infinite end, so that $D^{+}$ lies to the left. Its finite point is $i\alpha$, not zero. On the contour, $\mathrm{Re}\omega =0$; the factor $e^{ikx}$ decays because $x>0$.

For the boundary term, write $e^{-\omega t}{G}_j={\int }_0^t{e}^{-\omega (t-s)}{g}_j(s)ds$. Between the real axis and $\partial D^{+}$, $\mathrm{Re}\omega \ge 0$, so this factor has no exponential growth. Contour deformation and Jordan lemma, with the usual cutoffs for oscillatory integrals, therefore giveThis is the requested complex-plane representation with the unknown finite-time spectral boundary transform still present. The decay domain must agree with the sign of the drift in the dispersion relation.

The dispersion symmetry elimination of a boundary trace uses the symmetry of the dispersion relationFor $k\in \overline{D^{+}}$, $\mathrm{Im}\nu (k)\le 0$, so the global relation is valid at $\nu (k)$. Since $i\nu +\alpha =-ik$, it givesSubstitution into the contour integral representation produces an unwanted integral ${\int }_{\partial D^{+}}{e}^{ikx}\hat{u}(\nu (k),t)dk$. Its integrand is analytic in $D^{+}$ and decays on closing the contour upwards; Jordan lemma makes this integral zero for $x>0$. Hence the Fokas method eliminates the unknown normal derivative:All quantities here are determined by the prescribed initial and Dirichlet boundary data up to time $t$.

For verification and for numerical evaluation it is useful to evaluate the spectral contour integrals, giving a half-line drift reflection kernel. Putand define the half-line drift boundary kernelFubini's theorem, the Gaussian Fourier transform and contour deformation give the equivalent causal formulaFor the reflected initial term, ${\hat{u}}_0(i\alpha -k)={\int }_0^{\infty }{e}^{\alpha y+iky}{u}_0(y)dy$ on the contour; the evaluated Gaussian supplies the necessary large-$y$ decay. If $e^{\alpha y}{u}_0(y)$ is not integrable, truncate the initial conditions first, evaluate, and pass to the limit using Gaussian bounds. No extra exponential-decay assumption on the original data is needed for this kernel formula.

The boundary kernel follows particularly simply fromThis calculation independently checks both the sign and the coefficient of the boundary forcing.

The partial differential equation. Each term in $K_{\alpha }$ is a translated heat kernel and satisfies $({\partial }_t-{\partial }_x^2-\alpha {\partial }_x)K_{\alpha }=0$ for $t>0$. The same holds for $P_{\alpha }$, either by direct differentiation or because it is $-(2{\partial }_x+\alpha )H(x+\alpha t,t)$. For $x>0$, $P_{\alpha }(x,t)$ and its derivatives vanish faster than any power as $t\downarrow 0$. Therefore differentiating the boundary convolution creates no extra upper-endpoint term. Gaussian domination justifies differentiating both integrals on compact subsets of $x,t>0$, proving the advection-diffusion equation.

The initial condition. For fixed $x>0$, the first Gaussian in $K_{\alpha }$ is an approximate identity centered at $y=x+\alpha t\to x$. Its integral tends to $u_0(x)$. The reflected Gaussian is exponentially small as $t\downarrow 0$, because its center lies outside the half-line. The boundary convolution also tends to zero for $x>0$. Thus $u(x,t)\to u_0(x)$.

The Dirichlet boundary condition. The identityshows that $K_{\alpha }(0,y,t)=0$, so the initial-data contribution vanishes at zero. The boundary contribution must be evaluated as a limit, not by substituting $x=0$ inside its singular integral. The positive kernel satisfiesfor every $ϵ>0$. Hence it is a one-sided approximate identity at zero time, and continuity of $g_0$ givesThe mass formula follows from the Gaussian Laplace integral, or from the decaying solution of the corresponding constant-coefficient ordinary differential equation. Compatibility $u_0(0)=g_0(0)$ gives the continuous corner value. These arguments also verify the equivalent contour solution in (iii). In the decaying energy class the solution is unique: the difference $v$ of two solutions has zero data and $\frac{d}{dt}{\int }_0^{\infty }{v}^{2}dx=-2{\int }_0^{\infty }{v}_x^{2}dx$ for real solutions, with the analogous modulus identity for complex solutions.

The unheaded sine transform question. The direct classical Fourier sine transform does not close on $u$. Ifthen integration by parts givesThe drift introduces an unknown cosine transform; the usual scalar sine-transform solution of the heat equation is therefore unavailable directly.

A Dirichlet gauge transform for constant drift does provide a qualified alternative. Set $v=e^{\alpha x/2+{\alpha }^{2}t/4}u$. Then $v_t=v_{xx}$, with $v_0=e^{\alpha x/2}{u}_0$ and $v(0,t)=e^{{\alpha }^{2}t/4}{g}_0(t)$. If these weighted data have the decay needed for an ordinary Fourier sine transform, it solves the transformed problem and produces exactly the kernel above. Mere decay of $u_0$ does not ensure this weighted integrability. Thus not directly by the classical sine transform of $u$; yes after a gauge transformation when the required weighted-transform hypotheses hold, or after a justified cutoff/limit argument.

Use the Wirtinger derivatives ${\partial }_z=({\partial }_x-i{\partial }_y)/2$ and ${\partial }_{\bar{z}}=({\partial }_x+i{\partial }_y)/2$, and area measure $dA=dxdy=(i/2)dz\wedge d\bar{z}$. The supplied boundary-integral identity is the planar Generalized Stokes theorem; the usual Poincare lemma is a different local exactness result.

Let $f\in C^1(\overline{D})$ and $z\in D$. Remove a disk of radius $ϵ$ around $z$, and apply Generalized Stokes theorem to the one-form $f(\zeta )d\zeta /(\zeta -z)$ on the punctured domain. Away from the puncture,The outer boundary is counterclockwise and the small inner circle clockwise. Its counterclockwise integral tends to $2\pi if(z)$. Passing to the limit gives the Cauchy-Pompeiu formulaThe weak $1/|\zeta -z|$ singularity is locally integrable. If the boundary term vanishes on expanding $D$ to the whole plane, the formula becomes $f(z)={\pi }^{-1}{\int }_{\mathbb{C}}{f}_{\bar{\zeta }}(\zeta )/(z-\zeta )dA(\zeta )$. In particular it yields the distributional normalization ${\partial }_{\bar{z}}[1/(\pi z)]={\delta }_0$. For holomorphic functions the area term vanishes and one recovers the Cauchy integral formula.

Work first with $q$ in the Schwartz space, so that all Fourier manipulations and spectral contour integrals are justified; weaker classes follow by the usual density or distribution arguments. DefineThe phase is purely imaginary. Since ${\partial }_{\bar{z}}E=kE$, setting $F=E\Phi$ reduces the spectral equation to ${\partial }_{\bar{z}}\Phi =E^{-1}q$. The whole-plane Cauchy-Pompeiu formula therefore constructs the solution decaying spatially at infinity:The freedom to add $E$ times an entire function is removed by this decay condition. For every fixed $k$, the integral is a spatial Cauchy-Green operator applied to a modulated source.

Now differentiate in the conjugate spectral parameter. The two exponential derivatives produce $\zeta -z$, canceling the Cauchy denominator, soThis is the spectral dbar equation: its right-hand side is the forward transform of $q$ multiplied by a known plane wave. Apply the whole-plane Cauchy-Pompeiu formula again, now in $k$:The spatial spectral equation also gives $F(z,k)=-q(z)/k+o(k^{-1})$ as $|k|\to \infty$. One way to justify this is to integrate by parts in the first Cauchy integral: $F=-q/k+T_k(q_{\bar{z}})/k$, where the modulated Cauchy integral $T_k(q_{\bar{z}})$ tends to zero by the Riemann-Lebesgue lemma. Comparing the $1/k$ coefficient of the spectral contour integral therefore givesThis derives the transform pair from two uses of the Cauchy-Pompeiu formula, not from an assumed inversion formula.

To identify the usual normalization, write $z=x+iy$ and $k=k_1+i{k}_2$. Then $\bar{k}z-k\bar{z}=2i(k_{1}y-k_{2}x)$. Set ${\xi }_1=2k_2$, ${\xi }_2=-2k_1$, so $d{\xi }_{1}d{\xi }_2=4dA(k)$. The result is exactly the two-dimensional Fourier transform pairThe factor four in the real-frequency change of variables is essential.

Assume for this spectral construction that the source and known attenuation are sufficiently smooth and decaying, for instance compactly supported, and that $\mu \ge 0$ for the physical interpretation. The following steps identify both the forward attenuated Radon transform and the route to its inversion.

First write the complex transport operator asFor $k=e^{i\theta }$ this is the real directional derivative ${\mathbf{d}}_{\theta }\cdot \nabla$, with ${\mathbf{d}}_{\theta }=(\cos \theta ,\sin \theta )$. Put ${\mathbf{n}}_{\theta }=(-\sin \theta ,\cos \theta )$ and $x=p{\mathbf{n}}_{\theta }+s{\mathbf{d}}_{\theta }$. The spectral equation on the unit circle becomesIts minus sign fixes the appropriate endpoint condition: use $F(p,+\infty )=0$. The integrating factor givesThus the measured quantity at the opposite end of the line isThe exponent is the attenuation accumulated between the source point and the detector at the negative end. The common convention with detector at the positive end is the same transform after reversing the direction, ${\mathcal{R}}_{\mu }^{+}f(p,\theta )={\mathcal{R}}_{\mu }^{-}f(-p,\theta +\pi )$. Using an incoming zero condition at the negative end while retaining the printed minus sign would instead give a growing integrating factor, not physical attenuation.

Next, for $|k|\ne 1$, the operator $L_k$ is a complex elliptic first-order operator. Its decaying whole-plane Green function isThis is obtained by a real-linear change of variables in the Cauchy-Green operator; its change of orientation explains the sign. Write $G_{k}a$ for convolution with this kernel. Solve $L_{k}h=\mu$, and set $F=e^h\Psi$. This removes the attenuation:The normalized spectral solution is analytic separately inside and outside the unit circle. As $k$ approaches that circle, the Green function's characteristic singularity produces two limiting values. Their relation is computed from the weighted line integrals above together with transverse Hilbert transforms; the known attenuation determines the integrating factor weights. The two spectral limits are not individually just the incoming and outgoing real characteristic solutions: the singular-kernel prescription matters.

Finally formulate the resulting additive Riemann-Hilbert problem on the unit circle. Its jump is determined by the measured attenuated Radon transform and known $\mu$. With the unit circle oriented counterclockwise and jump $J=F_{\mathrm{inside}}-F_{\mathrm{outside}}$, a Cauchy integral formula reconstructs the normalized spectral solution:The normalization at zero supplies the compatibility condition $\int J/\zeta d\zeta =0$. Recover $f$ from $L_{k}F-\mu F$, or from its small-$k$ coefficient: if $F=k{F}_1+O(k^2)$, then $f={\partial }_{\bar{z}}{F}_1$. Equivalently the large-$k$ coefficient gives $f={\partial }_z{F}_{-1}$ when $F=k^{-1}{F}_{-1}+O(k^{-2})$. This spectral reconstruction of an attenuated Radon transform is the analogue of recovering $q$ from a coefficient in (ii). Known attenuation and full directed line data are inputs; one does not determine an arbitrary unknown attenuation and source simultaneously from this argument.

Put $E_{\lambda }=e^{-i\beta (\lambda z-\bar{z}/\lambda )}$. The Wirtinger derivatives give $(E_{\lambda }{)}_z=-i\beta \lambda E_{\lambda }$ and $(E_{\lambda }{)}_{\bar{z}}=i\beta E_{\lambda }/\lambda$. Write $W=Adz+Bd\bar{z}$, whereUsing the exterior derivative, $dW=(B_z-A_{\bar{z}})dz\wedge d\bar{z}$. The terms involving $u_z$ and $u_{\bar{z}}$ cancel, leavingBecause $E_{\lambda }$ never vanishes for $\lambda \ne 0$, $dW=0$ if and only if $u_{z\bar{z}}-{\beta }^{2}u=0$. In real coordinates this is the modified Helmholtz equation with mass parameter $2\beta$, not $\beta$. In fact any one nonzero spectral parameter suffices for the equivalence; the full family supplies many independent boundary tests.

Use a consistent number $n$ of vertices, with $z_{n+1}=z_1$, and set $m_j=(z_j+z_{j+1})/2$, $h_j=(z_{j+1}-z_j)/2$. Pull back the one-form to $z_j(s)=m_j+s{h}_j$. Since $dz=h_{j}ds$ and $d\bar{z}={\bar{h}}_{j}ds$, the side integrand isFor counterclockwise traversal, put ${\ell }_j=|h_j|$ and let $q_j={\partial }_{n}u$ be the outward normal derivative, while $g_j=u{|}_{S_j}$ is the Dirichlet boundary data. The outward unit normal is $-i{h}_j/{\ell }_j$ in complex notation. Therefore $h_j{u}_z-{\bar{h}}_j{u}_{\bar{z}}=i{\ell }_j{q}_j$, andThere is no tangential-derivative term: it cancels in this particular one-form. The Generalized Stokes theorem and $dW=0$ now give the polygonal global relationThe same zero identity holds with every side traversed clockwise, but then $h_j{u}_z-{\bar{h}}_j{u}_{\bar{z}}=-i{\ell }_j{q}_j$ for outward normals. One must change this sign consistently rather than mix the two orientations.

Interpret the printed coordinate notation as the square corners $1+i$, $1-i$, $-1-i$, $-1+i$. This listed order is clockwise, contrary to the counterclockwise convention in (ii). Keep the printed first side directed from top to bottom. Then $m_1=1$, $h_1=-i$, and $z(s)=1-is$. Put $g_1(s)=u(1,-s)$ and $q_1(s)=u_x(1,-s)$, where $q_1$ is the outward normal derivative on the right side. The pullback formula, without any orientation shortcut, givesIndeed $h_1{u}_z-{\bar{h}}_1{u}_{\bar{z}}=-i(u_z+u_{\bar{z}})=-i{u}_x$, and $i\beta (\lambda h_1+{\bar{h}}_1/\lambda )=\beta (\lambda -{\lambda }^{-1})$. If one reverses the side to match the counterclockwise convention, its parameter is $1+is$ and its integrand is $-W_1(-s,\lambda )$. Reversing every side multiplies the global relation by minus one, leaving its zero value unchanged.

The unheaded numerical reconstruction request. The four unknown Neumann boundary conditions are four functions on the sides, not four scalar values. The polygonal global relation is linear in their outward normal derivatives, and its remaining terms depend only on the prescribed Dirichlet boundary data.

Choose a finite approximation on each side, for example an expansion of $q_j(s)$ in Legendre polynomials or piecewise polynomials. Substitute those expansions into the consistently oriented global relation. At chosen nonzero spectral parameters for a linear boundary value problem, integrate the exponential kernels against each basis function to assemble a complex linear system; the known right-hand side is obtained by integrating the given Dirichlet boundary data. Use enough independent samples to resolve all side coefficients, and preferably oversample. The conjugate global relations for the modified Helmholtz equation provide useful companion tests; for complex data use the corresponding independent adjoint relation rather than assuming the traces real.

Solve the scaled system by least-squares solution using a stable factorization such as a singular value decomposition. Sampling directions should probe all sides, and exponential row scaling avoids overflow and poor conditioning. Refine the side approximation and spectral samples until the recovered traces and unused global relation residuals stabilize. Corner incompatibilities or limited corner regularity call for mesh refinement or enriched endpoint basis functions. This realizes a numerical Dirichlet-to-Neumann map without first discretizing the whole interior.

The underlying Dirichlet problem is uniquely solvable in the usual trace class for $\Delta u-4{\beta }^{2}u=0$ with $\beta >0$: the homogeneous problem has ${\int }_{\Omega }(|\nabla u{|}^2+4{\beta }^2|u{|}^2)dA=0$. This supports the boundary reconstruction, although uniqueness of the continuous problem alone does not guarantee that an arbitrary finite set of spectral samples is well conditioned.