### High-Order Summation-by-Parts Implicit Difference Operators for Wave Problems

CHIA-CHI CHANG, MIN-HUNGCHEN, CHUN-HAO TENG Department of Mathematics

National Cheng Kung University Tainan City 701, Taiwan

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### Motivation

**1** Why high-order methods?

Ans: Preferred for long time simulation.

**2** Why high-order finite difference methods?

Ans: Efficiency. For the same order of accuracy, the time step is larger than the Pseudo-spectral methods.

**3** What is SBP and why SBP?

Ans: To construct schemes leading to a bounded energy estimate for bounded domains with non-periodic boundary conditions.

Example

Consider the modal wave problem

∂u

∂t + ∂u

∂x = 0, x∈ [0, 1], t≥ 0, u(x, 0) = f (x), x∈ [0, 1], t= 0, u(0, t) = g(t), x= 0, t≥ 0.

Energy Method 1

2 dE(t)

dt = 1 2

Z 1 0

∂u^{2}

∂t dx= −1 2

Z 1 0

∂u^{2}

∂x dx

= 1

2(u^{2}(0, t) − u^{2}(1, t)) = 1

2(g^{2}(t) − u^{2}(1, t))

### Central Difference Scheme

Strongly Enforced Boundary Condition

dv_{i}

dt +v_{i+1}− v_{i−1}

2h = 0 i= 1, 2, . . . , N − 1 dvN

dt +vN− vN−1

h = 0

v0(t) = g(t)
v_{i}(0) = f (xi)

Accuracy: interior point (2nd order); boundary point (1st order) Energy estimate: Multiply civihto the scheme where cN =1

2, ci= 1 1

2 dED

dt = 1 2

N

X

i=1

civi

dvi

dth= 1

2(v1v0− v^{2}_{N})
Recall the continuous system

1 2

dE(t) dt = 1

2(g^{2}(t) − u^{2}(1, t))

### Central Difference Scheme

Weakly Enforced Boundary Condition- Penalty Method

Semi-discrete scheme (penalty method):

dv_{i}

dt +v_{i+1}− v_{i−1}

2h = 0 i= 1, 2, . . . , N − 1 dvN

dt +vN− vN−1

h = 0

dv0

dt +v1− v_{0}

h = −τ (v_{0}− g(t))
Energy estimate:

dED

dt = 1 2v0

dv0

dth+

N−1

X

i=1

vi

dvi

dth+1 2vN

dvN

dt h

### Central Difference Scheme

Weakly Enforced Boundary Condition

d dt

N

X

i=0

civ^{2}_{i}h= v1v0− v^{2}_{N}− v1v0+ v^{2}_{0}− τ hv^{2}_{0}+ τ hv0g(t)

= −v^{2}_{N}+ (1 − τ h)v^{2}_{0}+ τ hv_{0}g(t)

= −v^{2}_{N}+ (1 − τ h)

"

v^{2}_{0}+ τ h

1 − τ hv0g(t) +

τ h 2(1 − τ h)

^{2}
g^{2}(t)

#

− (1 − τ h)

τ h 2(1 − τ h)

2

g^{2}(t)

If τ h = 2

dE_{D}

dt = −v^{2}_{N}− (v_{0}− g(t))^{2}+ g^{2}(t)

### Central Difference Scheme

Matrix Vector Representation Let v = [v0(t) v1(t) · · · vN(t)]

Scheme

d dt

2 6 6 6 6 6 6 6 6 6 6 4

v_{0}
v_{1}
v_{2}
..
.
v_{N−1}

v_{N}
3
7
7
7
7
7
7
7
7
7
7
5

+ 1 2h

2 6 6 6 6 6 6 6 6 6 6 4

−2 2 0

−1 0 1 . ..

0 −1 0 1 . ..

. .. . .. . .. . ..

−1 0 1

0 −2 2

3 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4

v_{0}
v_{1}
v_{2}
..
.
v_{N−1}

v_{N}
3
7
7
7
7
7
7
7
7
7
7
5

= 2 6 6 6 6 6 6 6 6 6 6 4

−τ (v_{0}− g(t))
0
0
..
.
0
0

3 7 7 7 7 7 7 7 7 7 7 5

or d

dtv(t) + Dv(t) = −τ (v_{0}− g(t))e_{0}, e_{0}= [1 0 0 · · · 0 0]^{T}
Define H = diag(1

2, 1, 1, · · · , 1, 1
2)
v^{T}Hdv

