Hybrid BSQIWENO Based Numerical Scheme for Hyperbolic Conservation Laws
Abstract
In this paper, we intend to use a Bspline quasiinterpolation (BSQI) technique to develop higher order hybrid schemes for conservation laws. As a first step, we develop cubic and quintic Bspline quasiinterpolation based numerical methods for hyperbolic conservation laws in 1 space dimension, and show that they achieve the rate of convergence 4 and 6, respectively. Although the BSQI schemes that we develop are shown to be stable, they produce spurious oscillations in the vicinity of shocks, as expected. In order to suppress the oscillations, we conjugate the BSQI schemes with the fifth order weighted essentially nonoscillatory (WENO5) scheme. We use a weak local truncation based estimate to detect the high gradient regions of the numerical solution. We use this information to capture shocks using WENO scheme, whereas the BSQI based scheme is used in the smooth regions. For the time discretization, we consider a strong stability preserving (SSP) RungeKutta method of order three. At the end, we demonstrate the accuracy and the efficiency of the proposed schemes over the WENO5 scheme through numerical experiments.
Keywords:
Weighted Essentially NonOscillatory Scheme Bspline Quasiinterpolation∎
1 Introduction
In this work, we intend to use BSpline QuasiInterpolation (BSQI) to develop numerical schemes for the hyperbolic conservation laws
(1.1) 
where and is the flux function, along with the initial condition
(1.2) 
It is wellknown that the classical solution of (1.1) may ceases to exist in finite time, even the initial data is sufficiently smooth. The higher order finite difference, spectral, and finite element methods are well suitable in approximating the solution in the smooth regions, but fail to provide a nonoscillatory solution in the presence of jump discontinuities. Thus, it is challenging to devise an higher order numerical scheme, which preserves the higher order accuracy in smooth regions and produce nonoscillatory solution in a vicinity of shocks.
The BSQI approximations are known (see Sablonnière [29]) to be very accurate in approximating smooth functions and their derivatives. This motivated many researchers to use this technique to develop numerical methods for some partial differential equations. For instance, Li et al. [23] developed a meshfree method for the viscous Burgers’ equation, and Zhu and Kang [36] developed a finite difference method for the BurgersFisher equation. Recently, we (Kumar and Baskar [16]) have developed higher order schemes for some Sobolev type equations with special reference to equal width and BBMB equations. But these equations involve dissipative and dispersive effects that lead to smeared shocks and dispersive shocks, respectively. Zhu and Kang [35] made an initial effort in devising cubic Bspline quasiinterpolation (CBSQI) based numerical scheme for hyperbolic conservation laws. But their study is restricted only to scalar conservation laws and also they used quasilinear form of the equation to develop the numerical scheme and presented numerical results for Burgers’ equation up to the critical time at which the shock forms. To our knowledge, there is no work on developing a BSQI based conservative numerical method for hyperbolic conservation laws.
In the present work, our interest is to develop higher order methods for (1.1)(1.2) using cubic and quintic BSQI (denoted by CBSQI and QnBSQI, respectively) for approximating space derivative and the 3 order RungeKutta method for time derivative. Although the resulting schemes are conditionally stable, as they are higher order methods, they develop spurious oscillations near the shocks. It is known that such a problem is common for any higher order methods. There are at least two approaches to minimize the nonphysical oscillations in the numerical solution, namely, using limiters and using essentially nonoscillatory (ENO) reconstructions. Out of these, the ENO approach can minimize the oscillations without much compromise in the order of accuracy.
An improved version of the ENO approach is the well known Weighted ENO (WENO) reconstruction introduced by Liu et al. [24] (also see Jiang and Shu [13]). However, this reconstruction algorithm is expensive. In order to minimize the expense, there are many hybrid schemes introduced in the literature for the hyperbolic conservation laws. For instance, a nonconservative hybrid method is introduced by Harabetian and Pego [11], the multidomain methods in one and two space dimensions are developed using Fourier continuation and WENO scheme by Shahbazi et al. [32, 33]. Costa and Don [4, 5] obtained a high order hybrid method by combining the spectral and the central schemes with WENO. Cheng and Liu [3, 2] introduced multidomain RKDG methods with proper boundary treatments. The characteristics hybrid compactWENO scheme is proposed by Ren et al. [27].
In developing efficient hybrid schemes, choosing smooth indicator is an important component of the algorithm. Karni et al. [15], used a Weak Local Truncation Error (WLTE) based smooth indicator in developing an adaptive algorithm for hyperbolic systems. Further Kurganov and Liu [17] developed adaptive artificial viscosity methods for the system of Euler equations of motion using Weak Local Residual (WLR) (similar to WLTE). The WLR is used in many other problems, for instance, Mungkasi et al. [25] used WLR for the system of shallow water equations and Chen et al. [1] developed the adaptive artificial viscosity method for the SaintVenant system using this indicator. To reduce the computational cost in the case of system of Euler equations, Dewar et al. [8] introduced WLR pressure based smooth indicator. A comparison study of many switch indicators for hyperbolic conservation laws is performed by Li and Qiu [21, 22] for uniform and curvilinear meshes.
In this paper, we propose higher order hybrid BSQIWENO schemes for hyperbolic conservation laws. In a BSQIWENO scheme, we conjugate the BSQI based numerical scheme with a WENO scheme. To detect the high gradient region of the numerical solution, we use a WLR based estimate of Kurganov and Liu [17]. We use this information to capture shocks using WENO scheme, whereas the BSQI based scheme is used in the smooth region. For the time discretization, we consider a SSP RungeKutta method of order three. To validate the numerical scheme, we perform the numerical experiments using Burgers and BuckleyLeverett equations in the case of scalar conservation laws and Euler equations in the case of system of equations. It is also interesting and important to study the efficiency and the accuracy of the hybrid methods constructed using the BSQI based numerical schemes and WENO schemes.
The outline of the paper is as follows. In section 2, we discuss the BSQI and its application in approximating the derivative of a function. In section 3, we develop the semidiscrete form of hyperbolic conservation laws using CBSQI and QnBSQI approximations to the space derivative. We further perform the von Neumann stability analysis and provide the numerical evidence for the order of accuracy. We develop two hybrid numerical schemes, namely, the Hybrid4 and Hybrid6 by combining the CBSQI with WENO3 and QnBSQI with WENO5, respectively, in section 4. We perform some numerical experiments in section 5, to validate the schemes for scalar as well as system of equations. Finally, we test the accuracy and the efficiency of the proposed algorithms over the pure WENO in section 6.
2 BSpline QuasiInterpolation
In this section, we briefly discuss the construction of BSQI of degree upto 5 and give the numerical rate of convergence in each cases. For more details on the basic concepts of spline spaces and Bsplines we refer to Schumaker [31], De boor [6], Nürnberger [26], and Kvasov [19].
Let the space interval be divide into subintervals , of equal length . The boundary points of the subintervals are called the knots and let
denotes the set of knots. The spline space of degree is defined by (see Schumaker [31])
(2.3) 
where is the polynomial space of degree . The set forms a vector space of dimension and the set of Bsplines of degree forms a basis for this space. Bsplines can be constructed recursively by starting with the Bspline of degree 0 given by
(2.4) 
The recursive formula for the higher degree Bsplines takes the form
(2.5) 
The general form of the Bspline quasiinterpolant (BSQI) of degree for a function is given by (see Sablonnière [29])
(2.6) 
where the coefficients can be obtained by imposing the condition that the quasiinterpolant is exact if is a polynomial of degree less than or equal to . The coefficients are the linear combination of function (approximate) values at the nodes, where

