2 SBP Difference Operators

Chia-Chi Chang Min-Hung Chen C. H. Teng

Department of Mathematics

National Cheng Kung University, Tainan 701, Taiwan November 7, 2007

Correspondence: Min-Hung Chen

Email: mhchen@math.ncku.edu.tw Mailing address: Department of Mathematics

National Cheng Kung University

No.1 University Road, Tainan City 701, Taiwan Phone: +886-6-2757575 Ext 65127

Fax number: +886-6-2743191

1This work is partially supported by National Science Council grant No. NSC 95-2120-M- 001-003 and grant No. NSC 95-2115-M-006-010.

tors for approximating first order differential operators. We present these difference operators are global fifth and seventh order accurate schemes composed of eighth and tenth order stencils accurate at the interior points and enclosed by fourth and sixth order accurate numerical boundary sten- cils. The difference operators are specifically constructed such that these operators satisfy summation-by-parts rules. Numerical schemes based on the new difference operators incorporate with simultaneous approxima- tion term boundary conditions for solving model scalar waves problems were constructed. For temporal discretization, to achieve full high-order convergence in time and space, we march numerical solutions in time by implementing the m-stage mth-order strong stability-preserving Runge Kutta algorithms. Several applications of the difference operators are illustrated and their expected convergent patterns observed.

Key Words: summation-by-parts, high-order finite difference, simul- taneous approximation term, strong stability-preserving Runge Kutta, wave equations

1

1 Introduction

In the pioneering work contributed by Kreiss and Oliger [9], it was shown that high-order schemes, due to small accumulations of numerical phase error, are fa- vorable for simulating wave problems requiring long-time integrations. Serving as an important piece of evidence, this conclusion indeed has inspired researchers in the scientific computing community to construct high-order accurate simu- lation algorithms for solving time dependent problems. However, designing a high-order scheme remains a challenging task. A theoretically important issue is the stability of a high-order scheme.

For linear periodic problems, constructing a high-order accurate spatial ap- proximation based on central difference methods is a straight forward task, and stability of schemes can be examined by the classical von-Nuemann analysis or by energy methods (see [21]). However, to devise a scheme for a non-periodic problem, one needs to construct one-sided boundary stencils to enclose the do- main of interest. Without exercising great care, such an approach may cause a scheme to become unstable, even though the corresponding periodic scheme is stable in the von-Nuemann or in the energy sense.

For non-periodic problems, based on the concept of generalized summation- by-parts (SBP) property for differential equations introduced by Kreiss and Scherer [23], Strand [17] constructed high-order accurate SBP difference opera- tors. Such a methodology elucidated a way of constructing a high-order finite difference scheme which is stable in a discrete energy sense. However, it was also observed that numerical solutions computed by energy-stable SBP algorithms may still grow exponentially in time. Thus, this unwanted exponential growth may restrict the practical usage of such SBP computational methods in solving wave problems requiring long-time integrations.

Funaro and Gottlieb [5, 6] imposed a stable numerical boundary penalty term for solving hyperbolic problems by pseudospectral methods. With a set of free parameters introduced in a scheme, such a penalty formulation was in fact very constructive in the sense that the scheme can be made stable by properly choos- ing the values of the penalty parameters via a discrete energy estimate. More- over, it was shown that numerical solutions computed by these model spectral penalty schemes were time-stable. In addition to these conceptual schemes for model problems, advanced multidimensional spectral/pseudospectral penalty methods have also been developed for solving complex physical systems. For further details in the spectral penalty methods, we refer readers to the work by Hesthaven [8] and the references within.

There is an important similarity between the spectral penalty computa- tional framework and the SBP finite difference method. Notice that conducting an energy-based stability analysis for a pseudospectral penalty scheme hinges upon an important tool: Gauss’ quadrature integration rules. These quadra- ture integration rules, in fact, play a similar role to the SBP rules in high-order finite difference methods. Because of this similarity, one can construct an accu- rate time-stable scheme by adopting the penalized boundary formulation into a high-order SBP finite difference computational framework. Indeed, based on the compact finite difference methods, Carpenter et. al. [3] introduced such a boundary penalty formulation, called the simultaneous approximation term (SAT) procedure, and constructed SBP schemes for solving model scalar wave equations and systems of wave equations. It was shown that with the use of

the SAT boundary procedure, numerical solutions were stable in both an energy sense and in time.

Based on the SAT procedure and the theory established by Olsson [15, 16], Strand [18] further constructed SBP-SAT finite difference schemes and demon- strated time-stable wave computations. As demonstrated by Abarbanel, Cher- tock and Yefet [1, 2, 20] time-stable high-order compact difference schemes were successfully constructed for solving multidimensional scalar wave equations and systems of wave equations. Hence, we observe that the general framework pre- sented by Carpenter, Gottlieb and Abarbanel [3] is a road map toward to suc- cessfully and systematically constructing a stable and accurate scheme for wave problems. In this paper we follow a similar approach and present two high-order SBP schemes for solving wave problems. The features of the proposed schemes are as follows.

Spatial approximations are achieved by constructing one dimensional im- plicit difference operators which are of global fifth and seventh order accuracy.

