-H C ,C -H T TimeIntegrationSchemeM High-OrderFiniteDifferenceMethodsforWaveEquations

(1)

High-Order Finite Difference Methods for Wave Equations

Time Integration Scheme

MIN-HUNGCHEN, CHUN-HAO TENG Department of Mathematics

National Cheng Kung University Tainan City 701, Taiwan

(2)

Outline

1 Time Integration Scheme (I) Preliminaries

Stability of Runge-Kutta Methods High-Order Runge-Kutta methods Numerical Experiment

Reference

2 Time Integration Scheme (II) RK for inhomogeneous system Numerical experiment

Reference

(3)

Outline

Reference

(4)

1D Advection Equation

Example

Consider the 1D advection equation







u_t+ ux = 0 x∈ [0, 1] t≥ 0, u(x, 0) = f (x) x∈ [0, 1],

u(0, t) = u(1, t) t≥ 0, which has the exact solution

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

Refer to Numerical Scheme

(5)

Method of Lines (MOL) Discretization

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

vxj(t) = v_j+1(t) − v_j−1(t) 2hx

where hx= 1/L and

vj(t) ≈ u(xj, t), vxj(t) ≈ du dx(xj, t),

denote the numerical approximations of u and ^du_dx evaluated at the grid points, xj = jhx.

(6)

Method of Lines (MOL) Discretization (Contd.)

Let v(t) = [v₀(t), ..., vL(t)]^T, vx(t) = [v_x0(t), ..., v_xL(t)]^T. We can discretize the PDE by

dv

dt = Lv(t), where

L = −1 2hx







0 1 −1

−1 0 1

. .. ... ...

−1 0 1

1 −1 0







Then, we can solve the system numerically using numerical methods for ODE systems. Refer to PDE

(7)

Approximation of the Linear Operator

In general, we can write the method in a vector form dv

dt = Lv(t), which has the exact solution

v(t) = e^Ltv(0),

where e^Ltis the matrix exponential. The exponential e^zcan be approximated by an linear operator, for example, a finite Taylor series

e^Lh^t ' K(L, ht) =

m

X

i=0

(Lht)ⁱ i!

(8)

Stability of Linear Operators

The finite difference approximation is

v(t + ht) = vⁿ⁺¹= K(L, ht)vⁿ,

where t = nh_t, with step size h_t, and the matrix K(L, h_t)is the approximation to the matrix exponential e^Lh^t. The fully discrete method is strongly stable provided that

K(L, ht)n

≤ C(T).

A sufficient condition for strong stability is that

K(L, ht)

≤ 1 + κht,

for some bounded κ and all sufficiently small values of ht. In practice, condition is the above with κ = 0. Refer to ROS

(9)

Criteria to Choosing the Time Integration Scheme

Two classes of popular methods:

Linear Multi-step Methods.

Runge-Kutta Methods Criteria:

Accuracy Stability Efficiency Storage

(10)

Linear Multi-step Methods (LMM)

In general, an s-step LMM has the form

s

X

j=0

αjv^n+j = ht s

X

j=0

βjL(v^n+j, tn+j).

The value v^n+sis computed from this equation in terms of the previous values v^n+s−1, v^n+s−2,..., vⁿ and L values at these points (which can be stored and re-used if L is expensive to evaluate).

(11)

Runge-Kutta Methods (RK)

An s-stage explicit Runge-Kutta method is k1 = L(vⁿ, nht)

ki = L



vⁿ+ ht i−1

X

j=1

a_ijki, (n + ci)ht





vⁿ⁺¹= vⁿ+ ht s

X

i=1

b_iki,

where the choice of the constants aij, ci and bi determines the accuracy and efficiency of the overall method.

(12)

RK v.s. LMM

RK methods have several advantages over LMM:

The methods are self-starting.

The time step ht can be changed at any point.

RK requires less storage (for high order methods.) and some disadvantages:

Runge-Kutta methods require evaluating L several times each time step.

(Butcher barriers) Higher-order RK methods require more stages. (for u⁰(t) = f (u, t))

order p 1 2 3 4 5 6 7 8

Stages 1 2 3 4 5 6 7 8 9 10 11

Table:Number of stages to achieve a specified order

(13)

Outline

