# A New Sixth-Order Finite Difference Compact Reconstruction Unequal-Sized WENO Scheme for Nonlinear Degenerate Parabolic Equations

Authors: Liang Li; Yan Zhang; Jun Zhu
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this paper, a new sixth-order finite difference compact reconstruction unequal-sized weighted essentially nonoscillatory (CRUS-WENO) scheme is designed for solving the nonlinear degenerate parabolic equations on structured meshes. This new CRUS-WENO scheme only uses the information defined on three unequal-sized spatial stencils, obtains the optimal sixth-order accuracy in smooth regions, and preserves the second-order accuracy near strong discontinuities. This scheme can be applied to dispose the emergence of the negative linear weights and avoid the application of the mapped function. The corresponding linear weights can be artificially set to be any random positive numbers as well as their summation is one. The construction process of this scheme is very simple and can be easily extended to higher dimensions, since the new compact conservative formulation is applied to approximate the second-order derivatives. The new CRUS-WENO scheme uses narrower large stencil than that of the same order finite difference classical WENO schemes. Therefore, it has a better compactness and smaller truncation errors. Some benchmark numerical tests including the porous medium equation and the degenerate parabolic convection-diffusion equation are performed to illustrate the advantages of this new CRUS-WENO scheme.


---

## Body

## 1. Introduction

In this paper, a new sixth-order finite difference compact reconstruction unequal-sized WENO (CRUS-WENO) scheme is proposed for solving the nonlinear degenerate parabolic equations on structured meshes. We first mention a few superiorities of this new CRUS-WENO scheme. A new spatial compact WENO reconstruction methodology is adopted for the second-order derivative, which is different to many high-order WENO schemes [1] that need to dispose the appearance of the negative linear weights and the mapped nonlinear weights. The CRUS-WENO scheme uses the narrower large stencil and only three unequal-sized spatial stencils are applied in the spatial reconstruction, which is distinct from many classical WENO schemes [1–3] on structured meshes. Finally, the new CRUS-WENO scheme could obtain its optimal sixth-order accuracy in smooth regions, get sharp transitions, and simultaneously preserve the second-order accuracy close to strong discontinuities when solving for the nonlinear degenerate parabolic equations.The nonlinear degenerate parabolic equations arise in many applications, such as porous medium flow, collisional transport models in plasmas, radiative transport, and so on. Associated solutions of the nonlinear degradation equations often do not exist, even if their initial condition is smooth enough. So the research for its weak solutions has become the main task [4–6]. Many numerical schemes are used to study such nonlinear degradation equations, due to the nature of the nonlinear degradation equations is similar to that of the hyperbolic equations. Many numerical schemes were designed for solving the hyperbolic equations, such as the linear schemes [7–10], kinetic schemes [4], and some high-order unconditionally stable schemes which were used to solve the nonlinear degenerate advection-diffusion equations [12], and so on. Therefore, we will construct a new finite difference compact reconstruction unequal-sized WENO scheme to approximate second-order derivative.We briefly retrospect the developing history of the ENO and WENO schemes [13–22]. In 1985, Harten and Osher [23] designed a total variation diminishing standard [24] and a high-order ENO schemes. Two years later, Harten et al. [25] extended the scheme to finite volume framework. In 1994, Liu et al. [17] first proposed the construction idea of finite volume WENO scheme. In 1996, Jiang and Shu greatly improved this scheme in [2], which was the most widely used WENO scheme. In 2011, [1] was used to solve the degenerate parabolic equations for the first time. They directly used the conserved flux difference scheme to approximate the second-order derivative. In order to deal with the negative linear weights, a special treatment [18] was introduced to avoid the oscillation near strong discontinuities. To avoid loss of precision near critical points, the mapped function was also introduced. Abedian [26, 27] used another method to overcome this trouble. In 2021, Jiang [28] explored an alternative formulation to solve the degenerate parabolic equations.Although those WENO schemes could achieve arbitrarily high-order accuracy, they suffered from poor spectral resolution and increasingly wider spatial stencils. The compactness of the scheme has become a standard for measuring the quality of the scheme. In 1992, Lele [29] proposed new compact schemes to achieve high-order accuracy with narrower spatial stencils by using implicit interpolations for the first-order and second-order derivatives. These schemes could have better bandwidth resolution and obtain smaller truncation errors in smooth regions. Therefore, a lot of efforts have been done to enhance the compactness of WENO schemes. Hermite WENO (HWENO) scheme [30] can reduce the size of the spatial stencils by adding a control equation. New compact scheme was proposed in [31], which was based on linear compact schemes [29]. Ghosh and Baeder [32] developed a fifth-order finite difference compact WENO scheme, in which three compact smaller spatial stencils were used to form a larger stencil through the nonlinear weights. In 2002, Pirozzoli [33] proposed a hybrid compact WENO scheme. In 2019, Shi and Guo [34] proposed a fifth-order compact WENO scheme in an alternative formulation for hyperbolic conservation laws. In 2018, Ma and Wu [35] also constructed compact Hermite WENO scheme for solving different problems. However, these compact schemes are difficult to generalize from the approximate of the first-order derivative to the second-order derivative. So we provide a new method to overcome this difficulty.In this paper, a new sixth-order compact reconstruction unequal-sized WENO (CRUS-WENO) scheme is presented for solving the degenerate parabolic equations. The scheme using the information defined on one four-point stencil and two two-point stencils is designed to approximate the second-order derivatives. If there is a discontinuity lies in the big spatial stencil, the information of this big stencil is effectively abandoned and either of the two-point stencils is applicable. After that, the numerical flux value at half points can be directly computed by solving a tridiagonal system. Together with the application of a third-order TVD Runge-Kutta time discretization [36], the new WENO scheme could obtain high-order accuracy both in space and time.The organization of the paper is as follows. In Section2, the new sixth-order finite difference CRUS-WENO scheme is constructed in detail. In Section 3, some classical numerical tests are provided to demonstrate the high-order accuracy and high resolution of this new CRUS-WENO scheme in one and two dimensions, respectively. Concluding remarks are given in Section 4.


