Explicit unconditionally stable methods for the heat equation via potential theory

TL;DR

This paper demonstrates that explicit marching schemes for the heat equation, based on potential theory, can be unconditionally stable under certain boundary conditions, contrasting with traditional stability constraints.

Contribution

The authors establish unconditional stability of the forward Euler scheme for specific boundary value problems of the heat equation using spectral radius bounds and spectral analysis of Toeplitz matrices.

Findings

Unconditional stability for Dirichlet and Neumann problems on the unit ball.

Stability of the Robin problem with a timestep bound independent of spatial discretization.

Unconditional stability in the $L^ abla$-norm for Dirichlet problems on smooth convex domains.

Abstract

We study the stability properties of explicit marching schemes for second-kind Volterra integral equations that arise when solving boundary value problems for the heat equation by means of potential theory. It is well known that explicit finite difference or finite element schemes for the heat equation are stable only if the time step $\Delta t$ is of the order $O(\Delta x^2)$, where $\Delta x$ is the finest spatial grid spacing. In contrast, for the Dirichlet and Neumann problems on the unit ball in all dimensions $d\ge 1$, we show that the simplest Volterra marching scheme, i.e., the forward Euler scheme, is unconditionally stable. Our proof is based on an explicit spectral radius bound of the marching matrix, leading to an estimate that an $L^2$-norm of the solution to the integral equation is bounded by $c_dT^{d/2}$ times the norm of the right hand side. For the Robin problem on the half space in any dimension, with constant Robin (heat transfer) coefficient $\kappa$, we exhibit a constant $C$ such that the forward Euler scheme is stable if $\Delta t < C/\kappa^2$, independent of any spatial discretization. This relies on new lower bounds on the spectrum of real symmetric Toeplitz matrices defined by convex sequences. Finally, we show that the forward Euler scheme is unconditionally stable for the Dirichlet problem on any smooth convex domain in any dimension, in $L^\infty$-norm.

Equations

$$\frac{\partial u}{\partial t}(\mathbf{x},t)-\alpha \Delta u(\mathbf{x},t)$$

$$u(\mathbf{x},0)$$

$$u(\mathbf{x},t)=f(\mathbf{x},t){|}_{\mathbf{x}\in \Gamma ,t>0},\mathbf{x}\in \Gamma =\partial D,$$

$$\frac{\partial u(\mathbf{x},t)}{\partial {\mathbf{\nu }}_{\mathbf{x}}}=g(\mathbf{x},t){|}_{\mathbf{x}\in \Gamma ,t>0},\mathbf{x}\in \Gamma ,$$

$$\frac{\partial u(\mathbf{x},t)}{\partial {\mathbf{\nu }}_{\mathbf{x}}}+\kappa u(\mathbf{x},t)=h(\mathbf{x},t){|}_{\mathbf{x}\in \Gamma ,t>0},\mathbf{x}\in \Gamma ,$$

$$\frac{u_n^{m+1}-u_n^m}{\Delta t}={\Delta }_h[u{]}_n^m+F(x_n,t_m)$$

$$\frac{u_n^{m+1}-u_n^m}{\Delta t}={\Delta }_h[u{]}_n^{m+1}+F(x_n,t_m),$$

$$G(\mathbf{x},t)=\frac{1}{(4\pi t{)}^{d/2}}{e}^{-\frac{|\mathbf{x}{|}^2}{4t}},\mathbf{x}\in {\mathbb{R}}^d.$$

$$\mathcal{S}[\sigma ](\mathbf{x},t)={\int }_0^t{\int }_{\Gamma }G(\mathbf{x}-\mathbf{y},t-\tau )\sigma (\mathbf{y},\tau )ds(\mathbf{y})d\tau$$

$$\mathcal{D}[\sigma ](\mathbf{x},t)={\int }_0^t{\int }_{\Gamma }\frac{\partial G(\mathbf{x}-\mathbf{y},t-\tau )}{\partial \mathbf{\nu }(\mathbf{y})}\sigma (\mathbf{y},\tau )ds(\mathbf{y})d\tau ,$$

$$\mathcal{I}[u_0](\mathbf{x},t)={\int }_{D}G(\mathbf{x}-\mathbf{y},t-\tau )u_0(\mathbf{y})d\mathbf{y}$$

$$\mathcal{V}[F](\mathbf{x},t)={\int }_0^t{\int }_{D}G(\mathbf{x}-\mathbf{y},t-\tau )F(\mathbf{y},\tau )d\mathbf{y}d\tau .$$

$$u(\mathbf{x},t)=u^{(F)}(\mathbf{x},t)+u^{(B)}(\mathbf{x},t),$$

$$u^{(F)}(\mathbf{x},t)=\mathcal{I}[u_0](\mathbf{x},t)+\mathcal{V}[F](\mathbf{x},t),$$

$$\left(-\mbox{\small$\frac{1}{2}$}+\mathcal{D}\right)[\sigma](\mathbf{x},t)=\tilde{f}(\mathbf{x},t),\qquad(\mathbf{x},t)\in\Gamma\times[0,T],$$

$$\tilde{f}(\mathbf{x},t):=f(\mathbf{x},t)-u^{(F)}(\mathbf{x},t),\mathbf{x}\in \Gamma .$$

$$\sigma (\mathbf{x},n\Delta t)=$$

$$-2\tilde{f}(\mathbf{x},n\Delta t),n=1,2,\dots$$

$$\sigma (\mathbf{x},t)=-2\tilde{f}(\mathbf{x},t).$$

$$\Vert \sigma {\Vert }_2\le c_d{T}^{d/2}\Vert \tilde{f}{\Vert }_2$$

$$\Vert \sigma {\Vert }_2\le 7\Vert \tilde{f}{\Vert }_2$$

$$\lim_{\epsilon\rightarrow 0+}\frac{\partial\mathcal{S}[\sigma](\mathbf{x}_{0}\pm\epsilon{\bm{\nu}}(\mathbf{x}_{0}),t)}{\partial{\bm{\nu}}(\mathbf{x}_{0})}=\mp\mbox{\small$\frac{1}{2}$}\sigma(\mathbf{x}_{0},t)+\mathcal{S}_{{\bm{\nu}}}[\sigma](\mathbf{x}_{0},t),$$

$$\lim_{\epsilon\rightarrow 0+}\mathcal{D}[\sigma](\mathbf{x}_{0}\pm\epsilon{\bm{\nu}}(\mathbf{x}_{0}),t)=\pm\mbox{\small$\frac{1}{2}$}\sigma(\mathbf{x}_{0},t)+\mathcal{D}[\sigma](\mathbf{x}_{0},t),$$

$$\frac{\partial G(\mathbf{x}-\mathbf{y},t-\tau )}{\partial \mathbf{\nu }(\mathbf{y})}=\frac{(\mathbf{x}-\mathbf{y})\cdot \mathbf{\nu }(\mathbf{y})}{2^{d+1}{\pi }^{d/2}(t-\tau {)}^{1+d/2}}{e}^{-\frac{|\mathbf{x}-\mathbf{y}{|}^2}{4(t-\tau )}}$$

$$\frac{\partial G(\mathbf{x}-\mathbf{y},t-\tau )}{\partial \mathbf{\nu }(\mathbf{x})}=-\frac{(\mathbf{x}-\mathbf{y})\cdot \mathbf{\nu }(\mathbf{x})}{2^{d+1}{\pi }^{d/2}(t-\tau {)}^{1+d/2}}{e}^{-\frac{|\mathbf{x}-\mathbf{y}{|}^2}{4(t-\tau )}}.$$

$$\frac{\partial u}{\partial t}(\mathbf{x},t)-\Delta u(\mathbf{x},t)=F(\mathbf{x},t)$$

$$u(\mathbf{x},t)=\mathcal{D}[\sigma ](\mathbf{x},t)+\mathcal{I}[u_0](\mathbf{x},t)+\mathcal{V}[F](\mathbf{x},t)$$

$$\frac{\partial G(\mathbf{x}-\mathbf{y},t-\tau )}{\partial \mathbf{\nu }(\mathbf{y})}=-\frac{1-\mathbf{x}\cdot \mathbf{y}}{2^{d+1}{\pi }^{d/2}(t-\tau {)}^{1+d/2}}{e}^{-\frac{1-\mathbf{x}\cdot \mathbf{y}}{2(t-\tau )}}.$$

$$f(\theta )=\sum _{n=-\infty }^{\infty }{f}_n{e}^{in\theta },$$

$$f_n=\frac{1}{2\pi }{\int }_0^{2\pi }f(\theta )e^{-in\theta }d\theta .$$

$$H^2=\{f\in L^2(S^1)|f_n=0,n<0\},$$

Full text

1 Introduction

In this paper, we study the stability of integral equation methods for the heat equation

[TABLE]

for $0\leq t\leq T$, subject to suitable boundary conditions, in a smooth domain $D\subset\mathbb{R}^{d}$. Without loss of generality, we will assume that the diffusion coefficient (thermal conductivity) $\alpha$ is one in most of our discussion. We consider three standard boundary conditions: the Dirichlet boundary condition

[TABLE]

the Neumann boundary condition

[TABLE]

and the Robin boundary condition

[TABLE]

Here, $\kappa>0$ is the heat transfer coefficient, and (4) models heat transfer via Newton’s law of cooling [4]. For all three boundary conditions, we assume that proper compatibility conditions are satisfied between the initial and boundary data.

Before turning to the integral equation framework, we briefly review the finite difference approach. For this, we assume we are given a spatial mesh discretizing the domain $D$ with grid points $x_{n}$ and seek to approximate the solution $u_{m}^{n}\approx u(x_{n},t_{m})$ at time steps $t_{0},t_{1},\dots,t_{N}$ with $t_{m}=m\Delta t$. Two of the simplest schemes for solving (1) are the forward and backward Euler methods:

