Convergence of Chandrashekar’s Second-Derivative Finite-Volume Approximation

We consider a slightly modified local finite-volume approximation of the Laplacian operator originally proposed by Chandrashekar (Int J Adv Eng Sci Appl Math 8(3):174–193, 2016, https://doi.org/10.1007/s12572-015-0160-z). The goal is to prove consistency and convergence of the approximation on unstructured grids. Consequently, we propose a semi-discrete scheme for the heat equation augmented with Dirichlet, Neumann and Robin boundary conditions. By deriving a priori estimates for the numerical solution, we prove that it converges weakly, and subsequently strongly, to a weak solution of the original problem. A numerical simulation demonstrates that the scheme converges with a second-order rate.

show that injected Dirichlet boundary conditions also yield a stable scheme. That is, both approaches are applicable, and we have chosen the strong imposition to provide an alternative.
The paper is further organised as follows. Section 2 defines the problem, whereas a priori estimates for the continuous solution are found in Sect. 3. In Sect. 4 we state the weak formulation of the problem. Section 5 concerns the spatial discretisation and provides the proof of the slightly altered Laplacian operator being SBP. In Sect. 6, the numerical scheme that approximates our problem is stated. Furthermore, the SBP properties of the Laplacian operator are utilised to obtain discrete a priori estimates similar to those found for the continuous solution. Using these estimates, we show in Sect. 7 that the approximate solution obtained by the proposed numerical scheme converges weakly to a weak solution of the original problem. Furthermore, we show in Sect. 8 that the solution indeed converges strongly by using Aubin-Lions' lemma. The solution is subsequently shown to be unique in Sect. 9. Finally, Sect. 10 provides a numerical example that demonstrates the convergence of the scheme.
Remark 2. 1 We could have made all boundary conditions homogeneous by defining w differently. However, we choose non-zero Neumann and Robin data to keep the regularity assumptions on the boundary data to a minimum.

A Priori Estimates for the Continuous Problem
To obtain a priori estimates on u, we use the energy method (see e.g. [17]). By inserting F given in (4) into (3a) and integrating by parts, we obtain Using Cauchy-Schwarz's and Young's inequality on the first integral on the right-hand side, we obtain 1 2 reads Inserting the boundary conditions for w and u given in (2) and (3b)-(3d), respectively, we obtain 1 2 Consider the underlined boundary terms above. We follow [19], and bound these terms by first using the Cauchy-Schwarz inequality: and then by using the trace theorem, which states that u L 2 (∂ ) ≤ C u H 1 ( ) , C > 0 (see e.g. [1]): Here, we have introduced the notation a b for a ≤ Cb, where C > 0 is a constant. By employing Young's inequality, the boundary terms (7) finally read The preliminary estimate (6), can then be stated as The last term on the right-hand side is negative semi-definite, since α ≥ 0. We neglect it in the remaining analysis. Hence we have Consider the three last terms on the left-hand side of the above inequality. By adding and subtracting μ 2 u 2 L 2 ( ) , they can be rewritten as By choosing β sufficiently small μ−β Employing Grönwall's inequality (see e.g. [10]), we obtain u(·, ·, t) 2 The inequality holds for all 0 ≤ t ≤ T . In Sect. 2 we defined w ∈ L 2 (0, T ; H 1 ( )), w t ∈ L 2 (0, T ; H 1 ( )) and g N ,R ∈ L 2 (0, T ; L 2 (∂ N ,R )). Using this, we find that the right-hand side of the above inequality is bounded. Thus, u ∈ L ∞ (0, T ; L 2 ( )). Lastly, we integrate (8) in time to obtain where we have used u| t=0 ≡ 0 in the last step. Since T 0 u(·, ·, t) 2 dt ≤ constant, we observe from this inequality that ∇u ∈ L 2 (0, T ; L 2 ( )), and thus we have u ∈ L 2 (0, T ; H 1 ( )).

Weak Formulation of the Heat Equation
Next, we derive the weak formulation of (3). Let Integrating by parts and inserting the boundary conditions given in (3b)-(3d), give where we have used φ| ∂ D = 0. Using φ| t=T = 0, u| t=0 = 0 and partially integrating the left-hand side in time further yields the weak form of (3): where F given by (4) satisfies Remark 4.1 Since the forcing function is not the main focus of this work, we use φ F dx as short-hand notation for (12) and make comments about it where necessary. (12), we see that w ∈ H 1 ( ) is sufficient to bound the two first integrals on the right-hand side. Furthermore, the regularity of w t is determined by the regularity of the boundary data (see e.g. [17]). Thus, for γ (w t ) = g D t to be satisfied (where γ is the trace function), we must have w t ∈ H 1 ( ), and that is why we assumed that g D ∈ H 1 (0, T ; H 1/2 (∂ D )) in Sect. 2.

Spatial Discretisation
Let¯ h be a discretisation of¯ = ∪ ∂ into non-overlapping triangles K n , n = 1, . . . , N such that¯ h = ∪ N n=1 K n , and such that there are no hanging nodes in¯ h . The grid functions are defined on the vertices of the triangles. Furthermore, subdivide¯ h into a dual grid consisting of dual cells, V i , i = 1, . . . , I , such that¯ h = ∪ I i=1 V i . The dual cells are polygons surrounding a vertex, i. A dual-volume boundary consists of straight lines drawn from the mid-point of an edge adjacent to grid point i to the centroid of the triangles adjacent to the grid point (see Fig. 1). (These are the dual volumes of the standard node-centred finite-volume method, see e.g. [22]). We introduce the notation   [7] for the interior scheme. For triangles having at least one edge along the Dirichlet boundary, the Dirichlet condition was incorporated weakly in the gradient operator in [7]. Here, we use the approximation for interior triangles for every triangle in the grid. The approximation is given by where |K n | is the area of triangle K n ; i, j, k are the vertices of the triangle, andn n i, j,k are the outward pointing normal vectors of the triangle, opposite of the particular node (see Fig. 2). The length of the normal vectors,n n i, j,k , is equal to the length of the adjacent edge. Next, we introduce the following notation. Then the approximation of the Laplacian on a dual volume is found by approximating (see [7]) by Here,b(e) denotes the outward pointing normal vector at boundary edge e (see Fig. 3a). The superscript n(e) signifies the triangle that has the boundary edge e. The components of the approximation (15) is depicted in Fig. 3b. a b Fig. 3 a Example of a vertex i belonging to three triangles (K 1 , K 4 , K 7 ) where two of them (K 1 , K 7 ) have an edge along the boundary, depicted with the corresponding boundary normalsb(e). b Example of a dual cell, V i , and the components of the Laplace approximation (15) The approximation of the Laplacian (15) with Dirichlet boundary conditions taken into account, was demonstrated to satisfy the Summation-by-Parts (SBP) property in Theorem 1 in [7]. Here, we state the analogous result without any boundary conditions.

Theorem 5.1 Let u h and v h be two grid functions defined on¯ h such that u h
. . , u I ), and correspondingly for v h . Then the discrete approximation of the Laplacian operator (15) where the subscripts {i, 1} and {i, 2} indicate the two edges adjacent to the boundary node i.
Proof Multiply Eq. (15) by v i V i and sum over all vertices in the grid.
In the second equality, we have used that the set E i is empty for interior nodes. For the first term in the above equation, we change the order of summation and move ∇ h u n outside the summation over the vertices of a triangle K n in (16), to obtain For the boundary nodes, we have With (17) and (18), (16) can be written as In the last equality we have used the approximation of the gradient (13).

