Virtual element method for nonlinear Sobolev equation on polygonal meshes

In this work, the virtual element method (VEM) on convex polygonal meshes for the nonlinear Sobolev equations is developed, where the semi-discrete and fully discrete formulations are presented and analyzed. To overcome the complexity of nonlinear terms, the nonlinear coefficient is approximated by employing the orthogonal L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L^{2}} $$\end{document} projection operator, which is directly computable from the degrees of freedom. Under some assumptions about the nonlinear coefficient, the existence and uniqueness of the semi-discrete solution are analyzed. Furthermore, a priori error estimate showing optimal order of convergence with respect to the H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{H^{1}}$$\end{document} semi-norm was derived. Finally, some numerical experiments are conducted to illustrate the theoretical convergence rate.


Introduction
In this paper, we apply the conforming virtual element method scheme to study a nonlinear Sobolev boundary value problems: find u = u(x, t) such that where ⊂ R 2 be a bounded convex polygonal domain with smooth boundary ∂ , and a(u), b(u) and f (x, t) are given functions.The nonlinear Sobolev equation arises in many physical phenomena such as the flow of fluids through fissures rock [1], two-temperature heat conduction [2], the migration of moisture in soil [3] and other important applications.For details about the existence, uniqueness, and the stability of the solution of the problem (1.1), we can refer to [4][5][6].A large number of papers have focused on the numerical approximations of nonlinear Sobolev equations of the form (1.1).For example, the Galerkin methods were discussed for one-dimension and two-dimension nonlinear Sobolev equations in [7] and [8].In [9], a semi-discrete finite element method for nonlinear Sobolev equations with nonlinear boundary condition is constructed and the optimal L 2 error estimates are proven.The optimal convergence rates in H 1 and L 2 -norm are demonstrated with the time stepping along the characteristic finite element method to the semilinear and nonlinear Sobolev equations in [10].There are also some non-standard finite element methods based on the nonlinear Sobolev equations that have been developed such as the two-grid finite element method [11], mixed finite element method [12], split least-squares mixed finite element method [13] and low order characteristic-nonconforming finite element method [14].
The superconvergence result of order O(h 2 + τ 3 ) in H 1 -norm was displayed in [15], which developed a linearized three-step backward differential formula Galerkin finite element method.In addition, fully discrete discontinuous Galerkin approximations for nonlinear Sobolev equations were introduced in [16] and the interior penalty discontinuous Galerkin method was applied in [17].Recently, in [18], the authors considered a new hybrid high-order method for the Sobolev equation with convection-dominated term.
The virtual element method (VEM), first introduced in [19], is a generalization of the finite element method to deal with very general polygonal/polyhedral meshes.The main idea behind VEM is defining the virtual element space and a set of degrees of freedom, and using a suitable projection operator to approximate bilinear forms which avoid integrating the non-polynomial functions on the element but without loss of accuracy.The salient features of the method are as follows: (a) it does not require an explicit form of shape functions; (b) it can be applied to arbitrary polygonal/polyhedral meshes including convex and concave elements; (c) it can be easily handled in the adapted meshes (hanging nodes and interfaces are allowed) and (d) the stiffness matrix and mass matrix can be computed only use degrees of freedom.Due to its flexibility, the method has been successful for a wide range of problems.Some of them include linear elasticity problems [20,21], Stokes equation [22,23], inelastic problem [24,25], electromagnetic problems [26][27][28], contact problem [29], linear parabolic and hyperbolic problem [30,31], Sobolev equations [32] and the nonlocal elliptic problem of Kirchhoff type [33].The VEM framework has been concurrently extended to a number of different problems and applications, and, in particular, the literature on VEM for nonlinear problems is growing.By exploiting the external and internal L 2 projections, [34] approximated the trilinear term in the Navier-Stokes equation.The quasilinear elliptic problems in [35] were evaluated with the element-wise polynomial projection of the virtual element ansatz.Furthermore, exploiting internal L 2 projection operator, [36,37] deduce the virtual element formulation for nonlinear parabolic and hyperbolic problems asserting optimal order of convergence in L 2 and H 1 norms.However, to the best of our knowledge, the application of the VEM to the nonlinear Sobolev equation is scarce in the literature.
The purpose of this paper is to establish semi-discrete and fully discrete approximations for the nonlinear Sobolev equation on arbitrary polygonal meshes using the VEM described in [35].Using the projection operators P k h , P k h and R k h , we build three discrete bilinear forms that may be computed, with the nonlinear term is appropriately approximated to ensure VEM polynomial consistency and stability.It is worthy to be emphasized that due to the discretization of the nonlinear term, additional error has to be evaluated in the convergence analysis.Specifically, with the help of the virtual interpolation estimate given in [38], a priori error estimate showing the optimal order of convergence with respect to H 1 semi-norm was derived.Finally, some numerical examples are conducted in order to demonstrate the theoretical convergence rate.
The paper is organized as follows.In Section 2, we give some notations and necessary assumptions on the data and the regularity of the solution of the problem (1.1).In Section 3, we present the virtual element space and construct computable discrete virtual element forms.In Section 4, we build the semi-discrete scheme, demonstrate its well-posedness, and present the optimal error estimates of the H 1 semi-norm.In Section 5, we propose the fully discrete scheme and prove the corresponding optimal error estimates in the H 1 semi-norm.In Section 6, we establish the convergence of fixed-point iterations for the solution of the nonlinear system and carry out three numerical examples.Finally, we provide a brief conclusion of the paper which is drawn in Section 7.