[TABLE]

and

[TABLE]

respectively. Here $\Delta_{h}[u]^{m}_{n}$ denotes the finite difference approximation of the Laplacian evaluated at the grid point $x_{n}$ at time $t_{m}$. It is well known that the backward Euler scheme is unconditionally stable, while, in $d$ dimensions, the forward Euler scheme requires that the time step satisfy the condition $\Delta t<\frac{1}{2d}\Delta x^{2}$, for the case of the standard 2nd-order finite difference Laplacian stencil on a uniform spatial grid with step size $\Delta x$ in each direction (see, for example, [40, p. 158]). The constraint changes when using less standard stencils. For nonuniform grids, the time step restriction is more complicated to analyze, but generally requires that $\Delta t=\mathcal{O}(h_{min}^{2})$ where $h_{min}$ is the finest mesh spacing in the discretization.

The backward Euler scheme is implicit and requires the solution of a large sparse linear system at each time step $t_{m}$. The forward Euler scheme, on the other hand, is explicit and inexpensive. The stability restriction, however, forces extremely small time steps to be taken, making long-time simulations impractical. This has spurred the development of a variety of alternative approaches, including locally one-dimensional schemes, alternating direction implicit methods, etc. [34].

When finite difference methods are used to solve general initial-boundary value problems, GKSO (Gustafsson–Kreiss–Sundström–Osher) theory plays a critical role  [11, 12, 33, 38, 41], and requires that the interior marching scheme be Cauchy stable (that is, beyond the stability condition above, the discrete boundary conditions must satisfy additional criteria). In short, stability imposes rather intricate constraints on the coupling between the interior marching scheme and the boundary conditions themselves. Similar considerations are involved when using finite element methods.

An alternative to direct discretization of the governing PDE is to recast the problem as a boundary integral equation using heat potentials [18, 36]. The Green’s function for the heat equation is

[TABLE]

We assume that the boundary $\Gamma$ of $D$ is at least $C^{2}$, and let $\sigma$ be a square integrable function on $\Gamma\times[0,T]$. Then the single layer heat potential $\mathcal{S}$ is defined by the formula

[TABLE]

and the double layer heat potential $\mathcal{D}$ is defined by

[TABLE]

where ${\bm{\nu}}(\mathbf{y})$ is the unit outward normal vector at $\mathbf{y}\in\Gamma$. The initial potential is defined by

[TABLE]

and the volume potential is defined by

[TABLE]

By the linearity of the problem, we may decompose the solution into

[TABLE]

where all initial and volume data is captured by the free-space volume forced term

[TABLE]

while $u^{(B)}$ is the solution to a pure boundary value problem with zero initial data, zero volume forcing, and modified boundary data. Note that $u^{(F)}$ needs only evaluation of initial and volume potentials; it requires no linear solve. Thus, there is no stability issue with $u^{(F)}$, and its error is simply the quadrature error in evaluating the integrals that appear in (8) and (9). In other words, unlike finite difference or finite element methods, the volume part is completely decoupled from the boundary part in integral equation methods from the perspective of stability analysis.

For the Dirichlet problem, we proceed by representing $u^{(B)}(\mathbf{x},t)$ as a double layer potential with unknown density $\sigma$. The jump relation (see section 2) then leads to the following second kind Volterra equation,

[TABLE]

where $\mathcal{D}$ is interpreted in a principal value sense, and the corrected data is

[TABLE]

The main objective of this paper is to demonstrate certain advantages of integral equation methods by giving, for several combinations of archetypal geometries and boundary conditions, rigorous stability bounds for the simplest explicit time marching scheme, namely the forward Euler scheme. This scheme is derived by assuming $\sigma(\mathbf{y},t)$ is piecewise constant over each time interval $[j\Delta t,(j+1)\Delta t)$, taking on the value $\sigma(\mathbf{y},j\Delta t)$. For (12), this leads to a marching scheme of the form

[TABLE]

This falls into the class of collocation schemes [18, Sec. 13.3], as well as convolution quadrature schemes [23]. It is explicit, since $\sigma(\mathbf{x},n\Delta t)$ does not appear on the right-hand side. It is also first-order accurate (e.g. see section (5.1)). For the Neumann and Robin problems, second kind Volterra equations are obtained by representing $u^{(B)}(\mathbf{x},t)$ instead as a single layer potential; other than a change of kernel, the forward Euler scheme remains the same.

The principal reasons that integral equation methods have received relatively little attention for solving the heat equation has been that direct evaluation of layer (or volume) potentials require quadratic work in the total number of unknowns as well as the design of suitable quadrature rules. Recent advances in fast algorithms for heat potentials, however, have removed this obstacle. We refer the reader to the papers [8, 9, 10, 15, 24, 25, 37, 39, 43, 44, 45] and the references therein for further discussion.

We now summarize the results in this paper. Perhaps the simplest geometry is the half-space $D=\mathbb{R}^{d}_{+}:=\{\mathbf{x}=(x_{1},x_{2},\cdots,x_{d})\in\mathbb{R}^{d}\,|\,x_{d}\geq 0\}$ with $\Gamma=\partial D=\mathbb{R}^{d-1}$. It is easy to check that the integral kernel of $\mathcal{D}$ is identically zero on $\Gamma$, so (12) reduces to

[TABLE]

This is an analytic solution, so that stability follows trivially. A similar trivial analytic solution arises when the single layer potential is used to solve the Neumann problem on the half-space. Thus, we consider the Dirichlet and Neumann problems on possibly the next-simplest domain, the unit ball $B^{d}\subset\mathbb{R}^{d}$ (i.e., $\Gamma$ is the unit sphere $S^{d-1}$). For both these latter cases, we show that the forward Euler scheme is unconditionally stable in all dimensions $d\geq 1$. Specifically, we show that for $T\geq 1$,

[TABLE]

for all $N$, $\Delta t$ such that $N\Delta t\leq T$. Here $N$ is the total number of time steps, $\Delta t$ is the time step size, $\|\cdot\|_{2}$ denotes a space-time $L^{2}$-norm,777Explicitly, $\|\sigma\|_{2}^{2}:=\sum_{j=0}^{N}\int_{\Gamma}\sigma(\mathbf{x},j\Delta t)^{2}ds(\mathbf{y})$, i.e. the norm is $l^{2}$ in time $[0,T]$ but $L^{2}$ over the surface $\Gamma$. and $c_{d}$ is a positive constant depending on $d$. The estimate (15) is obtained by a Gershgorin spectral radius bound of the marching matrix; we show that this is no longer tight for the Dirichlet problem if a fairly mild condition is imposed on $\Delta t$. Indeed, we are able to show the improved estimate in two dimensions,

[TABLE]

for $\Delta t\leq 1$ and any $N$.

Returning to the $d$-dimensional half-space, the simplest boundary condition for which the integral equation is non-trivial is the Robin condition. We show that here the forward Euler scheme has a time step restriction determined by the physical parameter $\kappa$, namely $\Delta t<\frac{\pi}{c^{2}\kappa^{2}}$ with $c=3-\sqrt{2}$. Finally, considering more general domains, we prove that the forward Euler scheme for the Dirichlet problem is unconditionally stable for smooth convex domains in all dimensions, in the $L^{\infty}$-norm.

Firstly, in section 2 we summarize the necessary properties of layer potentials. Then in section 3 we present a lower bound for the spectrum of a Toeplitz operator defined by a convex sequence; this will be needed later to handle cases where the sequences are not summable and thus Gershgorin is inapplicable. The Dirichlet and Neumann problems on the unit ball are then treated in section 4, the Robin problem on the half space in section 5, and the Dirichlet problem on $C^{1}$ convex domains in section 6. We conclude in section 7. Finally, an appendix covers estimates on special functions used in the body of the paper.

2 Properties of heat potentials

By construction, the single and double layer heat potentials (6) and (7) satisfy the heat equation. They also satisfy certain well-known jump conditions when the target point $\mathbf{x}$ approaches the boundary from either side [18, 36]. In particular, for $\mathbf{x}_{0}\in\Gamma$, the normal derivative of the single layer potential $\mathcal{S}[\sigma]$ satisfies the relation

[TABLE]

and the double layer potential $\mathcal{D}[\sigma]$ satisfies the relation

[TABLE]

where both $\mathcal{S}_{{\bm{\nu}}}[\sigma](\mathbf{x}_{0},t)$ and $\mathcal{D}[\sigma](\mathbf{x}_{0},t)$ are interpreted in the Cauchy principal value sense. If we represent the solution to the heat equation (1) via a double layer potential $u(\mathbf{x},t)=\mathcal{D}[\sigma](\mathbf{x},t)$, then the integral equation (12) follows immediately from the jump relation (18).

The kernel of the double layer potential is given explicitly by

[TABLE]

and the kernel of $\mathcal{S}_{{\bm{\nu}}}$ is given by

[TABLE]

Finally, the initial potential (8) is well known to satisfy the homogeneous heat equation with initial data $u_{0}(\mathbf{x})$, while the volume potential (9) satisfies the inhomogeneous heat equation

[TABLE]

with zero initial data.

Remark 1**.**

Using these properties, it is straightforward to see that representing the solution to the Dirichlet problem in the form

[TABLE]

leads to the integral equation (12), with the only unknown corresponding to the double layer density $\sigma$.

Remark 2**.**

On the unit sphere ${S^{d-1}}$, ${\bm{\nu}}(\mathbf{y})=\mathbf{y}$ and $|\mathbf{x}|=|\mathbf{y}|=1$. Thus, $(\mathbf{x}-\mathbf{y})\cdot{\bm{\nu}}(\mathbf{y})=-(1-\mathbf{x}\cdot\mathbf{y})$, $|\mathbf{x}-\mathbf{y}|^{2}=2(1-\mathbf{x}\cdot\mathbf{y})$, and (19) reduces to