Reference

(14)

Stability of Runge-Kutta Methods

To study the linear stability of these methods we consider a linear scalar equation

du dt = λ u,

for which the general s-stage s-order method can be expressed as a truncated Taylor expansion of the exponential functions

vⁿ⁺¹=

s

X

i=1

(λ ht)ⁱ i! vⁿ.

Thus, the method is stable for a region of λ ht for which

s

X

i=1

(λ ht)ⁱ i!

≤ 1.

Stability

(15)

Stability of Runge-Kutta Methods-Linear System

Consider the linear system ^du

dt = A uwhere A is a constant m× m matrix. Suppose that A is diagonalizable and we have

A= RΛR⁻¹ and Λ = R⁻¹AR,

where Λ = diag(λ₁, ..., λm)is the diagonal matrix of eigenvalues and R is the matrix of eigenvectors.

Let w = R⁻¹uand we have dw

dt = Λ w

or dwp

dt = λpw_p, 1 ≤ p ≤ m.

The overall method is stable if each of the scalar problems is stable.

(16)

Stability regions for some explicit RK methods

Figure:Stability regions for Runge-Kutta methods (

Ps i=1

zⁱ i!

≤ 1.)

(17)

Outline

Reference

(18)

High-Order Runge-Kutta Methods

Classical fourth-order accurate, four-stage Runge-Kutta method

Low-storage Runge-Kutta methods

Linear m-stage m-order Runge-Kutta methods

(19)

Classical fourth-order accurate Runge-Kutta Methods

The classical fourth-order accurate, four-stage method is k1= L(vⁿ, nht)

k2= L

vⁿ+1

2htk1, n +1 2ht

k₃= L

vⁿ+1

2h_tk₂, n +1 2h_t

k4= L (vⁿ+ htk3, (n + 1)ht) vⁿ⁺¹= vⁿ+h_t

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

(20)

Low-storage RK methods

Thes-stage low-storage method is v0= vⁿ

∀ j ∈ [1, . . . , s] :

