1. Background
1.1. Elliptic PDEs on manifolds
The results of this article are motivated by the study of elliptic PDEs, possibly nonlinear, of the form
(5) 
Definition 1 (Elliptic equation).
The equation (5) is elliptic if for all , then
(6) 
where denotes that is a positive definite matrix.
The specific focus of the present article is linear divergence structure operators of the form
(7) 
which are defined on a compact manifold . These equations are elliptic if is a symmetric positive definite matrix.
1.2. Approximation of elliptic operators
To build approximation schemes for the PDE (1), we begin with a point cloud discretizing the underlying manifold and let
(9) 
denote the characteristic (geodesic) distance between discretization nodes. In particular, this guarantees that any ball of radius on the manifold will contain at least one discretization point.
In this manuscript, we will consider finite difference discretizations of the PDE (1) of the form
(10) 
Definition 2 (Consistency error).
We say that the approximation of the PDE operator has consistency error if for every smooth there exists a constant such that
for every sufficiently small .
Remark 3.
In this article, we assume solutions are . It is also possible to design approximation schemes that depend on higherorder derivatives of the solution; indeed, this is assumption is typically needed for schemes with higherorder consistency error (). The schemes analyzed in this article are required to satisfy an additional monotonicity assumption, which limits the consistency error to at most secondorder (). See [35, Theorem 4].
Another concept that has proved important in the numerical analysis of fully nonlinear elliptic equations is monotonicity [2]. At its essence, monotone schemes reflect at the discrete level the elliptic structure of the underlying PDE. This allows one to establish key properties of the discretization including a discrete comparison principle. Even in the linear setting, monotonicity can play an important role in establishing wellposedness and stability of the approximation scheme (10).
Definition 4 (Monotonicity).
The approximation scheme is monotone if it is a nondecreasing functions of its final two arguments.
Closely related to monotonicity is the concept of a proper scheme.
Definition 5.
The finite difference scheme is proper if it is an increasing function of its second argument.
We note that any consistent, monotone scheme can be perturbed to a proper scheme by defining
where as .
Monotone, proper schemes satisfy a strong form of the discrete comparison principle [35, Theorem 5].
Theorem 6 (Discrete comparison principle).
Let be a proper, monotone finite difference scheme and suppose that
for every . Then .
Finally, we make a continuity assumption on the scheme in order to guarantee the existence of a discrete solution.
Definition 7 (Continuity).
The scheme is continuous if it is continuous in its final two arguments.
Remark 8.
We recall that the domain of the first argument of is the discrete set . Thus it is not meaningful to speak about continuity with respect to the first argument.
Critically, continuous, monotone, and proper schemes always admit a unique solution [35, Theorem 8]. Moreover, under mild additional assumptions, it is easy to show that the solution can be bounded uniformly independent of .
Lemma 9 (Solution bounds).
Suppose the PDE (5) has a unique solution. Let be continuous, monotone, proper, and have consistency error . Suppose also that there exists a constant , independent of , such that for every ,
Then for every sufficiently small , the scheme (10) has a unique solution that is uniformly bounded independent of .
Proof.
Since is continuous, monotone, and proper, a solution exists by [35]. Let be the exact solution of (1). By consistency, we know that there exists a constant , that does not depend on , such that
Now let be any constant and compute
Thus by choosing , we find that
Then by the Discrete Comparison Principle 6
By an identical argument, we obtain
We conclude that
and thus is uniformly bounded. ∎
1.3. Tangent plane approximations
A variety of approaches are available for discretizing PDEs on manifolds. Particularly simple are methods that allow the surface PDE to be approximated using schemes designed for PDEs in Euclidean space [31, 32]. Here we overview one particular approach, which can easily be used to design monotone approximation schemes that satisfy the assumptions of our main result on error bounds (Theorem 1).
Consider the PDE operator
(11) 
at a particular point . We relate this to an equivalent PDE posed on the local tangent plane through a careful choice of local coordinates. In general, local coordinates will introduce distortions to the differential operators. However, this problem was avoided in [26, 27] with the use of geodesic normal coordinates, which preserve distance from the reference point
. In these coordinates the metric tensor is an identity matrix and the Christoffel symbols vanish at the point
.Given some neighbourhood of the point , we let denote geodesic normal coordinates. Because they are chosen to preserve distances from the point , they satisfy
where represents the geodesic distance along and the usual Euclidean distance on the tangent plane.
We can now introduce a local projection of onto the relevant tangent plane in a neighborhood of as follows
(12) 
This allows us to reexpress the PDE (11) at the point as an equivalent PDE on the local tangent plane. We define
(13) 
where now is the usual Euclidean differential operator. Because the particular choice of coordinates does not introduce distortions, the PDE operator will preserve its original form. In particular,
1.4. Monotone approximation schemes
There is a growing body of literature on the design of monotone finite difference approximations for PDEs (both linear and nonlinear) in Euclidean space [3, 4, 5, 6, 9, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 30, 34, 35, 36]. Many of these are specifically constructed on Cartesian grids; however, recent work on generalized or meshfree finite difference methods demonstrate how these can be adapted to unstructured grids [16, 37].
We briefly review the procedure for designing monotone generalized finite difference methods for approximating linear divergence structure operators of the form
(14) 
in Euclidean space. Generalization to nondivergence structure operators, many nonlinear elliptic equations, and firstorder operators is straightforward [16, 27, 37].
Consider the problem of approximating (14) at a point . It is natural to want to use values of at the “nearest neighbors” to accomplish this. Surprisingly, though, given any fixed stencil width, it is always possible to find a linear elliptic PDE operator that does not admit a consistent, monotone discretization on that stencil [28, 33]. For general degenerate elliptic operators, it is sometimes necessary to allow the stencil to grow wider as the grid is refined in order to achieve both consistency and monotonicity.
We attempt to discretize (14) at using points within some search neighborhood. To this end, we associate a search radius to the point cloud and require as . We define by
the set of neighbors to the reference point . A common choice is , which ensures that as the point cloud is refined (), the number of neighboring points grows without bound. For uniformly elliptic operators, a choice of may be sufficient.
We notice that the PDE operator can be written in the form
(15) 
where
This motivates us to seek a finite difference approximation of the form
(16) 
The monotonicity condition requires that the coefficient of be nonpositive for each . This leads to the set of linear inequality constraints
To achieve consistency, we Taylor expand the terms in (16) about the reference point . We then compare the coefficients of each term with the desired operator (15), which leads to a system of linear equations that must be satisfied by the coefficients .
In typical implementations, one possibility is to exploit the structure of the underlying PDE to set many of the coefficients to zero a priori and obtain closed form expressions for the (small) number of nonzero coefficients [16]. Another option is to use simple analytical and computational optimization tools to numerically determine the values of the coefficients and establish bounds needed to ensure consistency [25, 37].
Example. As an illustrative example, which easily generalizes to higherdimensions, consider the problem of approximating the onedimensional elliptic divergencestructure operator
(17) 
at the point .
Given a set of discretization points
we define . Note that for , we have . Then we seek an approximation of the form
We first Taylor expand about :
Next we multiply this out and Taylor expand the resulting products about the point :
Finally, we compare this with the desired operator (17) to obtain the following system of linear equations for the coefficients :
(18) 
These are coupled to the monotonicity condition
(19) 
The constrained linear system (18)(19) can be analyzed and solved numerically using optimization techniques as in [25, 37] or solved exactly as in [16].
Consider the special case of equally spaced nodes for which
The consistency and monotonicity conditions lead to an underdetermined system. One particular solution is
with otherwise. This leads to the approximation
which has a simple interpretation in terms of standard centered differences.
2. Empirical Convergence Rates in One Dimension
This section will consider the very simple example of Laplace’s equation on the onedimensional torus :
(20) 
which has the trivial solution .
We use this toy problem to demonstrate several surprising properties of consistent and monotone approximations on compact manifolds, which motivate and validate the main results presented in the remainder of this article. In particular, we observe that:

