Application of C-Bézier and H-Bézier basis functions to numerical solution of convection-diffusion equations

Convection-diffusion equation is widely used to describe many engineering and physical problems. The finite element method is one of the most common tools for computing numerical solution. In 2003, Wang et al. proposed C-Bézier and H-Bézier basis functions which are not only a generalization of classical Bernstein basis functions but also have a free shape parameter bringing a lot of flexibility to geometrical modeling. In this paper, we adopt C-Bézier and H-Bézier basis functions to construct test and trial function spaces of finite element method to get numerical solution of convection-diffusion equations. Compared with Lagrange basis functions, numerical accuracy is improved by 1−3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$1-3$\end{document} order-of magnitudes which implies a much better approximation in simulating convection-diffusion problems. Several examples are presented to verify the feasibility and effectiveness of our method.

The finite element method is an important tool for solving PDEs in various fields of science and engineering. In 1943, Courant [12] firstly proposed the finite element method and used piecewise linear functions to construct the Galerkin projection space. In the 1960s, Feng [13] published the paper entitled "Based on The Difference Scheme of Variational Principle" regarded as the first proof of convergence of finite element method. In the 1970s, Babuska and Brezzi [14] presented the general theoretical framework of the mixed finite element method. Subsequently, Falk and Osborn [15] improved the mixed finite element method.
Based on the classical finite element method, some scholars have proposed several new methods to solve convection-diffusion equations numerically. In 2005, Adjerid and Klauser [16] constructed a local discontinuous Galerkin method for dimensional-one transient convection-diffusion problems and gave an asymptotically posterior error estimation. Maleknejad [17] gave stable numerical algorithms for heat convection and integral equations. In 2008, Liu et al. [18] studied a mixed time discontinuous space-time finite element method for convection-diffusion equation; this method lowered the order of equations. In 2016, Jin et al. [19] proposed a Petrov-Galerkin finite element method, which combined "shifted" fractional powers with continuous piecewise linear elements to construct the test and trial function spaces. To solve convection-diffusion problems in convective dominance, Baumann and Oden [20] presented a discontinuous hp finite element method; it was appropriate for unrestricted mesh composed of arbitrary convex polygons. This new method has the features of finite element and finite volume method; numerical simulations show that it is robust and capable of providing high accuracy. In 2020, Mirzaee proposed a new meshless method to solve 2D time fractional Tricomi equations and stochastic time-fractional sine-Gordon equations [21,22].
Scholars have focused considerable attention on spline functions with excellent geometrical properties [23][24][25]. For instance, Qin et al. [26] constructed a class of new polynomial basis functions with n local shape parameters in 2013. These spline functions satisfy the conditions required for G 0 , G 1 , G 2 continuity. Many kinds of spline basis functions can be used to construct function spaces of finite element method. In 1979, Shi [27] introduced spline finite element method to solve the equilibrium problem of the plate-beam composite elastic structures in regular regions and derived a unified computation scheme for various boundary conditions. In 2005, inspired by CAD, Hughes [28], Academician of The American Academy of Sciences, proposed a new spline finite element method which applies NURBS basis functions in finite element analysis for isometric analysis. In 2007, Bhatti and Bracken [29] suggested that the method combining Bernstein basis function and Galerkin method can achieve better results in solving partial differential equations.
The spline basis functions are widely applied in finite element method to get better numerical solutions, such as Bernstein, B-spline, and NURBS functions. Bernstein basis functions form bases of polynomial space, which can be used to construct function spaces of finite element method, but they fail to express transcendental curves, for instance, the cycloid and helix. In addition, NURBS model also has several weaknesses. On the one hand, the rational form is unstable. On the other hand, it increases the complexity for computing integrals and derivatives. For overcoming the disadvantages of rational Bernstein basis functions, scholars made a lot of efforts to construct basis functions. In 2003, Wang et al. [30,31] applied an integral approach to construct C-Bézier basis functions for the space T = span{1, t, . . . , t n-2 , sin t, cos t} that extended the spaces of mixed algebra and trigonometric polynomial. After that, Wang developed H-Bézier basis functions for the space T = span{sinh t, cosh t, t k-2 , t k-3 , . . . , t, 1}.
C-Bézier and H-Bézier basis functions are one kind of spline functions that have a nonrational form and are capable of expressing circular arc and polynomial curves of high order exactly. Furthermore, C-Bézier and H-Bézier basis functions are generalization of classical Bézier basis and keep their excellent properties. In addition, the free variable α can be utilized to change the shape and geometry of curves. In previous work, scholars obtained good approximate solutions by combining spline functions with the finite element method [32][33][34]. Therefore, we consider using spline functions C-Bézier and H-Bézier basis to construct function spaces in the finite element method and expect the approximate solutions to be more accurate than the classical Lagrange basis function.
In this paper, we combine the Galerkin finite element method with C-Bézier and H-Bézier basis functions to solve convection-diffusion equation. C-Bézier and H-Bézier basis functions are selected to construct trial and test function spaces. The error analysis and numerical examples are given, and the numerical results indicate that our method has much better precision in solving convection-diffusion equations. The numerical solutions in this paper are generated in Matlab 2018a.
The rest of the paper is organized as follows. In Sect. 2, we recall the definitions and properties of C-Bézier and H-Bézier basis functions. In Sect. 3, we use C-Bézier and H-Bézier basis functions to construct the test and trial function spaces. Furthermore, the framework of the finite element method to solve convection-diffusion equations is presented. In Sect. 4, we discuss the stability of the difference scheme. In addition, error estimation of the finite element method with C-Bézier and H-Bézier basis functions is provided. In Sect. 5, several numerical examples are given to illustrate the theoretical results. In the final section, we not only summarize some comments on the overall work but also touch upon some avenues of future research.

