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

(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

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

∂u²

∂t dx= −1 2

Z 1 0

∂u²

∂x dx

= 1

2(u²(0, t) − u²(1, t)) = 1

2(g²(t) − u²(1, t))

(6)

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²_N) Recall the continuous system

1 2

dE(t) dt = 1

2(g²(t) − u²(1, t))

(7)

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₀

h = −τ (v₀− 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

civ²_ih= v1v0− v²_N− v1v0+ v²₀− τ hv²₀+ τ hv0g(t)

= −v²_N+ (1 − τ h)v²₀+ τ hv₀g(t)

= −v²_N+ (1 − τ h)

"

v²₀+ τ h

1 − τ hv0g(t) +

τ h 2(1 − τ h)

² g²(t)

#

− (1 − τ h)

τ h 2(1 − τ h)

2

g²(t)

If τ h = 2

dE_D

dt = −v²_N− (v₀− g(t))²+ g²(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

v₀ v₁ v₂ .. . 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₀ v₁ v₂ .. . 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₀− 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₀− g(t))e₀, e₀= [1 0 0 · · · 0 0]^T Define H = diag(1

2, 1, 1, · · · , 1, 1 2) v^THdv

dth+ v^TH D h v(t) = −τ (v0− g(t))v^TH 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

=Q_S+ Q_A

So

v^THDhv = v^TQ_Sv + v^TQ_Av = 1

2(−v²₀+ v²_N)

(11)

Central Difference Scheme

Interesting Property: Summation-by-Parts Rule

Hence

v^THdv

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

2 d

dtv^THvh = 1 2v²₀−1

2v²_N−τ h

2 v₀(v₀− g(t)) The rule

2(−v²₀+ v²_N) in fact mimics the action of

Z 1 0

u(x)∂u(x)

∂x dx= 1

2u²(1) − 1 2u²(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

2(−v²₀+ v²_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.

(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

(15)

Literature Review

(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

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 L₂inner product and norm for functions over I as

(f , g) = Z

I

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

Consider I²= [0, 1] × [0, 1]. The continuous L₂inner product and norm for functions over I²are defined as

(f , g) = Z

I²

f g dx dy, ||f ||²_I2 = (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 ∈ V_L+1, explicitly given by

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

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

(f , g)h,H = hf^THg, ||f ||²_h,H = (f , f )h,H. If H is an identity matrix then

(f , g)h= hf^Tg, ||f ||²_h = (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⁻¹Qv, or vx= Dv = h⁻¹P⁻¹Qv, (1) where

v = [v₀, v₁, ..., 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 ∈ M_L+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

Q^S = Q + Q^T

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

(23)

Summation-by-Parts in 1D Space I

Lemma S1 (Summation-by-Parts)

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

(v, Dv)_h,P= (v, h⁻¹P⁻¹Qv)_h,P= q₀₀v²₀+ qLLv²_L, for v ∈ V_L+1.

Back to the proof of Theorem

(24)

Summation-by-Parts in 1D Space II

Proof. Refer to Weighted Discrete L₂- norm

First we rewrite the inner product as

(v, Dv)_h,P= hv^TPDv = v^TQv = v^TQ^Sv + v^TQ^Av

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

(v, Dv)_h,P = v^TQ^Sv = q00v²₀+ qLLv²_L,

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

proof.

(25)

Outline

5 Conclusion

(26)

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)

(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 ⁴⁴⁵³⁷⁹₂₃₃₂₈ ³⁶⁸⁶⁸⁹₂₃₃₂₈ ^2828713₂₃₃₂₈ ²⁹³⁷⁷¹₂₃₃₂₈ −6 1 1 ⁻²¹²⁶³₂₅₉₂ ⁻⁴³⁰⁶⁹₂₃₃₂₈ ²⁹³⁷⁷¹₂₃₃₂₈ ^2821361₂₃₃₂₈ 15 −6 1

1 −6 15 120 15 −6 1

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







Q=







−70 ²⁸⁶³₂₇ ⁻¹¹²⁶⁷₂₄₃ ²⁷⁵³₂₄₃ −1

−2863

27 0 ²⁸⁵⁶¹₂₄₃ ⁻²³⁹⁵₃₂₄ ⁻³⁹⁹¹₉₇₂

11267 243

−28561

243 0 ⁷⁷²⁴⁵₉₇₂ ⁻¹⁰³³⁷₉₇₂ ⁷₃

−2753 243

2395 324

−77245

972 0 ²⁷⁵⁶₂₇ −21 ⁷₃

1 ³⁹⁹¹₉₇₂ ¹⁰³³⁷₉₇₂ ⁻²⁷⁵⁶₂₇ 0 105 −21 ⁷₃

−⁷₃ 21 −105 0 105 −21 ⁷₃ . .. . .. . .. . .. ...







(28)

Outline

5 Conclusion

(29)

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||²_I = g²(t) − u²(1, t),

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. (5)

(30)

Numerical Scheme

Semi-discrete scheme (penalty method) dv

dt + h⁻¹P⁻¹Qv = h⁻¹τ q₀₀(v₀− g(t))P⁻¹e₀, (6a) v(0) = f = [f (x0), f (x₁), ..., f (x_L)]^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 dimensional 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)||²_h,P≤ ||f ||²_h,P.

(32)

Proof of the Theorem

Proof.

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

2 d

dt||v||²_h,P= −(q₀₀v²₀+ qLLv²_L) + τ v₀q₀₀(v₀− 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||²_h,P= q₀₀(τ − 1)v²₀− q_LLv²_L≤ q₀₀(τ − 1)v²₀. So, taking τ ≥ 1 immediately yields a non-increasing energy

rate 1

2 d

dt||v||²_h,P≤ 0.

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

(33)

Outline

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 =







a₁₁B . . . a_1qB ... . .. ... a_p1B . . . 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= 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

(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∈ V_M+1,

Note that the sets {exi}^L_i=0and {eyj}^M_j=0consist bases of the vector spaces V_L+1and 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).

(38)

Framework in 2D Space I

Grid Point

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. (7)

We define

v = [v00v₁₀ . . . 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)h_xy,Pxy = hxhyf^T(Py⊗ Px)g, ||f ||²_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_yf^Tg, ||f ||²_h

xy = (f , f )hxy.

(40)

2D Difference operators

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 Py satisfySBP1 and Q_xand Q_y satisfySBP2. 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⁻¹_x P⁻¹_x Q_x, (9a) D_y= Dy⊗ I_L+1= h⁻¹_y P⁻¹_y Q_y⊗ I_L+1, (9b) where I_L+1∈ M_L+1 and I_M+1 ∈ M_M+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 q^x_ii0 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)_h

xy,Pxy = q^x₀₀||v_0j||²_h

y,Py+ q^x_LL||vLj||²_h

y,Py, (v, Dyv)_h

xy,Pxy = q^y₀₀||v_i0||²_h

x,Px + q^y_MM||viM||²_h

x,Px, where

v0j=

M

X

j=0

v_0jeyj, vLj=

M

X

j=0

vLjeyj, vi0=

L

X

i=0

v_i0exi, viM =

L

X

i=0

v_i0exi.

Skip the proof

(42)

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Â_x = −(P_y⊗ QÂ_x)^T, v^T(P_y⊗ QÂ_x)v = 0, we have (v, D_xv)_h

xy,Pxy = h_yv^T(P_y⊗ Q^S_x)v + h_yv^T(P_y⊗ Q^A_x)v = h_yv^T(P_y⊗ 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).

Applying Lemma K1 we have

(e^T_yj0⊗ e^T_xi0)(P_y⊗ Q^S_x)(e_yj⊗ e_xi) = e^T_yj0P_ye_yj⊗ e^T_xi0Q^S_xe_xi= δ_ii⁰q^x_ii0(e^T_yj0P_ye_yj) 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.

(43)

Outline

5 Conclusion

(44)

Two Dimensional Advection Equation I

Example

Consider the 2D advection equation

∂u

∂t + ∂u

∂x +∂u

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

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

u(0, y, t) = g₁(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 + D_xv + D_yv = 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⁻¹_x τ₁

M

X

j=0

q^x₀₀(v_0j− g1(y_j, t))(e_yj⊗ P⁻¹_x e_x0)

+ h⁻¹_y τ2 L

X

i=0

q^y₀₀(vi0− g2(xi, t))(P⁻¹_y 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.

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

||v(t)||²_h_xy_,P_xy ≤ ||f ||²_h_xy_,P_xy. (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||²_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⁰j⁰ = e_yj⁰⊗ exi⁰

(v, r)_h_xy_,P_xy= h_yτ₁q^x₀₀

L

X

i⁰=0 M

X

j⁰,j=0

vi⁰j⁰v0j(e^T_i0j⁰)(P_y⊗ P_x)(e_yj⊗ P⁻¹_x e_x0)

+ hxτ2q^y₀₀

L

X

i⁰,i=0 M

X

j⁰=0

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

(48)

Applying Lemma K1 one can easily show that

(e^T_i0j⁰)(Py⊗ Px)(eyj⊗ P⁻¹_x ex0) = δi⁰0(e^T_yj0Pyeyj) (e^T_i0j⁰)(P_y⊗ P_x)(P⁻¹_y e_y0⊗ e_xi) = δ_j⁰₀(e^T_xi0P_xe_xi) 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)

= τ₁q^x₀₀||v0j||²_h_y_,P_y+ τ₂q^y₀₀||vi0||²_h_x_,P_x

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

1 2

d dt||v||²_h

xy,Pxy = (τ₁− 1)q^x₀₀||v_0j||²_P

y − q^x_LL||v_Lj||²_P

y

+ (τ2− 1)q^y₀₀||vi0||²_P_x− q^y_MM||viM||²_P_x Taking τ1, τ₂ ≥ 1 yields

1 2

d

dt||v||²_h_xy_,P_xy ≤ 0, which implies Eq.(14). Thus, the scheme is stable.

(49)

Outline

5 Conclusion

(50)

Time Discretization

Example

Numerical Scheme: (SBP operators) Refer to Penalty method

dv

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

(51)

How do we implement a RK method?

Example: Classical RK4

k1= L(vⁿ, nht) =Lvⁿ+ Φ(tn) k2= L

vⁿ+1

2htk1, n + 1 2ht

=L(vⁿ+1

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

vⁿ+1

2htk2, n + 1 2ht

=L(vⁿ+1

2htk2) + Φ(t_n+1/2) k₄= L (vⁿ+ h_tk₃, (n + 1)ht) =L(vⁿ+ h_tk₃) + Φ(t_n+1) vⁿ⁺¹= vⁿ+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(vⁿ, nht) =Lvⁿ+ Φ(tn) k2= L

vⁿ+1

2htk1, n + 1 2ht

=L(vⁿ+1

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

vⁿ+1

2htk2, n + 1 2ht

=L(vⁿ+1

2htk2) + Φ(t_n+1/2) k₄= L (vⁿ+ h_tk₃, (n + 1)ht) =L(vⁿ+ h_tk₃) + Φ(t_n+1) vⁿ⁺¹= vⁿ+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:

k₁= Dvⁿ+ g₀e₀ k₂= D(vⁿ+1

2htk₁) + g₁e₀ k3= D(vⁿ+1

2h_tk₂) + g2e0

k₄= D(vⁿ+ htk₃) + g₃e₀ vⁿ⁺¹= vⁿ+ht

6(k₁+ 2k₂+ 2k₃+ k₄).

(54)

Implementation of Classical RK4-III

Approximation Solution:

uapp(ht) =

4

X

l=0

(Dht)^l l! uⁿ+

4

X

l=1 l−1

X

i=0

D^l−i−1h^l_t

l! φ⁽ⁱ⁾(tn) e₀ Numerical solution:

vⁿ⁺¹=

4

X

l=0

(Dht)^l

l! vⁿ+ (1 6h_t+1

6h²_tD+ 1

12h³_tD²+ 1

24h⁴_tD³)g₀e₀ + (1

3h_t+1

6h²_tD+ 1

12h³_tD²)g1e0

+ (1 3ht+1

6h²_tD)g₂e₀+ (1

6ht)g₃e₀