Consistent, monotone, proper schemes need not converge to the true solution unless the solvability condition (4) is carefully taken into account.

Typical approaches for proving convergence rates for linear elliptic PDEs with Dirichlet boundary conditions fail on compact manifolds.

Actual error bounds achieved by convergent schemes can be asymptotically worse than the truncation error of the finite difference approximation.

A simple consistent scheme for the gradient need not produce a convergent approximation of the gradient when applied to a numerically obtained solution.
2.1. A nonconvergent scheme
We begin by describing a natural “textbook” approach to attempting to solve (20) numerically, which does not lead to a convergent scheme.
Consider the uniform discretization of the onedimensional torus
where . Let be a consistent, monotone approximation of the Laplacian and let be a consistent approximation of the righthand side (which is zero in this case). We would like to solve the discrete system
(21) 
However, this does not enforce the additional uniqueness constraint . Adding this as an additional equation leads to an overdetermined system. Instead, a natural approach is to replace the equation (21) at with this additional constraint. This leads to the system
(22) 
As a specific implementation, we consider a widestencil approximation of the Laplacian, which mimics the type of scheme that is often necessary for monotonicity in higher dimension [28, 33]. We also make the scheme proper, which ensures that the system (22) has a unique solution [35].
Let be a perfect square (where ). We will build schemes with stencil width . Define
(23) 
and
(24) 
The resulting approximation (21) is consistent with (20), monotone, and proper. The use of wide stencils degrades the truncation error of the usual centered scheme from to , which is of the same order as the consistency error introduced by the proper term and the approximation of the righthand side. Nevertheless, the discrete solution obtained by solving the system (22) does not converge to the true solution of (20).
An issue that arises in this approach is that even though and are consistent with the original equation, they are not designed in a way that attempts to mimic the solvability condition (4) at the discrete level. As a result, all the work of imposing this compatibility condition must be made up for at the single point where no approximation of the Laplacian is explicitly enforced in (22). This is evident in Figure 1, which plots the value of (the “effective” truncation error of the scheme). This does not converge to zero as the grid is refined. In other words, the failure to incorporate the solvability condition at the discrete level has led to a scheme that is effectively inconsistent.
Enforcing a solvability condition at the discrete level is not straightforward: the discrete condition may not be known explicitly, and in many problems even the continuous solvability condition is not known explicitly [24].
The solution we propose and analyze in this work is to automatically “spread out” the effects of the solvability condition by first solving a discrete system that is consistent with Laplace’s equation at every grid point, then enforcing the uniqueness constraint in a second step. The resulting procedure is
(25) 
We notice that the resulting discrete solution satisfies the system
(26) 
This is consistent at all grid points since the firststep solution is uniformly bounded (Lemma 9). Moreover, the resulting solution automatically satisfies the uniqueness condition by construction.
2.2. The Dirichlet Problem
We are interested in establishing error bounds for solutions of (26) (and, of course, generalizations to nontrivial higherdimensional problems). To gain intuition and inspiration, we first review a standard approach to establishing error bounds for monotone schemes approximating the Dirichlet problem.
Consider as an example Poisson’s equation with Dirichlet boundary conditions on a domain .
(27) 
Suppose, in addition, that we have a consistent, monotone discretization scheme
(28) 
with truncation error on the exact solution given by
Let denote the solution error. We notice that satisfies the discrete system
(29) 
If the discrete linear operator and the underlying grid are sufficiently structured, we may be able to explicitly determine its eigenvectors and eigenvalues. In this case, we immediately obtain error bounds via
(30) 
If the discrete problem does not have a simple enough structure, we can instead choose some bounded such that
This can always be accomplished for a consistent approximation of a wellposed PDE. For example, we may choose to be the solution of the homogeneous Dirichlet problem
(31) 
Now we define the auxiliary grid functions
We notice that
Applying the discrete maximum principle, we obtain
Thus we find that
(32) 
In other words, the solution error is proportional to the truncation error of the underlying approximation scheme.
2.3. Error bounds on the 1D torus
It is natural to try to adapt the techniques used for the Dirichlet problem to error bounds for PDEs on manifolds without boundary. Indeed, we may attempt to interpret (20) as the “onepoint” Dirichlet problem
However, this is not a wellposed PDE and attempting to solve an analog of (31) for the auxiliary function will not lead to a function that is smooth on the torus.
We might attempt to carry this argument through at the discrete level, noticing that the solution of (25) does satisfy the following discrete version of a onepoint Dirichlet problem
The resulting discrete linear system involves a strictly diagonally dominant matrix. However, standard bounds on the inverse of such a matrix [10]
yield the estimate
which cannot provide any convergence guarantees when substituted into (30). The degradation of this bound as is due to the fact that, the the scheme (26) is proper, it is not uniformly proper as .
Instead, we attempt to utilize the techniques outlined above, which requires us to construct a function satisfying the system
(33) 
As this is a proper scheme, it does admit a unique solution. However, the numerically obtained solution is not uniformly bounded as the grid is refined (Figure 2) and the resulting estimate in (32) does not provide a useful error bound.
The approach we will use in section 3 to obtain error bounds involves effectively expanding the “Dirichlet” condition onto a larger set, which shrinks to a point as the grid is refined. A downside to this approach is that it degrades the error bounds from the size of the truncation error to the asymptotically worse rate of . Surprisingly, though, our simple onedimensional example indicates that this may be the best we can hope for.
Consider again the discrete solution obtained by solving (26), which has a truncation error of at all grid points on the onedimensional torus and exactly satisfies the uniqueness condition . We solve this system numerically and present the error in Figure 3. This example does display numerical convergence to the true solution. However, the observed accuracy is only , precisely what is predicted by Theorem 1.
2.4. Convergence of gradients
Finally, we recall the important and wellknown fact in numerical analysis that pointwise convergence of an approximation does not imply convergence of gradients. In particular, if we consider the approximation obtained by solving (25), we might try to obtain information about the solution derivative by using the standard centered difference scheme
However, this fails to converge to the true solution derivative as the grid is refined; see Figure 4.
This nonconvergence is perhaps unsurprising given that is a lowaccuracy approximation to . Indeed, a closer look at the centered difference scheme reveals that
The theoretical error of this approximation is potentially as large as , which is unbounded as .
Nevertheless, the numerical solution does still contain information about the true solution derivative. In order to obtain this, we will require approximations of the gradient that utilize sufficiently wide stencils to overcome potential highfrequency components in the solution error. This has the effect of making the size of the denominator in the finite difference approximation larger than the solution error, which leads to a convergent approximation as . This idea will be developed in section 4.
3. Convergence Rate Bounds
We now establish error bounds for a class of consistent, monotone approximations schemes for (1), (3). The main result is presented in Theorem 1. The approach we take here is to construct barrier functions, which are shown to bound the error via the discrete comparison principle. Importantly, the error estimates we obtain are consistent with the empirical convergence rates observed in section 2.
3.1. Hypotheses on Geometry and PDE
We begin with the hypotheses on the geometry and PDE (1) that are required by our convergence result.
Hypothesis 10 (Conditions on PDE and manifold).
The Riemannian manifold and PDE (1) satisfy:

The manifold is a 2D compact and connected orientable surface without boundary.

The matrix is symmetric positive definite.

The function satisfies .
Remark 11.
The compactness of the 2D manifold implies that it is geodesically complete, has injectivity radius strictly bounded away from zero, and that the sectional curvature (equivalent to the Gaussian curvature in 2D) is bounded from above and below [29].
Remark 12.
The fact that is dimensional can quite easily be generalized, since nothing in our convergence proof fundamentally depends upon the dimension (though the convergence bound does vary with ). We emphasize, however, that the lack of boundary on our manifold is an essential point in this paper.
3.2. Approximation Scheme
Next, we describe the class of approximation schemes that are covered by our convergence result. The starting point of the scheme is the idea that the uniqueness constraint (3) should be posed at the point , with a reasonable discrete approximation of the PDE posed on other grid points. However, as discussed in section 2, this approach may not yield a convergent scheme. Instead, we will create a small cap around and fix the values of at all points in this cap.
To construct an appropriate scheme, we begin with any finite difference approximation of the PDE operator (2) that is defined for and that satisfies the following hypotheses.
Hypothesis 13 (Conditions on discretization scheme).
We require the scheme to satisfy the following conditions:

is linear in its final argument.

is monotone.

There exist constants such that for every smooth the consistency error is bounded by
Next we define some regions in the manifold that will be used to create “caps” where is fixed in this scheme, and where additional conditions will be posed on barrier functions. Choose any . Define the regions
See Figure 5.
We then define the modified scheme as follows:
(34) 
Remark 14.
The condition can be relaxed provided the resulting discrete solution has a uniformly bounded Lipschitz constant in this region and the values of are close to zero. Pinning the value to zero has the particularly strong effect of setting the local Lipschitz constant to zero.
Note that the discretization is automatically proper by construction. Therefore, this scheme has a uniformly bounded solution by Lemma 9.
3.3. Convergence Rates
The idea in this section is to establish the convergence of the discrete solution of a monotone (and proper) scheme to the unique solution of the underlying PDE. We accomplish this by constructing a barrier function such that
(36) 
and then by invoking the discrete comparison principle to conclude that
(37) 
The barrier function can be chosen to satisfy . In this article, we explicitly treat the case . In Section 2, we saw for that the empirical convergence rate was , which is consistent with our theoretical error bound when . The factor appears because there is a contribution of from the dimension of the underlying manifold (which arises due to the solvability condition (4)), and a contribution of from deriving a Lipschitz bound (also constrained by the solvability condition). Thus, we see that it is the solvability condition on the manifold without boundary that leads to the reduced convergence rate overall of a monotone and proper discretization.
We state the main convergence result:
Theorem 1 (Convergence Rate Bounds).
Remark 16.
This convergence rate can be extended to more general dimensional manifolds as follows:
(39) 
3.3.1. Construction of barrier functions
We now define the barrier functions by solving a linear PDE on the manifold with an appropriately chosen (small) righthand side that satisfies the solvability condition (4). In particular, given a fixed (which will be determined later), we let be the solution of the PDE
(40) 
We emphasize that while the barrier function depends on the grid parameter , it is the solution of the PDE on the continuous level.
Now we outline the construction of an appropriate function ; see Figures 5 and 6 for two complementary visualizations of the resulting function . Let be a fixed constant, to be determined later. We let denote the volume of a set and note that in two dimensions,
We define the following real numbers
We record the fact that and for some . Finally, we introduce a smooth cutoff function
3.3.2. Properties of the barrier function equation
Next we verify several key properties of the righthand side function , which will in turn be used to produce estimates on the barrier function .
Lemma 17 (Meanzero).
Proof.
We can directly compute