[TABLE]

3 Spectral bounds for real symmetric Toeplitz operators

Although for many of the later results we can use simple Gershgorin spectral bounds for matrices, for the tight bound for the zeroth mode of the $d=2$ Dirichlet disc (section 4.2.2), and the Robin case in the half-space (section 5), a more delicate spectral bound on Toeplitz matrices is needed.

Let $S^{1}$ be the unit circle in the complex plane, parametrized by polar angle $\theta$ with normalized arc length measure $d\lambda=\frac{1}{2\pi}d\theta$. For any $f$ in the Hilbert space $L^{2}(S^{1})$, we write

[TABLE]

in terms of the orthogonal basis $\{e^{in\theta}\}_{n\in\mathbb{Z}}$, where $f_{n}$ ($n\in\mathbb{Z}$) is the $n$th Fourier coefficient of $f$ defined by

[TABLE]

The Hardy space $H^{2}$ is defined by

[TABLE]

and we let $P$ denote the orthogonal projection of $L^{2}(S^{1})$ onto $H^{2}$. The Toeplitz operator $T_{f}:H^{2}\rightarrow H^{2}$ with symbol $f\in L^{\infty}(S^{1})$, is defined by

[TABLE]

The operator $T_{f}$ is closely related to an infinite-dimensional Toeplitz matrix with entries $t_{ij},\ i,j\in\mathbb{N}$ that satisfy $t_{ij}=t_{i+1,j+1}$ for all $i,j$. That is, the matrix is constant along diagonals and determined by a two-sided sequence $(t_{n})_{n\in\mathbb{Z}}$ with $t_{ij}=t_{i-j}$. The Fourier transform maps $T_{f}$ onto the class of Toeplitz matrices on $l^{2}(\mathbb{Z}_{+})$; that is, if $\left(T_{f}(u)\right)_{n}$ denotes the $n$th Fourier coefficient of $T_{f}(u)$, then

[TABLE]

where $u_{m}$ is the $m$th Fourier coefficient of $u$.

Definition 1**.**

A sequence $\{a_{n}\}_{n\in\mathbb{Z_{+}}}$ is said to be convex if $\delta^{2}a_{n}\geq 0$ for every $n>0$, where $\delta^{2}a_{n}:=a_{n-1}-2a_{n}+a_{n+1}$ is the central second difference.

Recall that for $n\in\mathbb{Z_{+}}$ the Fejér kernel $F_{n}(x)$ is defined to be

[TABLE]

The following theorem can be found in [17, Chapter 1, Theorem 4.1].

Theorem 1**.**

If $a_{n}\rightarrow 0$ and the sequence $\{a_{n}\}_{n\in\mathbb{Z}_{+}}$ is convex, then the series

[TABLE]

converges in $L^{1}([-\pi,\pi])$ to a non-negative function, which is continuous except at $0,$ such that $v_{n}=a_{n}$.

It is often the case that the function $v(\theta)$ blows up as $\theta\to 0.$ Using the elementary estimate on the Fejér kernel

[TABLE]

[17, Chapter 1, formula 3.10] and the fact that, for a convex sequence tending to zero, we have $\lim_{n\to\infty}n(a_{n}-a_{n+1})=0,$ one can show that

[TABLE]

Bounds on the spectrum of finite Toeplitz matrices are of interest in many applications [5, 14, 19, 26]. When a real symmetric Toeplitz operator (or matrix) is generated by a positive sequence, the Gershgorin circle theorem [40, §3.3] often gives a satisfactory upper bound on its spectral radius or the largest eigenvalue. Curiously, satisfactory lower bounds on the smallest eigenvalue do not seem to be available. The following theorem leads to a tight lower bound on the smallest eigenvalue of a real symmetric Toeplitz matrix, defined by a convex sequence, even when the operator it defines is unbounded.

Theorem 2**.**

Suppose that $\{v_{n}\}_{n\in\mathbb{N}}$ is a convex sequence and $\lim_{n\rightarrow\infty}v_{n}=0.$ Set $v_{0}=2v_{1}-v_{2},$ and let $v(\theta)$ be the non-negative function defined by the sequence $\{v_{n}\}_{n\in\mathbb{Z}_{+}}$ as in Theorem 1. Suppose that $V$ is the self-adjoint Toeplitz matrix defined by $V_{ii}=0$ and $V_{ij}=v_{|i-j|}$. Then, for any ${\bm{u}}\in\mathbb{C}^{N}$, we have the lower bound

[TABLE]

Proof.

For a finite length vector ${\bm{u}}=(u_{0},\dots,u_{N},0,0,\dots)_{n\in\mathbb{Z}_{+}},$ define the function

[TABLE]

Theorem 1 implies that

[TABLE]

∎

Remark 3**.**

If $V_{N}$ is the upper left $N\times N$ principal submatrix of $V$, then, by an application of the Rayleigh–Ritz theorem, its spectrum is bounded below by $(v_{2}-2v_{1}).$

Remark 4**.**

For certain applications the sequence, $\{v_{n}\}_{n\in\mathbb{N}}$, generating $T_{v}$ is not convex. In this case, one may consider an operator of the form, $cI+T_{v}+T_{a}$ with $c$ and $\{a_{n}\}_{n\in\mathbb{N}}$ chosen so that $(c,v_{1}+a_{1},v_{2}+a_{2},\ldots)$ is a convex sequence. If $T_{a}$ is a bounded operator, then the previous theorem implies a lower bound on the spectrum of $V$

[TABLE]

Remark 5**.**

If the function $v(\theta)$ defined in (22) is unbounded, then the Toeplitz operator, $T_{v},$ it defines is not a bounded operator, and is not defined on all of $H^{2}.$ The discussion above easily applies to show that this operator is defined on a dense subset, and its closure is self-adjoint: Equation (23) implies that if $u\in H^{2},$ then $v(1-e^{i\theta})u\in L^{2}.$ Thus, $T_{v}w=P(vw)\in H^{2},$ for $w$ in the subspace $(1-e^{i\theta})H^{2}.$ It is not difficult to see that this subspace is dense. If $u\in H^{2}$ and $r>1,$ then

[TABLE]

and

[TABLE]

Since $\langle T_{v}w,w\rangle\geq 0,$ for $w$ in this domain, the Friedrichs extension of $T_{v}$ is a closed self-adjoint, non-negative operator defined on a dense subspace $D_{v}\subset H^{2}.$

4 The Dirichlet and Neumann problems on the unit ball $B^{d}$

We consider first the forward Euler scheme (13) for the Dirichlet problem (12). For general $d\geq 1$, our approximation of the unknown density $\sigma$ is piecewise constant in time,

[TABLE]

where $t_{j}=j\Delta t$. We restate (13) in the form

[TABLE]

for $j=0,1,2,\ldots$, where the tilde has been dropped from $f$, and where the action of each spatial integral operator $V_{j-k}:C(\Gamma)\to C(\Gamma)$ is defined by

[TABLE]

Here the operator kernel is itself the integral of the heat kernel over one time-step,

[TABLE]

Due to time-shift invariance, a simpler way to write the spatial kernel is

[TABLE]

and $\mathcal{V}_{0}(\mathbf{x},\mathbf{y})$ is set identically to [math]. For initialization of time-stepping we set $\sigma_{0}\equiv f_{0}\equiv 0$.

4.1 The Dirichlet problem in one dimension

The boundary $\Gamma$ of the unit ball in one dimension consists of only two points $\mathbf{x}=\pm 1$. Let the time-stepped density at these two points be $\bm{\sigma}^{\pm}=\{\sigma_{j}^{\pm}\}_{j=0}^{N}$, and the data ${\bm{f}}^{\pm}=\{f^{\pm}_{j}\}_{j=0}^{N}$. We will stack each pair into a single column, e.g. $[\bm{\sigma}^{-}\,,\,\bm{\sigma}^{+}]^{T}$. Recalling (14), the density at each boundary point is trivially coupled to the data at that same point; however, the coupling to the other boundary point will involve the double layer kernel acting at a distance of 2. Thus, (25) becomes a $2\times 2$ system with trivial diagonal blocks and Toeplitz off-diagonal blocks. Namely, after $N$ time-steps the stacked vectors are related by,

[TABLE]

where $I$ is the size-$(N+1)$ identity matrix. Here the action of the lower-triangular Toeplitz matrix $V\in\mathbb{R}^{(N+1)\times(N+1)}$ is given by

[TABLE]

with the convolution coefficients $\{v_{l}\}$ given by

[TABLE]

Here the underlying kernel is the double layer acting at a distance of 2,

[TABLE]

We denote the symmetric part of $V$ by $W$, and make its dependence on $N$ and $\Delta t$ explicit, thus

[TABLE]

We have the following lemma.

Lemma 1**.**

Fix $T>0$. Then, for any $N$ and $\Delta t$ with $N\Delta t\leq T$, the spectral radius $\rho(N;\Delta t)$ of the matrix $W(N;\Delta t)$ has the bound

[TABLE]

where

[TABLE]

Proof.

Using the Gershgorin circle theorem [40, §3.3], and the fact that the diagonal entries of $W$ are all zero, we have

[TABLE]

Now setting $t=N\Delta t$, we may collapse this sum into a single integral

[TABLE]

according to the definition (31) of the function $C_{1}$. Combining the last two results we have $\rho(N;\Delta t)\leq C_{1}(N\Delta t)$. The expression in (31) follows from the change of variables $u=\frac{1}{T-\tau}$. A further change of variables $x=\sqrt{u}$ leads to

[TABLE]