The Numerical Scheme and Discrete A Priori Estimates
To approximate the problem (1) we use (15) for the Laplacian approximation at the interior nodes. The Dirichlet condition is imposed strongly by injection (see e.g. [15,16]). The Neumann and Robin conditions are imposed weakly in the same way as in [7]. That is, by replacing the last term of (14) with the boundary data, we approximate the Neumann and Robin boundaries by: Remark 6.1 Imposing the Dirichlet condition by injection means in practice that the Dirichlet nodes are overwritten by the exact boundary data after each time step. (Equivalently, no equation is solved at these nodes, since u is equal to the boundary data.)

Remark 6.2
A boundary node is either of Dirichlet, Neumann or Robin type. The entire dual-cell boundary coinciding with the physical boundary is subsequently approximated as the same type as the boundary node, see Fig. 4. This means that in the junction between two boundary types, part of the computational boundary may be approximated as something different than the actual physical boundary. However, this is an O(h) error of the boundary integral which tends to zero with decreasing mesh sizes. Note that this is only necessary for the Dirichlet nodes where the boundary conditions are injected. For Neumann and Robin nodes, we could split the outer dual-cell boundary into a Neumann and a Robin part since these boundary conditions are set weakly. However, in the scheme (19) below, we use the first approach to reduce notation.
The above choices lead to the following discrete approximation scheme of (1) For boundary nodes, the whole dual cell boundary is approximated as the type of the boundary node Remark 6.3 For readers familiar with the simultaneous approximation term (SAT) (see e.g. the review papers [8,30]) we remark that the schemes for the Neumann and Robin nodes are equivalent to where the SATs take the form Remark 6.4 To simplify following energy analysis, we have defined the Dirichlet nodes in (19) as We emphasise that when implementing the scheme, the Dirichlet nodes should take the form v i = g D i in order to avoid discretisation errors from the time-stepping algorithm.
As for the continuous problem, we transform the scheme (19) into one that imposes homogeneous Dirichlet boundary conditions. That is, we construct a function w as defined in Sect. 2 and introduce u = v − w (see again [1,17]). Inserting v = u + w into the scheme (19), we obtain where Remark 6.5 By the Picard-Lindelöf theorem (see e.g. [25]), the ordinary differential equation (21) has a solution if the scheme is stable.
To obtain a priori estimates for the approximate solution u h = (u 1 , u 2 , . . . , u I ), we use the discrete energy method (see e.g. [17] for more details on the energy method). That is, we multiply the scheme (21) in each node, i, by u i V i and sum over all grid points.
, and all the sets are disjoint, and since the scheme for the Dirichlet nodes is zero, the underlined terms amount to summing over all nodes in the grid. That is, the above is equivalent to Using Theorem 5.1, we obtain where we in the last inequality have used that i∈∂ R h − 1 2 αu 2 i (|b i,1 (e)| + |b i,2 (e)|) ≤ 0 since α ≥ 0. We can further manipulate the Neumann boundary terms as follows Using Young's inequality, we obtain The Robin boundary terms can be manipulated the same way. Thus, (23) reads We introduce the following notation for the discrete equivalents of the L 2 -norms: Using the definitions (25)-(27), we can recast (24) as To obtain an estimate analogous to (8), we must consider the forcing term i∈¯ V h u i V i F i . Except for the time-derivative term in (22), F i takes the same form as the right-hand side of the scheme (21). By using the SBP property from Theorem 5.1 and Young's inequality, we obtain Thus, we have Similarly as for the continuous problem, if we choose = 1, we obtain −μ ∇ h u h 2 Finally using the trace theorem, we arrive at a similar estimate as in (8): Note that we have arrived at a semi-discrete equivalent of (8). Thus, by using Grönnwall's inequality followed by integration in time, as done in Sect. 3, we obtain u h ∈ L ∞ (0, T ; L 2 V ( )) and ∇ h u h ∈ L 2 (0, T ; L 2 K ( )). We may extend the numerical solution, u h , to the entire domain by a linear interpolation on the triangles. Let u h c denote this continuous piecewise linear function. We have that ∇ h u h c = ∇u h c = ∇ h u h . Hence ∇u h c ∈ L 2 (0, T ; L 2 K ( )) (and also, ∇u h c ∈ L 2 (0, T ; L 2 ( )) since ∇u h c is piecewise constant). Furthermore, the norm u h