The fifth order one is composed of an eighth order stencil accurate at the interior points and enclosed by a fourth order stencil accurate near the domain bound- aries. The seventh order one consists of tenth order accurate approximations at the interior points and sixth order accurate boundary stencils near the domain boundaries. Basic one dimensional difference operators are further employed to construct multidimensional difference operators used to compute numerical derivatives. In addition, SBP rules for the proposed difference operators, in- cluding the multidimensional ones, are established. Numerical SAT schemes for solving model wave problems in space of one dimension and in space of mul- tiple dimensions are constructed. Stability analysis of the schemes are then performed at the semi-discrete level. It is shown that these algorithms can be made asymptotically stable under certain discrete weighted norms by provid- ing suitable values of the penalty parameters. For temporal discretization we adopt the m-stage m-order strong stability-preserving (SSP) Runge Kutta (RK) algorithms [7] to march numerical solutions in time. Numerical validations of the proposed methods are tested through various model problems and we ob- served expected accuracy. The computational results indicate that the schemes are long-time stable. Thus, the proposed methods are favorable for simulating wave problems requiring long-time integration.

We organize the rest of the paper as follows. In Section 2 we present gen- eral concepts related to SBP implicit difference operators for approximating derivatives of functions defined on one dimensional and multidimensional spaces.

Section 3 is devoted to constructing the SBP schemes for solving model wave problems. Numerical validations of the schemes are presented in Section 4.

Concluding remarks are given last.

2 SBP Difference Operators

2.1 Basic Concepts and Notations

Let I = [0, 1]. Consider two functions f (x) and g(x) defined on I. We define the continuous L2 inner product and norm for functions over I as

(f, g) = Z

I

f g dx, ||f ||²I = (f, f )

The above concept can be extended for functions defined on multidimen- sional space. Consider I²= [0, 1] × [0, 1]. Let f and g be functions defined on I². The continuous L2 inner product and norm for functions over I² are defined as

(f, g) = Z

I²

f g dx dy, ||f ||²I²= (f, f ).

Likewise, for functions defined on I³ = [0, 1] × [0, 1] × [0, 1] we denote the con- tinuous L2 inner product and norm for functions over I³as

(f, g) = Z

I³

f g dx dy dz, ||f ||²I³= (f, f ).

Let L be a positive integer. Denote VLand MLthe vector spaces for vectors of length L and square matrices of size L, respectively.

We introduce a set of uniformly spaced grid points:

xi= ih, i = 0, 1, 2, ..., L, h = 1/L,

where h is the grid distance. Consider two vectors, f , g ∈ VL+1, explicitly given by

f = [f0, f1, ..., fL]^T, g= [g0, g1, ..., gL]^T,

where the superscript T denotes the transpose. Let M ∈ ML+1be a symmetric positive definite matrix. We define a weighted discrete L2 inner product and norm, with respect to the step size h and the matrix M , for vectors as

(f , g)h,M = hf^TM g, ||f ||²_h,M = (f , f )h,M.

If M is an identity matrix then we recover the un-weighted discrete L2 inner product and norm with respect to the step size h for vectors as

(f , g)h= hf^Tg, ||f ||²_h= (f , f )h.

To numerically approximate a function u and its derivative du/dx, we con- sider the difference approximation of the form

P vx= h⁻¹Qv, or vx= Dv = h⁻¹P⁻¹Qv, (1) where

v= [v0, v1, ..., vL]^T, vx= [v_x0, v_x1, ..., v_xL]^T,

denote the numerical approximations of u and u⁰ evaluated at the grid points, and D, P , Q ∈ ML+1.

Let u and ux denote vectors with components being, respectively, the col- located values of the functions u and du/dx at the grid points, i.e.,

u= [u(x0), u(x1), ..., u(xL)]^T, ux= du(x₀)

dx , du(x1)

dx , ..., du(xL) dx

T

.

The truncation error teof the scheme Eq.(1) is defined by P te= P ux− h⁻¹Qu,

and |te| = O(h^α_x, h^β_x) where α and β are the convergence rates of the approxi- mation at interior and boundary grid points, respectively.

In this paper, we devise implicit difference methods for approximating the differential operator d/dx by constructing a special class of P and Q satisfying the following properties;

SBP1: The matrix P is symmetric positive definite.

SBP2: The matrix Q is nearly skew-symmetric and satisfies the constraint Q+ Q^T

2 = diag(q00, 0, ..., 0, qLL), q00< 0, qLL= −q00> 0.

where q00 and qLL are the upper most and lower most diagonal elements of Q.

For an operator D = h⁻¹P⁻¹Q with P and Q satisfying the constraints SBP1 and SBP2, we have a result that mimics the integration-by-parts rule.

Lemma 1 (Summation-by-Parts). Consider the difference operator D = h⁻¹P⁻¹Qwhere P and Q satisfy SBP1 and SBP2, respectively. We have

(v, Dv)h,P = (v, h⁻¹P⁻¹Qv)h,P = q00v²₀+ qLLv_L², for v ∈ VL+1.

Proof. First we rewrite the inner product as (v, Dv)_h,P = v, h⁻¹Qv

h= v^TQ^Sv+ v^TQ^Av