C-Bézier basis and H-Bézier basis functions
In this section, we recall C-Bézier and H-Bézier basis functions proposed by Wang et al. A free shape parameter is introduced in these basis functions to portray polynomial and circular curves of high degree precisely.

Definition 1
The C-Bézier basis function {C n i (t)} n i=0 for the space T = span{1, t, . . . , t n-2 , sin t, cos t} of degrees n is defined as follows where α is referred as shape parameter [31].
While n = 2, C-Bézier basis function is expressed as follows.

Definition 3 The H-Bézier basis function {H
. . , t n-2 , sinh t, cosh t} of n degrees is defined by While n = 2, the H-Bézier basis function is expressed by Eq. (4).
The H-Bézier basis functions have many beneficial properties in common with Bernstein basis functions. Compared to classical Bernstein basis, C-Bézier and H-Bézier basis functions have a free variable allowing adjustment of the local or global shape and geometrical properties of the curves. In addition, even though shape parameter may change, the position and tangent direction of the curve's endpoints do not alter [35]. Furthermore, in the splicing of multiple C-Bézier curves, the shape parameter will not transform the original continuity between curves [36].

Discrete model of convection-diffusion equation
In this section, we consider combining the Galerkin finite element approach with C-Bézier and H-Bézier basis functions to solve the following convection-diffusion equation [37].  Denote by an open, bounded, connected subset of R 1 and ∂ the boundary of . An initial condition u 0 and a boundary condition g(x) on ∂ for u are prescribed. a = a(x, t) and b = b(x, t) are the diffusion and convection coefficients; the right-hand side, f = f (x, t) represents the source function [38]. To solve the above model with the finite element method, C-Bézier and H-Bézier basis functions are selected to construct trial and test function spaces. We give the weak formulation of one-dimensional convection-diffusion in Eq. (5) and discrete the weak formulation with the Galerkin method. First of all, denoted by H 1 0 ( ) the subspace of the Sobolev space H 1 ( ): Multiplying the first equation of Eq. (5) by test function v yields Using Green formula and divergence theorem, one can obtain the following: If we define two operators as then the variational scheme of Eq. (5) may be read as, looking for an u(x, t) satisfies saying that u is the generalized solution of Eq. (5) [39]. Now, we divide into mesh elements uniformly, denoted by h , then define C-Bézier basis functions or H-Bézier basis functions on each element E of h . The finite dimensional subspace U h = span{φ j } n j=1 of H 1 0 ( ) are constructed taking all the basis functions together, where φ j are actually C-Bézier basis functions or H-Bézier basis functions.
Subsequently, the Galerkin formulation of the above problem is to find u h ∈ U h , such that As φ j , j = 1, . . . , n is base of finite dimensional space U h , u h (x, t) and is expressed as Taking v i = φ i , i = 1, . . . , n, we get a system of ordinary differential equations for u 1 (t), u 2 (t), System Eq. (11) is expressed in vector form as We assume = [0, L], then is the stiffness matrix, is the load vector, is the unknown vector. So far, we gain a semi-discretization Galerkin function Eq. (11). If we further discretize the space-time domain [0, T], let 0 = t 0 < t 1 < · · · < t N = T, t k = kτ , where τ is the time step. Then obtain a completely discretized Galerkin equation. For instance, using the forward difference quotient using the backward difference quotient using the Crank-Nicolson format [40] 1 In which superscript k is the approximation at t = t k = kτ . We use the Crank-Nicolson format to discretize the space-time domain; its stability will be proved in the next chapter.