dth+ v^{T}H D h v(t) = −τ (v0− g(t))v^{T}H e0h

### Central Difference Scheme

Interesting Property: Summation-by-Parts Rule Observe

H D h = 1 2 2 6 6 6 6 6 6 6 6 4

1/2 0 0 1 . ..

. .. . .. . .. . .. 1 0

0 1/2 3 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 4

−2 2 0

−1 0 1

. .. . .. . ..

−1 0 1

−2 2 3 7 7 7 7 7 7 7 7 5

=1 2 2 6 6 6 6 6 4

−1 1 0

−1 0 1

. .. . .. . ..

−1 0 1

−1 1 3 7 7 7 7 7 5

=1 2 2 6 6 6 6 6 4

−1 0

. .. 1

3 7 7 7 7 7 5

+1 2 2 6 6 6 6 6 4

0 1 0

−1 0 0

. .. . .. . ..

0 0 1

−1 0 3 7 7 7 7 7 5

=Q_{S}+ Q_{A}

So

v^{T}HDhv = v^{T}Q_{S}v + v^{T}Q_{A}v = 1

2(−v^{2}_{0}+ v^{2}_{N})

### Central Difference Scheme

Interesting Property: Summation-by-Parts Rule

Hence

v^{T}Hdv

dth= −v^{T}H D h v − τ (v0− g(t))v^{T}H e0h
1

2 d

dtv^{T}Hvh = 1
2v^{2}_{0}−1

2v^{2}_{N}−τ h

2 v_{0}(v_{0}− g(t))
The rule

v^{T}HDhv = v^{T}Q_{S}v + v^{T}Q_{A}v = 1

2(−v^{2}_{0}+ v^{2}_{N})
in fact mimics the action of

Z 1 0

u(x)∂u(x)

∂x dx= 1

2u^{2}(1) − 1
2u^{2}(0)

### Remarks

To construct a high-order finite difference scheme we basically seek a differentiation matrix D and a positive definite matrix H such that a rule similar to

v^{T}HDhv = v^{T}Q_{S}v + v^{T}Q_{A}v = 1

2(−v^{2}_{0}+ v^{2}_{N})
exists.

Notice that the summation-by-part rule is only for

estimating the energy of a system. To stably impose the boundary conditions we still need the penalty methodology.

### Literature Review

Bo Strand (1994) “Summation by parts for finite difference approximations for d/dx.”

M.H. Carpenter, D. Gottlieb and S. Abarbanel (1994):

“Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: methodology and application to high-order compact schemes”

P. Olsson (1995). “Summation by parts, projections, and stability”

B. Strand (1998). “ Numerical studies of hyperbolic IBVP with high-order finite difference operators satisfying a summation by parts rule.”

J. Nordstrom, M.H. Carpenter, (2003). “High-order finite difference methods, multidimensional linear problems and curvilinear coordinates.”

### Literature Review

Bo Strand (1994) “Summation by parts for finite difference approximations for d/dx.”

M.H. Carpenter, D. Gottlieb and S. Abarbanel (1994):

“Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: methodology and application to high-order compact schemes”

P. Olsson (1995). “Summation by parts, projections, and stability”

B. Strand (1998). “ Numerical studies of hyperbolic IBVP with high-order finite difference operators satisfying a summation by parts rule.”

J. Nordstrom, M.H. Carpenter, (2003). “High-order finite difference methods, multidimensional linear problems and curvilinear coordinates.”

### Literature Review

Bo Strand (1994) “Summation by parts for finite difference approximations for d/dx.”

M.H. Carpenter, D. Gottlieb and S. Abarbanel (1994):

“Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: methodology and application to high-order compact schemes”

P. Olsson (1995). “Summation by parts, projections, and stability”

B. Strand (1998). “ Numerical studies of hyperbolic IBVP with high-order finite difference operators satisfying a summation by parts rule.”

J. Nordstrom, M.H. Carpenter, (2003). “High-order finite difference methods, multidimensional linear problems and curvilinear coordinates.”

### Literature Review - Contd.

Applications:

Time-domain Maxwell equations: H. M.Jurgens and D. W.

Zingg (2000).

General Relativity: L Lehner, O Reula, M Tiglio (2005) Compressible Navier-Stokes equations: M. Svärda, M. H.

Carpenterb and J. Nordström (2007) Related Issue - Order Reduction:

M.H. Carpenter, D. Gottlieb, S. Abarbanel, W.-S. Don (1995):

“The theoretical accuracy of Runge-Kutta time discretizations for the initial boundary value problem: a study of the boundary error.”

### Literature Review - Contd.

Applications:

Time-domain Maxwell equations: H. M.Jurgens and D. W.

Zingg (2000).

General Relativity: L Lehner, O Reula, M Tiglio (2005) Compressible Navier-Stokes equations: M. Svärda, M. H.

Carpenterb and J. Nordström (2007) Related Issue - Order Reduction:

M.H. Carpenter, D. Gottlieb, S. Abarbanel, W.-S. Don (1995):

“The theoretical accuracy of Runge-Kutta time discretizations for the initial boundary value problem: a study of the boundary error.”

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### Basic Concepts and Notations I

Let I = [0, 1]. Consider two functions f (x) and g(x) defined on I.

We define the continuous L_{2}inner product and norm for
functions over I as

(f , g) = Z

I

f g dx, ||f ||^{2}_{I} = (f , f )

Consider I^{2}= [0, 1] × [0, 1]. The continuous L_{2}inner product
and norm for functions over I^{2}are defined as

(f , g) = Z

I^{2}

f g dx dy, ||f ||^{2}_{I}2 = (f , f ).

### Basic Concepts and Notations II

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 ∈ V_{L+1},
explicitly given by

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

We define a weighted discrete L_{2}inner product and norm, with
respect to the step size h and the matrix H, for vectors as

(f , g)h,H = hf^{T}Hg, ||f ||^{2}_{h,H} = (f , f )h,H.
If H is an identity matrix then

(f , g)h= hf^{T}g, ||f ||^{2}_{h} = (f , f )h.

Refer to 1D Lemma

### Implicit Finite Difference Operators

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

Pvx = h^{−1}Qv, or vx= Dv = h^{−1}P^{−1}Qv, (1)
where

v = [v_{0}, v_{1}, ..., vL]^{T}, vx= [v_{x0}, v_{x1}, ..., v_{xL}]^{T},
denote the numerical approximations of u and u^{0} evaluated at
the grid points, and D, P, Q ∈ M_{L+1}.

### Summation-By-Parts Properties

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;

SBP Properties

**SBP1:**The matrix P is symmetric positive definite.

**SBP2:** The matrix Q is nearly skew-symmetric and satisfies the
constraint

Q^{S} = Q + Q^{T}

2 =diag(q_{00}, 0, ..., 0, qLL), q_{00}< 0, qLL= −q00.
where q_{00}and qLL are the upper most and lower most diagonal
elements of Q.

Back to the proof of Lemma S1 Back to the proof of Theorem

### Summation-by-Parts in 1D Space I

Lemma S1 (Summation-by-Parts)

Consider the difference operator D = h^{−1}P^{−1}Qwhere P and Q
satisfy**SBP1 and SBP2, respectively. We have**

(v, Dv)_{h,P}= (v, h^{−1}P^{−1}Qv)_{h,P}= q_{00}v^{2}_{0}+ qLLv^{2}_{L},
for v ∈ V_{L+1}.

Back to the proof of Theorem

### Summation-by-Parts in 1D Space II

Proof. Refer to Weighted Discrete L_{2}- norm

First we rewrite the inner product as

(v, Dv)_{h,P}= hv^{T}PDv = v^{T}Qv = v^{T}Q^{S}v + v^{T}Q^{A}v

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^{T}Q^{A}v = 0 since Q^{A} is antisymmetric. Thus, we have

(v, Dv)_{h,P} = v^{T}Q^{S}v = q00v^{2}_{0}+ qLLv^{2}_{L},

where the last equality is due to ** ^{SBP2}** ,. This completes the

proof.

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### New Fifth Order Scheme

Accuracy: interior point (8th order);

boundary point (4th order) globally 5th order

Scheme at interior points:

v_{x,j−3}− 6vx,j−2+ 15vx,j−1+ 120vx,j+ 15vx,j+1− 6vx,j+2+ vx,j+3

=1 h(−7

3v_{j−3}+ 21v_{j−2}− 105v_{j−1}+ 105vj+1− 21vj+2+ 7
3v_{j+3})

P=

21547 729

20560 729

−12734

729 2 1

20560 729

3971857 23328

−300533 23328

445379 23328

−21263

−12734 2592 729