where Q^S = (Q+Q^T)/2 and Q^A= (Q−Q^T)/2 are, respectively, the symmetric and anti-symmetric parts of the matrix Q. Notice that v^TQ^Av = 0 since Q^A is antisymmetric. Thus, we have

(v, Dv)_h,P = v^TQ^Sv = q00v²₀+ qLLv²_L, where the last equality is due to SPB2. This completes the proof.

Lemma 1 is of importance. For a continuous function f , with f⁰ denoting its derivative, we have

(f, f⁰) = Z

I

f f⁰dx = 1

2(f²(1) − f²(0)), (2) followed by the integration-by-parts rule, an important tool for showing well- posedness of initial boundary value problems by energy methods (see [21]). In fact, the result presented in Lemma 1 is a discrete analogue of Eq.(2), which will be used to establish stability results for numerical partial differential equa- tions. Later in this context, we will establish extended versions of Lemma 1 for multidimensional grid functions.

In this paper we adopt the methods shown by Carpenter et. al. and Chertock [3, 20] and construct two new high order SBP implicit difference operators: one is composed of an eighth order stencil accurate at interior points and a fourth order accurate stencil near the domain boundaries which yield a global fifth order approximation scheme; the other one is consisting of a tenth order stencil accurate at interior points and a sixth order stencil accurate near the boundaries which yields a global seventh order approximation scheme. Elements of matrices P and Q for these two operators are given in Appendix A.

2.2 Difference Approximations beyond 1D

We now generalize the introduced concepts into a multidimensional framework.

To clearly illustrate multidimensional SBP concepts it is important to have a clear set of notations. For such a purpose, we first present a set of notations by adopting those used in the studies [10, 11, 12, 13, 14] with some modifications.

We proceed as follows.

Let us introduce the Kronecker Product.

Definition 1 (Kronecker Product). Let A be a p×q matrix with its elements denoted by aij, and let B be an m × n matrix, then

A⊗ B =







a11B . . . a1qB ... . .. ... ap1B . . . apqB





.

The p × q block matrix A ⊗ B is called Kronecker Product.

Properties of Kronecker Product used through out this paper are summarized as follows:

Lemma 2.

(A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD), (A ⊗ B)^T = A^T ⊗ B^T

assuming that both the ordinary matrix multiplications AC and BD are defined.

Lemma 3. The Associative Property:

(A ⊗ B) ⊗ C = A ⊗ (B ⊗ C) Lemma 4. The Bilinear Property:

(A + B) ⊗ C = A ⊗ C + B ⊗ C A⊗ (B + C) = A ⊗ B + A ⊗ C (kA) ⊗ B = A ⊗ (kB) = k(A ⊗ B) Lemma 5. There exists a permutation matrix R such that

A⊗ B = R^T(B ⊗ A)R

Lemma 6. If {λi} and {µj} are eigenvalues of matrix A and B, then A ⊗ B has eigenvalue as the form λiµj for all i, j.

Lemma 7. Let A ∈ Mn and B ∈ Mm. If {λi} and {µj} are eigenvalues of matrix A and B, then I ⊗ A + B ⊗ I has eigenvalue as the form λi+ µj for all i, j.

Readers who are interested in the details of the Kronecker product are re- ferred to [22].

Let L, M and N denote positive integers. Through out the rest of this paper, we use i, j, k and i⁰, j⁰, k⁰ as dummy indices, and

0 ≤ i, i⁰≤ L, 0 ≤ j, j⁰≤ M, 0 ≤ k, k⁰≤ N

unless specified otherwise. Moreover, for simplicity we use the following short- handed notations:

L

X

i,i⁰=0

=

L

X

i=0 L

X

i⁰=0

,

M

X

j,j⁰=0

=

M

X

j=0 M

X

j⁰=0

,

N

X

k,k⁰=0

=

N

X

k=0 N

X

k⁰=0

.

Denote δµν the Kronecker delta function. Define exi= [δi0, δi1, ..., δiL]^T, exi∈ VL+1, eyj = [δj0, δj1, ..., δjM]^T, eyj ∈ VM +1, ezk = [δk0, δk1, ..., δkN]^T, ezk∈ VN +1.

Note that the sets {exi}^L_i=0, {eyj}^M_j=0 and {ezk}^N_k=0consist bases of the vector spaces VL+1, VM +1 and VN +1, respectively. In addition, we define the sets:

{eij= eyj⊗ exi| 0 ≤ i ≤ L, 0 ≤ j ≤ M },

{eijk= ezk⊗ eyj⊗ exi| 0 ≤ j ≤ M, 0 ≤ k ≤ N },

which consist bases of the vector spaces V(M +1)×(L+1)and V(N +1)×(M +1)×(L+1), respectively.

For two dimensional space problems we consider a 2D grid mesh with grid points (xi, yj) ∈ I² defined by

(xi, yj) = (ihx, jhy), hx= 1

L, hy= 1

M. (3)

We define

v= [v00v10. . . vL−1MvLM]^T =

L

X

i=0 M

X

j=0

vijeij, (4) in which vij denote the grid-function values at the grid points.

Notice that Pxand Py are symmetric positive definite. Thus, by Lemma 2 and Lemma 6, the matrix Pxy = Py⊗ Px∈ M(M +1)×(L+1) is also symmetric positive definite. We define a weighted discrete inner product and norm for vectors f , g ∈ V(M +1)×(L+1)as