## 2. New Sixth-Order CRUS-WENO Scheme

In this section, we propose the sixth-order compact spatial reconstruction procedure with the application of three unequal-sized compact stencils to deal with the second-order derivative in the conservative form. For simplicity, a uniform mesh is applied for the computational domain [a, b] with Δx=b−a/N,xi=iΔx in one dimension. A semidiscrete conservative finite difference scheme for solving one-dimensional degenerate parabolic equation(1)ut=fuxx,(in which f′(u) ≥ 0) has the form(2)ddtuit=Lu=f^i+1/2−f^i−1/2Δx2.The numerical fluxf^i+1/2=f^ui−p,…,ui+p is designed to satisfy the condition that the flux difference approximates the second-order derivative f (u)xx with the k-th order accuracy(3)f^i+1/2−f^i−1/2Δx2=fuxx+OΔxk.The construction of the numerical flux is our main objective in the following.


### 2.1. Reconstruction of Numerical Flux

We first define a functionh (x), such that(4)fux=1Δx2∫Δx/2x+Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ,and obtain(5)fuxxx=hx+Δx−2hx+hx−ΔxΔx2.We only need to approximate the functionhxi+Δx−hxi with numerical fluxf^i+1/2 to achieve the required order of accuracy in smooth regions. According to [1], we obtain(6)f^i+1/2=−2fui−2+25fui−1−245fui+245fui+1−25fui+2+2fui+3180,(7)f^i+1/2=hxi+Δx−hxi+1560d7hxdx7|x=xiΔx7+OΔx8.Now a sixth-order compact interpolation is given by(8)215f^i−1/2+1115f^i+1/2+215f^i+3/2=−9fui−1−153fui+153fui+1+9fui+2180.Now, we prove that this formula has sixth-order accuracy and calculate its truncation error. From the definition of (4), we get(9)−9fui−1−153fui+153fui+1+9fui+2180=−9/Δx2∫x−3Δx/2x−Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ−153/Δx2∫x−Δx/2x−Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ+153/Δx2∫x−3Δx/2x−3Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ+9/Δx2∫x−3Δx/2x−5Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ.180To analyze the accuracy of the numerical fluxf^i+1/2, we assume hx is smooth enough. Performing the Taylor series expansion for hx at the grid point xi(10)hθ=hxi+∑k=17θ−xikk!hkxi+OΔx8.By substituting the integral of the above formula, we get(11)f^i+1/2=∑k=16Δxkk!hkxi+145Δx77!h7xi+OΔx8.It is obvious that the Taylor series of the exact solution could expand to(12)hxi+Δx−hxi=∑k=16Δxkk!hkxi+Δx77!h7xi+OΔx8.Comparing (11) with (12), we can obtain the truncation error as(13)f^i+1/2=hxi+Δx−hxi+12800d7hxdx7|x=xiΔx7+OΔx8.We can see the coefficients of leading error term in compact interpolations are smaller than that of the noncompact interpolations. It can also be verified from the numerical examples that specified in the following. When the nonlinear degenerate parabolic equations contain the discontinuous solutions, these interpolation schemes might suffer from oscillations. So the high-order WENO schemes have been proposed as limiters to deal with this difficulty. The general form of the numerical flux reconstructed in the WENO fashion is(14)f^i+1/2=∑m=1rωmf^i+1/2m,where ωm is the nonlinear weights. The new way specified in this paper is to choose the unequal-sized stencils [28, 37, 38, 39]. We choose a six-point spatial stencil and two threepoint spatial stencils T1=xi−2,xi−1,xi,xi+1,xi+2,xi+3,T2=xi−1,xi,xi+1,T3=xi,xi+1,xi+2 and design the unequal degree polynomials defined on them, and then get their interpolations at half point xi+1/2:(15)f^i+1/21=−2fui−2+25fui−1−245fui+245fui+1−25fui+2+2fui+3180,f^i+1/22=fui+1−fui,f^i+1/23=fui+1−fui.Coincidentally, the approximations are the same on stencilT2 and T3. By the above mentioned formulas, we obtain the corresponding compact interpolations(16)215f^i−1/2+1115f^i+1/2+215f^i+3/2=−9fui−1−153fui+153fui+1+9fui+2180,12f^i−1/2+12f^i+1/2=fui+1−fui−12,12f^i+1/2+12f^i+3/2=fui+2−fui2.Similar to [38], three compact interpolations can be combined by the linear weights γ1,γ2,γ3 and the nonlinear weightsω1,ω2,ω3 as(17)ω11γ1215f^i−1/2+1115f^i+1/2+215f^i+1/2−γ2γ112f^i−1/2+12f^i+1/2−γ3γ112f^i+1/2+12f^i+1/2+ω212f^i−1/2+12f^i+1/2+ω312f^i+1/2+12f^i+3/2=ω11γ1−9fui−1−153fui+153fui+1+9fui+2180−γ2γ1fui+1−fui−12−γ3γ1fui+2−fui2+ω2fui+1−fui−12+ω3fui+2−fui2,in which, the linear weights and the nonlinear weights will be specified in the following. We use the similar definition [37] to construct the smoothness indicators:(18)β1=16048083708fui−22+2002064fui−12+7847144fui2+7847144fui+12−15310492fuifui+1+7205468fuifui+2−7738322fui+1fui+2+2002064fui+22+fui−2−791699fui−1+1457188fui−1308130fui+1+574412fui+2−99187fui+3−1308130fuifui+3+1457188fui+1fui+3−791699fui+2fui+3+83708fui+32+fui−1−7738322fui+7205468fui+1−3253987fui+2+574412fui+3,(19)β2=fui−1−2fui+fui+12,(20)β3=fui−2fui+1+fui+22.According to [40, 41], we define τ=β1−β3+β2−β3/22 Based on the linear weights and smoothness indicators, we define the nonlinear weights as(21)ωn=ω¯n∑l=13ω¯l,ω¯n=γn1+τε+βn,n=1,2,3.in which ε is a small positive number to avoid the denominator to become zero. We take ε = 10−6 in all computations. By performing such new spatial reconstruction procedure, we can achieve the sixth-order accuracy in smooth regions and preserve the essentially nonoscillatory property nearby the strong discontinuities.Through (17), the numerical flux value at half points can be directly computed by solving the following tridiagonal system(22)βγ0.........0αβγ.........00.........αβγ0.........0αβf^1/2f^i+1/2f^N−1/2=R^1/2−αf^1/2R^i+1/2R^N−1/2−γf^N−1/2,where R^i+1/2 is the right-hand side of (17), the coefficients α,β, and γ represent the coefficients on the left-hand side of 17. To solve this system of equations, we need to find the two unknown variables f^−1/2 andf^N+1/2 on the right-hand of 22, respectively. The values on the first and last half points along each grid line are reconstructed by the sixth-order US-WENO spatial reconstruction procedure [39]:(23)f^−1/2=ω1180γ1−2fu−3+25fu−2−245fu−1+245fu0−25fu1+2fu2+ω2−ω1γ1γ2+ω3−ω1γ1γ3fu0−fu−1,f^N+1/2=ω1180γ1−2fuN−2+25fuN−1−245fuN+245fuN+1−25fuN+2+2fuN+3+ω2−ω1γ1γ2+ω3−ω1γ1γ3fuN+1−fuN.The linear weights and nonlinear weights are calculated according to (21). Thus we can directly solve the three-diagonal linear equation system (22) by using the LU decomposition method.Remark 1.
The sixth-order finite difference WENO-JS scheme [1] uses the information of three equal-sized spatial stencils to construct three quadratic polynomials. Then their approximations are obtained at half point xi+1/2:(24)f^i+1/21=fui−2−3fui−1−9fui+11fui+112,f^i+1/22=fui−1−15fui+15fui+1−fui+212,f^i+1/23=−11fui+9fui+1+3fui+2−fui+312.
When we piece it together into higher order polynomials by linear and nonlinear weights, we have to solve the problem of negative linear weights and apply the mapping functions to obtain mapped nonlinear weights. It is not trivial to extend this scheme to compact scheme. After changing these three polynomials to compact polynomials, the corresponding linear weights do not exist. These difficulties can be easily solved by using the unequal-sized stencils that specified in this paper.Remark 2.
Since the compact reconstruction method needs to solve equation (22), its efficiency is slower than that of the traditional reconstruction. There has been relevant research in reference [32], and the efficiency will be about two times slower.