Finally, the above expression shows that $C_{1}(T)$ is a monotonically non-decreasing function of $T$, so that $\rho(N;\Delta t)\leq C_{1}(N\Delta t)\leq C_{1}(T)$. ∎

It is clear from (26) that to get a stability bound we need to control the gap between $C_{1}(T)$ and $\frac{1}{2}$. For $T\geq 1$, this turns out to shrink only polynomially in $T$:

[TABLE]

This very easily allows us to prove the following.

Theorem 3**.**

Suppose that $T\geq 1$. Then, using $\|.\|$ for the $l^{2}$-norm in $\mathbb{R}^{2(N+1)}$,

[TABLE]

for all $N$, $\Delta t$ such that $N\Delta t\leq T$. That is, for the $d=1$ unit ball where $\Gamma=\{-1,1\}$, the forward Euler scheme for solving the second kind Volterra integral equation (12) is unconditionally stable on any finite time interval $[0,T]$.

Proof.

We use a technique that will recur throughout this paper: we take the inner product of (26) with $-[\bm{\sigma}^{-},\bm{\sigma}^{+}]^{T}$, giving

[TABLE]

Applying the Cauchy–Schwarz inequality and (30) to the second term on the left side, and Cauchy–Schwarz to the right hand side, we obtain

[TABLE]

Using the arithmetic-geometric mean inequality on the left hand side of this estimate gives

[TABLE]

Finally dividing by $\left(\frac{1}{2}-C_{1}(T)\right)\|[\bm{\sigma}^{+},\bm{\sigma}^{-}]\|$ and applying (34) gives

[TABLE]

which completes the proof. ∎

4.2 The Dirichlet problem in two dimensions

We now consider (25) when $\Gamma$ is the unit circle ${S^{1}}$. We decompose both $\sigma_{j}(\mathbf{y})$ and $f_{j}(\mathbf{x})$ into Fourier series:

[TABLE]

From (20), writing $s=\theta-\phi$, the $n$th Fourier mode of the kernel is

[TABLE]

where, noting that the imaginary part of $e^{-ins}$ cancels by symmetry, we have

[TABLE]

Since $\{e^{in\theta}\}$ are orthonormal, each Fourier mode evolves independently. The marching scheme (or recurrence) (25) for the $n$th mode is then

[TABLE]

where the convolution coefficient $v^{n}_{l}$ is given by the formula

[TABLE]

and we set $v^{n}_{0}=0$. The system (39) for $j=0,1,\cdots,N$ can be written in matrix-vector form

[TABLE]

where $I$ is the $(N+1)\times(N+1)$ identity matrix, $V^{n}\in\mathbb{R}^{(N+1)\times(N+1)}$ with entries $v^{n}_{j,k}=v^{n}_{j-k}$, $\bm{\sigma}^{n}=\{\sigma_{j}^{n}\}_{j=0}^{N}$, ${\bm{f}}^{n}=\{f^{n}_{j}\}_{j=0}^{N}$. As before, the symmetric part of $V^{n}$ we denote by $W^{n}$,

[TABLE]

4.2.1 Stability analysis

We now prove two key results. The first is that the forward Euler scheme is unconditionally stable for any fixed time interval $[0,T]$ (Theorem 4). The second is that, when $\Delta t<1$, we have the stronger result that the $L^{2}$-norm of the solution is bounded by a constant multiple of the $L^{2}$-norm of $f$. We require the following lemma.

Lemma 2**.**

Fix $T>0$. Then, for any $N$ and $\Delta t$ with $N\Delta t\leq T$, and all $n\in\mathbb{Z}$, the spectral radius $\rho_{n}(N;\Delta t)$ of the matrix $W^{n}(N;\Delta t)$ has the bound

[TABLE]

where, in terms of the definition (38),

[TABLE]

Proof.

Let $n\in\mathbb{Z}$. Since the integrand in (38), excluding the $\cos ns$ factor, is non-negative, we observe that $|\gamma_{n}(t)|\leq\gamma_{0}(t)$, so $|v^{n}_{l}|\leq v^{0}_{l}$ for all $l\geq 1$. Using this, the Gershgorin theorem, and the fact that the diagonal entries of $W^{n}$ are all zero, we have

[TABLE]

Now setting $t=N\Delta t$, we may collapse this sum into a single integral

[TABLE]

according to the definition (44) of the function $C_{2}$. Combining the last two results we have $\rho_{n}(N;\Delta t)\leq C_{2}(N\Delta t)$. To prove the expression in (44) we insert (38), interchange the order of integration and apply the change of variables $\lambda=\frac{1-\cos(s)}{2(T-\tau)}$, thus

[TABLE]

Finally, the above expression shows that $C_{2}(T)$ is a monotonically increasing function of $T$, so that $\rho_{n}(N;\Delta t)\leq C_{2}(N\Delta t)\leq C_{2}(T)$. ∎

Analogous to before, it is clear from (41) that to get a stability bound we need to bound from below the gap between $C_{2}(T)$ and $\frac{1}{2}$. This motivates the following.

Proposition 1**.**

[TABLE]

where $I_{n}(\cdot)$ is the modified regular Bessel function of order $n$ (see Appendix). For $T\geq 1$,

[TABLE]

Proof.

(46) follows from the integral representation of $I_{0}(x)$ (102). (47) follows from the facts that $I_{0}(x)\leq 1+\frac{x^{2}}{2}$ [32, §10.25.2] and $e^{-x}\leq 1-\frac{x}{2}$ for $x\leq 1$. ∎

Theorem 4**.**

Suppose that $T\geq 1$. Then for all $n\in\mathbb{Z}$,

[TABLE]

for all $N$, $\Delta t$ such that $N\Delta t\leq T$. That is, when $\Gamma$ is the unit circle ${S^{1}}$, the forward Euler scheme for solving the second kind Volterra integral equation (12) is unconditionally stable on any finite time interval $[0,T]$.

Proof.

Since we are working in the Fourier domain, $\bm{\sigma}^{n}$ and ${\bm{f}}^{n}$ are complex-valued. Thus, we split (41) into two independent real systems

[TABLE]

where $\bm{\sigma}_{r}^{n}$ and $\bm{\sigma}_{i}^{n}$ are the real and imaginary part of $\bm{\sigma}^{n}$, respectively.

Multiplying both sides of the first equation in (49) by $-(\bm{\sigma}_{r}^{n})^{T}$, we have

[TABLE]

Applying (43) on the left side (50) and the Cauchy–Schwartz inequality on the right side, we obtain

[TABLE]

That is, finally applying Proposition 47,

[TABLE]

Similar result holds for $\|\bm{\sigma}_{i}^{n}\|$. As $\|\bm{\sigma}^{n}\|=\sqrt{\|\bm{\sigma}_{r}^{n}\|^{2}+\|\bm{\sigma}_{i}^{n}\|^{2}},$ these two inequalities give (48). ∎

4.2.2 Tighter bounds

We now show that the dependence on $T$ in (48) can be removed when the time step satisfies $\Delta t\leq 1$. This is a physically reasonable requirement since we have assumed that the diffusion coefficient is one, and the domain has area of order one. We first provide a bound on $\rho_{n}(N;\Delta t)$ for $n\neq 0$ that is independent of the total time $N\Delta t$.

Lemma 3**.**

Let $N$ and $\Delta t>0$ be arbitrary, and let $\rho_{n}(N;\Delta t)$ be the spectral radius of $W^{n}(N;\Delta t)$ defined in (42). Then for all $n\neq 0$,

[TABLE]

Proof.

Clearly, it is sufficient to prove (51) for $n>0$. For this, let us note that substituting (38) into (40), exchanging the order of integration, and making the change of variables $\lambda=\frac{1-\cos(s)}{2(l\Delta t-\tau)}$, we obtain

[TABLE]

By the integral representation (102) of $I_{n}$, we have

[TABLE]

From (52), defining $i_{n}(x):=e^{-x}I_{n}(x)$ and $x_{l}=1/(2l\Delta t)$, we consider the sum

[TABLE]

By Lemma 11 (see Appendix), the function $i_{n}(x)$ assumes its unique maximum at $r_{n}>0$; it increases monotonically on $[0,r_{n}]$ and decreases monotonically on $[r_{n},+\infty).$ We now consider (53) on a case-by-case basis.

1. (a)

All $x_{l}$ lie on $[0,r_{n}]$: Since $x_{l}<x_{l-1}$ and $i_{n}(x)$ increases on $[0,r_{n}]$, we have

[TABLE]

where the last inequality follows from (106). 2. (b)

All $x_{l}$ lie on $[r_{n},\infty)$: In this case, we have

[TABLE] 3. (c)

$x_{1}>\cdots>x_{m}\geq r_{n}>x_{m+1}>\cdots>x_{N}$: In this case, we have

[TABLE]

By (45) we have

[TABLE]

completing the proof. ∎

Corollary 1**.**

For all $n\neq 0$,

[TABLE]

Thus, all non-zero modes are unconditionally stable. The zeroth Fourier mode is a bit more subtle, and requires the convex sequence results of section 3. It brings in a weak restriction on $\Delta t$, as follows.

Lemma 4**.**

Suppose that $a=0.05$ and $\Delta t\leq 1$. Then $c_{2}I+W^{0}+aW^{1}$ is a positive definite matrix if

[TABLE]

Proof.