(f , g)hxy,Pxy = hxhyf^T(Py⊗ Px)g, ||f ||²_hxy,Pxy = (f , f )hxy,Pxy. Similarly, we define an un-weighted discrete inner product and norm for vectors f, g ∈ V(M +1)×(L+1) as

(f , g)hxy = hxhyf^Tg, ||f ||²_hxy = (f , f )hxy.

For three dimensional space problems, we consider a set of three dimensional grid points, (xi, yj, zk), where

(xi, yj, zk) = (ihx, jhy, khz), hx= 1

L, hy= 1

M, hz= 1

N. (5)

We define the three dimensional grid-function v by

v= [v000v100 . . . vL−1M NvLM N]^T =

L

X

i=0 M

X

j=0 N

X

k=0

vijkeijk, (6)

where vijk denote the grid-function values at the grid points (xi, yj, zk).

Notice that Px, Pyand Pz are symmetric positive definite matrices. Thus, by Lemma 2 and Lemma 6, the matrix Pxyz= Pz⊗ Py⊗ Pxis also symmetric positive definite. We define a weighted discrete inner product and norm for vectors f , g ∈ V(N +1)×(M +1)×(L+1) as

(f , g)hxyz,Pxyz= hxhyhzf^T(Pz⊗ Py⊗ Px)g,

||f ||²_hxyz,Pxyz = (f , f )hxyz,Pxyz.

Similarly, an un-weighted discrete inner product and norm for vectors f , g ∈ V(N +1)×(M +1)×(L+1) are defined as

(f , g)hxyz= hxhyhzf^Tg, ||f ||²_hxyz = (f , f )hxyz. 2.2.1 Difference Operators in 2D

We now seek difference operators for calculating numerical derivatives of a vector v∈ V(M +1)×(L+1). Denote the one dimensional difference operators in x and y directions, respective, by

Dx= h⁻¹_x P⁻¹_x Q_x, Dx, Px, Q_x∈ ML+1

Dy = h⁻¹_y P⁻¹_y Q_y, Dy, Py, Q_y ∈ MM +1,

where Pxand Pysatisfy SBP1 and Q_xand Q_y satisfy SBP2. To numerically differentiate v, we define the two dimensional difference operators for approxi- mating ∂/∂x and ∂/∂y, respectively, as

D_x= IM +1⊗ Dx= IM +1⊗ h⁻¹_x P⁻¹_x Q_x, (7a) D_y= Dy⊗ IL+1= h⁻¹_y P⁻¹_y Q_y⊗ IL+1, (7b) where IL+1∈ ML+1 and IM +1∈ MM +1 are identity matrices.

Equations (7) lead to the following result.

Lemma 8 (Summation-by-Parts in 2D). Consider the 2D difference oper- ators, Eqs.(7). Let q_ii^x0 and q^y_jj0 denote the elements of the matrices, Q_x and Q_y, respectively. For v ∈ V(L+1)×(M +1), we have

(v, Dxv)_hxy,Pxy = q₀₀^x||v0j||²_hy,Py+ q^x_LL||vLj||²_hy,Py, (v, Dyv)_h

xy,Pxy = q₀₀^y ||vi0||²_hx,Px+ q^y_{M M}||viM||²_hx,Px, where

v0j=

M

X

j=0

v0jeyj, vLj=

M

X

j=0

vLjeyj,

vi0=

L

X

i=0

vi0exi, viM =

L

X

i=0

vi0exi.

Proof. Denote Q^S_x = (Q_x+ Q^T_x)/2 and Q^A_x = (Q_x− Q^T_x)/2 the symmetric and anti-symmetric parts of Q_x, respectively. We have

(v, Dxv)_hxy,Pxy = hyv^T(Py⊗ Q^S_x)v + hyv^T(Py⊗ Q^A_x)v.

Since Py⊗ Q^A_x is anti-symmetric, the second term vanishes. We now treat the first term

hyv^T(Py⊗ Q^S_x)v

=

L

X

i,i⁰=0 M

X

j,j⁰=0

hyvijvi⁰j⁰(e^T_yj0⊗ e^T_xi0)(Py⊗ Q^S_x)(eyj⊗ exi)

in details. Applying Lemma 2 we have

(e^T_yj0⊗ e^T_xi0)(Py⊗ Q^S_x)(eyj⊗ exi)

=e^T_yj0Pyeyj⊗ e^T_xi0Q^S_xexi

=δii⁰q_ii^x0(e^T_yj0Pyeyj) where the last equality results from

Q^S_x = diag(q^x₀₀0, . . . , 0, q_LL^x ) due to SBP2. Therefore, we obtain

hyv^T(Py⊗ Q^S_x)v

=hyq₀₀^x

M

X

j⁰,j=0

v0j⁰e^T_yj0Pyeyjv0j+ hyq^x_LL

M

X

j⁰,j=0

vLj⁰e^T_yj0PyeyjvLj,

and prove the first statement.

The second part of the proof is similar and we omit the details. This com- pletes the proof.

2.2.2 Difference Operators in 3D