### 2.2. Time Discretization Method

After getting the value of flux, we use the third-order TVD Runge–Kutta time discretization method [36](25)u1=un+ΔtLun,u2=34un+14u1+14ΔtLu1,un+1=13un+23u2+23ΔtLu2,to discretize (2). The new CRUS-WENO scheme can also be easily applied for two-dimensional cases. For simplicity, we will not list the specific steps in detail.Remark 3.
For example, such asut+f1ux+f2uy=f3uxx+f4uyy, the time step is set as(26)Δt=CFLbx/Δx+cy/Δy+dx/Δx2+ey/Δy2,where(27)bx=maxuf1′u,cy=maxuf2′u,dx=maxuf3′u,ey=maxuf4′u,CFL=0.2.


## 2.1. Reconstruction of Numerical Flux

We first define a functionh (x), such that(4)fux=1Δx2∫Δx/2x+Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ,and obtain(5)fuxxx=hx+Δx−2hx+hx−ΔxΔx2.We only need to approximate the functionhxi+Δx−hxi with numerical fluxf^i+1/2 to achieve the required order of accuracy in smooth regions. According to [1], we obtain(6)f^i+1/2=−2fui−2+25fui−1−245fui+245fui+1−25fui+2+2fui+3180,(7)f^i+1/2=hxi+Δx−hxi+1560d7hxdx7|x=xiΔx7+OΔx8.Now a sixth-order compact interpolation is given by(8)215f^i−1/2+1115f^i+1/2+215f^i+3/2=−9fui−1−153fui+153fui+1+9fui+2180.Now, we prove that this formula has sixth-order accuracy and calculate its truncation error. From the definition of (4), we get(9)−9fui−1−153fui+153fui+1+9fui+2180=−9/Δx2∫x−3Δx/2x−Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ−153/Δx2∫x−Δx/2x−Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ+153/Δx2∫x−3Δx/2x−3Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ+9/Δx2∫x−3Δx/2x−5Δx/2∫τ−Δx/2τ+Δx/2hθdθdτ.180To analyze the accuracy of the numerical fluxf^i+1/2, we assume hx is smooth enough. Performing the Taylor series expansion for hx at the grid point xi(10)hθ=hxi+∑k=17θ−xikk!hkxi+OΔx8.By substituting the integral of the above formula, we get(11)f^i+1/2=∑k=16Δxkk!hkxi+145Δx77!h7xi+OΔx8.It is obvious that the Taylor series of the exact solution could expand to(12)hxi+Δx−hxi=∑k=16Δxkk!hkxi+Δx77!h7xi+OΔx8.Comparing (11) with (12), we can obtain the truncation error as(13)f^i+1/2=hxi+Δx−hxi+12800d7hxdx7|x=xiΔx7+OΔx8.We can see the coefficients of leading error term in compact interpolations are smaller than that of the noncompact interpolations. It can also be verified from the numerical examples that specified in the following. When the nonlinear degenerate parabolic equations contain the discontinuous solutions, these interpolation schemes might suffer from oscillations. So the high-order WENO schemes have been proposed as limiters to deal with this difficulty. The general form of the numerical flux reconstructed in the WENO fashion is(14)f^i+1/2=∑m=1rωmf^i+1/2m,where ωm is the nonlinear weights. The new way specified in this paper is to choose the unequal-sized stencils [28, 37, 38, 39]. We choose a six-point spatial stencil and two threepoint spatial stencils T1=xi−2,xi−1,xi,xi+1,xi+2,xi+3,T2=xi−1,xi,xi+1,T3=xi,xi+1,xi+2 and design the unequal degree polynomials defined on them, and then get their interpolations at half point xi+1/2:(15)f^i+1/21=−2fui−2+25fui−1−245fui+245fui+1−25fui+2+2fui+3180,f^i+1/22=fui+1−fui,f^i+1/23=fui+1−fui.Coincidentally, the approximations are the same on stencilT2 and T3. By the above mentioned formulas, we obtain the corresponding compact interpolations(16)215f^i−1/2+1115f^i+1/2+215f^i+3/2=−9fui−1−153fui+153fui+1+9fui+2180,12f^i−1/2+12f^i+1/2=fui+1−fui−12,12f^i+1/2+12f^i+3/2=fui+2−fui2.Similar to [38], three compact interpolations can be combined by the linear weights γ1,γ2,γ3 and the nonlinear weightsω1,ω2,ω3 as(17)ω11γ1215f^i−1/2+1115f^i+1/2+215f^i+1/2−γ2γ112f^i−1/2+12f^i+1/2−γ3γ112f^i+1/2+12f^i+1/2+ω212f^i−1/2+12f^i+1/2+ω312f^i+1/2+12f^i+3/2=ω11γ1−9fui−1−153fui+153fui+1+9fui+2180−γ2γ1fui+1−fui−12−γ3γ1fui+2−fui2+ω2fui+1−fui−12+ω3fui+2−fui2,in which, the linear weights and the nonlinear weights will be specified in the following. We use the similar definition [37] to construct the smoothness indicators:(18)β1=16048083708fui−22+2002064fui−12+7847144fui2+7847144fui+12−15310492fuifui+1+7205468fuifui+2−7738322fui+1fui+2+2002064fui+22+fui−2−791699fui−1+1457188fui−1308130fui+1+574412fui+2−99187fui+3−1308130fuifui+3+1457188fui+1fui+3−791699fui+2fui+3+83708fui+32+fui−1−7738322fui+7205468fui+1−3253987fui+2+574412fui+3,(19)β2=fui−1−2fui+fui+12,(20)β3=fui−2fui+1+fui+22.According to [40, 41], we define τ=β1−β3+β2−β3/22 Based on the linear weights and smoothness indicators, we define the nonlinear weights as(21)ωn=ω¯n∑l=13ω¯l,ω¯n=γn1+τε+βn,n=1,2,3.in which ε is a small positive number to avoid the denominator to become zero. We take ε = 10−6 in all computations. By performing such new spatial reconstruction procedure, we can achieve the sixth-order accuracy in smooth regions and preserve the essentially nonoscillatory property nearby the strong discontinuities.Through (17), the numerical flux value at half points can be directly computed by solving the following tridiagonal system(22)βγ0.........0αβγ.........00.........αβγ0.........0αβf^1/2f^i+1/2f^N−1/2=R^1/2−αf^1/2R^i+1/2R^N−1/2−γf^N−1/2,where R^i+1/2 is the right-hand side of (17), the coefficients α,β, and γ represent the coefficients on the left-hand side of 17. To solve this system of equations, we need to find the two unknown variables f^−1/2 andf^N+1/2 on the right-hand of 22, respectively. The values on the first and last half points along each grid line are reconstructed by the sixth-order US-WENO spatial reconstruction procedure [39]:(23)f^−1/2=ω1180γ1−2fu−3+25fu−2−245fu−1+245fu0−25fu1+2fu2+ω2−ω1γ1γ2+ω3−ω1γ1γ3fu0−fu−1,f^N+1/2=ω1180γ1−2fuN−2+25fuN−1−245fuN+245fuN+1−25fuN+2+2fuN+3+ω2−ω1γ1γ2+ω3−ω1γ1γ3fuN+1−fuN.The linear weights and nonlinear weights are calculated according to (21). Thus we can directly solve the three-diagonal linear equation system (22) by using the LU decomposition method.Remark 1.
The sixth-order finite difference WENO-JS scheme [1] uses the information of three equal-sized spatial stencils to construct three quadratic polynomials. Then their approximations are obtained at half point xi+1/2:(24)f^i+1/21=fui−2−3fui−1−9fui+11fui+112,f^i+1/22=fui−1−15fui+15fui+1−fui+212,f^i+1/23=−11fui+9fui+1+3fui+2−fui+312.
When we piece it together into higher order polynomials by linear and nonlinear weights, we have to solve the problem of negative linear weights and apply the mapping functions to obtain mapped nonlinear weights. It is not trivial to extend this scheme to compact scheme. After changing these three polynomials to compact polynomials, the corresponding linear weights do not exist. These difficulties can be easily solved by using the unequal-sized stencils that specified in this paper.Remark 2.
Since the compact reconstruction method needs to solve equation (22), its efficiency is slower than that of the traditional reconstruction. There has been relevant research in reference [32], and the efficiency will be about two times slower.