Define the sequence y_{j}=\mbox{\small\frac{1}{2}}(v_{j}^{0}+av_{j}^{1}) for $j\geq 1$ and $y_{0}=c_{2}$. Theorem 2 then shows that a sufficient condition for the positive semi-definiteness of $c_{2}I+W^{0}+aW^{1}$ is that the sequence $\{y_{j}\}_{j\in\mathbb{Z}_{+}}$ is convex. But $y_{1}=\frac{1}{4}f(x_{1})$, $y_{j}=\frac{1}{4}(f(x_{j})-f(x_{j-1}))$ ($j>1$), where $f$ is the function defined in Lemma 15, and $x_{j}=2j\Delta t$. That is, $y_{j}$ is the first order difference of $f$. Furthermore, the convexity of $\{y_{j}\}_{j\in\mathbb{N}}$ is equivalent to the non-negativity of the third order difference of $f$, which follows from the fact that $f^{\prime\prime\prime}(x)>0$ for all $x>0$ as proved in Lemma 15. For $j=0$, the convexity of the sequence requires that one choose $c_{2}$ such that

[TABLE]

By the integral representation (102) of $I_{0}$, it is easy to see that $e^{-x}I_{0}(x)$ is strictly decreasing. Thus, we have e^{-1/2}I_{0}\left(\mbox{\small\frac{1}{2}}\right)\geq e^{-\frac{1}{2\Delta t}}I_{0}\left(\frac{1}{2\Delta t}\right) for $\Delta t\leq 1$. Furthermore,

[TABLE]

by (106). Hence, (55) is achieved by choosing

[TABLE]

for $\Delta t\leq 1$. ∎

Corollary 2**.**

Suppose that $\Delta t\leq 1$. Then, for arbitrary $N$,

[TABLE]

Proof.

Set $a=0.05$. By Lemma 4, the smallest eigenvalue of $W^{0}$ is bounded by

[TABLE]

Thus a simple bound using the value of $c_{2}$ from Lemma 4 is

[TABLE]

completing the proof. ∎

4.3 The Dirichlet problem in higher dimensions

In dimensions $d>2$, we consider the Dirichlet problem on the unit ball, with data specified on the unit sphere ${S^{d-1}}$. The unknown density $\sigma$ is decomposed using the corresponding spherical harmonics [29]

[TABLE]

where

[TABLE]

Here, $a_{n,d}$ is the dimension of $H_{n}({S^{d-1}})$, the space of homogeneous harmonic polynomials of degree $n$ on $\mathbb{R}^{d},$ whose restrictions to the unit sphere are spanned by $\{Y^{m}_{n}\},$ the spherical harmonics of degree $n$. When $d=3$, $a_{n,d}=2n+1$, the inner summation is usually written as $\sum_{m=-n}^{n}$, and the spherical harmonics $Y^{m}_{n}(\theta,\phi)$ are defined by

[TABLE]

where $P_{n}^{m}(\cos\theta)$ is the associated Legendre polynomial [32, §18.3] of degree $n$ and order $m$.

The spherical harmonics admit the following integral representation [29]

[TABLE]

where $\omega_{d}$ is the area of ${S^{d-1}}$ defined in (87), and the $P_{n,d-1}$ are Gegenbauer polynomials [29, Chapter 2] (also called ultraspherical polynomials), defined by the Rodrigues formula

[TABLE]

The Funk–Hecke formula [29, Chapter 2, Theorem 2.39] states that

[TABLE]

where

[TABLE]

and $f$ is any measurable function such that

[TABLE]

In $d=3$ this reduces to $f\in L^{1}[-1,1]$.

We compute the double layer heat potential $nm$th Fourier mode,

[TABLE]

where, by analogy with (38),

[TABLE]

The third equality makes use of (56), (58), and exchanging the order of integration. The last step follows again from (56). Notice that $\gamma_{n,d}$ does not depend on the order $m$.

Since the $\{Y^{m}_{n}\}$ form an orthonormal basis for functions in $L^{2}({S^{d-1}})$ and (59) shows that each spherical harmonic evolves independently under the action of the double layer heat potential operator, we may consider the time evolution for each mode $nm$ separately.

For the forward Euler scheme, we again assume that $\sigma(\mathbf{x},t)$ takes the constant value $\sigma_{j}(\mathbf{x})=\sigma(\mathbf{x},j\Delta t)$ over each interval $[j\Delta t,(j+1)\Delta t]$, $j=0,1,\dots$. Equivalently, each spherical harmonic mode $\sigma^{nm}(t)$ takes the constant value $\sigma^{nm}_{j}=\sigma^{nm}(j\Delta t)$ over the interval. A straightforward calculation leads to the following recurrence for the $nm$th spherical harmonic mode, analogous to (39):

[TABLE]

where we use the abbreviations $\mu_{j}:=\sigma^{nm}_{j}$, $g_{j}=f^{nm}_{j}$, and the matrix elements

[TABLE]

involve the kernel modes (60), and, as before, $v^{n}_{0}=0$.

4.3.1 Stability analysis

The normalization in (57) leads to [28, 30]

[TABLE]

As the other terms in (60) are non-negative, we have

[TABLE]

An almost identical proof as in Lemma 44 leads to the following lemma.

Lemma 5**.**

Fix $T>0$. Then, for any $N$ and $\Delta t$ with $N\Delta t\leq T$, and all $n\in\mathbb{Z}_{+}$, the spectral radius $\rho_{n,d}(N;\Delta t)$, of the symmetric Toeplitz matrix $W^{n}(N;\Delta t)$ as defined by (42) with $v_{l}^{n}$ given by (62), has the bound

[TABLE]

where

[TABLE]

As before, we are also able to bound from below the gap between $C_{d}(T)$ and $\frac{1}{2}$, given a weak condition on $T$. For this, we interchange the order of integration and apply the change of variable $\lambda=\frac{1-x}{2(T-\tau)}$, giving

[TABLE]

and

[TABLE]

Assume now $T\geq 1$. Then for $x\in[-1,1]$, $\frac{1-x}{2T}\leq 1$. Thus, $e^{-\lambda}\geq e^{-1}$ for $\lambda\in[0,\frac{1-x}{2T}]$ and

[TABLE]

where the last equality follows from an integral identity in [7, §3.62].

Armed with this polynomial control of the gap, and following the same reasoning as used to show (48), we obtain the following theorem regarding the stability of the forward Euler scheme in higher dimensions.

Theorem 5**.**

Fix $d>2$, and $T\geq 1$. For all $n=0,1,\ldots$ and $m=1,\ldots,a_{n,d}$,

[TABLE]

for all $N$, $\Delta t$ such that $N\Delta t\leq T$. That is, when $\Gamma$ is the unit sphere ${S^{d-1}}$, the forward Euler scheme for solving the second kind Volterra integral equation (12) is unconditionally stable on any finite time interval $[0,T]$.

Remark 6**.**

When $d=2$, $H_{n}$ is spanned by $e^{in\theta}$ and $e^{-in\theta}$. The decomposition of $L^{2}(S^{1})$ into spherical harmonics is the usual Fourier series expansion. And if we identify $P_{n,1}(x)$ with the Chebyshev polynomials $T_{n}(x)$, then all calculations in this subsection are valid for $d=2$. We instead presented the analysis in two dimensions using the usual Fourier series for the reader’s convenience.

Remark 7**.**

It is easy to see that the bound (64) actually also includes the cases of $d=1$ and $d=2$ proved earlier, thus holds for all $d\geq 1$.

4.4 The Neumann problem on the unit ball

For the Neumann condition (3), we represent $u^{(B)}$ as the single layer potential $\mathcal{S}[\sigma]$. The jump relation (17) leads to the second kind Volterra integral equation

[TABLE]

where $\mathcal{S}_{{\bm{\nu}}}$ indicates the normal derivative of the single layer with respect to the target point, restricted to $\Gamma$, interpreted in a principal value sense, as in section 2. In (65) the right hand side is the corrected data

[TABLE]

On the unit sphere ${S^{d-1}}$, straightforward calculation shows that the kernel of the double layer potential $\mathcal{D}$ is exactly the same as that of $\mathcal{S}_{\bm{\nu}}$. Thus, the forward Euler scheme for (65) leads to the identical marching matrix except a sign change in diagonal entries. Since we prove the bound (64) by bounding the spectral radius of the marching matrix excluding the diagonal part, we observe that (64) holds as well for (65) with $f$ replaced by $\tilde{g}$. This leads to the unconditional stability of the forward Euler scheme for the Neumann problem on the unit ball, for all $d\geq 1$.

5 The Robin problem on the half space

For the Robin boundary condition (4), we also represent $u^{(B)}$ via a single layer potential $\mathcal{S}[\sigma]$. The jump relation (17) leads to the 2nd-kind Volterra equation

[TABLE]

with corrected Robin data

[TABLE]

When $D=\mathbb{R}^{d}_{+}$, where $\Gamma$ is naturally identified as $\mathbb{R}^{d-1}\subset\mathbb{R}^{d}$, the kernel of $\mathcal{S}_{\bm{\nu}}$ is identically zero due to the fact that $(\mathbf{x}-\mathbf{y})\cdot{\bm{\nu}}_{\mathbf{x}}=0$. Thus, (66) reduces to

[TABLE]

Here we assume that $\tilde{h}$ is sufficiently smooth and decays sufficiently fast at infinity so that the problem is well posed.

5.1 The Robin problem in one dimension

In one dimension, the boundary $\Gamma$ of the half line consists of a single point $x=0$. The integral equation (68) reduces to the Abel integral equation (multiplying both sides by two, and denoting the right-hand side by $f$ instead):

[TABLE]

Before discretizing, we show stability of the continuous problem for $\kappa>0$. The Riemann–Liouville fractional integral operator $\mathcal{R}_{\alpha}$ is defined by the formula

[TABLE]

where $\Gamma(\alpha)$ is the gamma function (88). Thus, the integral operator on the left side of (69) is simply $\frac{\Gamma(1/2)\kappa}{\sqrt{\pi}}\mathcal{R}_{1/2}=\kappa\mathcal{R}_{1/2}$. For all real functions $g$, $\mathcal{R}_{\alpha}$ satisfies the positivity property [31, Lemma 3.1]