To conduct numerical differentiation of v defined by Eq.(6), we take a similar approach as shown in Sec.2.2.1. Define the one dimensional difference operators along x, y and z directions, respective, by

Dx= h⁻¹_x P⁻¹_x Q_x, Dx, Px, Q_x∈ ML+1, Dy = h⁻¹_y P⁻¹_y Q_y, Dy, Py, Q_y ∈ MM +1, Dz= h⁻¹_z P⁻¹_z Q_z, Dz, Pz, Q_z∈ MN +1,

where Px, Py and Pz satisfy SBP1 and Q_x, Q_y and Q_z satisfy SBP2. The three dimensional difference operators for approximating ∂/∂x, ∂/∂y and ∂/∂z are defined as

D_x= IN +1⊗ IM +1⊗ Dx= IN +1⊗ IM +1⊗ h⁻¹_x P⁻¹_x Q_x, (8a) D_y= IN +1⊗ Dy⊗ IL+1= IN +1⊗ h⁻¹_y P⁻¹_y Q_y⊗ IL+1, (8b) D_z= Dz⊗ IM +1⊗ IL+1= h⁻¹_z P⁻¹_z Q_z⊗ IM +1⊗ IL+1, (8c) where IL+1∈ ML+1, IM +1∈ MM +1 and IN +1∈ MN +1 are identity matrices.

Equations (8) lead to the following result.

Lemma 9 (Summation-by-Parts in 3D). Consider the 3D difference oper- ators, Eqs.(8). Let q_ij^x, q_ij^y and q^z_ij denote the elements of the matrices Q_x, Q_y and Q_z, respectively. We have

(v, Dxv)hxyz,Pxyz = q₀₀^x||v0jk||²_hyz,Pyz+ q_LL^x ||vLjk||²_hyz,Pyz (v, Dyv)hxyz,Pxyz = q₀₀^y ||vi0k||²_hxz,Pxz+ q_{M M}^y ||viM k||²_hxz,Pxz (v, Dzv)h_xyz,P_xyz = q₀₀^z ||vij0||²_hxy,Pxy+ q^z_{N N}||vijN||²_hxy,Pxy where

v0jk=

M

X

j=0 N

X

k=0

v0jkezk⊗ eyj, vLjk=

M

X

j=0 N

X

k=0

vLjkezk⊗ eyj

vi0k=

L

X

i=0 N

X

k=0

vi0kezk⊗ exi, viM k =

L

X

i=0 N

X

k=0

viM kezk⊗ exi,

vij0=

L

X

i=0 M

X

j=0

vij0eyj⊗ exi, vijN =

L

X

i=0 M

X

j=0

vijNeyj⊗ exi.

Proof. Denote Q^S_x = (Q_x+ Q^T_x)/2 and Q^A_x = (Q_x− Q^T_x)/2 the symmetric and anti-symmetric parts of Q_x, respectively. We have

(v, Dxv)hxyz,Pxyz

=hyhzv^T(Pz⊗ Py⊗ Q^S_x)v + hyhzv^T(Pz⊗ Py⊗ Q^A_x)v.

Notice that the second term vanishes since the matrix Pz⊗ Py⊗ Q^A_x is anti- symmetric. We now treat the first term

hyhzv^T(Pz⊗ Py⊗ Q^S_x)v

=hyhz L

X

i,i⁰=0 M

X

j,j⁰=0 N

X

k,k⁰=0

vijkvi⁰j⁰k⁰e^T_i0j⁰k⁰(Pz⊗ Py⊗ Q^S_x)e^T_ijk

in details. Notice that

e^T_i0j⁰k⁰(Pz⊗ Py⊗ Q^S_x)eijk

=(e^T_zk0⊗ e^T_yj0⊗ e^T_xi0)(Pz⊗ Py⊗ Q^S_x)(ezk⊗ eyj⊗ exi)

=(e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj) ⊗ e^T_xi0Q^S_xexi

=δii⁰q_ii^x0(e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj).

So

(v, Dxv)h_xyz,P_xyz

=q₀₀^x

M

X

j,j⁰=0 N

X

k,k⁰=0

v0jkv0j⁰k⁰hyhz(e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj)

+q_LL^x

M

X

j,j⁰=0 N

X

k,k⁰=0

vLjkvLj⁰k⁰hyhz(e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj)

=q₀₀^x||v0jk||²_hyz,Pyz+ q_LL^x ||vLjk||²_hyz,Pyz,

and we prove the first part of the lemma.

Since the proofs of the other two statements are similar, we omit the deriva- tion details. This completes the proof.

3 Schemes for Model Wave Equations

In this section, we construct schemes for model wave problems using the dif- ference methods introduced in the previous section. We start by devising a scheme for solving one dimensional scalar wave problems. Stability analysis will be discussed at the semi-discrete level. We then illustrate the multidimensional setting, followed by the one dimensional one.

3.1 One Dimensional Advection Equation

Consider the advection equation

∂u

∂t +∂u

∂x = 0, x ∈ I, t ≥ 0, (10)

with the initial condition

u(x, 0) = f (x), x ∈ I, (11)

and the boundary condition

u(0, t) = g(t), t ≥ 0. (12)

Equation (10) leads to an energy rate d