## 2.2. Time Discretization Method

After getting the value of flux, we use the third-order TVD Runge–Kutta time discretization method [36](25)u1=un+ΔtLun,u2=34un+14u1+14ΔtLu1,un+1=13un+23u2+23ΔtLu2,to discretize (2). The new CRUS-WENO scheme can also be easily applied for two-dimensional cases. For simplicity, we will not list the specific steps in detail.Remark 3.
For example, such asut+f1ux+f2uy=f3uxx+f4uyy, the time step is set as(26)Δt=CFLbx/Δx+cy/Δy+dx/Δx2+ey/Δy2,where(27)bx=maxuf1′u,cy=maxuf2′u,dx=maxuf3′u,ey=maxuf4′u,CFL=0.2.


## 3. Numerical Tests

In this section, we intend to compare the sixth-order finite difference CRUS-WENO scheme with the sixth-order classical WENO-JS scheme narrated in [1]. In order to test the effect of the linear weights on the accuracy of the new CRUS-WENO scheme, we take three different types of the linear weights in the numerical accuracy tests as: 1γ1=0.98,γ2=0.01,γ3=0.01; 2γ1=1/3,γ2=1/3,γ3=1/3; 3γ1=0.01,γ2=0.495,γ3=0.495. Then we recover γ1=0.98,γ2=0.01,γ3=0.01 in the latter numerical examples, unless specified otherwise.Example 1.
We compute the one-dimensional heat equation.(28)−π≤x≤π,ut=uxx,ux,0=sinx.
The exact solution isux,t=e−tsinx. We use the sixth-order finite difference CRUS-WENO scheme to simulate it with a periodic boundary condition. The errors are computed according to L1=1/N+1∑i=0Nui−uxi,t and L∞=max0≤i≤Nui−uxi,t, in which ui is the numerical approximation of uxi,t at the final time t = 2. We can find a fact in Table 1 that the CRUS-WENO scheme with different types of the linear weights and WENO-JS scheme could achieve their designed order of accuracy. For different types of the linear weights, the truncation errors of the CRUS-WENO scheme are similar on the final mesh level. Table 1 also shows a fact that the CRUS-WENO scheme produces smaller truncation errors than classical WENO-JS scheme does.Table 1
Example1. Heat equation. L1 and L∞ errors.
NL1 errorOrderL∞ errorOrderNL1 errorOrderL∞ errorOrderWENO-JSCRUS-WENO (1)403.66E − 095.91E − 09401.08E − 102.02E − 10603.24E − 105.985.17E − 106.00603.67E − 128.331.18E − 117.00805.79E − 115.989.22E − 116.00804.96E − 137.761.57E − 127.001001.52E − 115.982.42E − 116.001001.37E − 135.773.30E − 137.001205.12E − 125.978.12E − 125.981203.91E − 146.889.21E − 147.00CRUS-WENO (2)CRUS-WENO (3)406.35E − 109.61E − 10408.99E − 101.38E − 09605.76E − 1211.601.18E − 1110.85607.02E − 1211.971.78E − 1111.75804.80E − 138.631.57E − 127.00804.75E − 139.361.57E − 127.001001.35E − 135.673.30E − 137.001001.35E − 135.653.30E − 137.001203.89E − 146.849.21E − 147.001203.89E − 146.819.22E − 146.99Example 2.
We solve the two-dimensional linear heat equation.(29)x,y∈−π,π×−π,π,ut=uxx+uyy,ux,0=sinx.with the periodic boundary conditions in two directions. This equation has a global smooth solution ux,y,t=e−2tsinx+y. The L1 and L∞ errors and numerical orders of accuracy computed by the CRUS-WENO scheme and the WENO-JS scheme are shown in Table 2. We can find a fact that the CRUS-WENO scheme with different types of the linear weights and the WENO-JS scheme can achieve sixth-order accuracy. The CRUS-WENO scheme generates smaller absolute truncation errors than that of the WENO-JS scheme.Table 2
Example2. 2D heat equation. L1 and L∞ errors.
N × NL1 errorOrderL∞ errorOrderL1 errorOrderL∞ errorOrderWENO-JSCRUS-WENO (1)40 × 406.31E − 109.93E − 102.96E − 118.53E − 1160 × 605.48E − 116.038.62E − 116.031.39E − 127.543.33E − 128.0080 × 809.71E − 126.021.53E − 126.021.93E − 136.875.28E − 136.40100 × 1002.54E − 126.013.99E − 126.014.54E − 146.491.30E − 136.29120 × 1208.45E − 136.041.33E − 126.041.26E − 147.013.61E − 147.01CRUS-WENO (2)CRUS-WENO (3)40 × 408.91E − 111.87E − 101.26E − 102.46E − 1060 × 601.43E − 1210.203.68E − 129.691.46E − 1211.003.92E − 1210.2180 × 801.92E − 136.975.18E − 136.821.91E − 137.075.13EE − -137.07100 × 1004.52E − 146.471.29E − 136.224.52E − 146.461.29E − 136.18120 × 1201.26E − 147.003.61E − 147.011.26E − 146.993.60E − 147.00Example 3.
Now, we simulate the porous medium equation (PME) [42], which is a classic nonlinear degradation equation. This equation often has no exact solution, but there is a famous weak Barenblatt solution [43, 44]. We set the Barenblatt solution at t = 1 as the initial condition. The numerical solutions on 200 grid points for m = 2, 3, 5 and 8 at the final time t = 2 are shown in Figure 1, respectively. The solid line is the exact solution. It is observed that the numerical results of two different WENO schemes basically coincide with the exact solution and there are no obvious oscillations nearby the strong discontinuous in Figure 1.Figure 1
Example3. Barenblatt solution for PME. T = 2.200 grid points. (a) m = 2, (b) m = 3, (c) m = 5, and (d) m = 8.
(a)(b)(c)(d)Example 4.
We consider the PME with two discontinuities in the initial values. This equation is often used to describe the appearance of two hot spots which are suddenly placed in two domains and the variableu in this equation is usually expressed as the temperature. We set the initial condition as(30)ux,0=1,−4≤x≤−1,2,0≤x≤3,0,otherwise.
The boundary condition is set asu (±6, t) = 0 and the computational domain is chosen as [−6, 6] which is divided into 160 grid points. In Figure 2, we show the changes of the numerical solution for the porous medium equation with m = 6 at different times. We can see that the evolution of the numerical solution completely conforms to the results in the literature [45], and there are basically no oscillations near strong discontinuities. To see the advantages of our scheme more clearly, we have plotted the errors at t = 0.04 in Figure 3. We can see that the error of the CRUS-WENO scheme is almost everywhere smaller than that of the WENO-JS scheme.Figure 2
Example4. Collision of two-box solution for PME. CRUS-WENO scheme (a), WENO-JS scheme (b).
(a)(b)Figure 3
Example4. Collision of two-box solution for PME. Errors of WENO-JS scheme and CRUS-WENO scheme with 160 grid points. T = 0.04. The reference solution is computed by WENO-JS scheme with 1600 grid points.Example 5.
Next, we compute the convection-diffusion Buckley–Leverett equation [46].(31)ut+fux=ενuuxx,ενu≥0.
This equation is mainly used for the estimation of development indicators and two-phase flow simulation in the early stage of oilfield development. Here we takeε = 0.01 and(32)νu=4u1−u,0,0≤u≤1,otherwise.
We first consider the situation without the influence of gravity(33)fu=u2u2+1−u2.
And, the continuous initial condition is(34)ux,0=1−3x,0≤x≤13,0,13≤x≤1.
For boundary conditions, the left boundary is set asu0,t=1 and the right boundary is extrapolated. Then, we compute the discontinuous Riemann initial condition(35)ux,0=0,0≤x≤1−12,1,1−12≤x≤1,
For boundary conditions, the right boundary is set asu1,t=1 and u0,t=0, respectively. Finally, we compute the Buckley–Leverett equation with gravitation at the Riemann initial condition. Due to the influence of gravity, the convection term changes to(36)fu=u2u2+1−u21−51−u2.
The convection term in the equation is discretized by using the WENO scheme [47], and the diffusion term is discretized by using the CRUS-WENO and WENO-JS schemes, respectively. Figure 4 shows all numerical results with different grids at t = 0.2. It is observed that there are no oscillations and the solution matches the results in [46] and [28] well in all cases. We also draw the errors of the scheme in Figure 5, and the errors of the CRUS-WENO scheme is relatively small.Figure 4
Example5. (a) shows the initial boundary value problem without gravity for CRUS-WENO scheme. (b) and (c) show the results of the Riemann problem without gravity and with gravity for CRUS-WENO scheme, respectively. (d) shows the initial boundary value problem without gravity for WENO-JS scheme. (e) and (f) show the results of the Riemann problem without gravity and with gravity for WENO-JS scheme, respectively.
(a)(b)(c)(d)(e)(f)Figure 5
Example5. Errors of WENO-JS scheme and CRUS-WENO scheme with 100 grid points. (a) shows the initial boundary value problem without gravity. (b) and (c) show the Riemann problem without gravity and with gravity, respectively. The reference solution is computed by WENO-JS scheme with 800 grid points.
(a)(b)(c)Example 6.
We continue to consider the previous equation(37)ut+fux=ενuuxx.
Whenνu is a discontinuous function, the equation becomes a strongly degenerate parabolic equation. We take ε = 0.1, f(u) = u2, and(38)νu=0,u≤0.25,1,u>0.25.
It causes the equation to be hyperbolic whenu ∈ [−0.25, 0.25] and parabolic elsewhere. We apply the sixth-order finite difference CRUS-WENO scheme to compute it with the initial.(39)ux,01,−12−0.4<x<−12+0.4,−1,12−0.4<x<12+0.4,0,otherwise.
The results on two different grids with the zero boundary conditionu±2,t=0 at time t = 0.7 are shown in Figure 6. Numerical results show the good performance of the high resolution and the accurate transition between the hyperbolic domain and the parabolic domain. It can be seen from the error Figure 7 that the maximum error of the CRUS-WENO scheme is smaller than that of the WENO-JS scheme.Figure 6
Example6. Riemann problem. T = 0.7. CRUS-WENO scheme (a); WENO-JS scheme (b).
(a)(b)Figure 7
Example6. Riemann problem. T = 0.7. Errors of WENO-JS scheme and CRUS-WENO scheme with 200 grid points. The reference solution is computed by WENO-JS scheme with 800 grid points.Example 7.
We compute the two-dimensional PME(40)ut=u2xx+u2yy,x,y∈−10,10×−10,10,with initial condition(41)ux,y,0=e−1/6−x−22−y+22x−22+y+22<6,e−1/6−x−22−y+22x+22+y−22<6,0otherwise.
And the periodic boundary conditions are applied in two directions. The numerical results of the CRUS-WENO and WENO-JS schemes at timet = 0.5, 1.0, and 4.0 are shown in Figure 8. It is observed that the sixth-order CRUS-WENO scheme captures the sharp interface very well without introducing noticeable oscillations near strong discontinuities. And associated numerical results coincide with that specified in [47]. We plot the errors of the two schemes at t = 0.5 along x + y = 0 in Figure 9, it is clear that the errors of the CRUS-WENO format are much smaller.Figure 8
Example7. 2D PME. CRUS-WENO scheme (top); WENO-JS scheme (bottom). 80 × 80 grid points.
(a)(b)(c)(d)(e)(f)Figure 9
Example7. Riemann problem. Errors of WENO-JS scheme and CRUS-WENO scheme with 200 × 200 grid points along x + y = 0. T = 0.5. The reference solution is computed by WENO-JS scheme with 800 × 800 grid points.Example 8.
We consider the two-dimensional Buckely–Leverett equation [48].(42)ut+f1ux+f2uy=εuxx+uyy,in which ε = 0.01. The flux functions are defined as f1u=u2/u2+1−u2 and f2u=u2−51−u2u2/u2+1−u2. The computation domain is [−1.5, 1.5] × [−1.5, 1.5], which is divided into 200 × 200 uniform grid points. We show the results computed by the CRUS-WENO scheme in Figure 10 with the initial condition(43)ux,0=0,x2+y2<0.5,1,otherwise.
This numerical results are in line with the results that specified in [47]. It is shown that the new sixth-order finite difference CRUS-WENO scheme is also very effective for this two-dimensional problem. We plot the errors of the two schemes at t = 0.5 along x = 0 in Figure 11, it is clear that the errors of the CRUS-WENO scheme are much smaller almost everywhere.Figure 10
Example8. 2D PME. T = 0.5. CRUS-WENO scheme (a, c); WENO-JS scheme (b, d). 200 × 200 grid points.
(a)(b)(c)(d)Figure 11
Example8. 2D PME. T = 0.5. Errors of WENO-JS scheme and CRUS-WENO scheme with 200 × 200 grid points along x = 0. The reference solution is computed by WENO-JS scheme with 800 × 800 grid points.Example 9.
Finally, we consider the two-dimensional strongly degenerate parabolic convection-diffusion equation(44)ut+fux+fuy=ενuuxx+ενuuyy,ενu≥0.with ε = 0.1, f (u) = u2, and vu are given in (30). The initial value is(45)ux,y,0=1,x−0.52+y+0.52<0.16,−1,x+0.52+y−0.52<0.16,0,otherwise.
The numerical solutions are shown in Figure12 at t = 0.5. And the computation domain [−1.5, 1.5] × [−1.5, 1.5] is divided into 200 × 200 uniform grid points. It is obviously observed that the new sixth-order finite difference CRUS-WENO scheme could work well for this two-dimensional test case.Figure 12
Example9. 2D strong parabolic equation. T = 0.5. CRUS-WENO scheme (a, c); WENO-JS scheme (b, d). 200×200 grid points.
(a)(b)(c)(d)


## 4. Concluding Remarks

In this paper, a new sixth-order finite difference compact reconstruction unequal-sized WENO (CRUS-WENO) scheme is designed for solving the nonlinear degenerate parabolic equations on structured meshes. Comparing with the same order finite difference classical WENO schemes, its major advantage is the simplicity and compactness in designing the highorder spatial reconstruction methodology. Its merits are specified in the following. The first is that we only perform a new spatial WENO reconstruction procedure for the second-order derivative values, which is very different to many high-order WENO schemes. The second is that the linear weights of this CRUS-WENO scheme can be any positive numbers as long as their summation is one. The third is that the new finite difference CRUS-WENO scheme uses the narrower large spatial stencil than the same order classical WENO schemes do. The fourth is that only three unequal-sized spatial stencils are applied in spatial reconstruction. It has a better compactness and smaller truncation errors in smooth regions. The fifth is that the new scheme will not introduce the negative linear weights and the mapping functions in spatial reconstructions in comparison with that specified in [1]. The sixth is that the new scheme could obtain the optimal sixth-order accuracy in smooth regions, and preserve the second-order accuracy nearby the strong discontinuities.



---

*Source: 1627069-2022-08-25.xml*