Error Estimates on Hybridizable Discontinuous Galerkin Methods for Parabolic Equations with Nonlinear Coefficients

HDG method has been widely used as an effective numerical technique to obtain physically relevant solutions for PDE. In a practical setting, PDE comes with nonlinear coefficients. Hence, it is inevitable to consider how to obtain an approximate solution for PDEwith nonlinear coefficients. Research on usingHDGmethod for PDEwith nonlinear coefficients has been conducted along with results obtained from computer simulations. However, error analysis on HDG method for such settings has been limited. In this research, we give error estimations of the hybridizable discontinuous Galerkin (HDG) method for parabolic equations with nonlinear coefficients. We first review the classical HDGmethod and define notions that will be used throughout the paper. Then, we will give bounds for our estimates when nonlinear coefficients obey “Lipschitz” condition. We will then prove our main result that the errors for our estimations are bounded.


Introduction
In this paper, we obtain uniform-in-time convergence error estimates for the semidiscretization by hybridizable discontinuous Galerkin (HDG) methods for the parabolic equation with nonlinear coefficient.
Parabolic equation (or sometimes denoted as "heat equation") describes the distribution of heat in a certain region over time with given boundary conditions. Hence, physically, can be interpreted as a time-dependent function that describes the temperature at a given location. Despite the importance, both practical and purely theoretical, of finding exact solutions for parabolic equations, it is very difficult to obtain closed-form solutions for the parabolic equation with nonlinear coefficients. Details regarding parabolic equations and the way of applying numerical approximations in various contexts can be found in diverse sources such as [1][2][3] or [4].
As said earlier, parabolic equations are particularly interesting as the equations inherently contain physical meanings, that is, information regarding heat transfer. And in many contexts, in order to explain naturally observable phenomenon, parabolic equations with linear coefficients are insufficient: we inevitably face those equations with nonlinear coefficients. One of the most fundamental difficulties of dealing with parabolic equations (or any PDE) with nonlinear coefficients is that finding an exact solution is extremely difficult, if impossible (whereas existence of the solution can be easily shown when we require certain conditions on coefficients). Thus, we tend to focus on finding numerical approximations.
Various methods have been developed to find approximate solutions for given parabolic equations having nonlinear coefficients (or more generally, any PDE). Some examples of popular methods include but are not limited to finite 2 Advances in Mathematical Physics volume method [5], finite difference method [6], continuous Galerkin (CG) method [7], discontinuous Galerkin (DG) method [3,[7][8][9], HDG method [10,11] (the method that we will mainly focus on), and mixed methods [12,13]. Among the listed methods, DG has been studied the most.
Since Douglas Jr. et al. [14,15] and Arnold [16] first introduced DG method for parabolic equations (and, at the same time, elliptic equations), it has been classically used to give physically meaningful approximate solutions [17], primarily because the method is notably advantageous over the continuous Galerkin (CG) method in many contexts. First, DG can be used in a much broader context as it can be used on general meshes. Furthermore, degrees of polynomials can be arbitrary. Also, DG can produce highly accurate discretization for convection-diffusion equations. It can also be applied to solve problems with unambiguous boundary conditions. Since DG can provide approximate solutions with high accuracy and it can be used on general meshes, the method has been widely used to solve nonlinear problems. Detailed explanations on DG methods can be found in [18].
However, there are some disadvantages on using DG method, and this inevitably gave a rise to HDG method, which tends to alleviate disadvantageous facets of DG. As [19] shows, HDG method, when compared to DG analyses (see [18,20] and/or [21], e.g.), DG method produces larger globally coupled degrees of freedom than CG when the mesh size is invariant, since the degree of freedom is not shared by the boundary elements. In short, DG method is computationally ineffective when it is compared to CG method [18], and this shortcoming is the very reason the standard hybridizable discontinuous Galerkin (HDG) method was first introduced in [17] and developed. In HDG setting, the degrees of freedom associated with the numerical trances of the field variables solely matter in algebraic system. Degrees of freedom are substantially reduced as the numerical traces are only defined on the interelement boundaries. Consequentially, HDG method is very efficient [19].
As we are dealing with numerical approximations, one may ask how accurate our approximations are, and this question is what we would like to answer in this study. Previously, optimal convergence order for convection-diffusion equations was studied in 2 norm of + 1 polynomials of degree (see [19,22]). Then, the choices of stabilization parameter were presented in numerical values and analyzed [23]. Then, based on optimal convergence and superconvergence of HDG methods, local postprocessing was developed in [24] for linear convection-diffusion equations and in [25] for nonlinear case to increase the convergence order of numerical solutions. Despite the importance of solving parabolic equations and the substantial number of researches on applying various techniques to give approximate solutions that tend to hint that, at least for now, HDG method seems to be the best tools we have in our pocket, error analysis has not been conducted yet. To have a sense on how accurate the approximate solution is, we obtained through HDG a method which is indispensable, primarily because we need to check whether the result obtained from the method is meaningful or not. We limit ourselves to parabolic equations with nonlinear coefficients and conduct error estimation to investigate how accurate our approximation is.
In this paper, we will conduct an error analysis within HDG frame. We will use HDG projections of exact solution, satisfying certain properties. Then, by computing the magnitude of the norm of the difference between the projected solution and the approximate solution, we compute how far HDG approximations can deviate from the exact solution. The main idea is to derive HDG error equations (that will be introduced in the third section), and, by using the error equations, we will derive several identities that will give an upper bound for the error.
The paper is organized in the following way: after this introduction, we will give some introduction on HDG method, and along with that, we will define notions regarding HDG method that will be used throughout the paper. Then, in Section 3, we will give a priori estimation while assuming certain bounds for the nonlinear coefficient. We will then give error estimations for the proposed HDG method.

Notations and Norms.
Let T ℎ be a conforming, shaperegular simplicial triangulation of our domain Ω [26]. For any element ∈ T ℎ , is defined to be the set of edges of when dim( ) = 2. When dim( ) = 3, it is defined as the set of the faces of and is denoted by . Let T ℎ fl ∪ ∈T ℎ . Now, let ℎ denote the set of all edges/faces of the triangulation T ℎ . 0 ℎ is the set of all interior faces of the triangulation. Now, for any element ∈ T ℎ , ℎ is the diameter of , and define ℎ fl max ∈T ℎ ℎ . Call this number the mesh size.
Throughout the paper, we will use the standard notations for Sobolev spaces and their norms on the domain Ω and the boundary. For example, ‖V‖ ,Ω , |V| ,Ω , ‖V‖ , Ω , |V| , Ω , > 0, denote the Sobolev norms and seminorms on Ω and its boundary Ω. For an integer , the Sobolev spaces are Hilbert spaces and the norms are defined by the 2 norms on their weak derivatives up to the order . When is not an integer, the spaces are defined by the interpolation [27].
Advances in Mathematical Physics 3
For each time on the interval [0, ], the method yields a scalar approximation ℎ ( ) to ( ), a vector approximation q ℎ ( ) to q( ), and a scalar approximation̂ℎ( ) to the trace of ( ) on element boundaries, in spaces of the form respectively, where Here, P ( ) fl [P ( )] and P ( ) is the space of polynomial of total degree at most . With the spaces, HDG method provides approximations ( ℎ , q ℎ ,̂ℎ) ∈ ℎ × V ℎ × ℎ , determined by the following five restrictions: For any ( , k, ) ∈ ℎ × V ℎ × ℎ , we require with a numerical trace for the flux defined for some nonnegative stabilization parameter defined on T ℎ , which we assume to be piecewise constant on T ℎ . Here, 0 is an 2 projection of 0 onto P . Now, for vector-valued functions u, k ∈ ( 2 ( )) , define Then, we introduce the following notations: [11], is defined as follows. Given for all faces of the simplex . Also, denotes the 2 orthogonal projection onto ℎ . Lemma 1. If the local spaces ( ( ), V( )) = (P ( ), P ( )) for ≥ 0 and are nonnegative, then the systems (9a), (9b), and (9c) are uniquely solvable for Π q and Π . Furthermore, there is a constant independent of the choice of and such that, for all 1 ≤ ≤ + 1, Proof. See [11].
The existence and the uniqueness (which depend on the Lipschitz condition that will be presented in the next section and the choices of the approximation spaces and the stabilization parameter ) on the solution for the system of equations (7a) to (7f) can be found in [28,29].
We will later use the above lemma to figure out the convergence orders between our estimated solution and the exact solution.

A Priori Estimate
Let us first give an assumption on ( ). Remember that ( ) = ( ) −1 . Hence, we are giving restrictions to our nonlinear coefficient.
is chosen in such a way that there exist positive constants 1 , 2 , 3 and 4 such that, for all ∈ R, we have the following inequalities: We further assume that both ( ) and ( ) are Lipschitz continuous.

Advances in Mathematical Physics
With this condition on the nonlinear coefficient, we have the following result:

Lemma 3. If Assumption 2 holds, then one has
where is a constant that does not depend on mesh size ℎ.
Proof. Take k = q ℎ in (7a), = ℎ in (7b), = −q ℎ ⋅ n in (7c), and = −̂ℎ in (7d). Adding the resulting four equations, we get where Using integration by parts and the definition ofq ℎ ⋅ n, (7f), we get Now use Cauchy-Schwarz and Young's inequality and observe that Integrating with respect to time over the interval (0, ) and using ℎ (0) = Π 0 along with the above results, we get the following inequality: Using Gronwall's lemma, we get This completes the proof.

Error Estimations
Let us now derive error estimates for the proposed method.

Error Equations.
For the remaining sections, we define Remark 4. Note that, by triangle inequality, The first term (i.e., ‖ −Π ‖) that appears in RHS is bounded by Lemma 1. Hence, we only need to find an upper bound for Π − ℎ =: . Similarly, we only need to bound e to give an upper bound for ‖q ℎ − q‖.

Lemma 5. One has
Proof. For all ( , k, ) ∈ ℎ × V ℎ × ℎ , the exact solution ( , q) obviously satisfies the following four equations: Since is the 2 projection into ℎ and satisfies the orthogonal property, we have Advances in Mathematical Physics
This completes the proof.
Remark 6. Note that, by definition, ‖ (0) = Π 0 − 0 ‖ ≤ ‖Π 0 − 0 ‖ + ‖ 0 − 0 ‖. Using Lemma 1, as in Remark 4, the two terms in the right hand side both vanish when ℎ is small enough. Hence, from now on, to ease the computation, we assume that (0) = 0; that is, although in reality, to take account for (0), we have two extra nonzero terms ‖Π 0 − 0 ‖ and ‖ 0 − 0 ‖ in equalities that we would show, we will ignore the two terms from now on for simplicity.

Estimations for in ∞ ( 2 )
Lemma 7. For any ∈ (0, ], one has where Proof. Take k = e in (21a), = in (21b), = −ê⋅ n in (21c), and = −̂in (21d). Add the four equations that we get. We then get where Using integration by parts and the definition of ê⋅ n, (21f), we have The identity we wanted to prove follows after integrating in time over the interval (0, ) and using the fact that (0) = 0 by (21e).

Theorem 8. If Assumption 2 is satisfied and
where is independent of mesh size ℎ.
Proof of Theorem 8. Consider S 's we defined in the previous lemma. Applying Cauchy-Schwarz and Young's inequalities to the equation that defines S 1 , we get the following inequality: Similarly, we get the following inequality for S 3 .
Since is Lipschitz continuous and q ∈ H div , we have 6

Advances in Mathematical Physics
Since is bounded below by 1 , we get Therefore, we have Take < 1 /2 and apply Gronwall's lemma. We have Therefore, (38) Corollary 9. If Assumption 2 is satisfied and ∈ 1 ( 2 ), one has where is independent of mesh size ℎ.
Proof. We obtained bounds for S 's from the proof of Theorem 8. Using Lemma 7, we obtain the desired upper bound.
Remark 10. Note that the above result gives (together with triangle inequality) an upper bound for ‖ − ℎ ‖. Depending on how we discretize the time interval (0, ], the convergence order with respect to time is determined. Note that the convergence order of ‖ − Π ‖ with respect to the given space was already determined by Lemma 1. Hence, we only need to sum (or integrate if our discretization was continuous) with respect to time to determine the convergence order.

Estimations for e in ∞ ( 2 )
Lemma 11. If Assumption 2 is satisfied and ∈ 1 ( 2 ), one has where is independent of a mesh size ℎ.
Proof. We differentiate all of error equations (21a)-(21e), with respect to time, and obtain Take k = e , = , = − ê⋅ n, and = −̂in the above four equations in that order. Adding the resulting four equations and following the steps of calculation as in the previous lemma, we get where Note first that Lemma 1 tells us that Thus, whenever we have suitably small ℎ, we will be able to make ‖ ‖ as small as we want. Keeping this in mind, using Cauchy-Schwarz' and Young's inequalities, we get

Advances in Mathematical Physics
To estimate T 2 , observe that (use the Lipschitz assumption on and the fact that is bounded) We then get the following estimations for T 2 : Apply Cauchy-Schwarz and Young's inequalities and get the following inequality for T 3 : Now let us derive an inequality for T. Remember that it is a well-known that when we have a finite dimensional vector space endowed with two norms ‖ ⋅ ‖ 1 and ‖ ⋅ ‖ 2 then the two norms are equivalent. That is, for any vector k in that vector space, there exists two constants 1 and 2 , independent of the choice of k, such that 1 ‖k‖ 1 ≤ ‖k‖ 2 ≤ 2 ‖k‖ 1 . Now, observe that Note that (using Lipschitz condition for and the fact that is bounded) the second term For the first term, observe that (using Hölder's inequality and the fact that all norms are equivalent in a finite dimensional vector space) Using Corollary 9, we can take (by decreasing the mesh size of ℎ) ‖e ‖ as small as we want. Proceeding (similarly to how we derived an inequality for T 1 ), we get Applying the Gronwall's lemma, we have 8