−300533 23328

273185 2592

368689 23328

−43069

23328 1

2 ^{445379}_{23328} ^{368689}_{23328} ^{2828713}_{23328} ^{293771}_{23328} −6 1
1 ^{−21263}_{2592} ^{−43069}_{23328} ^{293771}_{23328} ^{2821361}_{23328} 15 −6 1

1 −6 15 120 15 −6 1

. .. . .. ... ... ...

Q=

−70 ^{2863}_{27} ^{−11267}_{243} ^{2753}_{243} −1

−2863

27 0 ^{28561}_{243} ^{−2395}_{324} ^{−3991}_{972}

11267 243

−28561

243 0 ^{77245}_{972} ^{−10337}_{972} ^{7}_{3}

−2753 243

2395 324

−77245

972 0 ^{2756}_{27} −21 ^{7}_{3}

1 ^{3991}_{972} ^{10337}_{972} ^{−2756}_{27} 0 105 −21 ^{7}_{3}

−^{7}_{3} 21 −105 0 105 −21 ^{7}_{3}
. .. . .. . .. . .. ...

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### One Dimensional Advection Equation

Example

Consider the advection equation

∂u

∂t +∂u

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

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

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

Equation (2) leads to an energy rate d

dt||u||^{2}_{I} = g^{2}(t) − u^{2}(1, t),

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

||u(x, t)||^{2}_{I} ≤ ||u(x, 0)||^{2}_{I} = ||f (x)||^{2}_{I}. (5)

### Numerical Scheme

Semi-discrete scheme (penalty method) dv

dt + h^{−1}P^{−1}Qv = h^{−1}τ q_{00}(v_{0}− g(t))P^{−1}e_{0}, (6a)
v(0) = f = [f (x0), f (x_{1}), ..., f (x_{L})]^{T}, (6b)
where e0= [1, 0, ..., 0]^{T} and τ is the penalty parameter.

Back to Inhomogeneous ODE System

### Stability of The Numerical Scheme

Theorem

Assume that there exists a smooth solution to the one dimen- sional wave problem described by Eqs.(2-4). Then Eq.(6a) is stable at the semi-discrete level provided that

τ ≥ 1.

Moreover, v(t) satisfies the estimate

||v(t)||^{2}_{h,P}≤ ||f ||^{2}_{h,P}.

### Proof of the Theorem

Proof.

Multiplying hv^{T}Pto the scheme and invoking ^{Lemma S1} we obtain
1

2 d

dt||v||^{2}_{h,P}= −(q_{00}v^{2}_{0}+ qLLv^{2}_{L}) + τ v_{0}q_{00}(v_{0}− g(t)).

For the stability analysis, it is sufficient enough to consider the
scheme subject to g(t) = 0. Hence, by ** ^{SBP2}** ,

1 2

d

dt||v||^{2}_{h,P}= q_{00}(τ − 1)v^{2}_{0}− q_{LL}v^{2}_{L}≤ q_{00}(τ − 1)v^{2}_{0}.
So, taking τ ≥ 1 immediately yields a non-increasing energy

rate 1

2 d

dt||v||^{2}_{h,P}≤ 0.

Thus, the scheme is stable. This completes the proof.

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### Tools for 2D Scheme Construction

Definition (Kronecker Product)

Let A be a p × q matrix with its elements denoted by aij, and let Bbe an m × n matrix, then

A ⊗ B =

a_{11}B . . . a_{1q}B
... . .. ...
a_{p1}B . . . apqB

.

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

### Properties of Kronecker Product

Lemma K1

(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 K2

The Associative Property:

(A ⊗ B) ⊗ C = A ⊗ (B ⊗ C) Lemma K3

There exists a permutation matrix R such that
A ⊗ B = R^{T}(B ⊗ A)R

Lemma K4

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 K5

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 K6

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.

### Notations II

Denote δ_{µν} the Kronecker delta function. Define
exi= [δ_{i0}, δ_{i1}, ..., δiL]^{T}, exi∈ VL+1,
eyj= [δ_{j0}, δ_{j1}, ..., δjM]^{T}, eyj∈ V_{M+1},

Note that the sets {exi}^{L}_{i=0}and {eyj}^{M}_{j=0}consist bases of the
vector spaces V_{L+1}and V_{M+1}, respectively. In addition, we
define the sets:

{eij= eyj⊗ exi| 0 ≤ i ≤ L, 0 ≤ j ≤ M}, which consist bases of the vector space V(M+1)×(L+1).

### Framework in 2D Space I

Grid Point

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

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

L, hy = 1

M. (7)

We define

v = [v00v_{10} . . . vL−1MvLM]^{T} =

L

X

i=0 M

X

j=0

vijeij, (8)

in which vijdenote the grid-function values at the grid points.

### Framework in 2D Space II

Inner Product and Norm

Given that Px and Pyare symmetric positive definite, by Lemma K1 and Lemma K5, 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)h_{xy},Pxy = hxhyf^{T}(Py⊗ Px)g, ||f ||^{2}_{h}

