On the Reduction in Accuracy of Finite Difference Schemes on Manifolds without Boundary

We investigate error bounds for numerical solutions of divergence structure linear elliptic PDEs on compact manifolds without boundary. Our focus is on a class of monotone finite difference approximations, which provide a strong form of stability that guarantees the existence of a bounded solution. In many settings including the Dirichlet problem, it is easy to show that the resulting solution error is proportional to the formal consistency error of the scheme. We make the surprising observation that this need not be true for PDEs posed on compact manifolds without boundary. By carefully constructing barrier functions, we prove that the solution error achieved by a scheme with consistency error 𝒪(h^α) is bounded by 𝒪(h^α/(d+1)) in dimension d. We also provide a specific example where this predicted convergence rate is observed numerically. Using these error bounds, we further design a family of provably convergent approximations to the solution gradient.


page 1

page 2

page 3

page 4


Error Bounds for a Least Squares Meshless Finite Difference Method on Closed Manifolds

We present an error bound for a least squares version of the kernel base...

Can 4th-order compact schemes exist for flux type BCs?

In this paper new innovative fourth order compact schemes for Robin and ...

A formal proof of the Lax equivalence theorem for finite difference schemes

The behavior of physical systems is typically modeled using differential...

Parametric Complexity Bounds for Approximating PDEs with Neural Networks

Recent empirical results show that deep networks can approximate solutio...

The effect of time discretization on the solution of parabolic PDEs with ANNs

We investigate the resolution of parabolic PDEs via Extreme Learning Mac...

On compact 4th order finite-difference schemes for the wave equation

We consider compact finite-difference schemes of the 4th approximation o...

A Convergent Quadrature Based Method For The Monge-Ampère Equation

We introduce an integral representation of the Monge-Ampère equation, wh...

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

Definition 1 (Elliptic equation).

The equation (5) is elliptic if for all , then


where denotes that is a positive definite matrix.

The specific focus of the present article is linear divergence structure operators of the form


which are defined on a compact manifold . These equations are elliptic if is a symmetric positive definite matrix.

Given sets , and local coordinates we can locally recast this as the following linear divergence structure operator in Euclidean space:



is the metric tensor 


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


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


Critically, the approximation scheme (10) needs to be consistent with the underlying PDE (1).

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 higher-order derivatives of the solution; indeed, this is assumption is typically needed for schemes with higher-order consistency error (). The schemes analyzed in this article are required to satisfy an additional monotonicity assumption, which limits the consistency error to at most second-order (). 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 well-posedness and stability of the approximation scheme (10).

Definition 4 (Monotonicity).

The approximation scheme is monotone if it is a non-decreasing 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 .


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


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


This allows us to re-express the PDE (11) at the point as an equivalent PDE on the local tangent plane. We define


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,

The problem of approximating the PDE operator (11) at a point is now reduced to the problem of approximating the operator (13) at the point in the local tangent plane. This immediately allows one to make use of any existing method for designing monotone approximation of PDEs in Euclidean space.

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


in Euclidean space. Generalization to non-divergence structure operators, many nonlinear elliptic equations, and first-order 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



This motivates us to seek a finite difference approximation of the form


The monotonicity condition requires that the coefficient of be non-positive 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 non-zero 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 higher-dimensions, consider the problem of approximating the one-dimensional elliptic divergence-structure operator


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 :


These are coupled to the monotonicity condition


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 one-dimensional torus :


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:

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

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

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

  4. 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 non-convergent 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 one-dimensional torus

where . Let be a consistent, monotone approximation of the Laplacian and let be a consistent approximation of the right-hand side (which is zero in this case). We would like to solve the discrete system


However, this does not enforce the additional uniqueness constraint . Adding this as an additional equation leads to an over-determined system. Instead, a natural approach is to replace the equation (21) at with this additional constraint. This leads to the system


As a specific implementation, we consider a wide-stencil 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




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 right-hand side. Nevertheless, the discrete solution obtained by solving the system (22) does not converge to the true solution of (20).

Figure 1. LABEL: Maximum error and LABEL: effective maximum truncation error in the solution of (22).

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


We notice that the resulting discrete solution satisfies the system


This is consistent at all grid points since the first-step 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 non-trivial higher-dimensional 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 .


Suppose, in addition, that we have a consistent, monotone discretization scheme


with truncation error on the exact solution given by

Let denote the solution error. We notice that satisfies the discrete system


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


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 well-posed PDE. For example, we may choose to be the solution of the homogeneous Dirichlet problem


Now we define the auxiliary grid functions

We notice that

Applying the discrete maximum principle, we obtain

Thus we find that


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 “one-point” Dirichlet problem

However, this is not a well-posed 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 one-point 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


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.

Figure 2. The maximum norm of the auxiliary function obtained from (33).

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 one-dimensional 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 one-dimensional 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.

Figure 3. Maximum error in the solution of (26) on .

2.4. Convergence of gradients

Finally, we recall the important and well-known 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.

Figure 4. Maximum error in a centered difference approximation of obtained from the solution of the scheme (25).

This non-convergence is perhaps unsurprising given that is a low-accuracy 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 high-frequency 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:

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

  2. The matrix is symmetric positive definite.

  3. 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:

  1. is linear in its final argument.

  2. is monotone.

  3. 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.

Figure 5. The construction of a small cap about on the manifold

We then define the modified scheme as follows:

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.

Lemma 15.

Under the assumptions of Hypothesis 10,13, the discrete scheme


has a unique solution that is bounded uniformly independent of for sufficiently small .

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


and then by invoking the discrete comparison principle to conclude that


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).

Under the assumptions of Hypotheses 10 and 13, let be the solution of (1),  (3). Then the discrete solution solving (35) satisfies


where is a constant independent of .

Remark 16.

This convergence rate can be extended to more general -dimensional manifolds as follows:


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) right-hand 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


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

Now we define the right-hand side function by


In particular, this is chosen to be on the order of the local truncation error of (34) throughout most of the domain, but is allowed to take on larger values in the small cap in order to ensure the solvability condition is satisfied. See Figure 6.

Figure 6. The construction of the function from a “side profile” parametrized by distance from the point

3.3.2. Properties of the barrier function equation

Next we verify several key properties of the right-hand side function , which will in turn be used to produce estimates on the barrier function .

Lemma 17 (Mean-zero).

For every sufficiently small , the function defined in (41) satisfies the solvability condition (4)


We can directly compute