Weak Convergence to a Weak Solution
(¯ )). Since φ is smooth (in space), φ| V i , which is the restriction of φ to a dual cell, can be written as where h is a characteristic mesh size and p i (x i , y i , t) is a function of size O(1). The gradient approximation (13) is . This can easily be checked for equilateral triangles. Thereafter, one can prove the relation for a general triangle by transforming it to an equilateral one using a linear transformation. We denote the right-hand side of the scheme (21) by L h u h . To prove convergence to a weak solution, we test the numerical scheme (21) against φ. That is, we calculate We now use that φ| V i = φ i + hp i to obtain where we have used u h t = L h u h in the last step. Thus Inserting the specific form of L h u h (that is, the right-hand side of the scheme (21)) yields Since φ i is constant on each dual cell, V i and the Laplacian approximation is a scalar constant, the right-hand side above can be integrated exactly, leading to (34) As in the discrete analysis in Sect. 6, the underlined terms can be written as the sum over all grid points in¯ h as follows Using the SBP properties from Theorem 5.1 yields (36) Since ∇ h φ n and ∇ h u n are constant on each triangle K , we have that Thus, (36) can be written as

Partial integration in time yields
where we have used u h | t=0 = 0 and φ h | t=T = 0.

Remark 7.1
Here, φ h F h dx is the short-hand for the semi-discrete form of (12). By using the SBP property from Theorem 5.1, it can be written as We keep the symbolic expression to reduce notation. Since We utilise the following functional analysis theorem (see e.g. [10], and [5] for a proof). Here, we take T = × [0, T ]. Consider the O(1) term on the left-hand side of (39). Using Theorem 7.2, we have that The other O(1) terms can be treated in a similar way. Turning to the second term in (39), we have  (38)) are known and bounded in L 2 ( T ) (see the assumptions in Sect. 2), the weak convergence of the symbolic expression (38) follows trivially.