Stability and error analysis
Our first goal is to present an error estimate in L 2 when equation (5) is semi-discrete. (5) and equation (11), suppose u 0 = 0 on ∂ . Then,

Lemma 1 Let u and u h be the solutions of equation
Here, we require the solution of the continuous problem has the regularity implicitly assumed by the presence of the norms on the right [37].
In this paper, the Crank-Nicolson format is applied to discrete time t. Here, the semidiscrete form is discretized in a symmetric fashion around the point t k+ 1 2 = (k + 1 2 )τ . This will generate a second order accurate in time. Below, we prove that this format is stable with initial values. Let v h = u n+1 hu n h , f = 0, we have where · expresses L 2 norm. Taking advantage of the symmetry of a(·, ·), we obtain Therefore, the stability of format Eq. (15) is proved. Next, the error estimate reads as follows.
Theorem 1 Let U k be the solution of problem (15). U k represents the value of u h at t = t k . u is the solution of problem (5), and let u 0 hu 0 ≤ Ch r u 0 r and u 0 = 0 on ∂ . Then, Specific proof can be found in reference [41].

Numerical examples
In this section, several numerical examples are given to solve convection-diffusion functions with C-Bézier and H-Bézier basis functions, which are used to construct test and trial function spaces of finite element method. In contrast to Lagrange basis functions, the accuracy of numerical solution is improved by 1 -3 order-of magnitudes, which implies a much better approximation in simulating convection-diffusion problems. The errors for the finite element method shall be measured in three norms defined as follows.