dt||u||²I = g²(t) − u²(1, t),

followed by multiplying u to Eq.(10), integrating the resultant equation over the domain I, conducting integration-by-parts and applying the boundary condition.

For well-posed analysis it is sufficient to consider g = 0, and we obtain an energy estimate

||u(x, t)||²I ≤ ||u(x, 0)||²I = ||f (x)||²I. (13) We now devise a scheme, which mimics Eq.(13), for numerically solving the one dimensional scalar wave problem.

Consider a equally spaced partition:

xi= ih, h = 1/L.

With videnoting the approximation of u(xi), we seek a numerical solution v of the form

v(t) = [v0(t), v1(t), ..., vL(t)]^T, which satisfies the semidiscrete scheme

dv

dt + h⁻¹P⁻¹Qv= h⁻¹τ q00(v0− g(t))P⁻¹e0, (14a) v(0) = f = [f (x0), f (x1), ..., f (xL)]^T, (14b)

where P and Q are defined by Eq.(1), and e0= [1, 0, ..., 0]^T.

In Eq.(14a) it is shown that the boundary condition is incorporated into the scheme through the SAT method. The introduced variable τ is a free parameter, and it’s value is determined by conducting a discrete energy estimate for which Eq.(14a) is stable at the semi-discrete level. The stability analysis is summarized as a theorem.

Theorem 1. Assume that there exists a smooth solution to the one dimensional wave problem described by Eqs.(10-12). Then Eq.(14a) is stable at the semi- discrete level provided that

τ ≥ 1.

Moreover, v(t) satisfies the estimate

||v(t)||²_h,P ≤ ||f ||²_h,P.

Proof. Similar to the well-posed analysis yielding an energy estimate for the continuous problem, we now conduct a stability analysis for the scheme based on a discrete energy estimate in || · ||h,P norm.

Multiplying hv^TP to the scheme and invoking Lemma 1 we obtain 1

2 d

dt||v||²_h,P = −(q00v₀²+ qLLv²_L) + τ v0q00(v0− g(t)).

For the stability analysis, it is sufficient enough to consider the scheme subject to g(t) = 0. Hence,

1 2

d

dt||v||²_h,P = q00(τ − 1)v²₀− qLLv_L² ≤ q00(τ − 1)v₀²,

where the last inequality results from qLL> 0 demanded by SBP2. Recall that q00< 0. So, taking τ ≥ 1 immediately yields a non-increasing energy rate

1 2

d

dt||v||²_h,P ≤ 0, which leads to the estimate

||v(t)||²_h,P ≤ ||v(0)||²_h,P = ||f ||²_h,P. Thus, the scheme is stable. This completes the proof.

3.2 Multidimensional Wave Equations

3.2.1 Scheme for 2D Wave Problem Consider the 2D advection equation

∂u

∂t +∂u

∂x+∂u

∂y = 0, (x, y) ∈ I², t ≥ 0, (15) with the initial condition at t = 0

u(x, y, 0) = f (x, y), (x, y) ∈ I², (16)

and the boundary condition

u(0, y, t) = g1(y, t), 0 ≤ y ≤ 1, t ≥ 0, (17a) u(x, 0, t) = g2(x, t), 0 ≤ x ≤ 1, t ≥ 0. (17b)

We first give an energy estimate of the two dimensional problem. Multiplying u to Eq.(15), integrating the resultant equation over the domain I², conducting integration-by-parts and invoking Eqs.(17) with g1 = g2 = 0, we obtain an energy rate

d

dt||u||²I² ≤ 0, (18)

which implies

||u||²I² ≤ ||f ||²I², (19) followed by integrating Eq.(18) and further applying Eq.(16).

Similar to the 1D case, we now find an approximation that mimics Eq.(19).

Consider the 2D grid mesh with grid points (xi, yj) defined by Eq.(3). Let vij denote the grid-function values at the grid points (xi, yj). We define the numerical solution v ∈ V(L+1)×(M +1)of the form

v(t) =

L

X

i=0 M

X

j=0

vij(t)eij

satisfying the scheme dv

dt + Dxv+ Dyv= r, (20a)

v(t = 0) = f =

L

X

i=0 M

X

j=0

f (xi, yj)eij. (20b)

where

r= h⁻¹_x τ1 M

X

j=0

q₀₀^x(v0j− g1(yj, t))(eyj⊗ P⁻¹_x ex0)

+ h⁻¹_y τ2 L

X

i=0

q^y₀₀(vi0− g2(xi, t))(P⁻¹_y ey0⊗ exi).

As in the one dimensional case, the boundary condition is weakly imposed upon the scheme through the SAT methodology. The introduced variable τ1

and τ2 are free parameters, and their values are determined by conducting a discrete energy estimate for which Eq.(20) is stable at the semi-discrete level in

|| · ||Py⊗Px norm. The stability analysis is summarized as a theorem.

Theorem 2. Assume that there exists a smooth solution to the two dimensional wave problem described by Eqs.(15-17). Then Eq.(20a) is stable at the semi- discrete level provided that

τ1, τ2≥ 1.

Moreover, v(t) satisfies the estimate

||v(t)||²_hxy,Pxy ≤ ||f ||²_hxy,Pxy. (21)