Remark 7.4
Note that the boundary integrals over the computational boundaries converge to the integrals over the physical boundaries as h → 0. That is, as h → 0.

Strong Convergence to a Weak Solution
Next, we prove strong convergence to the weak solution.
We also need the following definition.
As we have assumed the spatial domain to be Lipschitz, the following result applies. Theorem 8.3 (see e.g. [1] or [27]) Any Lipschitz domain has the k-extension property.
For a bounded domain, , with the k-extension property, we have that H 1 ( ) is compactly embedded in L 2 ( ) (see e.g. [25]), which in turn is continuously embedded in H −1 ( ). To prove strong convergence, we need the Aubin-Lions Lemma: A Banach space X is compactly embedded in another Banach space Y , if the following two conditions hold (see [10]): Herein, we use to establish the strong convergence. That is, we need to show that the norm (see e.g. [25]) is bounded. To this end, we test the scheme (21) with a function φ ∈ C ∞ 0 (¯ ).
Note the resemblance to (29) (the only difference being the function φ that is now vanishing on the whole boundary ∂ ). From derivations analogous to (30)-(35), we can recast the above equation to By using φ| ∂ = 0 and the SBP property (see Theorem. 5.1) we have Since ∇ h u h ∈ L 2 (0, T ; L 2 K (¯ h )) and all terms of F h are properly bounded (see the assumptions in Sect. 2), letting h → 0 yields as lim h→0 (φ − hp) = φ. By inserting the specific form of φF dx and using the Cauchy-Schwarz inequality, we obtain This holds for all φ ∈ C ∞ 0 (¯ ), and by density, it follows that the inequality holds for all φ ∈ H 1 0 ( ). Integration in time finally yields Hence u h t ∈ L 1 (0, T ; H −1 ( )), and since (u h c ) t is u h t extended to the entire domain using a linear interpolant on the triangles, we also have (u h c ) t ∈ L 1 (0, T ; H −1 ( )). Thus, by Aubin-Lions' lemma 8.4, the family of functions, U = {u h c ∈ L 2 (0, T ; H 1 ( )) | (u h c ) t ∈ L 1 (0, T ; H −1 ( ))}, is compactly embedded in L 2 (0, T ; L 2 ( )), meaning that u h c converges strongly to the weak solution.

