Past exam of the mathematics course of the University of Cambridge / 2016 / iii / Paper 328 (Boundary Value Problems for Linear PDEs)

Put $a=\alpha /2$ and apply the Dirichlet gauge transform for constant drift $w(x,t)=e^{ax}u(x,t)$. Taking derivatives directly removes the advection term and gives the damped heat equation $w_t=w_{xx}-a^{2}w$. Thus $w_0(y)=e^{ay}{u}_0(y)$, while the two Dirichlet boundary conditions become $w(0,t)=g_0(t)$ and $w(L,t)=e^{aL}h(t)$.

The Dirichlet heat kernel on an interval and its damped version areThe method of images gives an alternative, often better at short times:To determine the boundary signs, multiply the equation for $w$ by the backward heat kernel $K_a(x,y,t-s)$ and use integration by parts in $y$. Since the kernel vanishes at $y=0,L$, the surviving boundary expression is $K_{a,y}(x,0,t-s)g_0(s)-K_{a,y}(x,L,t-s)e^{aL}h(s)$. Consequently an integral representation containing only the given data isThis is a Dirichlet boundary-forcing heat-kernel formula. For $0<x<L$ the short-time Gaussian decay makes the forcing integral well defined. The endpoint values are interior limits of the complete formula: evaluating a termwise sine series at an endpoint before taking the time integral loses the nonzero boundary values. At the initial corners a continuous classical solution requires $u_0(0)=g_0(0)$ and $u_0(L)=h(0)$; otherwise the same formula describes the solution away from those corners.

Here is a second integral representation, useful when the integral transforms are explicit. Extend $g_0,h$ by zero after $T$; the values of the extension cannot affect a solution at $t<T$. Write their Laplace transforms as $G(p),H(p)$ and set $q=\sqrt{p+a^2}$ using the principal square root. The transformed Dirichlet Green function for $-{\partial }_x^2+q^2$ isIt vanishes at both endpoints and its first derivative in $x$ jumps by $-1$, which verifies the sign of its unit source. The resolvent kernel for Dirichlet advection-diffusion on an interval gives the transformed solutionThe Bromwich inversion formula therefore yieldsThe square-root notation creates no genuine branch singularity here: all three kernels are even in $q$. Their only spatial resolvent operator poles are $p=-a^2-q_n^2$. This transformed formula and the causal heat kernel formula represent the same solution in the usual smooth-data class.

Use the Bromwich inversion formula in the real variable $\omega$, with $p=\gamma +i\omega$:When the relevant Laplace transforms and spatial integral transforms are explicit, every sample of this integrand can be evaluated directly, without a time-stepping approximation to the partial differential equation. Split the initial condition integral at $y=x$: the two pieces of $R_q$ are linear combinations of $e^{\pm qy}$, so truncated spatial exponential transforms suffice. For real data the negative-frequency half is the complex conjugate of the positive-frequency half.

A practical numerical integration is to truncate to $|\omega |\le \Omega$, apply an adaptive quadrature rule, and independently increase $\Omega$ and refine the quadrature. Choose $\gamma >0$ large enough to remain to the right of all singularities, but avoid an unnecessarily large $e^{\gamma t}$. Evaluate hyperbolic function ratios in scaled form; for example,This prevents overflow when $\mathrm{Re}qL$ is large. The boundary contributions at an interior point are damped by $e^{-x\mathrm{Re}q}$ or $e^{-(L-x)\mathrm{Re}q}$, with $\mathrm{Re}q\sim \sqrt{|\omega |/2}$. Near an endpoint more frequencies are needed; at the endpoint itself use the prescribed Dirichlet boundary condition. The initial condition resolvent has a generally only $O(1/\omega )$ tail. Treat its oscillatory tail accurately, subtract a known transform with the same leading term, or evaluate that contribution separately with the heat kernel.

An equally useful check is the eigenfunction expansion. DefineThe causal modal formula isIt follows either from the heat kernel formula or directly from integration by parts against a sine eigenfunction. Exponential-transform formulas evaluate the time integral explicitly when available. The initial condition contribution has a Gaussian function cutoff in $n$ for each positive time. Nonzero endpoint forcing produces a more slowly convergent sine series tail. A boundary lifting $\ell (x,t)=(1-x/L)g_0(t)+(x/L)e^{aL}h(t)$ gives $w=\ell +z$, where $z$ has zero Dirichlet boundary conditions and forcing $-{\ell }_t-a^2\ell$; expanding $z$ gives better convergence and imposes the endpoint values exactly. At very small times the method of images is often more efficient than retaining many modes.

Check convergence by increasing both the frequency cutoff and quadrature resolution, and compare with a separately truncated causal modal or image-kernel evaluation. Causality in finite-time Laplace contour inversion means a leftward contour deformation must respect both resolvent operator poles and growth of the transformed data. In particular, the transform of data extended by zero after $T$ can grow like $e^{-pT}$ in the left half-plane: for $t<T$ one cannot close that contour indiscriminately and discard its large arc. The modal forcing integral with upper limit $t$ avoids this causality error.