Proof. For stability analysis, it is sufficient to consider the problem subject to homogeneous boundary conditions. We have

1 2

d

dt||v||²_hxy,Pxy =

 v,dv

dt



hxy,Pxy

= − (v, Dxv)hxy,Pxy− (v, Dyv)hxy,Pxy+ (v, r)hxy,Pxy. (22) Notice that with ei⁰j⁰ = eyj⁰⊗ exi⁰

(v, r)h_xy,P_xy = hyτ1q₀₀^x

L

X

i⁰=0 M

X

j⁰,j=0

vi⁰j⁰v0j(e^T_i0j⁰)(Py⊗ Px)(eyj⊗ P⁻¹_x ex0)

+ hxτ2q₀₀^y

L

X

i⁰,i=0 M

X

j⁰=0

vi⁰j⁰vi0(e^T_i0j⁰)(Py⊗ Px)(P⁻¹_y ey0⊗ exi)

Applying Lemma 2 one can easily show that

(e^T_i0j⁰)(Py⊗ Px)(eyj⊗ P⁻¹_x ex0) = δi⁰0(e^T_yj0Pyeyj) (e^T_i0j⁰)(Py⊗ Px)(P⁻¹_y ey0⊗ exi) = δj⁰0(e^T_xi0Pxexi) Therefore,

(v, r)hxy,Pxy = τ1q^x₀₀

M

X

j⁰,j=0

hyv0j⁰v0j(e^T_yj0Pyeyj)

+ τ2q^y₀₀

L

X

i⁰,i=0

hxvi⁰0vi0(e^T_xi0Pxexi)

= τ1q^x₀₀||v0j||²_hy,Py+ τ2q₀₀^y ||vi0||²_hx,Px

Substituting the above expression into Eq.(22) and invoking Lemma 8, we have 1

2 d

dt||v||²_hxy,Pxy = (τ1− 1)q₀₀^x||v0j||²_Py− q^x_LL||vLj||²_Py + (τ2− 1)q^y₀₀||vi0||²_Px− q^y_{M M}||viM||²_Px Taking τ1, τ2≥ 1 yields

1 2

d

dt||v||²_hxy,Pxy ≤ 0,

which implies Eq.(21). Thus, the scheme is stable. This completes the proof.

3.2.2 Scheme for 3D Wave Problem Consider the advection equation

∂u

∂t +∂u

∂x+∂u

∂y +∂u

∂z = 0, (x, y, z) ∈ I³, t ≥ 0, (23) with the initial condition

u(x, y, z, 0) = f (x, y, z), (x, y, z) ∈ I³, (24)

and the boundary condition

u(0, y, z, t) = g1(y, z, t), 0 ≤ y, z ≤ 1, t ≥ 0, (25a) u(x, 0, z, t) = g2(x, z, t), 0 ≤ x, z ≤ 1, t ≥ 0, (25b) u(x, y, 0, t) = g3(x, y, t), 0 ≤ x, y ≤ 1, t ≥ 0. (25c)

We first give an energy estimate of the three dimensional problem. Integrat- ing the equation, resulting from multiplying u to Eq.(23), over the domain I³, performing integration-by-parts and applying Eqs.(25) with g1 = g2= g3 = 0, we yield an energy rate:

d

dt||u||²I³ ≤ 0, leading to

||u||²I³ ≤ ||f ||²I³. (26) Similar to 1D and 2D cases, we now seek an approximation possessing a discrete energy estimate mimicking Eq.(26).

Consider the 3D grid mesh defined by Eq.(5). Let vijk denote the grid- function values at the grid points (xi, yj, zk). We define the numerical solution of the form

v(t) =

L

X

i=0 M

X

j=0 N

X

k=0

vijk(t)eijk

satisfying the scheme dv

dt + Dxv+ Dyv+ Dzv = r, (27a) v(0) = f =

L

X

i=0 M

X

j=0 N

X

k=0

f (xi, yj, zk)eijk. (27b)

where

r= h⁻¹_x τ1 M

X

j=0 N

X

k=0

q₀₀^x(v0jk− g1(yj, zk, t))(ezk⊗ eyj⊗ P⁻¹_x ex0)

+ h⁻¹_y τ2 L

X

i=0 N

X

k=0

q^y₀₀(vi0k− g2(xi, zk, t))(ezk⊗ P⁻¹_y ey0⊗ exi)

+ h⁻¹_z τ3 L

X

i=0 M

X

j=0

q^z₀₀(vij0− g3(xi, yj, t))(P⁻¹_z ez0⊗ eyj⊗ exi)

As in the reduced dimensional cases, the boundary condition is weakly im- posed upon the scheme through the SAT methodology. The introduced variables τ1τ2and τ3are free parameters, and their values are determined by conducting a discrete energy estimate for which Eq.(27a) is stable at the semi-discrete level in || · ||Px⊗Py⊗Pz norm. The stability analysis is summarized as a theorem as follows.

We first prove the result.

Lemma 10. Let g1= g2= g3= 0. Then (v, r)hxyz,Pxyz

=τ1q₀₀^x||v0jk||²_hyz,Pyz+ τ2q₀₀^y ||vi0k||²_hxz,Pxz+ τ3q₀₀^z ||vij0||²_hxy,Pxy Proof. Notice that