Notations
In this paper, we denote ⊂ R 2 as a bounded convex polygonal domain with smooth boundary ∂ and T as a positive time that defines the time interval I := (0, T ].We adopt standard notations for the Lebesgue spaces L p ( ) and the Sobolev spaces H s ( ) and their associated norms • s, and semi-norms | • | s, .Also, for simplicity, we drop whenever it is possible.In addition, we use the notation L p (0, T ; H s ( )), 1 ≤ p ≤ ∞, s ≥ 0, for the Banach space of all L p integrable func- ) and standard modification at p = ∞.Throughout this paper, C denotes a generic positive constant, which does not depend on the spatial and temporal discretization parameters, and may be different at its different occurrences.

Assumptions
The domain is decomposed into non-overlapping polygonal elements {T h }, 0 < h ≤ 1, satisfying the following mesh-regularity assumptions: Assumption 2.1 [19] For any fixed h > 0 and k ∈ N, we have: where L ∈ N is a positive constant.
Furthermore, for each element K ∈ {T h }, we denote by P k (K ), the set of polynomials of degree k on K .Conventionally, P −1 (K ) = {0}.Meanwhile, we use M k (K ) and M k (e) to denote the basis sets of P k (K ) and P k (e), that is, x K is the centroid of K , x e is the midpoint of e, and h e is the length of e. Hereafter, we denote the L 2 norm by • .
Furthermore, we assume that the following regularity conditions on a, b, f and u are satisfied: Assumption 2.2 [17] For t ∈ [0, T ] and (x, p) ∈ × R, we have: • There exist constants a 0 and a * , s.t.0 ≤ a 0 ≤ a(x, p), b(x, p) ≤ a * .• Let a(x, p), b(x, p) and f (x, p) be continuously differentiable with respect to each variable and assume uniform bounds for ) is a unique solution to (1.1), and u ∈ L 2 (0, t; H s ( )), u t ∈ L 2 (0, T ; H s ( )), for s ≥ 2. • ∇u and ∇u t are bounded in L ∞ ( × (0, T ]).

Variational formulation
The standard variational formulation of (1.1) is to find u(•, t) ∈ L 2 (0, T ; where Henceforth, it will be assumed that the problem (2.1) has a unique solution u, and in the appropriate places to follow, additional conditions on the regularity of u which guarantee the convergence results will be imposed.

Virtual element spaces
In this subsection, the local and global virtual element spaces are displayed for any order k and arbitrary polygonal meshes K .To compute the discrete form, it is sufficient to define the following useful element-wise projections on T h : with, in order to uniquely determine R k h , the addition of the following condition: 123 From now on, we follow the construction of the simplest local C 0 conforming virtual element discrete space on general meshes [39]: where for each edge e of ∂ E .
Then, we define the following set (DoF) of linear operators from where Therefore, the global virtual element space V h is obtained from the local spaces The global degrees of freedom for V h are the direct extension of the local ones (DoF).Moreover, the virtual element space has three f undamental properties: • The above projection operators and the standard L 2 projection operator 0,K k : V h (K ) → P k (K ) are computable just by accessing the local degrees of freedom given by DoF, • The global virtual element space V h ⊂ H 1 0 ( ) is a finite-dimensional subspace.

Construction of the computable discrete virtual element forms
We now define discrete versions of the bilinear form From Section 3.1, we see that V h (K ) contains not only the polynomials in P k (K ) but also non-polynomial functions satisfying the constraints.By making use, for all K ∈ T h , of the L 2 projection operator R h := R K h , P h := P k h and where S K m and S K E are symmetric bilinear functions, which are detailed below.Parameters a and b in (3.4) and (3.5) are chosen to guarantee two positive constants λ * and λ * independent of h and K such that Consider local symmetric bilinear S K m and A variety of computable choices for S m and S E can be found in the literature.The global discrete forms

as the sum of local contributions as
) 4 The semi-discrete scheme In this section, we establish the virtual element discretization of problem (2.1) for using the discrete bilinear forms proposed in Section 3, prove its well-posedness, and analyze the error estimates in the H 1 semi-norm.The semi-discrete scheme for (2.1) is to find where u I ,0 is the interpolation of u 0 onto V h .We proceed to derive the following theorem which shows the existence and uniqueness of the virtual element problem (4.1).