Uniqueness of the Weak Solution
Assume that there are two weak solutions u, v to the problem (1) satisfying the boundary and initial data. Then w = u − v is also a weak solution with homogenous data (F = g D = g N = g R ≡ 0). Take φ = w in (10) to obtain Integrating the right-hand side by parts, and using the fact that the boundary data is zero, we obtain Hence, w L 2 (0,T ;L 2 ( )) = u − v L 2 (0,T ;L 2 ( )) = 0 and thus the weak solution is unique in L 2 (0, T ; L 2 ( )).

Numerical Simulations
We implement the scheme (1) and consider the manufactured solution used in [7]. That is, the exact solution is given by u(x, y, t) = e −8π 2 t sin(2π x) sin(2π y) + e −32π 2 t sin(4π x) sin(4π y), which yields a zero forcing function. Furthermore, we let μ = α = 1. We consider a square domain = [0, 1] × [0, 1] containing a hole. The hole is located at (x, y) = (0.5, 0.5), and has radius r = 1 8 . We pose Dirichlet boundary conditions on the boundary of the hole, Neumann boundary conditions on y = 0, y = 1 and Robin boundary conditions on x = 0, x = 1. The boundary data is given by (44). t = 0.05 is used as the final time. The scheme was run on grids containing 398, 1394, 5097, 19457 and 76166 nodes. A typical grid is depicted in Fig. 5a. All grids were generated using Gmsh (see [14]). The scheme was implemented using the Julia programming language (see [4]). Figure 5b shows the convergence rate together with a reference line representing secondorder convergence. We conclude that the scheme converges at approximately a rate of two.

Conclusion
Herein, we have considered a slightly modified local finite-volume approximation of the Laplacian operator proposed by Chandrashekar in [7] for discretising the heat equation in a b Fig. 5 a A typical mesh. b Convergence rate obtained for simulations using N = 398, 1394, 5097, 19457, 76166 grid points two spatial dimensions on general triangular grids. The equation was augmented with Dirichlet, Neumann and Robin boundary conditions. The Dirichlet boundary condition was imposed strongly by injection, while the Neumann and Robin conditions were imposed weakly. We demonstrated that this modification satisfies the SBP property proved in [7]. By using the energy method, a priori estimates for the numerical solution were derived. From these estimates, we were able to prove the weak convergence of the numerical solution to a weak solution of the heat equation. Thus, consistency, in a weak sense, of the Laplacian operator was established. Subsequently, we demonstrated that the numerical solution converges strongly to a weak solution by using Aubin-Lions' lemma. Finally, the weak solution was shown to be unique. To the best of our knowledge, this is the first proof of convergence for a local finite-volume method for the Laplacian on general triangular grids. The theory presented here is straightforwardly applicable to three spatial dimensions, provided that the Laplacian approximation can be generalised to such domains. A numerical simulation, which included Dirichlet, Neumann and Robin conditions was run on an unstructured triangulated grid containing a hole. By using the method of manufactured solutions, we demonstrated that the numerical solution converged with a second-order rate.
Funding Open access funding provided by University of Bergen (incl Haukeland University Hospital) The authors have not disclosed any funding.
Data availability Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Declarations
Competing interests The authors declare that they have no known competing interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.