(v, r)hxyz,Pxyz = τ1q^x₀₀r1+ τ2q^y₀₀r2+ τ3q₀₀^z r3, where with Pxyz= Pz⊗ Py⊗ Px

r1=

L

X

i⁰=0 M

X

j⁰,j=0 N

X

k⁰,k=0

hyhzvi⁰j⁰k⁰v0jk(e^T_i0j⁰k⁰)(Pxyz)(ezk⊗ eyj⊗ P⁻¹_x ex0)

r2=

M

X

j⁰=0 L

X

i⁰,i=0 N

X

k⁰,k=0

hxhzvi⁰j⁰k⁰vi0k(e^T_i0j⁰k⁰)(Pxyz)(ezk⊗ P⁻¹_y ey0⊗ exi)

r3=

N

X

k⁰=0 L

X

i⁰,i=0 M

X

j⁰,j=0

hxhyvi⁰j⁰k⁰vij0(e^T_i0j⁰k⁰)(Pxyz)(P⁻¹_z ez0⊗ eyj⊗ exi)

Observe that

(e^T_i0j⁰k⁰)(Pxyz)(ezk⊗ eyj⊗ P⁻¹_x ex0)

=(e^T_zk0 ⊗ e^T_yj0⊗ e^T_xi0)(Pz⊗ Py⊗ Px)(ezk⊗ eyj⊗ P⁻¹_x ex0)

= (e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj)
⊗ (e^T_xi0ex0)

=δi⁰0(e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj) Hence

r1=

M

X

j⁰,j=0 N

X

k⁰,k=0

hyhzv0j⁰k⁰v0jk(e^T_zk0⊗ e^T_yj0)(Pz⊗ Py)(ezk⊗ eyj)

= ||v0jk||²_hyz,Pyz

Following similar approaches we obtain

r2= ||vi0k||²_hxz,Pxz, r3= ||vij0||²_hxy,Pxy. This completes the proof.

Theorem 3. Assume that there exists a smooth solution to the three dimen- sional wave problem described by Eqs.(23-25). Then Eq.(27a) is stable at the semi-discrete level provided that

τ1, τ2, τ3≥ 1.

Moreover, v(t) satisfies the estimate

||v(t)||²_hxyz,Pxyz ≤ ||f ||²_hxyz,Pxyz.

Proof. For stability analysis, it is sufficient to consider the problem subject to homogeneous boundary conditions, i.e., g1= g2= g3= 0

We have

1 2

d

dt||v||²_hxyz,Pxyz =

 v,dv

dt



hxyz,Pxyz

= − (v, Dxv)hxyz,Pxyz− (v, Dyv)hxyz,Pxyz

− (v, Dzv)hxyz,Pxyz+ (v, r)hxyz,Pxyz. Invoking Lemma 9 and Lemma 10, we have

1 2

d

dt||v||²_hxyz,Pxyz = (τ1− 1)q₀₀^x||v0jk||²_hyz,Pyz− q_LL^x ||vLjk||²_hyz,Pyz + (τ2− 1)q^y₀₀||vi0k||²_hxz,Pxz− q_{M M}^y ||viM k||²_hxz,Pxz + (τ3− 1)q^z₀₀||vij0||²_hxy,Pxy− q^z_{N N}||vijN||²_hxy,Pxy

Recall that q^x₀₀, q^y₀₀, q^z₀₀< 0, and q_LL^x , q^y_{M M}, q^z_{N N} > 0 due to SBP2. So, taking τ1, τ2, τ3≥ 1 yields

1 2

d

dt||v||²_hxyz,Pxyz≤ 0, which implies

||v(t)||²_hxyz,Pxyz ≤ ||f ||²_hxyz,Pxyz. Therefore, the scheme is stable. This completes the proof.

3.3 Time Discretization

Let us write the semi-discrete schemes as d

dtC = LC + S(t) (28)

where C(t) = (cj(t))_1≤j≤K, L consists the difference operator modified by the SAT term, and S(t) contains the time explicit boundary conditions. To approx- imate the solution to Eq (28), we implement a scheme which is an extension to the SSP-RK scheme introduced in Gottlieb et al. [7] for autonomous systems.

Denoting Cⁿ≡ C(tn), we seek an m-th order, m-stage scheme in the form







C⁽⁰⁾ = Cⁿ,

C⁽ⁱ⁾ = C⁽ⁱ⁻¹⁾+ htLC⁽ⁱ⁻¹⁾+ htS⁽ⁱ⁾, i = 1, ..., m, Cⁿ⁺¹ =Pm

k=0αm,kC^(k),

(29)

where the ht denotes the time step and the coefficients αm,k are those corre- sponding to the scheme in Ref.[7], namely







α1,0 = 1

αm,k =¹_kαm−1,k−1, k = 1, ..., m − 2, αm,m =_m!¹ , αm,m−1= 0, αm,0= 1 −Pm

k=1αm,k.

(30)

and the source terms S⁽ⁱ⁾are derived in Chen et. al. [4]

S⁽ⁱ⁾= (I + ht∂t)ⁱ⁻¹S(tn), (31) where I denotes the identity operator.