Theorem 4.1 If the functions a, b and f in (1.1) satisfy Assumption 2.2 then there exists a unique semi-discrete VEM approximation u h (x, t) satisfying (4.1) and u h (t)
1 is continuous with respect to t.
Then (4.1) deduces a system of nonlinear ordinary differential equations such that which implies the following initial value problem where u = (u T is a vector.Two matrices G and S and a vector F are defined by From the stability of A h and m h , and the Assumption (2.2), we can obtain which implies G is positive definite.By the regularity conditions on a, b and f , and the theory of ODE, u(t) exists and is unique for t > 0. This result suffices to show that u i (t) is continuous with respect to t so that u h (t) 1 is also continuous.

Error analysis
First, we recall some valid results that will help us prove our estimate.Consider the local polynomial interpolation estimates: Lemma 4.1 [40] There exists a positive constant C depending on k and ρ such that, for all K ∈ T h and any v ∈ H s+1 (K ) and w ∈ (H s+1 (K )) 2 , it holds: The approximation properties of the virtual element space are characterized by the following interpolation error bound, whose proof can be found in [38].Lemma 4.2 Suppose that Assumption 2.1 is satisfied and let s be a positive integer such that 2 ≤ s ≤ k + 1.Then for any w ∈ H s ( ) there exists an element where C is positive constants depending only on k and on ρ.
Remark 1 If K is a bounded domain and v ∞,K < C, v 0,K < C, then using the inverse inequality on polynomials, property of P h and standard estimation, we have Lemma 4.3 [35] For f ∈ L 2 (0, T ; H s−1 ( )), 0 ≤ s ≤ k, there exists a positive constant C, independent of h and f , such that The previous Lemmas allow us to derive the following error in the H 1 semi-norm.
Theorem 4.2 Let u be the solution of (2.1) and u h be its approximation defined by where r = min{s, k + 1} and C a positive constant independent of h.
Proof Let u I ∈ V h be the interpolant of u.Decompose the error as follows: which are then estimated separately.Using Lemma 4.2, for all t ∈ (0, T ), we can easily find (4.12) Combining (1.1) and (4.1), adding and subtracting the term

13)
Letting v h = θ and rearranging, we have We start with I 1 .Using Lemma 4.3 and the fact that r ≤ s, we have Making use of the definition of m h (•, •) and the property of P h , we obtain Again, using property (3.9),Cauchy-Schwarz inequality, generalized Hölder inequality, property of P h and Lemma 4.1, we get ) By the definition of projection, and using the fact that ∇ R h q is a (vector) polynomial of degree ≤ k − 1, we have which implies Thus, in view of the stability of Using the same justification as before, we can estimate Using (3.8), (3.9), Assumption 2.2, Cauchy-Schwarz inequality and Lemma 4.1, we get On the other hand, using (3.7), (3.8), property of P h and P h , we can estimate Integrating (4.27) with respect to t from θ(x, 0) = 0, we get By using Gronwall's inequality, it yields Combining (4.12) and (4.29), we complete the proof of the theorem.

The fully discrete scheme
In this section, we consider a backward Euler scheme for (1.1) and present the corresponding optimal error analysis.
Let τ and u h (t n ) ∈ V h be the time step and the approximation of u(t) at time t = t n respectively.Now we define the fully discrete scheme for (1.1) as follows: Find where N is a positive integer, The error estimates for | U n −u n | 1 with u n = u(x, t n ) are derived in Theorem 5.1.
Theorem 5.1 Let u n and U n be the solutions of (2.1) and (5.1) at t = t n , respectively.Assume that (5.2) where r = min{k + 1, s} and C a positive constant independent of h.