Use the printed coordinate derivative $u_x(0,t)$, not an outward normal derivative; at the left endpoint these have opposite signs. Assume the usual bounded or admissibly growing solution at infinity, with an initial trace $b_0=u_0(0)$ and forcing possessing a Laplace transform. Set $a=\alpha /2$, so $\beta =a^2$, and use $z=\sqrt{p}$ on the principal square root branch. DefineThe half-line Dirichlet Green function for $p-{\partial }_x^2$ givesWrite the transformed solution as $U=U_D+e^{-zx}B(p)$. The Laplace transform of a derivative in the dynamic boundary condition is $pB-b_0$, soThe special choice of $\beta$ gives the perfect squareIn particular, the initial boundary trace cannot be omitted. Withan integral representation isChoose the Bromwich contour to the right of the forcing's growth abscissa and every pole; $\gamma >\max (0,a^2)$ and sufficiently large for the data is a safe choice for real $a$. No sign of $\alpha$ is explicitly imposed in this question.

The inverse transform can also be performed explicitly. For the repeated-root dynamic-boundary heat kernel, setFirst integrate the Heat Poisson kernel against $e^{ar}$, or differentiate the result with respect to $a$:Since differentiation of $(\sqrt{p}-a{)}^{-1}$ with respect to $a$ produces $(\sqrt{p}-a{)}^{-2}$, and ${\partial }_a{E}_a=(-x+2at)E_a+2\sqrt{t/\pi }{e}^{-x^2/(4t)}$, the desired kernel isHere $\mathrm{erfc}$ is the complementary error function. A direct integral characterization, also proving the sign and the transform, isThe convolution theorem for Laplace transforms now gives the fully real time-domain representationAll quantities here are known from the prescribed data. For $x>0$ it tends to $u_0(x)$ as $t\downarrow 0$. At $x=0$, $J_a(0,t)\to 1$ and the spatial correction tends to zero, recovering $b_0$. A solution classical through the initial corner additionally needs $u_0^{\prime \prime }(0)+\alpha u_0^{\prime }(0)+a^2{b}_0=f(0)$; weaker corner regularity does not invalidate the formula for positive times.

The repeated-root unstable heat boundary mode gives a useful sign check on the dynamic boundary condition for the heat equation. If $a>0$, $(\sqrt{p}-a{)}^2$ has a genuine double pole at $p=a^2$ on the physical branch, and the exact homogeneous mode $e^{-ax+a^{2}t}$ satisfies both the heat equation and the printed boundary condition. Generic data can also excite a $t{e}^{a^{2}t}$ contribution. A contour deformation must retain this repeated-pole contribution. If $a<0$, the putative root $\sqrt{p}=a$ is outside the chosen branch and is not a physical pole. At $a=0$, $J_0=\mathrm{erfc}(x/(2\sqrt{t}))$, as expected when the boundary trace satisfies $b^{\prime }(t)=f(t)$.

The formal adjoint of $\mathcal{L}=\Delta -k^2$ is itself. With the bilinear convention used in Green second identity, take any complex-valued $v$ satisfying $\mathcal{L}v=0$. The product rule for divergence gives the familyEquivalently, the differential one-formis a closed differential form. Arbitrary linear combinations of adjoint solutions give further members of this family. One can also write it as $\nabla \cdot [v^2\nabla (u/v)]=0$ wherever $v\ne 0$, but the first formula has no division and remains valid at zeros of $v$. For a sesquilinear convention replace $v$ by $\overline{v}$; no such conjugation is built into the bilinear identity used below.

Write $z=x+iy$ and $\overline{z}=x-iy$. A convenient spectral parameter for a linear boundary value problem is any $\lambda \in \mathbb{C}\setminus \{0\}$. Define the modified Helmholtz adjoint plane wavewhereIndeed $A^2+B^2=k^2$, so $(\Delta -k^2)v_{\lambda }=0$. Equivalently, $4{\partial }_z{\partial }_{\overline{z}}{v}_{\lambda }=k^2{v}_{\lambda }$. The exclusion of zero is only the coordinate singularity of this spectral parametrization; it still supplies a whole complex one-parameter family. The conjugate companion isIt also solves the adjoint equation.

Let $\Omega$ be a bounded piecewise smooth domain, $f=u{|}_{\partial \Omega }$ and $q={\partial }_{n}u$, with $n$ the outward normal vector. Integrating the divergence form from part (i) and using the divergence theorem gives two global relations for a linear boundary value problem:These conjugate global relations for the modified Helmholtz equation are identities between the prescribed Dirichlet boundary condition and its unknown normal derivative; neither is a pointwise boundary equation. Both hold for every nonzero boundary spectral parameter under the usual trace regularity assumptions.

If $u$ is real, its two boundary traces $f,q$ are real. Since ${\tilde{v}}_{\lambda }=\overline{v_{\overline{\lambda }}}$ and the outward normal vector is real,Thus apply complex conjugation to the first global relation for a linear boundary value problem, and replace its parameter by $\overline{\lambda }$, to get the second. For complex $u$, conjugation would also replace its traces by their conjugates, so that argument requires the stated reality assumption. There is also an algebraic redundancy in this particular parametrization: ${\tilde{v}}_{\lambda }=v_{1/\lambda }$, hence $\tilde{\mathcal{G}}(\lambda )=\mathcal{G}(1/\lambda )$ even for complex data. This does not change the requested conjugation identity.