[TABLE]

Taking the inner product of (69) with $\sigma$ over a fixed interval $[0,T]$, and using (70), gives

[TABLE]

where Cauchy–Schwartz was used in the last step. So on any finite interval $[0,T]$ this gives the continuous version of the $L^{2}$ stability bound

[TABLE]

We now proceed to discretization. Recall that the forward Euler scheme uses a piecewise constant approximation $\sigma(t)\approx\sigma_{m}:=\sigma(t_{m})$ on $[t_{m},t_{m+1})$ on the uniform grid $t_{m}=m\Delta t$. Then performing the integrals exactly in (69) gives the explicit marching rule

[TABLE]

with the lower-triangular Toeplitz matrix weights

[TABLE]

and where $f_{n}:=f(t_{n})$, and $h:=\kappa^{2}\,\Delta t/\pi$. For smooth solutions $\sigma\in C^{1}([0,T])$, this rule can be proved to be first-order accurate by combining compactness of the integral operator, Céa’s lemma, and noting that the piecewise constant approximant has error ${\mathcal{O}}(\Delta t)$ (see [18, Sec. 13.1–3]).

For initialization, as before we set $\sigma_{0}=f_{0}=0$. Let us define the vectors $\bm{\sigma}$ and ${\bm{f}}$ by $\{\sigma_{n}\}_{n=0}^{N},\{f_{n}\}_{n=0}^{N}\in\mathbb{R}^{N+1}$, respectively. Using this notation, (71) takes the form of the lower-triangular Toeplitz linear system

[TABLE]

where $V\in\mathbb{R}^{(N+1)\times(N+1)}$ has elements $v_{n,m}=v_{n-m}$ for $n>m$, and $v_{n,m}=0$ otherwise. Here, $v_{n}$ is defined in (72) with $h=\kappa^{2}\,\Delta t/\pi$.

There is a substantial literature on the numerical analysis and stability of Volterra equations in the one-dimensional setting. For a discussion of convergence theory and step-size control, see [1, 16] and the monograph [3]. Much work on stability has been devoted to an analysis of the model problem

[TABLE]

or to problems with a continuous kernel [16, 27]. In [20], a more relevant stability result is obtained for systems of the form (73), but assuming that the sequence $\{v_{j}\}$ is in $l^{1}$, which is not the case here.

For previous work on Abel-type equations with singular kernels, we refer the reader to [6, 21, 22, 42]. These papers, however, are mostly concerned with implicit marching schemes. An exception is Lubich’s 1986 paper [23], which does a careful stability analysis for a variety of schemes and makes clear the connection between completely monotonic sequences and stability. An interesting result from that paper is Corollary 2.2, which states that “the stability region of an explicit convolution quadrature $\dots$ is bounded.” Theorem 6 below, which is consistent with Lubich’s result, gives a precise value for the time step restriction. It also guarantees that $\sigma$ decays once the right-hand side $f$ has switched off.

Theorem 6**.**

There is a constant $0<c<3-\sqrt{2}$ such that, for any $N$ and any ${\bm{f}}\in\mathbb{R}^{N+1}$, the solution to (73) obeys

[TABLE]

where $\|.\|$ denotes the $l^{2}$-norm. That is, the marching scheme (71) is stable for $h<0.39<(1/c)^{2}$ or $\Delta t<\pi/(c^{2}\kappa^{2})$, where $\kappa$ is the heat transfer coefficient.

Proof.

We first show that there exists a constant $c>0$ such that

[TABLE]

*i.e. *that the smallest eigenvalue of $V$ is bounded from below. Writing \sqrt{h}W_{N+1}:=\mbox{\small\frac{1}{2}}(V+V^{T}) as the scaled symmetric part of $V$, note that $\bm{\sigma}^{T}V\bm{\sigma}=\sqrt{h}\bm{\sigma}^{T}W_{N+1}\bm{\sigma}$, and that $W_{N+1}$ is independent of the time-step. Note that $W_{N+1}$ is the $(N+1)\times(N+1)$ upper left principal submatrix of the infinite symmetric Toeplitz matrix $T_{v}$, defined by the sequence $0,v_{1},v_{2},\ldots$ with

[TABLE]

It is straightforward to check that the sequence $\{v_{j}\}_{j\in\mathbb{N}}$ is convex and that $\lim_{j\rightarrow\infty}v_{j}=0$. By Theorem 2 and Remark 3, we have

[TABLE]

That is, (75) holds if $c=2v_{1}-v_{2}=3-\sqrt{2}$. To complete the proof, take the inner product of (73) with $\bm{\sigma}$ to get

[TABLE]

Applying (75) to the left-hand side and the Cauchy–Schwartz inequality to the right-hand side, we have

[TABLE]

from which (74) follows for any $\bm{\sigma}\neq{\bm{0}}$. It holds trivially when $\bm{\sigma}={\bm{0}}$. ∎

Remark 8**.**

The above proof gives $c=3-\sqrt{2}\approx 1.5858$. By numerically computing the smallest eigenvalue of successively larger Toeplitz matrices $V$, or, better, by evaluating $v(\pi)=2\sum_{j>0}(-1)^{j-1}v_{j}$, one can obtain an optimal estimate of $c\approx 1.52041925043874$. We omit the details of this computation and mention it only to illustrate that the explicit bound is within about 4% of the optimal one.

Remark 9**.**

With unit diffusion constant, the transfer coefficient $\kappa$ has units $(\mbox{length})^{-1}$. Thus our time-step condition $\Delta t<\pi/(c\kappa)^{2}$ is proportional to the square of the physical length $1/\kappa$. Although reminiscent of the explicit finite-difference stability condition $\Delta t<c\Delta x^{2}$, our stability condition is, by contrast, independent of any spatial discretization. (Indeed, in practice the only spatial discretization needed would be quadrature to evaluate (11) to get $f(t)$, as in (67). With $f(t)$ computed, there is no spatial variable left to discretize.)

5.2 The Robin problem in higher dimensions

In higher dimensions $d\geq 2$, the boundary $\Gamma$ of the half space $\mathbb{R}^{d}_{+}$ can be identified as $\mathbb{R}^{d-1}$ by natural embedding. The integral equation (68) is rewritten as

[TABLE]

Again we have multiplied both sides by two and denote the right hand side by $f$. We observe that the kernel inside the spatial integral on the left side of (76) is exactly the heat kernel in $\mathbb{R}^{d-1}$. It is well known that the Fourier transform of the heat kernel $G(\mathbf{x},t)$ in $\mathbb{R}^{d-1}$ is simply $e^{-|{\bm{\xi}}|^{2}t}$, in terms of the Fourier variable ${\bm{\xi}}\in\mathbb{R}^{d-1}$. Using this fact and that the convolution in physical space becomes pointwise multiplication in frequency, and taking the Fourier transform in $\mathbb{R}^{d-1}$ of both sides of (76), we obtain

[TABLE]

Note that in the special case ${\bm{\xi}}=0$ this recovers (69).

Fixing ${\bm{\xi}}$, we proceed similarly as in the one-dimensional case. That is, we approximate $\hat{\sigma}({\bm{\xi}},t)$ by a constant $\hat{\sigma}_{m}({\bm{\xi}}):=\hat{\sigma}_{m}({\bm{\xi}},t_{m})$ on $[t_{m},t_{m+1})$ with $t_{m}=m\Delta t$, and perform the integrals exactly. Let us define the vectors $\hat{\bm{\sigma}}({\bm{\xi}})$ and $\hat{\bm{f}}({\bm{\xi}})$ by $\{\hat{\sigma}_{n}({\bm{\xi}})\}_{n=0}^{N},\{\hat{f}_{n}({\bm{\xi}})\}_{n=0}^{N}\in\mathbb{R}^{N+1}$, respectively. Using this notation, the forward Euler scheme for (77) takes the form of the lower-triangular Toeplitz linear system

[TABLE]

where $\hat{V}({\bm{\xi}})\in\mathbb{R}^{(N+1)\times(N+1)}$ has elements $v_{n,m}({\bm{\xi}})=v_{n-m}({\bm{\xi}})$ for $n>m$, and $v_{n,m}=0$ otherwise. Here, $v_{n}$ is defined by

[TABLE]

with, as before, $h=\kappa^{2}\,\Delta t/\pi$.

Lemma 6**.**

For any $\Delta t>0$ and any fixed ${\bm{\xi}}$, the sequence $\{v_{n}({\bm{\xi}})\}_{n\in\mathbb{N}}$ is convex.

Proof.

Let $x=|{\bm{\xi}}|^{2}\Delta t$. Applying the change of variables $u=n-\tau/\Delta t$ on the integral in (79) leads to

[TABLE]

Thus, in order to show that $\{v_{n}({\bm{\xi}})\}_{n\in\mathbb{N}}$ is a convex sequence, we only need to show that the function

[TABLE]

is convex for $t\geq 1$. Here $x\geq 0$ is a fixed parameter. Differentiating $g(t)$ twice leads to

[TABLE]

with

[TABLE]

Now

[TABLE]

which is positive for any $x\geq 0$ and $t>0$. This shows that $p^{\prime}(t)$ is monotonically increasing for $t>0$. Thus, $p^{\prime}(t)>p^{\prime}(t-1)$ for $t\geq 1$, and $g^{\prime\prime}(t)>0$ for $t\geq 1$, completing the proof. ∎

Lemma 7**.**

For any ${\bm{\xi}}\in\mathbb{R}^{d-1}$,

[TABLE]

Proof.

Using the expression (80), we only need to show that

[TABLE]

for $x\geq 0$. For this, we calculate

[TABLE]