Proof
In analogy with the previous theorem, let u n I ∈ V h be the interpolant of u n onto V h .We denote 3) The term η n can be bounded as in (4.12) (5.4) As usual, for the estimate for θ n , we subtract (5.2) from (1.1) to arrive at (5.5) Let v n = θ n , we can rewrite (5.5) as (5.6) Using (3.8), (3.9), Assumption 2.2, Cauchy-Schwarz inequality and Lemma 4.1, we get (5.7) Using the definition of m h (•, •), and add and subtract (u n − u n−1 , θ n ), and the property of P h , we get (5.8) Using the definition of A h (•; •, •) once again, we obtain (5.9) Using (3.7), property (3.9),Lemma 4.1 and generalized Hölder inequality, we have next adding and subtracting τ a(P h U n ) + a(u n ) P h ∇u n t + a(u n )∇u n t , P h ∇θ n K , and using the property of P h , we get (5.11) Summing this new bound of (5.10) and (5.11), it follows that (5.12) Similarly, we can rewrite I I 4 as (5.13) Using (3.7), property (3.9),Lemma 4.1 and generalized Hölder inequality, we deduce Adding and subtracting (b(P h U n + b(u n ))∇u n , P h , we have (5.15)Then, we can deduce the bound of I I 4 (5.16) Applying the Lemma 4.3, we can bound the term I I 5 as follows (5.17) It is easy to observe that and (5.19)

Implementation aspects
In all tests, we consider the time interval [0, 1], the domain = [0, 1] × [0, 1] with distortion square meshes Q h and distortion Voronoi meshes V h .These polygonal meshes are generated by using the data in VEM [41].An example of the adopted meshes is shown in Fig. 1 Furthermore, we use the fixed-point iterations to solve the nonlinear system.
Starting with a given U n 0 ∈ V h , we construct a sequence U n l , l ≥ 0, such that which implies The convergence in H 1 -norm of the sequence U n l as l → ∞ to a fixed-point of (6.1), and hence a solution of (5.1).
with c * = min{m * , a 0 α * + τ a 0 β * }.Then using (6.2), we obtain where ) .There, F(U n l ) is a decreasing sequence and, in view of the fact that

Table 1
The relative error E 1 h,τ for Example 1 on the meshes ) → 0, as n → ∞, which completes the proof.
Meanwhile, we evaluated the relative errors in the approximate H 1 semi-norm at the final time T , which are defined by We solve Example 1 by using the lowest-order virtual element (k = 1).The corresponding numerical results for E 1 h,τ relative meshes are shown in Tables 1  and 2 respectively.As the exact solutions are smooth, as expected that both meshes yield similar accuracy.In this test, we notice that the error is almost constant in h for large values of τ , recalling Theorem 5.1, the error due to spatial discretization dominates the error generated by the time discretization.The coefficients are same as Example 1, and the right-hand side f is calculated accordingly.
We also solve Example 2 by using the lowest-order virtual element ( k = 1).The corresponding numerical results for E 1 h,τ on relative meshes are shown in Tables 3  and 4 respectively.The convergence orders in the direction of the space with τ = 1/64 on V h and with τ = 1/256 on Q h for Example 2 are displayed in Fig. 2. We can observe indeed that, for small values of τ , the expected convergence order in h is O(h) for the approximate H 1 semi-norm.

Example 3
We consider the nonlinear Sobolev problem with a(u) = u +1, b(u) = u 2 .The exact solution is set to be u = sin(π x)(1 − y)ye 0.1t   and the load term f and the initial data u 0 are chosen with the exact solution.
We solve Example 3 with the VEM discretization coupled with the backward Euler method on V h and Q h .Errors and their rates of convergence are calculated for the orders of approximation k = 1, k = 2 and k = 3 respectively.The corresponding numerical results for the relative error E 1 h,τ are shown in Tables 5 and 6.From Tables 5  and 6, we can notice that the expected convergence order in h for the approximate H 1 semi-norm, that is, O(h k ) for k = 1, 2, 3.

Conclusions
In this paper, the VEM is employed to solve the nonlinear Sobolev equation, where the spatial variable is discretized by the conforming virtual element space, and the temporal variable is discretized by using the backward Euler method.In particular, for the computation of nonlinear term, which appears on the left-hand side of the equation, a L 2 -orthogonal projection operator is used, and P h is used for the computation of local bilinear forms involved in discrete formulation.Our main result in this paper is the optimal error estimates in H 1 semi-norm in Theorem 5.1, which shows that the VEM can achieve optimal convergence rate, and the convergence order due to the temporal variable discretization is one order.Furthermore, three numerical examples are given to verify the validity of the proposed method.The experimental results show a complete agreement with the theoretical analysis results.

Theorem 6 . 1 Fig. 1
Fig. 1 Example of polygonal meshes: on the left is distortion square meshes, on the right is distortion Voronoin meshes

Fig. 2
Fig. 2 Convergence orders in space direction for Example 2 on the mesh of V h (left) and Q h (right) for τ = 1/64 and τ = 1/256 respectively

Example 2
In order to investigate the convergence order in h, we weaken the influence of the temporal discretization and consider replacing the exact solution in Example 1 with u = x y(1 − x)(1 − y)e 0.1t .
satisfies star convexity with respect to a ball B ρ , where ρ is a positive constant satisfying ρ ≥ γ h K and γ is h independent positive constant.•For each side e ⊂ ∂ K , | e |> ch K , | • | denotes the measure and c is a positive constant independent of mesh-size h K .• For all elements K ∈ {T h }, the total number of edges n

Table 5
The relative error E 1 h,τ for Example 3 on the meshes

Table 6
The relative error E 1 h,τ for Example 3 on the meshes