Advances in Mathematical Physics
Applying the error estimates of Theorem 8 and Corollary 9, we obtain Therefore, and we are done.

Lemma 12.
For any ∈ (0, ], one has where Proof. We keep all the error equations except for (21a) and (21c); instead of the two, we use the equations obtained from differentiating the two equations with respect to time . That is, we have the following equations: for all ( , k, ) ∈ ℎ ×V ℎ × ℎ , where ê⋅n = e ⋅n+ ( −̂) on T ℎ .
Substitute k = e to the first equation, = to the second one, = −ê⋅ n to the third one, and = −̂to the fourth one. Adding the resulting four equations, we get where Integrating in time over the interval (0, ), we can get the desired identity.

Conclusion
In this research, we have presented error estimates on HDG methods for parabolic equations having nonlinear coefficients. Under Lipschitz assumption, that is, Assumption 2, we have found upper bounds for errors e , . Using Lemma 1 and "integrate with respect to time " (or, if time is discretized, summing), we can find how far our estimated solutions can deviate, in terms of a mesh size ℎ from the actual solution. As HDG method is very practical and efficient when compared to classical DG methods such as [16,18], or [30], the error estimation gives us how reliable our estimation is. There are some possible generalizations of our result. First, we can replace Lipschitz assumption by more general conditions on which the existence of solutions for PDE can be guaranteed. Other further research topics include changing the local spaces for the projection. Our main tool for the research was the projection. We projected the solutions to the space of polynomials, and we can possibly use 2 projection instead.
In our future work, we plan to propose and develop the method by implementing it to the multiscale settings, which will give how the method can be applied to find an approximate solution with satisfying error bounds, without being too time-costly.
As aforementioned, parabolic problems with heterogeneous coefficients inherently are practically applicable in many situations such as porous media. We hope our method would further enrich related studies.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.