Example 1 Consider following convection-diffusion equation with Dirichlet conditions
Its exact solution is u = t sin x cos x.
In this example, we construct test and trial function spaces consisting quadratic C-Bézier basis function. Divide the [0, 1] into n elements, the ith element is denoted by Finite element nodes include all the mesh nodes and the middle points of all the mesh elements. We firstly consider the reference C-Bézier basis function on the reference element [0, 1]. As shown in Eq. (1), the range of t is [0, α], let t = mα, m ∈ [0, 1]. Then the reference C-Bézier basis function on the reference element [0, 1] is 1-cos α , φ 2 = -1-cos(mα)+cos α-cos(α-mα) Here, m ∈ [0, 1], α ∈ (0, π]. Next, we use the affine mapping between an arbitrary interval [x i , x i+1 ] and the reference interval [0, 1] to construct the local basis functions from the reference. Equation (19) is the quadratic C-Bézier basis function on A finite dimensional subspace will be constructed by putting these local basis functions together.
The comparison of the errors in L ∞ norm with quadratic C-Bézier and quadratic Lagrange basis functions is indicated in Fig. 5(a) (h = 1 4 , t = 1). Figure 5(b) presents the corresponding colormap of convection diffusion when 0 ≤ x ≤ 1, 0 ≤ t ≤ 1. The numerical errors obtained from C-Bézier and Lagrange basis functions are set out in Table 1.  In Table 1, h is the step size of subdivision, and α is the shape parameter. It stands out that the precision of numerical solutions with C-Bézier basis function is one order-of magnitude higher than Lagrange basis function with the same step size h. Example 2 In this example, the convection-diffusion function with Dirichlet conditions is given as The exact solution of the above problem is u = t 2 sin(πx); we use quadratic C-Bézier and Lagrange basis functions to create trail space and test space. Figure 6(a) compares the errors in L ∞ norm of C-Bézier and Lagrange basis functions. The table below illustrates the specific data of errors.
Comparing the errors of numerical solution obtained by C-Bézier and Lagrange basis functions with h unchanged, it stands out in Table 2 that the numerical solutions of C-Bézier basis function are 4 -7 order-of magnitudes more accurate than Lagrange basis function.  The error graph and convection diffusion diagram Table 3 The numerical errors Example 3 Considering the following equation where the exact solution is u = t cosh x. We make quadratic H-Bézier and Lagrange basis functions construct trail space and test space. The errors in L ∞ norm got from Lagrange and H-Bézier basis functions can be compared in Fig. 7(a) (h = 1 4 , t = 1). Figure 7(b) provides the picture of convection diffusion when 0 ≤ x ≤ 1, 0 ≤ t ≤ 1. Numerical values for the errors are given in Table 3.
What is striking about the figures in this table is the precision of numerical solutions with H-Bézier basis function is 5 -6 order-of magnitudes higher than Lagrange basis function. Example 4 In this example, we consider The exact solution for above problem is u = 4t cosh 2 ( x 4 ). The quadratic H-Bézier and Lagrange basis functions are used to solve Example 4. Figure 8(a) shows the comparison  of errors in L ∞ norm, which obtained by Lagrange and H-Bézier basis functions (h = 1 4 , t = 1), numerical errors are given in Table 4.
The difference of errors between H-Bézier and Lagrange basis functions is highlighted in Table 4. Compared to the numerical results obtained by Lagrange basis function, the accuracy of numerical results got by H-Bézier basis function is improved by 2 -5 order-of magnitudes.
Example 5 In this example, we verify the effectiveness of H-Bézier basis function with shape parameter α ≥ 1. Consider the following equation with exact solution u(x, t) = t sinh x, quadratic H-Bézier and Lagrange basis functions are used to solve the above problem. Figure 9 shows the numerical errors in L ∞ norm (h = 1 8 , t = 1) and convection diffusion diagram. Table 5 provides the rate of convergence of H-Bézier and Lagrange basis functions in L ∞ , L 2 , and H 1 norm on uniform rectangular mesh with size h.
From the date in the Table 5, it can be seen that while Lagrange basis function and H-Bézier basis function have the same rate of convergence, the mesh size h for Lagrange basis function is much smaller, which implies that the computational complexity of Lagrange basis function is far greater than that of H-Bézier basis function.  From the above examples, it is surprising that the accuracy of numerical solutions with H-Bézier and C-Bézier functions is higher than Lagrange basis function. Consequently, we have verified the effectiveness and feasibility of our method.

Conclusions and future work
In this contribution, we adopt C-Bézier and H-Bézier basis functions to construct the test and trial function spaces of the finite element method. In contrast to classical Lagrange basis functions, our scheme gets much better approximation precision for convectiondiffusion equations. The limitation of this paper is that there is no systematic method to obtain the optimal shape parameters to minimize the error between the numerical solution and the exact solution.
There are two interesting aspects to be explored in future research. Firstly, how to get better shape parameter for solving PDEs numerically. Secondly, how to analyze the error of exact results and numerical solutions generated by C-Bézier and H-Bézier basis functions theoretically.