Hybrid Spectral Difference Methods for an Elliptic Equation (original) (raw)

1 Introduction

Over the last two decades, high-order methods have received considerable attention because of their potential in delivering higher accuracy with lower cost than low-order methods. Many types of high-order methods have been developed to deal with a diverse range of problems [26]. See also conference proceedings on spectral and high-order methods [18].

Finite difference (FD) methods are amongst the most popular methods solving differential equations in science and engineering. In an effort to develop high-order methods, compact finite difference schemes are proposed for various problems [19, 9, 24, 23]. FD methods have low geometric flexibility, in particular, for high-order schemes.

It is well known that finite element formulations allow high-order approximations with high geometric flexibility. But, their implementation is not so straightforward. Recently, nodal discontinuous Galerkin (DG) methods and hybridized DG methods have been proposed to enhance efficiency of the DG schemes. Nodal DG methods [8] presented by Hesthaven and Warburton provide an efficient way of implementing high-order DG methods for various problems including Poisson equations and Maxwell equations. Hybridized DG methods [3, 12, 13, 14, 15] exploit Lagrange multipliers for hybridization to relax constraints imposed on the cell-interfaces, and static condensation via Schur complement allows efficient implementation while keeping most of the advantages of DG methods. Recently, hybrid high-order methods [5, 4, 2] have been presented in the context of HDG methods.

On the other hand, in order to develop high-order conservative schemes, Wang and his collaborators proposed and studied so-called spectral volume (SV) methods [25, 22] and spectral difference (SD) methods [21]. These methods have been successfully applied to, but limited to, mostly hyperbolic conservation laws. The SV method is designed as a high-order finite volume scheme by considering some local reconstruction idea in the DG schemes. The SD method utilizes a local pseudospectral representation of the solution to obtain spectral-like resolutions. The SD method first defines two sets of points, i.e., the solution points and flux points in each cell. Then, the conservative variables defined at the solution points provide the DOFs. To evolve these DOFs, the differential form of the governing equation is applied at solution points and the divergence of flux is evaluated in terms of values at flux points, cf. [27]. Note that placements of these points are staggered to enhance stability of the approximation to hyperbolic conservation laws. Stability of SD methods is investigated in [1, 10].

In this article, we develop and analyze a locally conservative hybrid spectral difference method (HSD) for the Poisson equation. Similar to the standard SD method, we use two sets of points, i.e., the solution points and flux points in each cell. But, placements of these points are chosen to be axiparallel unlikely to the standard SD method and no grid staggering is needed. Moreover, in our scheme, the differential form of the governing equation is applied at interior solution points to provide the cell finite difference, while flux points are related to the interface finite difference and provide global DOFs. After static condensation built in the procedure, a global stiffness matrix system is obtained in the interface unknowns only. The HSD aims to provide an alternative for the finite difference method and the finite element method for partial differential equations of elliptic type.

Figure 1

A Q2-mesh: |R1|=h1×k1,|R2|=h2×k1.

For simplicity, we present our scheme for the Poisson model problem with the Dirichlet boundary condition:

{-Δ⁢u=fon ⁢Ω,u=0on ⁢Γ=∂⁡Ω,

where Ω is an axiparallel polygonal domain. Let 𝒯h be a partition of Ω by rectangular cells. The simplest HSD is briefly explained, based on Figure 1. The hybrid difference is composed of two kinds of finite difference approximations as follows. The 5-point stencil FD approximation of the Poisson equation, for example, at the cell centered point η4 is the cell finite difference:

(1.1)-Δh⁢u⁢(η4):=-u⁢(η3)+2⁢u⁢(η4)-u⁢(η5)(h1/2)2+-u⁢(η8)+2⁢u⁢(η4)-u⁢(η1)(k1/2)2=f⁢(η4),

and the FD approximation of the flux continuity at η5 is the interface finite difference:

[[∂νh⁡u]]η5=3⁢u⁢(η5)-4⁢u⁢(η4)+u⁢(η3)h1+3⁢u⁢(η5)-4⁢u⁢(η6)+u⁢(η7)h2=0.

Static condensation on cell-centered unknowns reduces the global couplings, resulting in the system of equations in the cell interface unknowns only.

Numerical experiments on high-order hybrid difference methods (hybrid spectral difference) for elliptic and flow problems, including the Poisson equation, Stokes and Navier–Stokes equations, can be found in [11, 16]. The HSD can be viewed as a finite difference version of the hybridized discontinuous Galerkin method [3, 12, 13, 15]. The HSD is comparable with the finite difference method (FDM). The main difference between the FDM and HSD is that the FD formula of a single type is deployed for all interior nodes in the FDM, while the cell finite difference and the interface finite difference are combined in the HSD. The finite difference method is simple to implement and it can solve many physical problems efficiently [6, 20]. The HSD is as easy to implement as the FDM, and it apparently seems to possess several advantages over the FDM [11]. Those advantages are listed below.

1. The method can be applied to nonuniform grids, retaining the optimal order of convergence, and numerical methods with an arbitrarily high-order convergence can be obtained easily.
2. Problems on a complicated geometry can be treated reasonably well, and the boundary condition can be imposed exactly on the exact boundary.
3. The inf-sup stability is obtained without introducing staggered grids [16].
4. The mass conservation property holds in each cell, and flux continuity holds across inter-cell boundaries.
5. The embedded static condensation property considerably reduces global degrees of freedom.

The aim of the current paper is to propose and analyze the HSD for the Poisson equation. Numerical analysis of the Poisson equation is a very first step toward theoretical understanding of the more complicated PDEs that arise from flow problems or elasticity problems. There can be a couple of different approaches in deriving the FD formulas; one based on the Lagrange polynomial interpolation and the other based on the polynomial degrees of precision. Both approaches yield the same FD formulas for the same derivatives. The latter approach was previously used in the study of pseudospectral (PS) methods. Indeed, the PS method can be seen either as the limit of increasing order FD methods, or as approximations by basis functions, such as Fourier or Chebyshev [7].

The paper is organized as follows. In Section 2 we derive one-dimensional finite difference formulas based on the polynomial degrees of precision. In Section 3, a one-dimensional coercivity result is derived in a semi-discrete energy norm. In Section 4 we extend the coercivity result in Section 3 to the two-dimensional case. In Section 5 numerical experiments are performed for the Poisson equations and numerical results confirm our numerical analysis. Moreover, the Poisson equation on a quarter disk is considered for numerical experiments, and the HSDs provide quite good numerical results. Section 6 concludes the paper with some remarks.

2 Finite Difference Formulas

We derive one-dimensional finite difference formulas in terms of degrees of precision.

Let us introduce an increasing sequence of abscissas {η0,η1,…,ηm} in [a,a+h]. For μ=1,2,…,m, consider the problem to find(wj⁢k(μ))j=0,…,m;k=0,…,m such that

(2.1)1hμ∑j=0m(ηj-η0)ℓwj⁢k(μ)=dμ⁢(x-η0)ℓd⁢xμ|x=ηk,ℓ=0,…,m

for k=0,…,m. For fixed μ and k, equation (2.1) has a unique solution(wj⁢k(μ))j=0,…,m as system (2.1) forms a Vandermonde matrix system. In particular, it follows from (2.1) that

(2.2)∑j=0m(ηj-η0)ℓwj⁢k(μ)={0,ℓ=0,…,μ-1,hμ⁢ℓ!(ℓ-μ)!⁢(ηk-η0)ℓ-μ,ℓ=μ,…,m,

for k=0,…,m.

Definition 2.1

Based on(wj⁢k(μ))j=0,…,m;k=0,…,mas in (2.1), define a class of general μ-th order finite difference operators(Dxh)μ, μ=1,…,m, as follows:

(Dxh)μ⁢f⁢(ηk)=1hμ⁢∑j=0mf⁢(ηj)⁢wj⁢k(μ),0≤k≤m,for all ⁢f∈C⁢[a,a+h].

Theorem 2.2

Let f∈Cm+1⁢[a,a+h] and (Dxh)μ, μ=1,2,…,m be given as in Definition 2.1. Then,

|(Dxh)μ⁢f⁢(ηk)-dμd⁢xμ⁢f⁢(ηk)|≲hm-μ+1⁢∥f(m+1)∥L∞⁢[a,a+h] for all ⁢k=0,…,m.

Throughout the paper, the notation L≲R means that there exists a constant C>0, independent of h, such that L≤C⁢R.

Proof.

Suppose that f∈Cm+1⁢[a,a+h] is given and 1≤μ≤m is fixed. Then, by the Taylor expansion,

(2.3)f⁢(ηj)=f⁢(η0)+f′⁢(η0)⁢(ηj-η0)+⋯+f(m)⁢(η0)m!⁢(ηj-η0)m+f(m+1)⁢(γj(0))(m+1)!⁢(ηj-η0)m+1

for some γj(0)∈[a,a+h] and all j=0,1,…,m. In the mean while,

(2.4)f(μ)⁢(ηk)=f(μ)⁢(η0)+f(μ+1)⁢(η0)⁢(ηk-η0)+⋯+f(m)⁢(η0)(m-μ)!⁢(ηk-η0)m-μ+f(m+1)⁢(γk(μ))(m-μ+1)!⁢(ηk-η0)m-μ+1,

for some γk(μ)∈[a,a+h] for all k=0,1,…,m. Multiplying by wj⁢k(μ) term by term in (2.3) and summing on j from 0 to m, with (2.2) and (2.4) invoked, one gets

∑j=0mf⁢(ηj)⁢wj⁢k(μ)=f⁢(η0)⁢∑j=0mwj⁢k(μ)+f′⁢(η0)⁢∑j=0m(ηj-η0)⁢wj⁢k(μ)+⋯+f(m)⁢(η0)m!⁢∑j=0m(ηj-η0)m⁢wj⁢k(μ)

+∑j=0mwj⁢k(μ)⁢f(m+1)⁢(γj(0))(m+1)!⁢(ηj-η0)m+1

=hμ⁢[f(μ)⁢(η0)+f(μ+1)⁢(η0)⁢(ηk-η0)+⋯+f(m)⁢(η0)(m-μ)!⁢(ηk-η0)m-μ]

+∑j=0mf(m+1)⁢(γj(0))(m+1)!⁢(ηj-η0)m+1⁢wj⁢k(μ)

=hμ⁢(f(μ)⁢(ηk)+Ek),

where

Ek=1hμ⁢∑j=0mf(m+1)⁢(γj(0))(m+1)!⁢(ηj-η0)m+1⁢wj⁢k(μ)-f(m+1)⁢(γk(μ))(m-μ+1)!⁢(ηk-η0)m-μ+1.

Since

|Ek|≲hm-μ+1⁢∥f(m+1)∥L∞⁢[a,a+h],

the theorem follows immediately. ∎

3 One-Dimensional Coercivity Analysis

Consider m+1 abscissas {ξ0,ξ1,…,ξm} with ξ0=-1 and ξm=1 on the reference interval [-1,1] and the corresponding abscissas {η0,η1,…,ηm} on an interval I=[a,a+h] that are given as ηj=ϕ⁢(ξj), where ϕ⁢(x)=h⁢(x+1)/2+a.

Let us introduce a quadrature on the function space C⁢([a,a+h]):

Qh⁢(f)=h2⁢∑j=1m-1σj⁢f⁢(ηj)=h2⁢∑j=1m-1σj⁢f~⁢(ξj).

Here, f~=f⁢(ϕ-1). From here and on, the weights {σj}j=1m-1 and the interior abscissas {ξj}j=1m-1 are the Gaussian weights and abscissas on the reference interval [-1,1], respectively. Therefore, the Gaussian quadrature has the polynomial degree of precision 2⁢m-3.

Let us introduce a bilinear form:

(3.1)𝒜h⁢(u,v)=-(Dx⁢xh⁢u,v)I,h+〈∂νh⁡u,v〉∂⁡I,u,v∈C⁢(I),

where

(u,v)I,h=Qh⁢(u⁢v),〈∂ν⁡f,g〉∂⁡I=fx⁢(ηm)⁢g⁢(ηm)-fx⁢(η0)⁢g⁢(η0).

Since the bilinear form (3.1) is defined by using the function values of u and v at the given m+1 abscissas on an interval I, we can regard u and v as interpolating polynomials of degree m, that is, u,v∈Pm⁢(I).

Theorem 3.1

Theorem 3.1 (Symmetry)

For u,v∈Pm⁢(I) the following symmetry holds:

𝒜h⁢(u,v)=𝒜h⁢(v,u).

Proof.

For u,v∈Pm⁢(I), we have u=pm-1+cm⁢xm and v=qm-1+dm⁢xm with pm-1,qm-1∈Pm-1⁢(I). Notice that ∂ν⁡u=∂νh⁡u at η0, ηm, andDx⁢x⁢u=Dx⁢xh⁢u at ηj (j=1,…,m-1) for u∈Pm⁢(I). The integration by parts yields

-(Dx⁢xh⁢u,v)I,h+〈∂νh⁡u,v〉∂⁡I=-(Dx⁢xh⁢u,v)I,h+(ux,vx)I+(Dx⁢x⁢u,v)I.

Since the quadrature has a degree of precision 2⁢m-3, we get

-(Dx⁢xh⁢u,v)I,h+(Dx⁢x⁢u,v)I=m⁢(m-1)⁢cm⁢dm⁢(∫Ix2⁢m-2⁢𝑑x-Qh⁢(x2⁢m-2)).

Therefore,

𝒜h⁢(u,v)=(ux,vx)I+m⁢(m-1)⁢cm⁢dm⁢(∫Ix2⁢m-2⁢𝑑x-Qh⁢(x2⁢m-2)).

By the same way,

𝒜h⁢(v,u)=(vx,ux)I+m⁢(m-1)⁢cm⁢dm⁢(∫Ix2⁢m-2⁢𝑑x-Qh⁢(x2⁢m-2)).∎

Theorem 3.2

Theorem 3.2 (Coercivity)

Suppose

∫Ix2⁢m-2⁢𝑑x-Qh⁢(x2⁢m-2)≥0.

Then, the following coercivity holds:

(3.2)𝒜h⁢(u,u)≥∥ux∥L2⁢(I)2,u∈Pm⁢(I).

Proof.

For u∈Pm⁢(I), u=pm-1+cm⁢xm we have

𝒜h⁢(u,u)=(ux,ux)I+m⁢(m-1)⁢cm2⁢(∫Ix2⁢m-2⁢𝑑x-Qh⁢(x2⁢m-2)).

Equation (3.2) is obtained immediately. ∎

Theorem 3.3

It holds that

∫abx2⁢m-2⁢𝑑x-Qh⁢(x2⁢m-2)=(b-a)2⁢m-1⁢[(m-1)!]4(2⁢m-1)⁢[(2⁢m-2)!]2>0,m≥2.

Proof.

The proof can be found in [17]. ∎

As a conclusion the operator 𝒜h is symmetric and elliptic for each order m as long as the interior abscissas are chosen so that the corresponding Gaussian quadrature has the degree of precision 2⁢m-3.

4 Numerical Analysis in Two Dimensions

Let Ω be a rectangular domain with the boundary Γ. In the hybrid difference method we need a rectangular partition 𝒯h of the domain Ω. The partition 𝒯h is shape regular and quasi uniform, and the parameter h represents the maximum diameter of the partition. The notation Kh represents the skeleton of a mesh generation 𝒯h. To define the hybrid difference method we need to generate nodes in each cell, and Figure 2 represents the reference cell configuration (a two-dimensional version of the five abscissa case). In Figure 2 the abscissas {ηi⁢j}i,j=0,…,m (here, m=4) are represented and the nodes {ηi⁢j}i,j=0,m are introduced as virtual abscissas (for convenience of presentation). Let 𝒩⁢(Ω) denote the set of nodes in the closure of a domain and 𝒩⁢(Kh) the set of nodes on its skeleton.

Let C⁢(Ω) be the space of continuous functions in Ω and introduce an equivalence relation u≡v for u,v∈C⁢(Ω) ifu⁢(η)=v⁢(η) for all η∈𝒩h⁢(Ω)∪𝒩h⁢(Γ). Denote by Ch⁢(Ω) the space of its equivalent classes. By using the equivalences of the function spaces we regard C~h⁢(R)=Qm⁢(R), where Qm⁢(R) is the space of polynomials in R of degree less than or equal to m in each variable. Then we define

Qm(Ω)={u∈C(Ω):u|R∈Qm(R)}.

Figure 2

A Q4-mesh: |R|=h1×h2. The point values at the four vertices {ηj⁢k:i,j=0,4} are not used in computation.

Let us introduce the FD approximations for the Laplacian and normal derivatives:

Δh⁢u=Dx⁢xh1⁢u+Dy⁢yh2⁢u,∂νh⁡u=(Dxh1⁢u,Dyh2)⋅ν.

Here, the superscript h in Δh and ∂νh represents the level of discretization as well. The Gaussian quadratures on the reference cell R and its boundary ∂⁡R are defined as follows:

(f,g)R,h=h1⁢h24⁢∑i=1m-1∑j=1m-1σi⁢σj⁢f⁢(ηi⁢j)⁢g⁢(ηi⁢j),

〈f,g〉∂⁡R,h=h22⁢∑j=1m-1σj⁢(f⁢(η0⁢j)⁢g⁢(η0⁢j)+f⁢(ηm⁢j)⁢g⁢(ηm⁢j))

+h12⁢∑j=1m-1σj⁢(f⁢(ηj⁢0)⁢g⁢(ηj⁢0)+f⁢(ηj⁢m)⁢g⁢(ηj⁢m⁢j)).

There exists a natural composite quadrature defined on the whole domain Ω as follows:

(f,g)h=∑R∈𝒯h(f,g)R,h.

Let us consider an elliptic problem with the Dirichlet boundary condition:

(4.1){-Δ⁢u=fon ⁢Ω,u=0on ⁢Γ=∂⁡Ω.

Then, the solution u satisfies the cell problem

-Δ⁢u=f on ⁢R

for each R∈𝒯h and it satisfies the interface condition

[[∂νu]]e=∂νu|e+∂ν′u|e=0 on e=∂R∩∂R′.

Here, ∂ν⁡u:=∂⁡u∂⁡ν and ν and ν′ are the outward unit normal vectors from R and R′, respectively.

The hybrid spectral difference method (HSD) is composed of two kinds of finite difference approximations. The cell finite difference solves the cell problem:

(4.2)-Δh⁢uh⁢(p)=f⁢(p),p∈𝒩⁢(Ω)∖𝒩⁢(Kh),

and the interface finite difference solves the interface condition

(4.3)[[∂νh⁡uh]]p=0,p∈𝒩⁢(Kh)∖Γ.

Of course, we require uh⁢(p)=0 for p∈𝒩⁢(Kh)∩Γ. The HSD, (4.2) and (4.3) can be unified in a variational form:

(4.4)-(Δh⁢uh,v)h+∑R∈𝒯h〈∂νh⁡uh,v〉∂⁡R,h=(f,v)h,v∈C0⁢(Ω).

It is easy to see that the exact solution satisfies

-(Δ⁢u,v)h+∑R∈𝒯h〈∂ν⁡u,v〉∂⁡R,h=(f,v)h,v∈C0⁢(Ω).

The continuous function space with zero trace is represented by C0⁢(Ω). Introducing eh=uh-u, we have

(4.5)-(Δh⁢eh,eh)h+∑R∈𝒯h〈∂νh⁡eh,eh〉∂⁡R,h=-(Δ⁢u-Δh⁢u,eh)h+∑R∈𝒯h〈∂ν⁡u-∂νh⁡u,eh〉∂⁡R,h.

The following theorem is essential for unique solvability of the HSD (4.4).

Theorem 4.1

For u∈Qm⁢(R),

-(Δ⁢u,u)R,h+〈∂νh⁡u,u〉∂⁡R,h≥h⁢∑j=1m-1(∥ux∥L2⁢(Ij)2+∥uy∥L2⁢(Jj)2),

whereIj=η0⁢j⁢ηm⁢j¯ and Jj=ηj⁢0⁢ηj⁢m¯. The notation x⁢y¯ represents the line from x to y.

Proof.

The quadrature 〈∂νh⁡u,u〉∂⁡R,h can be expressed as the sum of 2⁢(m-1) line components as follows:

〈∂νh⁡u,u〉∂⁡R,h=h22⁢∑j=1m-1σj⁢〈∂νh⁡u,u〉∂⁡Ij+h12⁢∑j=1m-1σj⁢〈∂νh⁡u,u〉∂⁡Jj.

By the one-dimensional ellipticity analysis in the previous section, we have

-(ux⁢x,u)Ij,h+〈∂νh⁡u,u〉∂⁡Ij≥(ux,ux)Ij,

-(uy⁢y,u)Jj,h+〈∂νh⁡u,u〉∂⁡Jj≥(uy,uy)Ij.

Note that

(Δ⁢u,u)R,h=h22⁢∑j=1m-1σj⁢(ux⁢x,u)Ij,h+h12⁢∑k=1m-1σk⁢(uy⁢y,u)Jk,h.

The theorem is immediate. ∎

Let us introduce the norms

∥u∥0,h2=∑R∈𝒯h∥u∥R,0,h2,∥u∥R,0,h2=(u,u)R,h

and the energy semi-norm

|u|E,h2=∑T∈𝒯h|u|E,R,h2,

where

(4.6)|u|E,R,h2=h⁢∑j=1m-1(∥ux∥L2⁢(Ij)2+∥uy∥L2⁢(Jj)2).

Lemma 4.2

Lemma 4.2 (Poincaré inequality)

Let Ω=⋃i,j=1NRi⁢j. For u∈Qm⁢(Ω),

∥u∥0,h≲|u|E,h.

Figure 3

_N_-stacks of rectangles Sj and lines Lk for Q3.

Proof.

Since Ω=⋃i,j=1NRi⁢j, there exist _N_-stacks of rectangles, and each stack is in the form Sj=⋃i=1NRi⁢j (see Figure 3). Each stack contains m-1 sets of m⁢N+1 aligned nodes, and suppose Nk={p0,k,p1,k,…,pm⁢N,k} is one of such an aligned node set. Note that u⁢(p0,k)=u⁢(pm⁢N,k)=0. Then,

|u⁢(pM,k)|≤∫p0,kpM,k|ux|⁢𝑑x≤∫Lk|ux|⁢𝑑x≲(∫Lk|ux|2⁢𝑑x)1/2,

where p0,k⁢pM,k¯⊂Lk=p0,k⁢pm⁢N,k¯. Therefore,

∥u∥R,0,h2≲h2⁢∑R∩Lk≠∅∥ux∥L2⁢(Lk)2.

Then,

∥u∥Sj,0,h2≲h⁢∑Sj∩Lk≠∅∥ux∥L2⁢(Lk)2.

Now, we have

∥u∥0,h2≲h⁢∑j=1N∑Sj∩Lk≠∅∥ux∥L2⁢(Lk)2≲|u|E,h2.∎

Theorem 4.3

Suppose u∈Cm+1⁢(Ω). Then, we have

(4.7)|(Δ⁢u-Δh⁢u,eh)h|≲hm⁢∥u(m+1)∥∞⁢|eh|E,h,

(4.8)|〈∂ν⁡u-∂νh⁡u,eh〉h|≲hm⁢∥u(m+1)∥∞⁢|eh|E,h.

Here, u(m+1) represents the sum of all partial derivatives of order m+1.

Proof.

For the error estimate, let us define the following norm:

∥f∥R,∞,h=max1≤i,j≤m-1⁡|f⁢(ηi⁢j)|.

Using the above definition, Theorem 2.2 implies that

∥Δ⁢u-Δh⁢u∥R,∞,h≤∥Dx⁢x⁢u-Dx⁢xh1⁢u∥R,∞,h+∥Dy⁢y⁢u-Dy⁢yh2⁢u∥R,∞,h

≲max1≤j≤m-1⁡h1m-1⁢∥u(m+1)∥L∞⁢(Ij)+max1≤j≤m-1⁡h2m-1⁢∥u(m+1)∥L∞⁢(Jj)

≲hm-1⁢∥u(m+1)∥L∞⁢(R).

Then simple calculation yields that

|(Δ⁢u-Δh⁢u,eh)R,h|≲∥Δ⁢u-Δh⁢u∥R,∞,h⁢(h1⁢h24⁢∑i,j=1m-1σi⁢σj⁢|eh⁢(ηi⁢j)|)

≲h⁢∥Δ⁢u-Δh⁢u∥R,∞,h⁢(h1⁢h24⁢∑i,j=1m-1σi⁢σj⁢|eh⁢(ηi⁢j)|2)1/2

≲hm⁢∥u(m+1)∥L∞⁢(R)⁢∥eh∥R,0,h.

Estimate (4.7) follows by the Poincaré inequality.

Now, we turn to the proof of estimate (4.8). Let w=∂ν⁡u-∂νh⁡u. Then,

〈w,eh〉∂⁡R,h=h22⁢∑j=1m-1σj⁢〈w,eh〉∂⁡Ij+h12⁢∑j=1m-1σj⁢〈w,eh〉∂⁡Jj.

We have

h2⁢〈w,eh〉∂⁡Ij=h2⁢(w⁢(ηm⁢j)⁢eh⁢(ηm⁢j)-w⁢(η0⁢j)⁢eh⁢(η0⁢j)),

and

|w⁢(ηm⁢j)⁢eh⁢(ηm⁢j)|≤|w⁢(ηm⁢j)⁢(eh⁢(ηm⁢j)-eh⁢(ηj0,j))|+|w⁢(ηm⁢j)⁢eh⁢(ηj0,j)|,

|w⁢(η0⁢j)⁢eh⁢(η0⁢j)|≤|w⁢(η0⁢j)⁢(eh⁢(η0⁢j)-eh⁢(ηj0,j))|+|w⁢(η0⁢j)⁢eh⁢(ηj0,j)|

for some 0<j0<m. Thus, we have

|h2⁢〈w,eh〉∂⁡Ij|≲maxi=0,m⁡|w⁢(ηi⁢j)|⁢(h3/2⁢∥Dx⁢eh∥L2⁢(Ij)+∥eh∥R,0,h).

Here, the first term of the right-hand side can be bounded:

maxi=0,m⁡|w⁢(ηi⁢j)|≤maxi=0,m⁡(|(Dx⁢u⁢(ηi⁢j)-Dxh⁢u⁢(ηi⁢j))⁢ν1|+|(Dy⁢u⁢(ηi⁢j)-Dyh⁢u⁢(ηi⁢j))⁢ν2|)

≤maxi=0,m⁡(|Dx⁢u⁢(ηi⁢j)-Dxh⁢u⁢(ηi⁢j)|+|Dy⁢u⁢(ηi⁢j)-Dyh⁢u⁢(ηi⁢j)|)

≲max1≤j≤m-1⁡h1m⁢∥u(m+1)∥L∞⁢(Ij)+maxj=0,m⁡h2m⁢∥u(m+1)∥L∞⁢(Jj)

≲hm⁢∥u(m+1)∥L∞⁢(R),

where ν1 and ν2 are the first and second component of ν. Then we have the following estimate by using (4.6) and Lemma 4.2:

|h2⁢〈w,eh〉∂⁡Ij|≲hm⁢∥u(m+1)∥L∞⁢(R)⁢(h⁢|eh|E,R,h+|eh|E,R,h).

By the same way, we have

|h1⁢〈w,eh〉∂⁡Jj|≲hm⁢∥u(m+1)∥L∞⁢(R)⁢(h⁢|eh|E,R,h+|eh|E,R,h).

Therefore, (4.8) is proved by combining the above two estimates:

|〈w,eh〉h|≲hm⁢∥u(m+1)∥∞⁢(h⁢|eh|E,h+|eh|E,h).∎

Theorems 4.1 and 4.3 with (4.5) result in the following convergence estimate.

Theorem 4.4

Let u∈Cm+1⁢(Ω) be the exact solution of (4.1) and uh the numerical solution of the HSD (4.4). Then,

|u-uh|E,h≲hm⁢∥u(m+1)∥∞.

5 Numerical Experiments

The uniform and geometric cell generations of the unit square Ω=(0,1)×(0,1) are considered for numerical experiments. A quarter disk domain is also considered to exhibit that the HSD can manage the curved boundary quite well. The geometric cell subdivision is obtained by considering a mesh transformation

(x~,y~)=(xτxτ+(1-x)τ,yτyτ+(1-y)τ),

where (x,y) is a uniform griding and τ is a positive integer, which is called the order of mesh transform.

Figure 4

Numerical results with respect to the h-refinement (left) and p-refinement (right) for the uniformgrid cell with Example 5.1.

Example 5.1

Example 5.1 (Smooth Solution)

Consider the Poisson equation

{-Δ⁢u=fin ⁢Ω,u=gon ⁢Γ,

where f and g are found from the exact solution u⁢(x,y)=exp⁡(x)⁢cos⁡(y)+(x2+y2). In this example the uniform grid on the unit square is considered.

Figure 4 shows the convergence history for Example 5.1 with respect to both the h-refinement and p-refinement. Noting that _degrees of freedom_≈h-2 for the uniform mesh, the orders of convergence concord well with the theoretical result in Theorem 4.4. For the 3-abscissa subgrid case the L2 and H1 convergence are of the same second order. For the other cases the order of the L2 convergence is one higher than that of the H1 convergence. The L2 error analysis will be a subject of future research. In contrast to the results with respect to the h-refinement, the convergence history for the p-refinement shows the possibility to achieve exponential convergence rates in terms of the number of degrees of freedom.

Figure 5

Numerical results with respect to the h-refinement (left) and p-refinement (right) for the uniformgrid cell with Example 5.2.

Figure 6

Numerical results with respect to the h-refinement (left) and p-refinement (right) for the geometricgrid (τ=3) cell with Example 5.2.

Example 5.2

Example 5.2 (Singular Solution)

Consider the Poisson equation

{-Δ⁢u=fin ⁢Ω,u=gon ⁢Γ,

where f and g are given to have the exact solution u⁢(x,y)=x3/2+y3/2. In this example the uniform and geometric grids on the unit square are considered.

Figures 5 and 6 represent numerical results for Example 5.2. By the low regularity of the exact solution the convergence order is limited, hence, similar convergence rates are observed from various abscissa subgrid methods for the uniform grid with respect to the h-refinement. In contrast to Example 5.1, the shape of the convergence history shows linear behavior with respect to the p-refinement. On the other hand the optimal convergence rates are recovered by introducing the geometric grid with a suitably ordered mesh transform as shown in Figure 6. Also exponential convergence behavior is observed for the p-refinement with the geometric grid.

Example 5.3

Example 5.3 (Curved Domain)

We consider the same Poisson equation as in Example 5.1 on a quarter disk.

In order to use the boundary information exactly for the curved domain, the following modification is proposed. If at least one vertex of a cell belongs to the curved boundary, abscissas’ locations are moved to the boundary so that the new abscissas lie on the axiparallel line (see Figure 7).

Figure 8 shows sample grids of a quarter disk, and Figure 9 represents the convergence history for the HSD on the quarter disk. The grid in Figure 8 is obtained by equally dividing the arc length. Therefore, cells tend to be more concentrated toward where the slope of the arc is 0 or ∞. It seems that it is not easy for the HSD to generate a quasi-uniform and shape-regular cell subdivision for domains like a quarter disk.

As shown in Figure 9 the convergence results are almost the same as those with the uniform grid on the unit square (Figure 4). Numerical experiments suggest that the HSD can be as well an effective numerical method for PDEs on the domain with a curved boundary.

Figure 7

Change of abscissas’ location in a boundary element.

Figure 8

Example grids for a quarter disk.

Figure 9

Numerical results with respect to the h-refinement (left) and p-refinement (right) for Example 5.3 on the quarter disk.

Table 1

Flux conservation property for Example 5.1:|∫∂⁡Dh∂νh⁡uh⁢d⁢s-2|,D=(0,1)×(0,12).

Remark 5.4

The HSD is designed to satisfy the mass conservation property. For a subdomain D⊂Ω consisting of a finite number of rectangular cells, the HSD satisfies the discrete mass conservation:

∫Dhf⁢𝑑x+∫∂⁡Dh∂νh⁡uh⁢d⁢s=0,

where ∫Dh and ∫∂⁡Dh denote the composite Gaussian quadratures on given abscissas.

Let Πh⁢f be the L2-orthogonal projection onto the standard binomial function space Qk-2⁢(R) for the _k_-abscissa case. If one substitutes f by Πh⁢f in (1.1), the exact flux conservation property holds:

Resflux=∫Df⁢𝑑x+∫∂⁡Dh∂νh⁡uh⁢d⁢s=0,

using ∫D(f-Πh⁢f)⁢𝑑x=0.

Table 2

Flux conservation property for Example 5.2 without the projection on the uniform and thegeometric grids: |∫∂⁡Dh∂νh⁡uh⁢d⁢s-34⁢(2⁢0.5+1)|,D=(0,1)×(0,12).

Table 3

Flux conservation property for Example 5.2 with the projection on the uniform and thegeometric grids: |∫∂⁡Dh∂νh⁡uh⁢d⁢s-34⁢(2⁢0.5+1)|,D=(0,1)×(0,12).

Table 4

Flux conservation property for Example 5.3:|∫∂⁡Dh∂νh⁡uh⁢d⁢s-1|,D=(0,12)×(0,12).

Tables 1–4 display the value |Resflux|. All the results confirm the flux conservation property except Table 2. Since our algorithm is based on the finite difference and collocation, we have

∫Dℐh⁢f⁢𝑑x+∫∂⁡Dh∂νh⁡uh⁢d⁢s=0

for the interpolation ℐh⁢f onto Qk-2⁢(R) for each R. The large errors in Table 2 reflect the interpolation error, which is |Resflux|=|∫D(f-ℐh⁢f)⁢𝑑x|. By replacing f with the L2 projection Πh⁢f, the desired conservation holds as in Table 3.

6 Conclusion

A hybrid spectral difference method is developed in this paper. The method is a finite difference version of the hybrid discontinuous Galerkin method [12, 13], resulting in an efficient conservative scheme. Numerical analysis and numerical experiments in the case of the Poisson equation in two space dimensions are presented. The optimal order of convergence in the discrete energy norm is proved. The L2 convergence theory can be a subject of future research. We observe that the numerical results concord well with our theoretical findings. Moreover, the HSD can manage the complicated geometry problems quite well without introducing any boundary transform.

For numerical and theoretical results of the Stokes and Navier–Stokes equations we refer to [11, 16], where some numerical experiments are performed, and the inf-sup condition is proved.

References

[1] Balan A., May G. and Schoberl J., A stable high-order spectral difference method for hyperbolic conservation laws on triangular elements, J. Comput. Phys. 231 (2012), 2359–2375.10.1016/j.jcp.2011.11.041Search in Google Scholar

[2] Cockburn B., Di Pietro D. and Ern A., Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods, ESAIM Math. Model. Numer. Anal. 50 (2016), no. 3, 635–650.10.1051/m2an/2015051Search in Google Scholar

[3] Cockburn B., Gopalakrishnan J. and Lazarov R., Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365.10.1137/070706616Search in Google Scholar

[4] Di Pietro D. and Ern A., A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Methods Appl. Mech. Engrg. 283 (2015), 1–21.10.1016/j.cma.2014.09.009Search in Google Scholar

[5] Di Pietro D., Ern A. and Lemaire S., An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators, Comput. Methods Appl. Math. 14 (2014), no. 4, 461–472.10.1515/cmam-2014-0018Search in Google Scholar

[6] Ferziger J. H. and Peric M., Computational Methods for Fluid Dynamics, Springer, Berlin, 2002.10.1007/978-3-642-56026-2Search in Google Scholar

[7] Fornberg B., A Practical Guide to Pseudospectral Methods. Vol. 1, Cambridge University Press, Cambridge, 1998.Search in Google Scholar

[8] Hesthaven J. and Warburton T., Nodal Discontinuous Galerkin Methods. Algorithms, Analysis, and Applications, Springer, New York, 2008.10.1007/978-0-387-72067-8Search in Google Scholar

[9] Ito K. and Qiao Z., A high order compact MAC finite difference scheme for the Stokes equations: Augmented variable approach, J. Comput. Phys. 227 (2008), 8177–8190.10.1016/j.jcp.2008.05.021Search in Google Scholar

[10] Jameson A., A proof of the stability of the spectral difference method for all orders of accuracy, J. Sci. Comput. 45 (2010), 348–358.10.1007/s10915-009-9339-4Search in Google Scholar

[11] Jeon Y., Hybrid difference methods for PDEs, J. Sci. Comput. 64 (2015), 508–521.10.1007/s10915-014-9941-ySearch in Google Scholar

[12] Jeon Y. and Park E.-J., A hybrid discontinuous Galerkin method for elliptic problems, SIAM J. Numer. Anal. 48 (2010), no. 5, 1968–1983.10.1137/090755102Search in Google Scholar

[13] Jeon Y. and Park E.-J., New locally conservative finite element methods on a rectangular mesh, Numer. Math. 123 (2013), 97–119.10.1007/s00211-012-0477-5Search in Google Scholar

[14] Jeon Y., Park E.-J. and Sheen D., A hybridized finite element method for the Stokes problem, Comput. Math. Appl. 68 (2014), no. 12B, 2222–2232.10.1016/j.camwa.2014.08.005Search in Google Scholar

[15] Jeon Y. and Sheen D., A locking-free locally conservative hybridized scheme for elasticity problems, Japanese J. Indust. Appl. Math. 30 (2013), 585–603.10.1007/s13160-013-0117-1Search in Google Scholar

[16] Jeon Y. and Sheen D., Upwind hybrid spectral difference methods for the steady-state Navier–Stokes equations, submitted.10.1007/978-3-319-72456-0_28Search in Google Scholar

[17] Kahaner D., Moler C. and Nash S., Numericl Methods and Software, Prentice-Hall, Upper Saddle River, 1989.Search in Google Scholar

[18] Kirby R. M., Berzins M. and Hesthaven J. S., Spectral and High Order Methods for Partial Differential Equations, Lect. Notes Comput. Sci. Eng. 106, Springer, Berlin, 2015.10.1007/978-3-319-19800-2Search in Google Scholar

[19] Lele S., Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992), 16–42.10.1016/0021-9991(92)90324-RSearch in Google Scholar

[20] LeVeque R., Finite Volume Methods for Hyperbolic Problems, Cambridge Texts Appl. Math., Cambridge University Press, Cambridge, 2002.10.1017/CBO9780511791253Search in Google Scholar

[21] Liu Y., Vinokur M. and Wang Z., Spectral difference method for unstructured grids I: Basic formulation, J. Comput. Phys. 216 (2006), 780–801.10.1016/j.jcp.2006.01.024Search in Google Scholar

[22] Sun Y., Wang Z. and Liu Y., Spectral (finite) volume method for conservation laws on unstructured grids VI: Extension to viscous flow, J. Comput. Phys. 215 (2006), 41–58.10.1016/j.jcp.2005.10.019Search in Google Scholar

[23] Turkel E., Gordon D., Gordon R. and Tsynkov S., Compact 2D and 3D sixth order schemes for the Helmholtz equation with variable wave number, J. Comput. Phys. 232 (2013), no. 1, 272–287.10.1016/j.jcp.2012.08.016Search in Google Scholar

[24] Wang Y. and Zhang J., Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D poisson equation, J. Comput. Phys. 228 (2009), 137–146.10.1016/j.jcp.2008.09.002Search in Google Scholar

[25] Wang Z., Spectral (finite) volume method for conservation laws on unstructured grids: Basic formulation, J. Comput. Phys. 178 (2002), 210–251.10.1006/jcph.2002.7041Search in Google Scholar

[26] Wang Z., Fidkowski K., Abgrall R., Bassi F., Caraeni D., Cary A., Deconinck H., Hartmann R., Hillewaert K., Huynh H., Kroll N., May G., Persson P.-O., van Leer B. and Visbal M., High-order CFD methods: Current status and perspective, Int. J. Numer. Methods Fluids 72 (2013), no. 8, 811–845.10.1002/fld.3767Search in Google Scholar

[27] Xie L., Xu M., Zhang B. and Qiu Z., A new spectral difference method using hierarchical polynomial bases for hyperbolic conservation laws, J. Comput. Phys. 284 (2015), 434–461.10.1016/j.jcp.2014.11.008Search in Google Scholar

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.