• 沒有找到結果。

# C -C C ,M -H C ,C -H T High-OrderSummation-by-PartsImplicitDifferenceOperatorsforWaveProblems

N/A
N/A
Protected

Share "C -C C ,M -H C ,C -H T High-OrderSummation-by-PartsImplicitDifferenceOperatorsforWaveProblems"

Copied!
73
0
0

(1)

### 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

(2)

### 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

(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

(4)

### 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.

(5)

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

∂u2

∂t dx= −1 2

Z 1 0

∂u2

∂x dx

= 1

2(u2(0, t) − u2(1, t)) = 1

2(g2(t) − u2(1, t))

(6)

### Central Difference Scheme

Strongly Enforced Boundary Condition

dvi

dt +vi+1− vi−1

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

dt +vN− vN−1

h = 0

v0(t) = g(t) vi(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− v2N) Recall the continuous system

1 2

dE(t) dt = 1

2(g2(t) − u2(1, t))

(7)

### Central Difference Scheme

Weakly Enforced Boundary Condition- Penalty Method

Semi-discrete scheme (penalty method):

dvi

dt +vi+1− vi−1

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

dt +vN− vN−1

h = 0

dv0

dt +v1− v0

h = −τ (v0− g(t)) Energy estimate:

dED

dt = 1 2v0

dv0

dth+

N−1

X

i=1

vi

dvi

dth+1 2vN

dvN

dt h

(8)

### Central Difference Scheme

Weakly Enforced Boundary Condition

d dt

N

X

i=0

civ2ih= v1v0− v2N− v1v0+ v20− τ hv20+ τ hv0g(t)

= −v2N+ (1 − τ h)v20+ τ hv0g(t)

= −v2N+ (1 − τ h)

"

v20+ τ h

1 − τ hv0g(t) +

 τ h 2(1 − τ h)

2 g2(t)

#

− (1 − τ h)

 τ h 2(1 − τ h)

2

g2(t)

If τ h = 2

dED

dt = −v2N− (v0− g(t))2+ g2(t)

(9)

### 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

v0 v1 v2 .. . vN−1

vN 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

v0 v1 v2 .. . vN−1

vN 3 7 7 7 7 7 7 7 7 7 7 5

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

−τ (v0− g(t)) 0 0 .. . 0 0

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

or d

dtv(t) + Dv(t) = −τ (v0− g(t))e0, e0= [1 0 0 · · · 0 0]T Define H = diag(1

2, 1, 1, · · · , 1, 1 2) vTHdv

dth+ vTH D h v(t) = −τ (v0− g(t))vTH e0h

(10)

### 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

=QS+ QA

So

vTHDhv = vTQSv + vTQAv = 1

2(−v20+ v2N)

(11)

### Central Difference Scheme

Interesting Property: Summation-by-Parts Rule

Hence

vTHdv

dth= −vTH D h v − τ (v0− g(t))vTH e0h 1

2 d

dtvTHvh = 1 2v20−1

2v2N−τ h

2 v0(v0− g(t)) The rule

vTHDhv = vTQSv + vTQAv = 1

2(−v20+ v2N) in fact mimics the action of

Z 1 0

u(x)∂u(x)

∂x dx= 1

2u2(1) − 1 2u2(0)

(12)

### 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

vTHDhv = vTQSv + vTQAv = 1

2(−v20+ v2N) 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.

(13)

### 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.”

(14)

### 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.”

(15)

### 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.”

(16)

### 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.”

(17)

### 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.”

(18)

### 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

(19)

### Basic Concepts and Notations I

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

We define the continuous L2inner product and norm for functions over I as

(f , g) = Z

I

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

Consider I2= [0, 1] × [0, 1]. The continuous L2inner product and norm for functions over I2are defined as

(f , g) = Z

I2

f g dx dy, ||f ||2I2 = (f , f ).

(20)

### 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 ∈ VL+1, explicitly given by

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

We define a weighted discrete L2inner product and norm, with respect to the step size h and the matrix H, for vectors as

(f , g)h,H = hfTHg, ||f ||2h,H = (f , f )h,H. If H is an identity matrix then

(f , g)h= hfTg, ||f ||2h = (f , f )h.

Refer to 1D Lemma

(21)

### 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−1Qv, or vx= Dv = h−1P−1Qv, (1) where

v = [v0, v1, ..., vL]T, vx= [vx0, vx1, ..., vxL]T, denote the numerical approximations of u and u0 evaluated at the grid points, and D, P, Q ∈ ML+1.

(22)

### 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

QS = Q + QT

2 =diag(q00, 0, ..., 0, qLL), q00< 0, qLL= −q00. where q00and 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

(23)

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

Lemma S1 (Summation-by-Parts)

Consider the difference operator D = h−1P−1Qwhere P and Q satisfySBP1 and SBP2, respectively. We have

(v, Dv)h,P= (v, h−1P−1Qv)h,P= q00v20+ qLLv2L, for v ∈ VL+1.

Back to the proof of Theorem

(24)

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

Proof. Refer to Weighted Discrete L2- norm

First we rewrite the inner product as

(v, Dv)h,P= hvTPDv = vTQv = vTQSv + vTQAv

where QS = (Q + QT)/2 and QA = (Q − QT)/2 are, respectively, the symmetric and anti-symmetric parts of the matrix Q. Notice that vTQAv = 0 since QA is antisymmetric. Thus, we have

(v, Dv)h,P = vTQSv = q00v20+ qLLv2L,

where the last equality is due to SBP2 ,. This completes the

proof. 

(25)

### 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

(26)

### New Fifth Order Scheme

Accuracy: interior point (8th order);

boundary point (4th order) globally 5th order

Scheme at interior points:

vx,j−3− 6vx,j−2+ 15vx,j−1+ 120vx,j+ 15vx,j+1− 6vx,j+2+ vx,j+3

=1 h(−7

3vj−3+ 21vj−2− 105vj−1+ 105vj+1− 21vj+2+ 7 3vj+3)

(27)

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 44537923328 36868923328 282871323328 29377123328 −6 1 1 −212632592 −4306923328 29377123328 282136123328 15 −6 1

1 −6 15 120 15 −6 1

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

Q=

−70 286327 −11267243 2753243 −1

−2863

27 0 28561243 −2395324 −3991972

11267 243

−28561

243 0 77245972 −10337972 73

−2753 243

2395 324

−77245

972 0 275627 −21 73

1 3991972 10337972 −275627 0 105 −21 73

73 21 −105 0 105 −21 73 . .. . .. . .. . .. ...

(28)

### 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

(29)

Example

∂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||2I = g2(t) − u2(1, t),

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

||u(x, t)||2I ≤ ||u(x, 0)||2I = ||f (x)||2I. (5)

(30)

### Numerical Scheme

Semi-discrete scheme (penalty method) dv

dt + h−1P−1Qv = h−1τ q00(v0− g(t))P−1e0, (6a) v(0) = f = [f (x0), f (x1), ..., f (xL)]T, (6b) where e0= [1, 0, ..., 0]T and τ is the penalty parameter.

Back to Inhomogeneous ODE System

(31)

### 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)||2h,P≤ ||f ||2h,P.

(32)

### Proof of the Theorem

Proof.

Multiplying hvTPto the scheme and invoking Lemma S1 we obtain 1

2 d

dt||v||2h,P= −(q00v20+ qLLv2L) + τ v0q00(v0− 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||2h,P= q00(τ − 1)v20− qLLv2L≤ q00(τ − 1)v20. So, taking τ ≥ 1 immediately yields a non-increasing energy

rate 1

2 d

dt||v||2h,P≤ 0.

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

(33)

### 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

(34)

### 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 =

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

.

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

(35)

### Properties of Kronecker Product

Lemma K1

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

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 = RT(B ⊗ A)R

(36)

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.

(37)

### Notations II

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

Note that the sets {exi}Li=0and {eyj}Mj=0consist bases of the vector spaces VL+1and VM+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).

(38)

### Framework in 2D Space I

Grid Point

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

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

L, hy = 1

M. (7)

We define

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

L

X

i=0 M

X

j=0

vijeij, (8)

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

(39)

### 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)hxy,Pxy = hxhyfT(Py⊗ Px)g, ||f ||2h

xy,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 = hxhyfTg, ||f ||2h

xy = (f , f )hxy.

(40)

### 2D Difference operators

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

Dx= h−1x P−1x Qx, Dx, Px, Qx ∈ ML+1

Dy= h−1y P−1y Qy, Dy, Py, Qy ∈ MM+1,

where Pxand Py satisfySBP1 and Qxand Qy satisfySBP2. To numerically differentiate v, we define the two dimensional difference operators for approximating ∂/∂x and ∂/∂y, respectively, as

Dx = IM+1⊗ Dx = IM+1⊗ h−1x P−1x Qx, (9a) Dy= Dy⊗ IL+1= h−1y P−1y Qy⊗ IL+1, (9b) where IL+1∈ ML+1 and IM+1 ∈ MM+1are identity matrices.

(41)

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

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

Consider the 2D difference operators, Eqs.(9). Let qxii0 and qyjj0

denote the elements of the matrices, Qx and Qy, respectively.

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

xy,Pxy = qx00||v0j||2h

y,Py+ qxLL||vLj||2h

y,Py, (v, Dyv)h

xy,Pxy = qy00||vi0||2h

x,Px + qyMM||viM||2h

x,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.

Skip the proof

(42)

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

Proof.

Denote QSx= (Qx+ QTx)/2 and QAx = (Qx− QTx)/2.

Since Py⊗ QAx = −(Py⊗ QAx)T, vT(Py⊗ QAx)v = 0, we have (v, Dxv)h

xy,Pxy = hyvT(Py⊗ QSx)v + hyvT(Py⊗ QAx)v = hyvT(Py⊗ QSx)v

=

L

X

i,i0=0 M

X

j,j0=0

hyvijvi0j0(eTyj0⊗ eTxi0)(Py⊗ QSx)(eyj⊗ exi).

Applying Lemma K1 we have

(eTyj0⊗ eTxi0)(Py⊗ QSx)(eyj⊗ exi) = eTyj0Pyeyj⊗ eTxi0QSxexi= δii0qxii0(eTyj0Pyeyj) Therefore, we obtain

hyvT(Py⊗ QSx)v = hyqx00

M

X

j0,j=0

v0j0eTyj0Pyeyjv0j+ hyqxLL

M

X

j0,j=0

vLj0eTyj0PyeyjvLj,

and prove the first statement.

The second part of the proof is similar.

(43)

### 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

(44)

### Two Dimensional Advection Equation I

Example

∂u

∂t + ∂u

∂x +∂u

∂y = 0, (x, y) ∈I2, t≥ 0, (10) with the initial condition at t = 0

u(x, y, 0) = f (x, y), (x, y) ∈I2, (11) and the boundary condition

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

(45)

### Two Dimensional Advection Equation III

Semi-discrete scheme dv

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

v(t = 0) = f =

L

X

i=0 M

X

j=0

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

where

r = h−1x τ1

M

X

j=0

qx00(v0j− g1(yj, t))(eyj⊗ P−1x ex0)

+ h−1y τ2 L

X

i=0

qy00(vi0− g2(xi, t))(P−1y ey0⊗ exi).

(46)

### 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)||2hxy,Pxy ≤ ||f ||2hxy,Pxy. (14)

Skip the proof

(47)

### 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||2h

xy,Pxy=

 v,dv

dt



hxy,Pxy

= − (v, Dxv)hxy,Pxy− (v, Dyv)hxy,Pxy+ (v, r)hxy,Pxy. (15) Notice that with ei0j0 = eyj0⊗ exi0

(v, r)hxy,Pxy= hyτ1qx00

L

X

i0=0 M

X

j0,j=0

vi0j0v0j(eTi0j0)(Py⊗ Px)(eyj⊗ P−1x ex0)

+ hxτ2qy00

L

X

i0,i=0 M

X

j0=0

vi0j0vi0(eTi0j0)(Py⊗ Px)(P−1y ey0⊗ exi)

(48)

Applying Lemma K1 one can easily show that

(eTi0j0)(Py⊗ Px)(eyj⊗ P−1x ex0) = δi00(eTyj0Pyeyj) (eTi0j0)(Py⊗ Px)(P−1y ey0⊗ exi) = δj00(eTxi0Pxexi) Therefore,

(v, r)hxy,Pxy= τ1qx00

M

X

j0,j=0

hyv0j0v0j(eTyj0Pyeyj) + τ2qy00

L

X

i0,i=0

hxvi00vi0(eTxi0Pxexi)

= τ1qx00||v0j||2hy,Py+ τ2qy00||vi0||2hx,Px

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

1 2

d dt||v||2h

xy,Pxy = (τ1− 1)qx00||v0j||2P

y − qxLL||vLj||2P

y

+ (τ2− 1)qy00||vi0||2Px− qyMM||viM||2Px Taking τ1, τ2 ≥ 1 yields

1 2

d

dt||v||2hxy,Pxy ≤ 0, which implies Eq.(14). Thus, the scheme is stable.

(49)

### 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

(50)

### Time Discretization

Example

Numerical Scheme: (SBP operators) Refer to Penalty method

dv

dt + h−1P−1Qv = h−1τ q00(v0− g(t))P−1e0, 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.

(51)

How do we implement a RK method?

Example: Classical RK4

k1= L(vn, nht) =Lvn+ Φ(tn) k2= L

 vn+1

2htk1, n + 1 2ht



=L(vn+1

2htk1) + Φ(tn+1/2) k3= L

 vn+1

2htk2, n + 1 2ht



=L(vn+1

2htk2) + Φ(tn+1/2) k4= L (vn+ htk3, (n + 1)ht) =L(vn+ htk3) + Φ(tn+1) vn+1= vn+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

(52)

How do we implement a RK method?

Example: Classical RK4

k1= L(vn, nht) =Lvn+ Φ(tn) k2= L

 vn+1

2htk1, n + 1 2ht



=L(vn+1

2htk1) + Φ(tn+1/2) k3= L

 vn+1

2htk2, n + 1 2ht



=L(vn+1

2htk2) + Φ(tn+1/2) k4= L (vn+ htk3, (n + 1)ht) =L(vn+ htk3) + Φ(tn+1) vn+1= vn+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

(53)

### Implementation of Classical RK4-II

ODE system: d

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

k1= Dvn+ g0e0 k2= D(vn+1

2htk1) + g1e0 k3= D(vn+1

2htk2) + g2e0

k4= D(vn+ htk3) + g3e0 vn+1= vn+ht

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

(54)

### Implementation of Classical RK4-III

Approximation Solution:

uapp(ht) =

4

X

l=0

(Dht)l l! un+

4

X

l=1 l−1

X

i=0

Dl−i−1hlt

l! φ(i)(tn) e0 Numerical solution:

vn+1=

4

X

l=0

(Dht)l

l! vn+ (1 6ht+1

6h2tD+ 1

12h3tD2+ 1

24h4tD3)g0e0 + (1

3ht+1

6h2tD+ 1

12h3tD2)g1e0

+ (1 3ht+1

6h2tD)g2e0+ (1

6ht)g3e0

Nicolas Scaringella, Giorgio Zoia, Daniel Mlynek, IEEE SIGNAL PROCESSING MAGAZINE.. IEEE SIGNAL PROCESSING MAGAZINE

– During each time unit, the backlight device at none(C), low light display (DL), and high light display (HL) consumes E C , E DL , and E DH of energy, respectively.. • Of course

• Theorem: boosting tends to increase margins boosting tends to increase margins boosting tends to increase margins boosting tends to increase margins of training boosting tends

(A)Universal Serial Bus(USB)(B)High Definition Multimedia Interface(HDMI) (C)Video Graphics Array(VGA)(D)Integrated Drive

We compare the results of analytical and numerical studies of lattice 2D quantum gravity, where the internal quantum metric is described by random (dynamical)

(C)John’s love for graffiti began in junior high s chool when his school held a graffiti contest.. (D)John’s love for crows start s in junior high school when he joined

Methods involving finite differences for solving boundary-value problems replace each of the derivatives in the differential equation by an appropriate