That is, $f$ is monotonically increasing for $x\geq 0$, completing the proof. ∎

Lemma 6 together with Theorem 2 and Remark 3 leads to

[TABLE]

Combining the above estimate with (81), we obtain

[TABLE]

where

[TABLE]

An argument similar to that in the proof of Theorem 4 then gives

[TABLE]

Taking the $L^{2}$-norm in Fourier space and then applying the Plancherel theorem, we have

[TABLE]

That is, we obtain exactly the same bound (74) as in one dimension, which shows that the forward Euler scheme is stable for (76) if $\Delta t<\pi/(c^{2}\kappa^{2})$, where $\kappa$ is the heat transfer coefficient.

Remark 10**.**

In the limit $\kappa\rightarrow 0$, the scheme is unconditionally stable. This is to be expected, since when $\kappa=0$, the Robin boundary condition becomes a Neumann condition and the integral equation (68) yields the analytic solution $\sigma(\mathbf{x},t)=2\tilde{h}(\mathbf{x},t)$.

6 The Dirichlet problem on an arbitrary smooth convex domain

We now study the stability property of the forward Euler scheme (13) for the Dirichlet problem on an arbitrary $C^{1}$ convex domain, i.e., the boundary integral equation (12).

We first establish a connection between the heat kernel and the Laplace kernel. The Green’s function for the Laplace equation in $\mathbb{R}^{d}$ is

[TABLE]

where

[TABLE]

is the area of the unit sphere ${S^{d-1}}\subset\mathbb{R}^{d}$. Here $\Gamma$ is the gamma function defined by the formula

[TABLE]

The kernel of the Laplace double layer potential operator is given by

[TABLE]

It is well known to satisfy Gauss’ Lemma [18]:

[TABLE]

Lemma 8**.**

[TABLE]

Proof.

By (19), we have

[TABLE]

The change of variables $\lambda=\frac{|\mathbf{x}-\mathbf{y}|^{2}}{4(t-\tau)}$ leads to

[TABLE]

Taking the limit $t\rightarrow\infty$ and using the definition of the gamma function (88), we obtain (91). ∎

The following provides the key ingredient for the stability of the forward Euler scheme in an arbitrary smooth convex domain.

Lemma 9**.**

Suppose that $D$ is a $C^{1}$ convex domain. Then

[TABLE]

and

[TABLE]

For $t\in(0,\infty)$, define

[TABLE]

Then $C(t)$ is a monotonic increasing function of $t$ and

[TABLE]

Proof.

(93) follows from the expressions (19) and (89) and the fact that $\mathbf{x}\cdot\mathbf{y}\leq 0$ for $\mathbf{x},\mathbf{y}\in\Gamma$ when $D$ is convex due to the convex separation theorem [2]. (94) follows from (90) and (91). The monotonic increasing property of $C(t)$ follows from (92) and the fact that the integrand is of the same sign everywhere by (93). Finally, (96) is a simple consequence of (93) and (94). ∎

Recall that the forward Euler scheme (13) for the Dirichlet problem is

[TABLE]

Here we have dropped the tilde from $f$ again.

Theorem 7**.**

Let $D\subset\mathbb{R}^{d}$ be a bounded, convex domain with $C^{1}$-boundary. Fix $T>0$. The solution $\sigma$ to (97) satisfies

[TABLE]

for any $N$, $\Delta t$ such that $N\Delta t\leq T$. Here $C(T)$ is defined in (95), $\|\cdot\|_{\infty}$ denotes the $L^{\infty}$ norm in space and the $l^{\infty}$ norm in the discrete temporal variable. In other words, the forward scheme (13) is unconditionally stable on $[0,T]$ for any $T>0$.

Proof.

Taking the absolute value on both sides of (97), we have

[TABLE]

where the first inequality follows from the triangle inequality, the second one following from taking the $L^{\infty}$ norm in the spatial variable for both $\sigma$ and $f$, and the third one follows from taking the maximum norm in the discrete temporal variable. We continue our calculation

[TABLE]

Since the above inequality is valid for any $\mathbf{x}\in\Gamma$ and any $n$ such that $n\Delta t\leq T$, its left hand side can be replaced by \mbox{\small\frac{1}{2}}\|\sigma\|_{\infty}, completing the proof. ∎

7 Conclusions and further remarks

We have analyzed the stability of the forward Euler scheme for solving the Dirichlet and Neumann problems for the heat equation in the unit ball, with data specified on the unit sphere ${S^{d-1}}\subset\mathbb{R}^{d}$, using second-kind Volterra time-domain boundary integral equations. While finite difference methods require that the Courant number $\Delta t/(\Delta x)^{2}$ be $\mathcal{O}(1)$, we have shown that integral equation methods can be both explicit and unconditionally stable for any fixed final time $T$.

We have also studied the Robin problem on the half space in all dimensions and shown that stability of the forward Euler scheme follows if $\Delta t<\frac{\pi}{c^{2}\kappa^{2}}$, where $c=3-\sqrt{2}$ and $\kappa$ is the heat transfer coefficient. As pointed out in Remark 8, this bound is very close to the optimal bound where $c\approx 1.52041925043874$.

A critical element in the proof of unconditional stability of the forward Euler scheme is the pointwise non-positivity of the double layer heat kernel on the unit sphere ${S^{d-1}}$, a property which extends to any convex domain. Combining this with the elementary fact that a unit double-layer density generates a surface potential approaching -\mbox{\small\frac{1}{2}} enabled us to extend this stability result to arbitrary smooth convex domains, in the Dirichlet case and the $L^{\infty}$-norm.

A key ingredient in the Robin proofs was a bound on the smallest eigenvalue of real symmetric Toeplitz matrices via the convexity of the associated sequence. This may be of independent interest in signal processing applications. Another ingredient for the proofs was a tight rational function bound for the ratio of modified Bessel functions of the first kind with large positive real argument, which may be of interest in its own right. A detailed analysis combining these ingredients showed that in the Dirichlet disc ($d=2$), the density is bounded in norm by the data, uniformly in time, so long as $\Delta t\leq 1$.

While this paper is purely analytic, we note that the numerical experiments in [43] are consistent with the theory presented here. More detailed experiments will be reported in a forthcoming paper [45] that considers the full initial-boundary value problem including forcing terms.

Some other questions arise naturally from our study. First, for the Dirichlet problem on the unit ball in higher dimensions, one may ask whether the scheme is stable for all time given some mild constraint on $\Delta t$. Second, one may ask about the stability analysis of the Robin problem on the unit ball in all dimensions. Third, it is natural to inquire about the stability of other explicit time marching schemes such as Adams–Bashforth multistep methods or explicit Runge–Kutta methods. Fourth, it would be interesting to see if the convexity assumption could be relaxed, and stability proved for arbitrary, sufficiently smooth domains. Integral equation methods become difficult to analyze when the boundary of the domain is not at least $C^{1}.$ We are currently investigating these issues and will report our findings in the future.

Acknowledgments

We are grateful for a discussion with Marcus Webb of KU Leuven on the Fourier series approach. S. Jiang was supported by NSF under grant DMS-1720405 and by the Flatiron Institute, a division of the Simons Foundation. We are also grateful to the anonymous referees for a careful reading of our paper and many useful suggestions for improvement.

Appendix A Properties of the modified Bessel functions of the first kind

The modified Bessel function of the first kind $I_{\nu}(x)$ is defined by the formula [32, Chapter 10]

[TABLE]

It satisfies the recurrence relations [32, §10.29.2]

[TABLE]

When $\nu$ is fixed and $x\rightarrow\infty$ [32, §10.30.4],

[TABLE]

When $\nu$ is an integer $n$, $I_{n}$ admits the integral representation  [32, §10.32.3]

[TABLE]

The following results can be found in [46].

Lemma 10**.**

Let $W_{\nu}(x)=\frac{xI_{\nu}(x)}{I_{\nu+1}(x)}$ and $S_{p,\nu}=W^{2}_{\nu}(x)-2pW_{\nu}(x)-x^{2}$. Then $S_{\nu,\nu-1}(x)$ is monotonically decreasing from [math] to $-\infty$ on $(0,\infty)$ for $\nu>1/2$,

[TABLE]

and

[TABLE]

for \nu\geq\mbox{\small\frac{1}{2}}, with $x\in(0,\infty)$.

Lemma 11**.**

Let $n$ be a positive integer. Then

1. (a)

There is only one zero $r_{n}$ for the equation

[TABLE]

on $(0,+\infty)$. Furthermore,

[TABLE] 2. (b)

The function $e^{-x}I_{n}(x)$ increases monotonically on $[0,r_{n}]$ and decreases monotonically on $[r_{n},+\infty)$. 3. (c)

The maximum value of $e^{-x}I_{n}(x)$ on $[0,\infty)$ satisfies

[TABLE]

Proof.

1. (a)

Using the recurrence (101), we have

[TABLE]

Thus,

[TABLE]

When $\frac{I^{\prime}_{n}(x)}{I_{n}(x)}=1$, $S_{n,n-1}=-n^{2}$. By the monotonicity and the range of $S_{n,n-1}(x)$, $S_{n,n-1}$ takes the value $-n^{2}$ at only one point and we denote that point by $r_{n}$.

Substituting $x=r_{n}$ into (103) and (104) with $\frac{I^{\prime}_{n}(r_{n})}{I_{n}(r_{n})}=1$ and simplifying the resulting expressions, we obtain (105). 2. (b)

We have

[TABLE]