xy,Pxy = (f , f )h_{xy},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 = hxh_{y}f^{T}g, ||f ||^{2}_{h}

xy = (f , f )hxy.

### 2D Difference operators

Denote the one dimensional difference operators in x and y directions, respective, by

Dx= h^{−1}_{x} P^{−1}_{x} Q_{x}, Dx, Px, Q_{x} ∈ ML+1

Dy= h^{−1}_{y} P^{−1}_{y} Q_{y}, Dy, Py, Q_{y} ∈ MM+1,

where Pxand Py satisfy**SBP1 and Q**_{x}and Q_{y} satisfy**SBP2. To**
numerically differentiate v, we define the two dimensional
difference operators for approximating ∂/∂x and ∂/∂y,
respectively, as

D_{x} = I_{M+1}⊗ Dx = I_{M+1}⊗ h^{−1}_{x} P^{−1}_{x} Q_{x}, (9a)
D_{y}= Dy⊗ I_{L+1}= h^{−1}_{y} P^{−1}_{y} Q_{y}⊗ I_{L+1}, (9b)
where I_{L+1}∈ M_{L+1} and I_{M+1} ∈ M_{M+1}are identity matrices.

### Summation-by-Parts in 2D Space I

Lemma S2 (Summation-by-Parts in 2D)

Consider the 2D difference operators, Eqs.(9). Let q^{x}_{ii}0 and q^{y}_{jj}0

denote the elements of the matrices, Q_{x} and Q_{y}, respectively.

For v ∈ V(L+1)×(M+1), we have
(v, Dxv)_{h}

xy,Pxy = q^{x}_{00}||v_{0j}||^{2}_{h}

y,Py+ q^{x}_{LL}||vLj||^{2}_{h}

y,Py,
(v, Dyv)_{h}

xy,Pxy = q^{y}_{00}||v_{i0}||^{2}_{h}

x,Px + q^{y}_{MM}||viM||^{2}_{h}

x,Px, where

v0j=

M

X

j=0

v_{0j}eyj, vLj=

M

X

j=0

vLjeyj, vi0=

L

X

i=0

v_{i0}exi, viM =

L

X

i=0

v_{i0}exi.

Skip the proof

### Summation-by-Parts in 2D Space II

Proof.

Denote Q^{S}_{x}= (Q_{x}+ Q^{T}_{x})/2 and Q^{A}_{x} = (Q_{x}− Q^{T}_{x})/2.

Since Py⊗ Q^{A}_{x} = −(P_{y}⊗ Q^{A}_{x})^{T}, v^{T}(P_{y}⊗ Q^{A}_{x})v = 0, we have
(v, D_{x}v)_{h}

xy,Pxy = h_{y}v^{T}(P_{y}⊗ Q^{S}_{x})v + h_{y}v^{T}(P_{y}⊗ Q^{A}_{x})v = h_{y}v^{T}(P_{y}⊗ Q^{S}_{x})v

=

L

X

i,i^{0}=0
M

X

j,j^{0}=0

hyvijvi^{0}j^{0}(e^{T}_{yj}0⊗ e^{T}_{xi}0)(Py⊗ Q^{S}_{x})(eyj⊗ exi).

Applying Lemma K1 we have

(e^{T}_{yj}0⊗ e^{T}_{xi}0)(P_{y}⊗ Q^{S}_{x})(e_{yj}⊗ e_{xi}) = e^{T}_{yj}0P_{y}e_{yj}⊗ e^{T}_{xi}0Q^{S}_{x}e_{xi}= δ_{ii}^{0}q^{x}_{ii}0(e^{T}_{yj}0P_{y}e_{yj})
Therefore, we obtain

hyv^{T}(Py⊗ Q^{S}_{x})v = hyq^{x}_{00}

M