if is odd, then the nodes are simply the knots, i.e. , and the set of nodes involved in calculating is
(2.7) and

if is even, then the nodes are taken as and the set of nodes involved in calculating is
(2.8)
For more details, we refer to Sablonnière [29].
In case when support of Bspline
contains distinct knots, the coefficients for quadratic (QBSQI), cubic (CBSQI), quadric (QdBSQI), and quintic (QnBSQI) are given as follows:
Quadratic:
(2.9) 
Cubic:
(2.10) 
Quadric:
(2.11) 
Quintic:
(2.12) 
We have the following theorem for general BSQI in approximating a function:
Theorem 2.1
There exists a constant such that for all and for any set of knots on with step size , we have
From the construction we can see that a Bsplines of degree is times continuously differentiable. Thus, we can also approximate the derivatives of a function with the corresponding derivatives of BSQI as long as it is possible, i.e.
(2.13) 
Approximation  BSpline  
Function 
Quadratic  0  0  0  0  
()  Cubic  0  0  0  0  
Quadric  
Quintic  
1 derivative 
Quadratic  0  0  0  0  0  
()  Cubic  0  0  0  0  0  
Quadric  0  
Quintic  0 
The approximation of function and its first derivative at point are as follows:
Case I. When is even, we have
(2.14) 
and
Case II. when is odd, we have
(2.15) 
where coefficients are to be obtained from (2.13) and the derivatives of the Bsplines. For these coefficients are listed in Table 1.
In the following example, we demonstrate the accuracy and the rate of convergence in approximating the first derivative of a smooth function.
Example 1
Consider the function
(2.16) 
whose first derivative is given by In order to make the comparison among the BSQI’s in terms of their accuracy, we consider the error norm given by
(2.17) 
where the superscript denotes the derivative and denotes the error associated with the derivative. In Figure 1(a), we depict the loglog plot of the error and the number of mesh points, in the case of the first derivative approximation of the model problem (2.16). The QBSQI converges to the first derivative with rate two, whereas CBSQI and QdBSQI achieves the fourth order of convergence, which is clear from Figure 1(b). The QnBSQI has a higher order accuracy among the considered BSQI, which has the rate of convergence six in the case of the first derivative approximation.
3 BSQI Scheme
In this section, we discuss the semidiscrete form of (1.1) obtained using the BSQI based spatial approximation. For the sake of notation simplicity, we restrict to , and consider the scalar conservation laws of the form
(3.18) 
The BSQI formulation for a general 1D system (1.1) can be done componentwise in a similar way.
We treat the flux function in (3.18) as and use the BSQI approximation (2.6) in (3.18) to get
(3.19) 
Let us use the notation , for , to denote the approximation of the solution obtained using (3.19) and write the semidiscrete form of (3.18) as
(3.20) 
The above methods can work for any degree of BSQI. But in this work, we intend to use only the CBSQI and QnBSQI formulations, which are given, respectively, by
(3.21) 
for and
(3.22) 
for , where we have used the notation . The above semidiscrete formulations can be written in the conservation form as
(3.23) 
where denotes the numerical flux given by
(3.24) 
for the CBSQI scheme (3.21) and
(3.25) 
for the QnBSQI scheme (3.22).
Observe that (3.23) forms a system of ODEs. When the flux is linear, then (3.23) leads to a system of linear ODEs. In this case, by assuming that the solution of (3.23) in the form
(3.26) 
where and with being the wavenumber, we get
where is given by
For the time discretization, we use the Strong Stability Preserving (SSP) RungeKutta method of order three given by (see [10, 28] for more details)
Note that the spectrum of the RungeKutta method of order three is known to contain some part of the imaginary axis. Thus, this method along with the CBSQI and QnBSQI space discretizations can be stable under the condition
where and (see Hirsch [12] for more details). This implies that, we can find the condition on such that the scaled spectrum of CBSQI and QnBSQI schemes are contained in the stability spectrum of the RungeKutta method of order three.
In the following example, we demonstrate the order of accuracy of the CBSQI and QnBSQI schemes in the case of linear advection equation.
Example 2
Consider the linear advection equation
(3.27) 
with the initial data
(3.28) 
The rate of convergence of a numerical scheme, denoted by , for the error , is given by
(3.29) 
where and denote the error and the number of mesh points, respectively. From Example 1, we expect the rate of convergence, with respect to the space variable, of the CBSQI and the QnBSQI schemes to be four and six, respectively. In order to achieve the rate numerically, we take the time step . The numerical solution is computed over the domain at time . Table 2 and Table 3, clearly shows that, as the number of mesh points increase, the numerical solution obtained with CBSQI and QnBSQI schemes converge to the exact solution with rate four and sixth, respectively.
error  order  error  order  error  order  
20  3.188811e04  –  0.001283  –  1.283344e03  – 
40  2.020435e05  3.980281  8.097906e05  3.986215  8.097906e05  3.986215 
80  1.267336e06  3.994795  5.068834e06  3.997823  5.068834e06  3.997823 
160  7.925669e08  3.999122  3.170066e07  3.999069  3.170066e07  3.999069 
320  4.954097e09  3.999839  1.981728e08  3.999682  1.981728e08  3.999682 
error  order  error  order  error  order  
20  3.998450e05  –  0.000026  –  2.840256e05  – 
40  6.376781e07  5.970469  4.075987e07  5.974427  4.522632e07  5.972714 
80  1.013043e08  5.976060  6.453003e09  5.981035  7.166100e09  5.979831 
160  1.585622e10  5.997504  1.009394e10  5.998409  1.121209e10  5.998061 
320  2.633005e12  5.912195  1.673563e12  5.914423  1.859182e12  5.914243 
We now show numerically that the CBSQI and QnBSQI schemes achieve the expected order of accuracy in the case of Burgers equation as long as the solution remains smooth.
Example 3
Consider the Burgers’ equation
(3.30) 
In order to show accuracy of the BSQI schemes, we consider (3.30) with the following initial data
(3.31) 
error  order  error  order  error  order  
40  1.211681e03  –  0.000161  –  3.518324e04  – 
80  9.369829e05  3.692843  1.048627e05  3.937758  2.465823e05  3.834747 
160  6.419103e06  3.867579  6.624985e07  3.984440  1.586015e06  3.958591 
320  4.077475e07  3.976624  4.155074e08  3.994971  9.982610e08  3.989845 
640  2.553719e08  3.997004  2.598439e09  3.999157  6.250058e09  3.997476 
error  order  error  order  error  order  
40  3.802114e04  –  0.000047  –  1.137541e04  – 
80  1.388433e05  4.775273  9.976366e07  5.547309  2.911152e06  5.288185 
160  2.604317e07  5.736409  1.771954e08  5.815101  5.314697e08  5.775459 
320  4.307657e09  5.917858  2.849299e10  5.958591  8.664328e10  5.938756 
640  6.843121e11  5.976105  4.613797e12  5.948509  1.370955e11  5.981835 
In Table 4 and Table 5, we compare the , , and error and their rate of convergence obtained using the CBSQI and the QnBSQI schemes, respectively. The numerical solution is computed over the at time , at which the solution remains smooth. In order to achieve the expected rate of convergence, we choose the time step significantly small to make time discretization error small. In both cases, we choose the time step . From the tables, it is clear that the CBSQI scheme converges to the exact solution with a rate more then four, whereas the QnBSQI scheme initial has order nearly four which increased further to order more then six. ∎
4 Hybrid Scheme
The CBSQI and the QnBSQI schemes developed in the previous section are stable with 3 order RungeKutta time discretization and also they achieve the expected order of accuracy when the solution remains smooth, as demonstrated in Example 2 and 3. But the main disadvantage of these schemes is that the numerical solution develops spurious oscillations in the case of nonsmooth solution. In order to minimize the oscillations without much compromise in the order of accuracy, we devise an hybrid scheme, where we use the BSQI scheme in the region where the solution is smooth, whereas in the nonsmooth region we propose to use WENO scheme. Thus, we propose the hybrid semidiscrete scheme for (3.18) to be
(4.32) 
where
(4.33) 
with being a smooth indicator, which is a constant over the cell i.e . Using this in (4.33), the above scheme can be written as
(4.34) 
Here, denotes either CBSQI or QnBSQI flux given by (3.24) and (3.25), respectively. Whereas denotes the WENO flux of appropriate order. More precisely, we conjugate WENO flux of order 3 with CBSQI and form a hybrid scheme with numerical flux
(4.35) 
and we combine QnBSQI with WENO flux of order 5 to get
(4.36) 
To be selfcontent, we quickly recall the WENO fluxes and for more details, we refer to Jiang and Shu [13].
The flux can be decomposed into positive and negative parts as
(4.37) 
such that and . For instance, these can be taken as
where . We recall that the stencil of WENO3 is and , and the corresponding numerical flux is written as
where the positive part is given by
with
and and . The negative part can be written in a similar way.
The stencil of WENO5 is
and the positive part of the numerical flux is written as
where
and
(4.38) 
Here the smooth indicators are given by
and , and .
Thus, the hybrid scheme switches between WENO and BSQI scheme, depending upon the smooth indicator. Once the shock region is detected, the hybrid scheme switches to the expensive shock capturing WENO scheme, otherwise the solution is computed using the BSQI scheme.
For the smooth indicator, we use the weak local truncation error (WLTE) introduced by Karni and Kurganov [14]. The WLTE is derived on the fact that a weak solution of (3.18) satisfies the integral form given by (see Godlewski and Raviart [9] )
(4.39) 
for all . The variation in the value of from smooth to discontinuous regions and vice versa for computed solution can be taken as a measure of smooth indicator. The is referred as weak truncation error for with respect to . In practice, the calculation of weak truncation error is seem to be a difficult task, since is a general test function. To overcome this difficulty Kurganov and his coworkers (see [17, 14]) used test function with Bsplines as
(4.40) 
where
(4.41) 
and
(4.42) 
are the quadratic and the linear Bsplines with the localized supports of size and , respectively. Putting (4.40) in (4.39), we can arrive at
(4.43) 
Karni and Kurganov [14] established an estimate for the WLTE (4.43) and remarked that the WLTE bound can be converted into error bound in the smooth region. Although their result is valid for one dimension scalar problem, they showed numerical that the WLTE behavior is similar in the case of 1D system of equations. With the aid of these estimates and further numerical experiments, Kurganov and Liu [17] proposed the following estimate and used it as a smooth indicator:
where , and is the accuracy of the numerical scheme. In our hybrid scheme (4.34), we propose to use the smooth indicator defined by
(4.44) 
where is some empirically chosen nonnegative real number. In general, we take in our numerical experiments. To make smooth transition between smooth and rough regions, we extend the rough region to the neighboring node points of . That is, once the discontinuity is detected in a vicinity of , which is the case if , then we take
(4.45) 
where is again an empirically chosen positive integer. In our numerical computation, we take . The smooth transition can also be made using a smaller value of , but that lead to over estimate of discontinuous regions.
5 Numerical Experiments
This section comprises the numerical implementation of hybrid BSQIWENO schemes for linear and nonlinear problems. The BSQI schemes are oscillatory in nature at shock positions, whereas the hybrid BSQIWENO schemes developed in the previous section are nonoscillatory because of the usage of the WENO scheme in the nonsmooth regions and the high order accuracy is maintained in the smooth regions using BSQI approximation. The developed algorithm works for any degree of BSQI, but to demonstrate the idea we use CBSQI and QnBSQI for space approximation in conjugation with WENO3 and WENO5, respectively. We study the accuracy of the CBSQIWENO3 and the QnBSQIWENO5 schemes in terms of rate of convergence. Along with the accuracy, we also demonstrate the efficiency of the proposed hybrid BSQIWENO.
We start the numerical experiment with smooth initial data.
Example 4 (Smooth Initial Data)
In order to measure the accuracy and convergence rate of the hybrid algorithm, we consider the Example 3.1 and Example 3.2. In both cases, we take in smooth indicator. We fixed the time step and increase the mesh points in space. In Figure 2(a), we compare the accuracy of the Hybrid4 and the Hybrid6 schemes with the WENO5 scheme in case of linear advection equation at time . Both the Hybrid4 and the Hybrid6 schemes achieve their expected convergence rate. In terms of accuracy, we observe that the Hybrid6 scheme is better than both the WENO5 and the Hybrid4 schemes. Similar behavior is observed from Figure 2(b) in the case of Burgers’ equation at time .
As a next step, we perform the numerical experiment in the case of the linear advection equation with discontinuous initial data.
Example 5 (Linear advection equation)
Consider the linear advection equation (3.27) with initial data
(5.46) 
which has jump discontinuities at points and moves towards right as time progress. The CBSQI and QnBSQI schemes capture the correct position of the discontinuities but with oscillations in the profile (not shown here). We now compute the solution using the Hybrid6 scheme whose numerical flux is given by (4.36). In Figure 3(a), we compare the numerical solution with the exact solution at time along with the initial data. The number of mesh points is taken to be and the CFL number is . In the figure, the symbol ‘’ represents the region where the WENO5 is used whereas the symbol ‘’ represents the region where QnCBSQI is used. We observe that the discontinuities are captured without oscillations. In Figure 3(b), the plot of smooth indicator over the plane is shown. It is clear from the figure that the smooth indicator moves along the characteristic direction and detect the discontinuity positions correctly. Also, we observe that the solution is computed with WENO5 only at discontinuity positions, whereas in the smooth regions QnBSQI is used to compute the solution.∎
We now test the Hybrid6 scheme in the case of some nonlinear scalar conservation laws. We first start with the inviscid Burgers equation, which is a common test case for convex flux.
Example 6 (Nonlinear equation with convex flux  Burgers equation)
To show the shock capturing ability of the Hybrid6 scheme (4.36), we consider the Burgers’ equation (3.30) with discontinuous initial data (5.46). The numerical solution is computed at time over the domain with mesh points. The CFL number is taken to be . The solution contains a rarefaction wave and a shock moving towards right. As expected the QnBSQI scheme is oscillatory at the shock position, which spreads over the smooth region of the solution as the time progresses, which can be observed from Figure 4(a).
Figure 4(b) depicts that the numerical solution at time obtained using the Hybrid6 scheme (4.36). Here, we observe that the scheme captures shock and rarefaction accurately without any oscillation. The smooth indicator over the plane is shown in Figure 4(c), where we observe that the indicator moves along the discontinuity. The left state of the initial discontinuity at the point is connected with the right state by the rarefaction fan, where the singularity at the lower part of the rarefaction wave is smoothen after time due to numerical dissipation, which is well illustrated by the smooth indicator.
In Figure 4(d), we show the plot of the percentage of the computational domain in which WENO5 is used verses the number of grid points. We kept the time fixed at and increased the number of mesh points and observed that the percentage of the region where WENO5 scheme is used has been reduced drastically from to when the number of mesh points is increased from to . The low percentage of usage of the WENO scheme for large number of mesh points, makes the hybrid BSQIWENO schemes efficient.∎
We now turn our attention to (3.18) with nonconvex fluxs. Since vanishes at one or more points, (3.18) fails to be genuinely nonlinear which makes the solution more complicated as it may involve the formation of the composite waves. The composite waves are generally the combination of joined rarefaction and shock waves, for more detail see LeVeque [20]. Many popular schemes fail to approximate the entropy solution of (3.18) with nonconvex flux, for instance semidiscrete centralupwind schemes (see Kurganov et al. [18]). Here we test the Hybrid6 scheme (4.36) in case of the BuckleyLeverett problem.
Example 7 (Nonconvex flux  BuckleyLeverett equation)
Consider the conservation law (3.18) with nonconvex flux, called the BuckleyLeverett equation,
(5.47) 
along with the initial data (5.46). In Figure 5(a), the numerical solution obtained using QnBSQI scheme at time is compared with the reference solution over the domain with and CFL is taken to be . In this figure, we observe that the QnBSQI scheme fails to produce nonoscillatory result. The rarefaction fan is missing and the numerical solution is oscillatory at the rarefaction position, which makes the solution to moves with different speed. As a result, the QnBSQI scheme fails to capture the correct positions of the shocks. In comparison to QnBSQI, the Hybrid6 scheme (4.36) produces nonoscillatory solution with correct positions of shocks as is apparent from Figure 5(b) and thus converges to the entropy solution. In Figure 5(c), we show the smooth indicator on plane and the percentage of the usage of WENO5 scheme is depicted in Figure 5(d). As observed in the case of Burgers equation, here also we see that the usage of WENO5 is approximately when the number of mesh points approaches . This shows that, the Hybrid6 scheme captures the entropy solution with less computational complexity. ∎
We continue with the scalar conservation law (3.18) with another nonconvex flux.
Example 8 (Another nonconvex flux)
Another test case for nonconvex flux is (see for instance Kurganov et al. [18])
(5.48) 
which belongs to , whose convexity changes at . We consider the following initial data
(5.49) 
The exact solution is given by