Label the sides by $L,R,B,T$, and parametrize the two vertical sides by $s=y\in [-1,1]$, the horizontal sides by $s=x\in [-1,1]$. Their outward normal derivatives are $q_L=-u_x(-1,s)$, $q_R=u_x(1,s)$, $q_B=-u_y(s,-1)$ and $q_T=u_y(s,1)$. Let $f_j$ denote the prescribed side values. Assume compatible, sufficiently regular boundary traces; the square's corners have zero arclength measure and do not require separate normal values.

For sine collocation of square modified Helmholtz global relations, use Legendre polynomials for the known Dirichlet boundary condition and a Fourier sine series for the unknown normal derivative:The known coefficients are $d_{j\ell }=(2\ell +1){\int }_{-1}^1{f}_j{P}_{\ell }ds/2$. The sine functions have ${\int }_{-1}^1{\varphi }_m{\varphi }_{n}ds={\delta }_{mn}$ and form a complete basis in $L^2(-1,1)$. Choosing them for the normal derivative does not impose zero flux at a corner: the expansion is an $L^2$ representation, and endpoint values are not determined by it. If preserving corner values of the approximated Dirichlet trace is necessary, subtract its endpoint-interpolating line before polynomial approximation, and add that line back. All subsequent known-data integrals can alternatively be evaluated with the exact $f_j$.

Define the entire function basis transformsAt $\tau =\pm i{p}_m$ the quotient is evaluated by its removable limit or by the defining integral. For example, ${\mathsf{P}}_0=2\sinh \tau /\tau$, and polynomial expansion of $P_{\ell }$ expresses every ${\mathsf{P}}_{\ell }$ in derivatives of ${\mathsf{P}}_0$. PutFor compactness write $A=A(\lambda )$ and $B=B(\lambda )$. Substituting these expansions into the first global relation for a linear boundary value problem givesThe signs are fixed by the outward normal vectors, rather than by a choice of traversal direction. The second approximate global relation for a linear boundary value problem replaces $A$ by $-A$, leaving $B$ unchanged:Here $\simeq 0$ records omitted boundary-expansion tails. In a finite spectral method the selected equations are set equal to zero to solve for the $4N$ unknown real coefficients $c_{jm}$. For real data the two families obey the same complex conjugation relation as their exact counterparts.

A particularly simple constructive choice uses paired collocation points for a global relation. For $n=1,\dots ,N$, putFour suitable sets, one for each side, areThe braces in each set range over all $n$. For large $p_n/k$, compute $r_n^{-1}=k/({\omega }_n+p_n)$ to avoid cancellation. The corresponding $(A,B)$ pairs are, respectively,Each pair has the same exponentially growing normal direction and opposite tangential frequencies. Take the phase-weighted difference of the two collocated global relations for a linear boundary value problem. For example, the bottom pair gives the adjoint testThe other three phase-weighted differences giveThus the four spectral sets supply four real sine-tested equations per mode. We retain these $4N$ equations; requiring every uncombined complex equation as well would instead form an overdetermined finite approximation. For real traces each retained equation is the appropriate real combination of its conjugate pair.

The key advantage is that $v_{Bn},v_{Tn}$ vanish on the vertical sides, while $v_{Ln},v_{Rn}$ vanish on the horizontal sides. Orthogonality then eliminates all off-mode unknowns. Define the known right-hand sidesThe retained global relations for a linear boundary value problem areFor example, with $f_{jn}={\int }_{-1}^1{f}_j^{(M)}(s){\varphi }_n(s)ds$, the bottom and top right-hand sides are explicitlyFor completeness the left and right right-hand sides areAll these integrals use known polynomial approximations or the exact Dirichlet boundary condition.

Scale each row by $e^{-{\omega }_n}$. Every block becomesThere is just one off-diagonal entry in each row. This proves diagonal dominance of paired square global-relation collocation. Therefore the assembled $4N\times 4N$ system is strictly diagonally dominant, with row margin $1-{\delta }_n>1-e^{-\pi }$. This proof applies to the explicitly paired and scaled collocation system, rather than assuming arbitrary raw spectral rows are diagonally dominant. Each block is a positive-definite matrix and have spectral condition number of a positive-definite matrixFor instance, setting ${\rho }_{jn}=e^{-{\omega }_n}{R}_{jn}$ gives the concise answerwith the identical formula for the opposite vertical sides. If the known-data integrals use exact $f_j$, these are exactly the first $N$ sine coefficients of the true normal traces: no discarded unknown mode contributes to a retained row. With a finite Legendre polynomial approximation, the error is solely that of the known-data approximation and the final normal-trace truncation. These square modified Helmholtz Dirichlet-to-Neumann coefficients determine the unknown boundary normal derivative; Green's third identity or the separated-variable construction can then recover the interior solution.