Using (103), it follows that $\frac{I^{\prime}_{n}(x)}{I_{n}(x)}>1$ for x<n^{2}-\mbox{\small\frac{1}{2}} and $\frac{I^{\prime}_{n}(x)}{I_{n}(x)}<1$ for $x>n^{2}+n$. Combing these facts with (a), we have $\frac{I^{\prime}_{n}(x)}{I_{n}(x)}>1$ for $x<r_{n}$ and $\frac{I^{\prime}_{n}(x)}{I_{n}(x)}<1$ for $x>r_{n}$. That is, $\frac{d}{dx}(e^{-x}I_{n}(x))>0$ for $x<r_{n}$ and $\frac{d}{dx}(e^{-x}I_{n}(x))<0$ for $x>r_{n}$, which completes the proof of (b). 3. (c)

By the identity §10.35.5 in [32], we have

[TABLE]

Section 10.37 of [32] states that for fixed $x>0$, $I_{\nu}(x)$ is positive and decreasing for $0<\nu<\infty$. Hence,

[TABLE]

which completes the proof.

∎

The following lemma about differential inequalities can be found in [13, Chapter III, §4]. See also [35].

Lemma 12** (Petrovitsch 1901).**

Suppose that $f(y,t)$ is continuous in an open domain $D$. Suppose further that $y$ is the solution to the Cauchy problem

[TABLE]

1. (a)

(Increasing $t$). Suppose that $u$ satisfies the inequalities

[TABLE]

Then

[TABLE]

The inequality in (108) is reversed if both inequalities in (107) are reversed. 2. (b)

(Decreasing $t$). Suppose that $u$ satisfies the inequalities

[TABLE]

Then

[TABLE]

The inequality in (110) is reversed if both inequalities in (109) are reversed.

Lemma 13**.**

Let

[TABLE]

Then $g_{0}(x)$ has a unique zero, denoted as $x^{\ast}$, on $(\frac{3}{4},\infty)$. Furthermore, $g_{0}(x)<0$ on $[\frac{3}{4},x^{\ast})$ and $g_{0}(x)>0$ on $(x^{\ast},\infty)$.

Proof.

Let $r_{\nu}(x)=\frac{I_{\nu}(x)}{I_{\nu+1}(x)}$. In particular,

[TABLE]

From §10.37 of [32], we know that $I_{\nu}(x)$ is positive and increasing on $(0,\infty)$ for fixed $\nu(\geq 0)$ and $I_{\nu}(x)$ is decreasing on $0<\nu<\infty$ for fixed $x$. Thus, $r_{\nu}(x)>1$ on $(0,\infty)$ for $\nu\geq 0$. Let

[TABLE]

Then it is clear that the sign of $g_{0}(x)$ is determined by comparing $r_{0}(x)$ with $l_{0}(x)$. First, $\lim_{x\rightarrow\frac{3}{4}^{+}}l_{0}(x)=+\infty$ and thus $l_{0}(x)>r_{0}(x)$ as $x\rightarrow\frac{3}{4}^{+}$. Second, the series expansion of $l_{0}(x)$ and the asymptotic expansion of $r_{0}(x)$ are as follows:

[TABLE]

[TABLE]

Hence, $r_{0}(x)>l_{0}(x)$ as $x\rightarrow\infty$. Combining these two facts, there is at least one point $x^{\ast}\in(\frac{3}{4},\infty)$ where $r_{0}(x^{\ast})=l_{0}(x^{\ast})$. Or equivalently,

[TABLE]

By the recurrence relations (101), $r_{\nu}$ satisfies the following Riccati equation

[TABLE]

In particular, for $\nu=0$,

[TABLE]

We now calculate

[TABLE]

By Lemma 12, we have

[TABLE]

Equivalently,

[TABLE]

completing the proof. ∎

Remark 11**.**

Numerical computation shows that $x^{\ast}\approxeq 1.452165365078841\ldots$.

Corollary 3**.**

Let

[TABLE]

where $g_{0}(x)$ is defined in (111). Then $h_{0}(x)\geq 0$ on $[\frac{3}{4},x^{\ast}]$ and $[2,\infty)$; $h_{0}(x)\leq 0$ on $[x^{\ast},2)$.

Lemma 14**.**

Let

[TABLE]

Then $h_{1}(x)>0$ on $[\frac{3}{4},\infty)$.

Proof.

Let $x_{1}^{\ast}=\frac{\sqrt{33}+9}{8}=1.843\ldots$ be the larger root of $4x^{2}-9x+3$. Then $4x^{2}-9x+3>0$ for $x>x_{1}^{\ast}$ and $4x^{2}-9x+3<0$ for $x\in[\frac{3}{4},x_{1}^{\ast})$. We break $[\frac{3}{4},\infty)$ into several subintervals and show the positivity of $h_{1}(x)$ on each subinterval.

1. (a)

$x\in[x_{1}^{\ast},\infty)$. Let

[TABLE]

Then

[TABLE]

which is greater than zero if $x>3$ and less than zero if $x<3$. At $x=3$, $u_{0}(3)=\frac{5}{4}=1.25$ and $r_{0}(3)=1.23459\ldots<1.25=u_{0}(3)$. Thus, Using Lemma 12 in the increasing direction we have $r_{0}(x)<u_{0}(x)$ on $[3,\infty)$; and using Lemma 12 in the decreasing direction, we still have $r_{0}(x)<u_{0}(x)$ on $[x_{1}^{\ast},3)$. Equivalently, $h_{1}(x)>0$ on $[x_{1}^{\ast},\infty)$. 2. (b)

$x\in[\frac{7}{4},x_{1}^{\ast}]$. On this subinterval, we have $4x^{2}-7x\geq 0$ and $-4x^{2}+9x-3\geq 0$. Hence, $h_{1}(x)>0$, since $I_{1}(x)$ and $I_{0}(x)$ are always positive on $[0,\infty)$. 3. (c)

$x\in[\frac{3}{4},\frac{7}{4}]$. By (114), we have $u^{\prime}_{0}(x)-(1+\frac{1}{x}u_{0}(x)-u^{2}_{0}(x))\leq 0$ on $[\frac{3}{4},\frac{7}{4}]$. Also, $u_{0}(\frac{3}{4})=2<r_{0}(\frac{3}{4})=2.8\ldots$. Using Lemma 12, we have $r_{0}(x)>u_{0}(x)$, or equivalently $h_{1}(x)>0$ on $[\frac{3}{4},\frac{7}{4}]$.

∎

Lemma 15**.**

Let $f_{0}(x)=e^{-\frac{1}{x}}I_{0}(\frac{1}{x})$, $f_{1}(x)=e^{-\frac{1}{x}}I_{1}(\frac{1}{x})$, $f(x)=f_{0}(x)+af_{1}(x)$ with $a=0.05$. Then $f^{\prime\prime\prime}(x)>0$ on $(0,\infty)$.

Proof.

Using the recurrence relation (101), we obtain

[TABLE]

where $h_{0}(x)$ is defined in (112). Similarly,

[TABLE]

where $h_{1}(x)$ is defined in (113). Thus, in order to show that $f^{\prime\prime\prime}(x)>0$ on $(0,\infty)$, we only need to show that $h_{0}(x)+ah_{1}(x)>0$ on $(0,\infty)$.

We break it into several steps.

1. (a)

$x\in[0,1/4]$. On this interval, $3-4x\geq 2$, $0\leq 1-4x\leq 1$, $2-x\geq 1.75$, thus $h_{0}(x)\geq 1.75(2I_{0}(x)-I_{1}(x))>1.75I_{0}(x)$. And $4x^{2}-7x\geq-1.5$, $4x^{2}-9x+3\leq 3$, thus $h_{1}(x)\geq-1.5I_{1}(x)-3I_{0}(x)>-4.5I_{0}(x)$. Combinging these results, we have

[TABLE] 2. (b)

$\frac{1}{4}\leq x\leq\frac{9-\sqrt{33}}{8}<0.5$. On this interval, $3-4x>1$, $4x-1\geq 0$, $2-x>1.5$, thus $h_{0}(x)>1.5I_{0}(x)$. And $4x^{2}-7x>-2.5$, $0\leq 4x^{2}-9x+3\leq 1$, thus $h_{1}(x)>-2.5I_{1}(x)-I_{0}(x)>-3.5I_{0}(x)$. Combining these results, we have

[TABLE] 3. (c)

$\frac{9-\sqrt{33}}{8}\leq x\leq 3/4$. On this interval, $3-4x\geq 0$, $4x-1>0.6$, $2-x>1$, thus $h_{0}(x)>0.6I_{1}(x)$. And $4x^{2}-7x\geq-3$, $-(4x^{2}-9x+3)\geq 0$, thus $h_{1}(x)\geq-3I_{1}(x)$. Combining these results, we have

[TABLE] 4. (d)

$x\in[\frac{3}{4},x^{\ast}]\cup[2,\infty)$. On these two subintervals, both $h_{0}(x)$ and $h_{1}(x)$ are positive by Corollary 3 and Lemma 14. Thus $h_{0}(x)+ah_{1}(x)>0$. 5. (e)

$x\in(x^{\ast},2)$. We calculate

[TABLE]

where $g_{0}(x)$ is defined in (111). By Lemma 13, $g_{0}(x)>0$ on $(x^{\ast},\infty)$. Thus, $h^{\prime}_{1}(x)<0$ on $(x^{\ast},2)$. This shows that $h_{1}(x)>h_{1}(2)\approx 0.901688$ on $(x^{\ast},2)$. On the other hand, it is straightforward to show that $g^{\prime}_{0}(x)>0$ and $g^{\prime\prime}_{0}(x)<0$ on $(x^{\ast},2)$. Hence, $h^{\prime\prime}_{0}(x)=g^{\prime\prime}_{0}(x)(x-2)+2g^{\prime}_{0}(x)>0$ on $(x^{\ast},2)$, indicating that $h_{0}(x)$ achieves its minimum at exactly one point. Numerical calculation shows that

[TABLE]

Hence,

[TABLE]

∎