X

j^{0},j=0

v0j^{0}e^{T}_{yj}0Pyeyjv0j+ hyq^{x}_{LL}

M

X

j^{0},j=0

vLj^{0}e^{T}_{yj}0PyeyjvLj,

and prove the first statement.

The second part of the proof is similar.

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### Two Dimensional Advection Equation I

Example

Consider the 2D advection equation

∂u

∂t + ∂u

∂x +∂u

∂y = 0, (x, y) ∈I^{2}, t≥ 0, (10)
with the initial condition at t = 0

u(x, y, 0) = f (x, y), (x, y) ∈I^{2}, (11)
and the boundary condition

u(0, y, t) = g_{1}(y, t), 0 ≤ y ≤ 1, t≥ 0, (12a)
u(x, 0, t) = g2(x, t), 0 ≤ x ≤ 1, t≥ 0. (12b)

### Two Dimensional Advection Equation III

Semi-discrete scheme dv

dt + D_{x}v + D_{y}v = r, (13a)

v(t = 0) = f =

L

X

i=0 M

X

j=0

f(x_{i}, y_{j})e_{ij}. (13b)

where

r = h^{−1}_{x} τ_{1}

M

X

j=0

q^{x}_{00}(v_{0j}− g1(y_{j}, t))(e_{yj}⊗ P^{−1}_{x} e_{x0})

+ h^{−1}_{y} τ2
L

X

i=0

q^{y}_{00}(vi0− g2(xi, t))(P^{−1}_{y} ey0⊗ exi).

### Two Dimensional Advection Equation IV

Theorem

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

τ_{1}, τ_{2}≥ 1.

Moreover, if g = 0, v(t) satisfies the estimate

||v(t)||^{2}_{h}_{xy}_{,P}_{xy} ≤ ||f ||^{2}_{h}_{xy}_{,P}_{xy}. (14)

Skip the proof

### Two Dimensional Advection Equation IV

Proof.

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

1 2

d
dt||v||^{2}_{h}

xy,Pxy=

v,dv

dt

hxy,Pxy

= − (v, Dxv)hxy,Pxy− (v, Dyv)hxy,Pxy+ (v, r)hxy,Pxy. (15)
Notice that with ei^{0}j^{0} = e_{yj}^{0}⊗ exi^{0}

(v, r)_{h}_{xy}_{,P}_{xy}= h_{y}τ_{1}q^{x}_{00}

L

X

i^{0}=0
M

X

j^{0},j=0

vi^{0}j^{0}v0j(e^{T}_{i}0j^{0})(P_{y}⊗ P_{x})(e_{yj}⊗ P^{−1}_{x} e_{x0})

+ hxτ2q^{y}_{00}

L

X

i^{0},i=0
M

X

j^{0}=0

vi^{0}j^{0}vi0(e^{T}_{i}0j^{0})(Py⊗ Px)(P^{−1}_{y} ey0⊗ exi)

Applying Lemma K1 one can easily show that

(e^{T}_{i}0j^{0})(Py⊗ Px)(eyj⊗ P^{−1}_{x} ex0) = δi^{0}0(e^{T}_{yj}0Pyeyj)
(e^{T}_{i}0j^{0})(P_{y}⊗ P_{x})(P^{−1}_{y} e_{y0}⊗ e_{xi}) = δ_{j}^{0}_{0}(e^{T}_{xi}0P_{x}e_{xi})
Therefore,

(v, r)hxy,Pxy= τ1q^{x}_{00}

M

X

j^{0},j=0

hyv0j^{0}v0j(e^{T}_{yj}0Pyeyj) + τ2q^{y}_{00}

L

X

i^{0},i=0

hxvi^{0}0vi0(e^{T}_{xi}0Pxexi)

= τ_{1}q^{x}_{00}||v0j||^{2}_{h}_{y}_{,P}_{y}+ τ_{2}q^{y}_{00}||vi0||^{2}_{h}_{x}_{,P}_{x}

Substituting the above expression into Eq.(15) and invoking Lemma S2, we have

1 2

d
dt||v||^{2}_{h}

xy,Pxy = (τ_{1}− 1)q^{x}_{00}||v_{0j}||^{2}_{P}

y − q^{x}_{LL}||v_{Lj}||^{2}_{P}

y

+ (τ2− 1)q^{y}_{00}||vi0||^{2}_{P}_{x}− q^{y}_{MM}||viM||^{2}_{P}_{x}
Taking τ1, τ_{2} ≥ 1 yields