( kj = ajkj−1+ htL vj, (n + cj)ht

vj = v_j−1+ bjkj

vⁿ⁺¹= vs

where the constants aj, bjand cjare chosen to yield the desired order, s − 1, of the method. For the method to be self-starting we require that a₁ = 0. We need only two storage levels (2N) containing kjand vjto advance the solution.

(21)

Five-stage Fourth-order Low Storage RK Methods

The constants for a five-stage fourth-order Runge-Kutta method (RK4S5) are

a₁ = 0 b₁= 1432997174477

9575080441755 c₁= 0 a₂ = −567301805773

1357537059087 b₂= 5161836677717

13612068292357 c₂= 1432997174477 9575080441755

a₃ = −2404267990393

2016746695238 b₃= 1720146321549

2090206949498 c₃= 2526269341429 6820363962896

a₄ = −3550918686646

2091501179385 b₄= 3134564353537

4481467310338 c₄= 2006345519317 3224310063776

a₅ = −1275806237668

842570457699 b₅= 2277821191437

14882151754819 c₅= 2802321613138 2924317926251

See Carpenter and Kennedy (1994)

(22)

Stability regions for RK4S5

Figure:Stability regions for RK4S5

(23)

Linear m-stage m-order RK Methods

SSPRK linear (m,m) The class of m stage Methods given by:

v⁽ⁱ⁾ = v⁽ⁱ⁻¹⁾+ htLv⁽ⁱ⁻¹⁾, i= 1, . . . , m − 1 v^(m)=

m−2

X

k=0

α_m,kv^(k)+ α_m,m−1

v^(m−1)+ htLv^(m−1) ,

where α_1,0 = 1 and α_m,k= 1

kα_m−1,k−1, k= 1, . . . , m − 2 αm,m−1= 1

m!, αm,0= 1 −

m−1

X

m=1

αm,k

See Gottlieb, Shu and Tadmor (2001) Modified Scheme

(24)

Stability regions for linear SSPRK Methods

Figure:Stability regions for Runge-Kutta methods

Compact Scheme

(25)

RK4S5 v.s. SSPRK

Comparison:

Low Storage: RK4S5:2N; SSPRK:3N Efficiency: ?

(26)

RK4S5 v.s. SSPRK:Stability regions

Figure:Stability regions for Runge-Kutta methods

Compact Scheme

(27)

RK4S5 v.s. SSPRK: Efficiency

3D test problem

N method Error dt CFL T/dt ops Ratio

80 RK4S5 4.00e − 9 1.25e − 3 .3/3 800 4000 SSPRK7 1.48e − 9 3.55e − 3 .85/3 282 1974 2.03 40 RK4S5 3.39e − 7 3.75e − 3 .45/3 267 1335

SSPRK7 2.14e − 7 7.09e − 3 .85/3 141 987 1.35 2D test problem

80 RK4S5 3.33e − 9 1.88e − 3 .3/2 533 2665 SSPRK7 1.34e − 9 5.32e − 3 .85/2 188 1316 2.02 40 RK4S5 1.78e − 7 3.75e − 3 .3/2 267 1335

SSPRK7 1.94e − 7 1.06e − 2 .85/2 94 658 2.03 Table: RK4S5 v.s. SSPRK: (7th order SBP operator)

Note: ops = T/dt · # of stages.

(28)

RK4S5 v.s. SSPRK: Efficiency (Contd.)

3D test problem

80 RK4S5 3.32e − 7 3.745e − 3 .9/3 267 1335 SSPRK5 2.96e − 7 7.194e − 3 1.73/3 139 695 1.92 40 RK4S5 8.20e − 6 8.3333e − 3 1/3 120 600

SSPRK5 9.90e − 6 1.449e − 2 1.73/3 69 345 1.74 2D test problem

80 RK4S5 4.31e − 7 6.25e − 3 1/2 160 800

SSPRK5 2.68e − 7 1.09e − 2 1.73/2 92 460 1.74 40 RK4S5 6.85e − 6 1.25e − 2 1/2 80 400

SSPRK5 8.64e − 6 2.17e − 2 1.73/2 46 230 1.74 Table: RK4S5 v.s. SSPRK: (5th order SBP operator)

Note: ops = T/dt · # of stages.

(29)

Outline

Reference

(30)

Sample Code

Location:

http://www.math.ncku.edu.tw/∼mhchen/HOFD/sample.zip advection_pbc.m: Solve the 1D advection equation with periodic boundary condition.

Diff_cpt4.m: 4th order compact difference operator.

Spectrum_ROS_pbc.m: Plot the ROS and the spectrum.

f1.m-f8.m: Describe the boundary of ROS of RK.

(31)

Example







ut+ ux= 0 x∈ [0, 1] t≥ 0, u(x, 0) = sin(2πx) x∈ [0, 1],

u(0, t) = u(1, t) t≥ 0, which has the exact solution

u(x, t) = sin(2π(x − t)), x∈ [0, 1], t≥ 0.

(32)

Compact Fourth-Order Scheme

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

1

4v_xj−1+ vxj+1

4v_xj+1= 3 2

v_j+1− v_j−1 2h

where

v_j ≈ u(x_j) vx = vxj≈ du dx(xj).

(33)

Eigenvalue Spectrum of the Compact Scheme

Figure:Eigenvalue Spectrum:(Left) N=40; (Right) N=80

ROS_RK ROS_RK4S5

(34)

Figure:(left) RK4 with CFL= 1.5; (right) RK4S5 with CFL= 1.8.

(35)

Summary

Steps to solve a PDE

1 Approximate the spatial differential operators. (FDM, FVM, FEM, Spectral Methods)

2 Study the stability of the resulting semi-discrete scheme

3 Study the eigenvalue spectrum of the resulting ODE system

4 Choose a time integrating scheme (RK, LMM, etc) and determine the CFL number

5 Solve the ODE system by an appropriate time integration scheme

(36)

Outline

Reference

(37)

Reference

J. S. Hesthaven, S. Gottlieb and D. Gottlieb, “Spectral Methods for Time-Dependent Problems.”

Randall LeVeque, “Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems.”

(38)

Outline

Reference

(39)

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, and Φ contains the time explicit boundary conditions.

(40)

Question

How do we implement a RK method to solve the semi-discrete scheme?

Classical RK4

k₁= L(vⁿ, nht) =Lvⁿ+ Φ(tn) k₂= L

vⁿ+1

2h_tk₁, n + 1 2h_t

=L(vⁿ+1

2h_tk₁) + Φ(t_n+1/2) k₃= L

vⁿ+1

2htk₂, n + 1 2ht

=L(vⁿ+1

2htk₂) + Φ(t_n+1/2) k₄= L (vⁿ+ htk₃, (n + 1)ht) =L(vⁿ+ htk₃) + Φ(tn+1) vⁿ⁺¹= vⁿ+ht

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

(41)

Question

How do we implement a RK method to solve the semi-discrete scheme?

Classical RK4

k₁= L(vⁿ, nht) =Lvⁿ+ Φ(tn) k₂= L

vⁿ+1

2h_tk₁, n + 1 2h_t

=L(vⁿ+1

2h_tk₁) + Φ(t_n+1/2) k₃= L

vⁿ+1

2htk₂, n + 1 2ht

=L(vⁿ+1

2htk₂) + Φ(t_n+1/2) k₄= L (vⁿ+ htk₃, (n + 1)ht) =L(vⁿ+ htk₃) + Φ(tn+1) vⁿ⁺¹= vⁿ+ht

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

(42)

Implementation of Classical RK4-I

Run Wave_RK4C_1D.m

N Error Ratio

20 1.544e − 2 - 40 2.013e − 3 7.67 80 2.5428e − 4 7.92 Table: Classical RK4

(43)

Implementation of Classical RK4-II

ODE system: d

dtu = Du + φ(t)e0, D = Q/hx

Numerical solution:

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

2htk₁) + g₁e0

k₃= D(vⁿ+1

2h_tk₂) + g₂e0

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

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

(44)

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) e0

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²)g₁e0

+ (1 3ht+1

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

6ht)g₃e₀

(45)

Implementation of Classical RK4-IV

D³: 1

24D³h⁴_tg0= 1

24D³h⁴_tφ(tn)

⇒ g0= φ(tn) D²: 1

12D²h³_tg0+ 1

12D²h³_tg1= 1

6D²h³_tφ(t_n) + 1

24D²h⁴_tφ⁰(t_n)

⇒ g1= φ(tn) +1 2htφ⁰(tn) D¹: . . .

⇒ g₂= φ(t_n) + φ(t_n)h_t

2 +φ⁰⁰(t_n) 2! (h_t

2)² D⁰: . . .

⇒ g₃= φ(t_n) + φ⁰(t_n)h_t+φ⁰⁰(t_n)

2! h²_t +φ⁽³⁾(t_n) 3! h³_t

(46)

Implementation of Classical RK4-V

Run Wave_RK4CM_1D.m

N Error Ratio Order

20 4.677e − 4 - -

40 3.111e − 5 15.034 3.91 80 1.977e − 6 15.734 3.98 80 1.242e − 7 15.922 3.99 Table: Modified Classical RK4

(47)

Modified SSP-RK Methods







v⁽⁰⁾ = vⁿ,

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

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

where the ht denotes the time step and the coefficients α_m,k, 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. and

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

(48)

Outline

Reference

(49)

Sample Code-SBP 1D

Location:

http://www.math.ncku.edu.tw/∼mhchen/HOFD/wave.zip Wave_SSPRK_1D.m: Solve the 1D problem with SSPRK.

Wave_RK4S5_1D.m: Solve the 1D problem with RK4S5.

Wave_SSPRK_2D.m: Solve the 2D problem with SSPRK.

Mx_2D.m, My_2D.m:

Diff_CompFD77451_e1.m: 5th order SBP operator.

coefb.m and coefbin.m: Coefficients for SSPRK.

Spectrum_ROS.m: Plot the ROS and the spectrum.

f1.m-f8.m: Describe the boundary of ROS of RK.

(50)

1D Example

Example







u_t+ ux= 0 x∈ [0, 1] t≥ 0, u(x, 0) = sin(ωx) x∈ [0, 1],

u(0, t) = sin(−ωt) t≥ 0, which has the exact solution

u(x, t) = sin(ω(x − t)), x∈ [0, 1], t≥ 0.

Numerical Scheme: (SBP operators) dv

dt + h⁻¹P⁻¹Qv = h⁻¹τ q00(v0− g(t))P⁻¹e0,