1 2

d

dt||v||^{2}_{h}_{xy}_{,P}_{xy} ≤ 0,
which implies Eq.(14). Thus, the scheme is stable.

### Outline

**1** SBP Operators
Introduction

Basic Concepts and Notations New Fifth Order SBP Operator

**2** SBP Operators for Wave Problems

One Dimensional Advection Equation - Penalty Methods 2D Difference Operators

Scheme for 2D wave problem

**3** Time Integration Scheme

RK for inhomogeneous system - Order Reduction SSP-RK for inhomogeneous system

**4** Numerical experiment
1D Numerical experiment
2D Numerical experiment

**5** Conclusion

### Time Discretization

Example

Numerical Scheme: (SBP operators) Refer to Penalty method

dv

dt + h^{−1}P^{−1}Qv = h^{−1}τ q00(v0− g(t))P^{−1}e0,
Let us write the semi-discrete schemes as

d

dtv = Lv + Φ(t)

where v(t) = (vj(t))_{1≤j≤K}, L consists the difference operator
modified by the SAT term (penalty method), and Φ contains the
time explicit boundary conditions.

How do we implement a RK method?

Example: Classical RK4

k1= L(v^{n}, nht) =Lv^{n}+ Φ(tn)
k2= L

v^{n}+1

2htk1, n + 1 2ht

=L(v^{n}+1

2htk1) + Φ(t_{n+1/2})
k3= L

v^{n}+1

2htk2, n + 1 2ht

=L(v^{n}+1

2htk2) + Φ(t_{n+1/2})
k_{4}= L (v^{n}+ h_{t}k_{3}, (n + 1)ht) =L(v^{n}+ h_{t}k_{3}) + Φ(t_{n+1})
v^{n+1}= v^{n}+ht

6(k1+ 2k2+ 2k3+ k4).

Numerical Result

N Error Ratio Order

20 1.544e − 2 - -

40 2.013e − 3 7.67 2.94 80 2.543e − 4 7.92 2.99

Table: Classical RK4

How do we implement a RK method?

Example: Classical RK4

k1= L(v^{n}, nht) =Lv^{n}+ Φ(tn)
k2= L

v^{n}+1

2htk1, n + 1 2ht

=L(v^{n}+1

2htk1) + Φ(t_{n+1/2})
k3= L

v^{n}+1

2htk2, n + 1 2ht

=L(v^{n}+1

2htk2) + Φ(t_{n+1/2})
k_{4}= L (v^{n}+ h_{t}k_{3}, (n + 1)ht) =L(v^{n}+ h_{t}k_{3}) + Φ(t_{n+1})
v^{n+1}= v^{n}+ht

6(k1+ 2k2+ 2k3+ k4).

Numerical Result

N Error Ratio Order

20 1.544e − 2 - -

40 2.013e − 3 7.67 2.94 80 2.543e − 4 7.92 2.99

Table: Classical RK4

### Implementation of Classical RK4-II

ODE system: d

dtu = Du + φ(t)e0. Numerical solution:

k_{1}= Dv^{n}+ g_{0}e_{0}
k_{2}= D(v^{n}+1

2htk_{1}) + g_{1}e_{0}
k3= D(v^{n}+1

2h_{t}k_{2}) + g2e0

k_{4}= D(v^{n}+ htk_{3}) + g_{3}e_{0}
v^{n+1}= v^{n}+ht

6(k_{1}+ 2k_{2}+ 2k_{3}+ k_{4}).

### Implementation of Classical RK4-III

Approximation Solution:

uapp(ht) =

4

X

l=0

(Dht)^{l}
l! u^{n}+

4

X

l=1 l−1

X

i=0

D^{l−i−1}h^{l}_{t}

l! φ^{(i)}(tn) e_{0}
Numerical solution:

v^{n+1}=

4

X

l=0

(Dht)^{l}

l! v^{n}+ (1
6h_{t}+1

6h^{2}_{t}D+ 1

12h^{3}_{t}D^{2}+ 1

24h^{4}_{t}D^{3})g_{0}e_{0}
+ (1

3h_{t}+1

6h^{2}_{t}D+ 1

12h^{3}_{t}D^{2})g1e0

+ (1 3ht+1

6h^{2}_{t}D)g_{2}e_{0}+ (1

6ht)g_{3}e_{0}