Refer to Inhomogeneous ODE System

(51)

Eigenvalue Spectrum

Numerical Scheme: (SBP operators) dv

dt + h⁻¹P⁻¹Qv = h⁻¹τ q₀₀(v₀− g(t))P⁻¹e₀,

To investigate the stability of the one dimensional scheme, we need to find the eigenvalue spectra of the operators:

L = −h⁻¹P⁻¹Q + h⁻¹τ q₀₀P⁻¹E₀

= h⁻¹P⁻¹(−Q + τ q₀₀E₀) where

E₀ =diag(1, 0, 0, ..., 0) ∈ M_L+1.

(52)

Figure:Eigenvalue spectra of the fifth order accurate operator

−h⁻¹P⁻¹Q + h⁻¹τ q00P⁻¹E0, for various values of N and τ .

(53)

Eigenvalue Spectrum with ROS

Figure:Stability envelope of fifth order accurate SSP-RK and the spectra of the 5th order scheme with CFL = 1.73. Computational parameters: N = 40 and τ = 2

(54)

Figure:L₂-error history of the fifth order (left) and the seventh order (right) 1D schemes with CFL = 1.73 .

(55)

2D Example

Example

Let us consider the two dimensional advection problem











ut+ ux+ uy = 0 x, y ∈I²= [0, 1]² t≥ 0 u(x, y, 0) = sin(ω(x + y)) x, y ∈I²

u(0, y, t) = sin(ω(y − 2t)) y∈ [0, 1] t≥ 0 u(x, 0, t) = sin(ω(x − 2t)) x∈ [0, 1] t≥ 0 which has the exact solution of the form

u(x, y, t) = sin(ω(x + y − 2t)), x, y ∈ [0, 1], t≥ 0.

(56)

Numerical Scheme

dv

dt + Dxv + Dyv = r, where

D_x = I_M+1⊗ h⁻¹_x P⁻¹_x Q_x, D_x = h⁻¹_y P⁻¹_y Q_y⊗ I_L+1,

r = h⁻¹_x τ₁

M

X

j=0

q^x₀₀(v_0j− g₁(yj, t))(eyj⊗ P⁻¹_x e_x0)

+ h⁻¹_y τ₂

L

X

i=0

q^y₀₀(v_i0− g₂(xi, t))(P⁻¹_y e_y0⊗ exi).

(57)

Eigenvalue Spectrum-2D

For the two dimensional scheme, the corresponding operator L takes the form

L = IM+1⊗ h⁻¹_x P⁻¹_x (−Q_x+ τ₁q^x₀₀E_x0) + h⁻¹_y P⁻¹_y (−Q_y+ τ₂q^y₀₀E_y0) ⊗ I_L+1

where E_x0∈ M_L+1and E_y0∈ M_M+1are explicitly given by Ex0=diag(1, 0, 0, ..., 0), Ey0=diag(1, 0, 0, ..., 0).

(58)

Eigenvalue Spectrum-2D (contd.)

Let λ^x_i, i = 0, 1, ..., L and λ^y_j, j = 0, 1, ..., M denote the eigenvalues of the one dimensional operators

h⁻¹_x P⁻¹_x (−Q_x+ τ₁q^x₀₀E_x0) and h⁻¹_y P⁻¹_y (−Q_y+ τ₂q^y₀₀E_y0), respectively. Then, denoted by λijthe eigenvalues of L are

λij= λ^x_i + λ^y_j,

followed by Lemma. Immediately, this gives a rough estimate on the largest eigenvalue in magnitude as

max |λij| ≤ 2 max(|λ^x_i|, |λ^y_j|).

(59)

Figure:L2-error history of the fifth order 2D schemes. Computational time steps are computed with CFL = 1.73/2 .

(60)

Outline

Reference

(61)

Reference

M. H. Carpenter, D. Gottlieb, S. Abarbanel and W.S. Don,

“The theoretical accuracy of Runge-Kutta time

discretizations for the initial boundary value problem: a study of the boundary error.” SIAM Journal on Scientific Computing, 1995.

M-H Chen, B. Cockburn and F. Reitich, “High-order RKDG Methods for Computational Electromagnetics.” Journal of Scientific Computing, 2005.