ERROR BOUNDED SCHEMES FOR TIME-DEPENDENT HYPERBOLIC PROBLEMS∗
JAN NORDSTR ¨OM†
Abstract. In this paper we address the error growth in time for hyperbolic problems in first
order form in bounded domains. The energy method is used to study when an error growth or a fixed error bound is obtained. It is shown that the choice of boundary procedure is a crucial point. Numerical experiments corroborate the theoretical findings.
Key words. temporal error, error bound, hyperbolic problem, boundary conditions, energy method, well posed problem, stability, weak formulation, penalty technique
AMS subject classifications. 35M40, 35M45, 35M50, 65M06, 65M12, 65M15 DOI. 10.1137/060654943
1. Introduction. Stable approximations of hyperbolic problems in first order
form often exhibit a linear (or near linear) error growth in time; see, for example, [16], [14], [7]. In many of those cases, the wave is restricted or trapped in the domain for long times. Typical examples include periodic problems or problems where the wave is trapped in cavities. However, cases with a definite bound on the error as time passes can also be observed. Typically, in such cases, the wave propagates through the domain for a limited time as in an inflow-outflow problem. This means that (i) the calculation is performed on a bounded domain and that (ii) boundary conditions are imposed. These are the two prerequisites in our investigation and we will focus on the latter.
To initiate our investigation we consider the problem
(1.1) ut+ ux= F (x, t), −∞ < x < ∞, t ≥ 0,
where F is a 1-periodic forcing function in x. The 1-periodic solution which we consider known in this example provides us with initial (f (x)) and boundary (g(t)) data. Assume that we want to compute the numerical solution to (1.1) on the domain x∈ [0, 1]. There are two different ways of imposing the periodic boundary conditions. The first one does it directly by demanding u(0, t) = u(1, t) (periodic boundary con-dition). A second possible choice reads u(0, t) = g(t) (Dirichlet boundary concon-dition). We compute (details of the numerical procedures are discussed later) the two solutions and subtract the exact solution. The error as a function of time for the two types of boundary conditions is shown in Figure 1.1. The computation using the first type leads to an essentially linear error growth in time, while the second type leads to an error bound.
A hint as to why the results are different can be obtained by studying (1.1) on x∈ [0, 1] with a perturbed forcing function F + δF and identical initial and boundary ∗Received by the editors March 23, 2006; accepted for publication (in revised form) April 10, 2007; published electronically November 7, 2007.
http://www.siam.org/journals/sisc/30-1/65494.html
†Department of Computational Physics, FOI, The Swedish Defense Research Agency, SE-164 90 Stockholm, Sweden, and Department of Information Technology, Scientific Computing, Uppsala University, SE-751 05 Uppsala, Sweden, and Department of Aeronautical and Vehicle Engineering, KTH, The Royal Institute of Technology, SE-100 44 Stockholm, Sweden (Jan.Nordstrom@foi.se).
0 2 4 6 8 10 −4.5 −4 −3.5 −3 −2.5 −2 time log10(Err) Dirichlet BC Periodic BC
Fig. 1.1. The error as a function of time.
conditions. Denote the solution to the perturbed problem by v and subtract the equations. That yields the error equation
(1.2) et+ ex= δF (x, t), 0≤ x ≤ 1, t ≥ 0,
where e = v− u and e(x, 0) = 0. The two corresponding boundary conditions are e(0, t) = e(1, t) and e(0, t) = 0, respectively.
The energy method applied to (1.2) yields
(1.3) ||e||t≤ −λ(t)||e|| + ||δF ||,
where||e||2=01e2dx. For the boundary condition e(0, t) = e(1, t) we have λ(t) = 0 which leads to the error growth shown in Figure 1.1. More interesting is the fact that we get λ(t) = (e(1, t)/||e||)2/2 by using the boundary condition e(0, t) = 0. A constant
positive λ leads directly to error boundedness. Actually, (as we will show) the weaker condition that0tλ(η)dη is a monotonically growing function of time suffices for error boundedness; see Figure 1.1.
The error growth in time is particularly problematic for hyperbolic problems. For parabolic problems (add uxxto the right-hand side of (1.1)), a temporal error bound
is obtained by using Poincar´e’s inequality (||ex||/||e|| ≥ δ > 0) even with vanishing
boundary terms. However, for vanishing , as the parabolic problem converges to a hyperbolic problem, the error growth appears since the coefficient λ = δ in (1.3) goes to zero.
Some of the phenomena above were discussed in [1], [3], [4]. In this paper we will both extend and simplify that analysis and focus on the error behavior in time for hyperbolic problems. Our theoretical tools will be the energy method, summation-by-parts operators, and penalty procedures for implementing boundary conditions. A related problem was discussed in [17]. The error bound was obtained by making use of the presence of so-called lacunae (pockets in space-time) where the solutions to the wave equation in odd space dimensions are identically zero.
The remainder of this paper will proceed as follows. In section 2, we consider different boundary procedures for a general continuous hyperbolic problem and derive an error estimate. The corresponding semidiscrete finite-difference approximation and
error estimate are considered in section 3. Numerical experiments are performed in section 4, and finally we draw conclusions in section 5.
2. The continuous problem. We will consider a general three-dimensional
hyperbolic problem in a half space of the form
Sut+ Aux+ Buy+ Cuz= ˜F , x≥ 0, t ≥ 0, Lu = ˜G, x = 0, t≥ 0, (2.1)
u = ˜H, x≥ 0, t = 0.
The subscripts t, x, y, z denote partial differentiation with respect to time (t) and space (x, y, z), respectively. The solution u has m components and we assume that all m×m matrices (S, A, B, C) are symmetric (or simultaneously symmetrizable; see [2]) and that S is a diagonal strictly positive matrix. In terms of applications this means that we consider aerodynamic, aeroacoustic (the Euler equations), and electromagnetic (Maxwell’s equations) problems. We also assume that u, ˜F , ˜H have compact support in x≥ 0.
Remark. Note that although the problem (2.1) is formulated on a half space, we consider bounded domains in this paper. The assumption of compact support simply enables us to treat one boundary at a time. In applications (see section 4 below) only bounded domains exist, and the whole boundary must be treated as we treat the boundary at x = 0.
Definition 2.1. Equation (2.1) is said to be strongly well posed if a unique solution exists and the estimate
(2.2) u2+ t 0 u2 Γdτ ≤ Kceηct ˜H2+ t 0 ( ˜F2+ ˜G2Γ)dτ
holds. Kc and ηc are constants and do not depend on ˜F , ˜G, or ˜H. · Γ and · are
suitable continuous norms.
Remark. If one can estimate the solution for zero boundary data only, the problem is called well posed; see [8] for more details.
2.1. Well posedness. The energy method applied to (2.1) yields
(2.3) ((u, u)S)t=
+∞
−∞
+∞
−∞
(uTAu)x=0 dydz + 2(u, ˜F )I,
where I is the identity matrix and we have introduced the scalar product
(2.4) (u, v)T = +∞ −∞ +∞ −∞ +∞ 0 uTT v dxdydz.
T is a symmetric positive definite matrix and (u, u)T = ||u||2T defines a norm. In
Definition 2.1 above, ||u||2 I =||u|| 2 and +∞ −∞ +∞ −∞ (uTu)x=0dydz =u2Γ were used.
We have the splitting
where A+= XTΛ+X and A− = XTΛ−X. The matrix X has the normalized
eigen-vectors of A as columns. The superscripts +,− indicate the positive and negative portions of the eigenvalue matrix Λ and A. Let ˜g denote the value of the solution u at the boundary x = 0. The (characteristic) boundary condition
(2.6) Lu = (XTu)+= ˜G = (XTg)˜ + or A+u = A+˜g yields the estimate
(2.7) ((u, u)S)t= +∞ −∞ +∞ −∞ (˜gTA+˜g + uTA−u)x=0dydz + 2(u, ˜F )I.
A direct integration of (2.7) yields an estimate of the form (2.2), which proves the following proposition.
Proposition 2.2. Equation (2.1) augmented with (2.6) is strongly well posed.
2.2. The continuous error equation. Consider (2.1) with a perturbed forcing
function ˜F + ˜δF . Denote the solution to the perturbed problem by v and subtract the equations. That yields the error equation
Set+ Aex+ Bey+ Cez= ˜δF , x≥ 0, t ≥ 0, Le = 0, x = 0, t≥ 0, (2.8)
e = 0, x≥ 0, t = 0. The energy method applied to (2.8) and the relations
(2.9) ((e, e)S)t= 2||e||S(||e||S)t, (e, ˜δF )I ≤ ||e||I|| ˜δF||I
yields (2.10) (||e||S)t≤ C0||e||S+ C1|| ˜δF||I, where (2.11) C0= +∞ −∞ +∞ −∞(eTAe)x=0dydz 2||e||2 S , C1= ||e||I ||e||S .
Remark. We call the boundary procedure dissipative if, by using the boundary conditions, we can make sure that C0 ≤ 0. Dissipative boundary treatments, for
example, using (2.6), lead to strong stability. The problem (2.10) is of the form
(2.12) ut≤ −λ(t)u + F (t), u(0) = f, λ ≥ 0,
where F has a maximum value|F |max.
We will need the following lemma.
Lemma 2.3. The solution u to (2.12) is bounded for t≥ 0 provided that λ is such that the estimate
(2.13) θ(ξ, t) =
t ξ
λ(η)dη≥ δ0(t− ξ), δ0> 0,
()
(0)
t
t
Fig. 2.1. The nonnegative λ and the corresponding integral θ(0, t).
Proof. A direct calculation yields
(2.14) u(t)≤ e−θ(0,t)f +
t
0
e−θ(ξ,t)F (ξ)dξ.
The estimate (2.13) and (2.14) lead to
(2.15) |u| ≤ e−δ0tf +|F |max(1− e−δ0t)/δ0,
which means that u remains bounded as time increases.
Remark. Note that Lemma 2.3 and the estimate (2.13) hold even if λ is zero oc-casionally; see Figure 2.1. Clearly θ(ξ, t) is a positive monotonically growing function of time if λ does not vanish completely. Roughly speaking, λ consists of the bound-ary error divided by the total error in the domain (see (2.11) with λ = −C0). The
assumption that λ is such that (2.13) holds is reasonable, since the solution to (2.1) and (2.8) consists of waves that travel in all directions and frequently (if not always) exist on the boundary.
We can now state one of the main results of this paper.
Proposition 2.4. If the boundary procedure in (2.1) is dissipative, i.e., C0≤ 0 in (2.10), and the estimate (2.13) holds with λ =−C0, then the error is bounded in
time. It is bounded by the size of the forcing function. Proof. Lemma 2.3 applied to (2.10) leads to the result.
Remark. We include a nonhomogeneous initial condition in (2.12) even though (2.10) has zero initial data. The reason is that C0can be zero initially since e(0) = 0.
If so, the error grows linearly for a short time and forms the nonzero initial error. Remark. The error bound in Proposition 2.4 is obtained if the calculation is on a bounded domain and the boundary procedure is dissipative.
A direct result of Propositions 2.2 and 2.4 is the following corollary. Corollary 2.5. Strong well posedness leads to a temporal error bound.
3. The semidiscrete problem. Let u, uξ be the numerical approximation of
the scalar quantities u, uξ, respectively. The approximation of the first derivative uξ
is introduced as
(3.1) uξ = P−1ξ Qξu,
where Pξ and Qξ are matrices. If a spatial operator is of the form (3.1) and the
conditions (i) and (ii) below are fullfilled, the operator is referred to as a summation by parts (SBP) operator (see [12]).
(i) Pξ is symmetric, positive definite, and bounded. ΔξpI ≤ Pξ ≤ ΔξqI, where p > 0 and q are bounded independent of 1/Δξ.
(ii) Qξ is almost skew-symmetric. Qξ+ QTξ = diag(−1, 0, . . . , 0, 1).
Details of the structure of the matrices Pξ and Qξ are given in [6]. The nth order
interior accuracy is complemented by n− 1 order accuracy at the boundaries. In this paper we use n = 2, 4, 6.
The boundary and interface conditions will be imposed by the weak penalty tech-nique referred to as the simultaneous approximation term (SAT) procedure. For de-tails on penalty procedures for finite difference, finite volume, and spectral methods, see [5], [15], and [9], respectively.
The unknowns are organized in a vector u where ui,j,k corresponds to the grid
point (xi, yj, zk). The semidiscrete version of problem (2.1) on the finite-difference
form and using the SBP operators above can be written Ix⊗ Iy⊗ Iz⊗ S ut+ P−1x Qx⊗ Iy⊗ Iz⊗ A u + Ix⊗ P−1y Qy⊗ Iz⊗ B u + Ix⊗ Iy⊗ P−1z Qz⊗ C u (3.2) = P−1x Bx⊗ Iy⊗ Iz⊗ Σ (u− ˜g) + ˜F , u(0) = ˜H.
The symbol ⊗ denotes the Kronecker product [13]. The matrices P−1x Qx, P−1y Qy,
P−1z Qz are one-dimensional SBP difference operators in the x, y, z directions, respec-tively. The data ˜g, ˜G, ˜F , and ˜H of the continuous problem are stored in the discrete vectors ˜g, ˜G, ˜F , and ˜H, respectively. Ix, Iy, Iz are identity matrices of appropriate
size. Bxis identically zero except that (Bx)11= 1 and Σ will be determined later.
The first term on the right-hand side of (3.2) imposes the boundary conditions at x = 0 weakly using the SAT technique. The characteristic boundary condition (2.6) can be imposed in this way by a special choice of Σ.
Remark. Also, in the semidiscrete problem (3.2) we simplify the analysis and treat one boundary only. See the related remark for the continuous problem in section 2.
Definition 3.1. Equation (3.2) is said to be strongly stable if, for a sufficiently fine mesh Δx, Δy, Δz < h, there is a unique solution that satisfies
(3.3) u2+ t 0 u2 Γdτ ≤ Kdeηdt ˜H2+ t 0 ( ˜F2+ ˜G2Γ)dτ .
Kdand ηdare constants and do not depend on ˜F , ˜H, or ˜G. ·Γand· are suitable
discrete norms.
Remark. If one can estimate the solution for zero boundary data only, the problem is called stable; see [8] for more details.
3.1. Stability. By applying the energy method to (3.2) (we multiply from the
left by uT(P
x⊗ Py⊗ Pz⊗ I) and add the transpose) we obtain
((u, u)P,S)t= uT0 Py⊗ Pz⊗ A + Σ + ΣT u0− 2uT0 Py⊗ Pz⊗ Σ ˜ g + 2(u, ˜F )P,I, (3.4)
where u0= (Bx⊗Iy⊗Iz⊗I)u denotes the restriction of u to the boundary at x = 0
and (3.5) (u, v)P,T = uT Px⊗ Py⊗ Pz⊗ T v
is the discrete analogue to the continuous scalar product (2.4). T is a symmetric positive definite matrix and (u, u)P,T = ||u||2P,T defines a norm. In Definition 3.1
above, ||u||2 P,I=||u|| 2 and uT 0 Py⊗ Pz⊗ I u0=u2Γ were used.
Clearly, we have a bounded growth if φT(A + Σ + ΣT)φ = 0 for an arbitrary
nonzero vector φ and ˜g = 0. We have proved the following proposition. Proposition 3.2. The approximation (3.2) with Σ chosen such that
φT(A + Σ + ΣT)φ = 0 for an arbitrary nonzero vector φ is stable.
A particularly useful and strong concept of stability is obtained in the following way. By introducing the previous notation A = A++ A− we find that the boundary terms in (3.4) can be written
BT = uT0 Py⊗ Pz⊗ A++ Σ + ΣT u0 + uT0 Py⊗ Pz⊗ A− uo− 2uT0 Py⊗ Pz⊗ Σ ˜ g. (3.6)
The choice Σ = σA+means that A++Σ+ΣT = (1+2σ)A+and σ≤ −1/2 means that
the quadratic terms in u0 give a negative contribution. The special choice σ =−1,
(3.4), and (3.6) yield the estimate ((u, u)P,S)t= ˜gT Py⊗ Pz⊗ A+ ˜ g + uT0 Py⊗ Pz⊗ A− u0
+ 2(u, ˜F )P,I− Rest,
(3.7) where (3.8) Rest = u0− ˜g T Py⊗ Pz⊗ A+ u0− ˜g .
The discrete estimate (3.7) is completely analogous to the continuous estimate (2.7) except for the small negative definite rest term−Rest.
A direct integration of (3.7) yields an estimate of the form (3.3) which proves the following proposition.
Proposition 3.3. The approximation (3.2) with Σ = σA+ and σ ≤ −1/2 is
Remark. Note that the energy method combined with the penalty technique more or less automatically leads to the right number of boundary conditions (equal to the number of positive eigenvalues) and the most appropriate form (the character-istic boundary conditions (2.6)). Note also that the rest term is proportional to the accuracy at the boundary since
(3.9) Rest = (XTu)0− ˜G T Py⊗ Pz⊗ Λ+ (XTu)0− ˜G .
3.2. The semidiscrete error equation. By inserting the continuous solution
u (and using the correct boundary data, forcing function, and initial data) into the scheme (3.2) and subtracting the discrete solution we obtain the error equation
Ix⊗ Iy⊗ Iz⊗ S et+ P−1x Qx⊗ Iy⊗ Iz⊗ A e + Ix⊗ P−1y Qy⊗ Iz⊗ B e + Ix⊗ Iy⊗ P−1z Qz⊗ C e (3.10) = P−1x Bx⊗ Iy⊗ Iz⊗ Σ e + T e, e(0) = 0.
Clearly the only source of error is the truncation error T e.
By following the recipe outlined above (the energy method) we obtain (3.11) ((e, e)P,S)t= eT0
Py⊗ Pz⊗ A + Σ + ΣT
e0+ 2(u, T e)P,I,
where e0 denotes the restriction of e to the boundary at x = 0, i = 0.
Remark. The error equation (3.11) completely determines the accuracy of the solution for long times. With a given truncation error (given by the scheme) we obtain a more or less favorable error development depending on how we choose our boundary conditions (choice of Σ).
The relations
(3.12) ((e, e)P,S)t= 2||e||P,S(||e||P,S)t, (e, T e)P,I≤ ||e||P,I||T e||P,I
inserted into (3.11) yield
(3.13) (||e||P,S)t≤ C0||e||P,S+ C1||T e||P,I,
where (3.14) C0= eT 0 Py⊗ Pz⊗ (A + Σ + ΣT) e0 2||e||2 P,S , C1= ||e|| P,I ||e||P,S .
Remark. If C0 vanishes, we have the usual (essentially) linear error growth in
time for hyperbolic problems (see, for example, [16], [14], [7]). We can now state the other main result of this paper.
Proposition 3.4. If the boundary procedure in (3.2) is dissipative, i.e., C0≤ 0 in (3.13), and the estimate (2.13) holds with λ =−C0, then the error is bounded in
time. It is bounded by the size of the truncation error. Proof. Lemma 2.3 applied to (3.13) leads to the result.
The assumption that (2.13) holds is equally reasonable in the discrete case, and the general requirements for the error bound are the same; see the two remarks just before Proposition 2.4.
Propositions 3.3 and 3.4 immediately lead to the following corollary. Corollary 3.5. Strong stability leads to a temporal error bound.
0 50 100 150 200 −4.5 −4 −3.5 −3 −2.5 −2 −1.5 −1 time log10(Err) Dirichlet BC Periodic BC
Fig. 4.1. The error for a very long time.
4. Examples and numerical experiments. By performing numerical
calcu-lations we will check the validity of Proposition 3.4.
4.1. The one-way wave equation. The numerical approximation of (1.1) on
the domain [0, 1] including the penalty terms for the boundary condition has the general form given in (3.2). The one-dimensional approximation is given by inserting Py = Pz = Iy= Iz= 1, B = C = R = 0, and A = 1. We get
(4.1) ut+ P−1Qu = P−1σ(u0− g(t))a0+ F ,
where a0 = (1, 0, . . . , 0, 0)T. We have used the stable choice σ =−1/2. Recall that F is a 1-periodic forcing function. The exact 1-periodic solution which we consider in this example provides us with the boundary data (g(t)).
The energy method on the corresponding error equation directly yields
(4.2) (||e||P)t≤ −λ(t)||e||P+||T e||P,
where ||e||2 P = e TP e and λ = (−(1 + 2σ)e2 0+ e2N)/2||e|| 2 P. By the choice 1 + 2σ ≤ 0
we get an error bound in time; see Proposition 3.4. To compare with the general formulation (see the description below (3.6)), note that A+ = 1, Σ = σ, and A++
Σ + ΣT = 1 + 2σ. The calculation using the scheme (4.1) and the calculation where
the periodicity is imposed directly,
(4.3) ut+ P−1Qu = P−1σL(u0− uN)a0+ P−1σR(uN − u0)aN + F ,
are shown in Figures 1.1 and 4.1. In (4.3) we have used the stable choices σL = −1/2, σR= 1/2, and aN = (0, 0, . . . , 0, 1)T. The calculation (4.1) using the Dirichlet
boundary condition and (4.3) using the periodic boundary condition are initiated with the exact same periodic function. The amplitude of the error is proportional to Δx4
since we use a fourth order method.
Remark. The result in Figure 4.1 is truly remarkable. By mesh refinement, we can obtain an arbitrarily high accuracy at any future time.
4.2. Maxwell’s equations. The distribution of electromagnetic fields are
de-scribed by Maxwell’s equations
(4.4) μ∂H
∂t =−∇ × E, ∂E
∂t =∇ × H − J,
(4.5) ∇ · E = ρ, ∇ · μH = 0,
combined with the equation of continuity ρt+∇ · J = 0. Here E is the electric field, H is the magnetic field, J is the electric current density, and ρ is the charge density. and μ are permittivity and permeability, respectively.
With J = 0 we can formulate the two-dimensional version of (4.4) in transverse electric (TE) mode as
(4.6) Sut+ Aux+ Buy= 0,
where we consider the domain (x, y)∈ [−1, 1] × [0, 1] and
A = ⎡ ⎣ 00 00 10 1 0 0 ⎤ ⎦ , B = − ⎡ ⎣ 01 10 00 0 0 0 ⎤ ⎦ , S = ⎡ ⎣ μ0 0 00 0 0 ⎤ ⎦ , u = ⎡ ⎣ HExz Ey ⎤ ⎦ .
4.2.1. The continuous problem. Using the energy method on (4.6) and using
that the matrices A and B are symmetric lead to
(4.7) d dtu 2 S =− 1 0 uTAu|x=1x=−1dy− 1 −1 uTBu|y=1y=0dx,
whereu2S = uTSu. The problem (4.6) augmented with the initial condition is well posed if we specify the Ey component at x =−1 and x = 1 and the Ex component
at y = 0 and y = 1.
Remark. In terms of an error equation where the boundary data and forcing function are zero, this means that the time derivative of the norm of the error is zero. With a nonzero forcing function included, we would get the classical linear error growth.
The problem (4.6) augmented with the initial condition is strongly well posed if we specify the ingoing characteristic variable at each boundary. This yields the estimate d dtu 2 S = 1 0 (uTA−u + gTA+g)x=−1dy− 1 0 (uTA+u + gTA−g)x=1dy + 1 −1 (uTB−u + gTB+g)y=0dx− 1 −1 (uTB+u + gTB−g)y=1dx, where A+,−= X AΛ+,A−XAT, B+,− = XBΛ+,B−XBT, and XA= ⎡ ⎣ −1/ √ 2 0 1/√2 0 1 0 1/√2 0 1/√2 ⎤ ⎦ , ΛA= ⎡ ⎣ −1 0 00 0 0 0 0 1 ⎤ ⎦ , XB= ⎡ ⎣ −1/ √ 2 1/√2 0 1/√2 1/√2 0 0 0 1 ⎤ ⎦ , ΛB = ⎡ ⎣ 10 −1 00 0 0 0 0 ⎤ ⎦ .
-1 0 1 0 1 -1 0 1 (a) -1 0 1 0 1 -1 0 1 (b) -1 0 1 0 1 -1 0 1 (c)
Fig. 4.2. (a) The wave propagating to the left. (b) The wave propagating to the right. (c) The
total wave, Θi=π6.
Remark. In terms of an error equation where the boundary data and forcing function are zero, this means the time derivative of the norm of the error is negative semidefinite. With a nonzero forcing function included, we will get an error bound; see Proposition 2.4 above.
4.2.2. Calculations. We will compute the wave propagation over a material
discontinuity. All details (especially concerning the treatment of the discontinuity itself) for this calculation are given in [16]. Here we will focus on the effect of the outer boundary conditions. The numerical approximation of (4.6) including the penalty terms for boundary conditions has the general form in (3.2). The more restricted two-dimensional approximation is obtained with Pz= Iz= 1 and C = 0 and becomes
Ix⊗ Iy⊗ S ut+ P−1x Qx⊗ Iy⊗ A u + Ix⊗ P−1y Qy⊗ B u = P−1x Bx⊗ Iy⊗ Σ (u− ˜g) + BT + IT . (4.8)
For clarity we have included only the explicit form of the outer boundary condition at x = −1. The other outer boundary conditions at x = 1, y = 0, and y = 1 are collected in the term BT . The interface conditions at x = 0 for the treatment of the discontinuity are included in the term IT .
The different types of boundary conditions are obtained by varying the form of the penalty matrix Σ. For the time integration we use the classical fourth order
0 50 100 150 200 250 10−4 10−3 10−2 10−1 100 t ||error E ||L2 n = 2 n = 4 n = 6
Fig. 4.3. The error using noncharacteristic boundary conditions.
0 50 100 150 200 250 10−6 10−5 10−4 10−3 10−2 10−1 t ||error E ||L2 n = 2 n = 4 n = 6
Fig. 4.4. The error using characteristic boundary conditions.
Runge–Kutta method with a very small time-step. We have an exact solution which is a plane wave that propagates over the discontinuity. That solution provides us with boundary and initial data, and we impose the two different sets of boundary conditions mentioned above.
By using the tangential boundary condition at the boundaries, i.e., the Ey
com-ponent at x =−1 and x = 1 and the Ex component at y = 0 and y = 1, we obtain
a stable approximation according to Proposition 3.2. For example, at the boundary x =−1 that means that we have
Σ = ⎡ ⎣ 00 00 −10 0 0 0 ⎤ ⎦ , and consequently A + Σ + ΣT = ⎡ ⎣ 00 00 00 0 0 0 ⎤ ⎦ , which yields C0= 0. An error bounded scheme cannot be expected; see Proposition
By specifying the ingoing characteristic variables at the boundaries, i.e., the (Ey+ Hz)/ √ 2 at x =−1, (Ey−Hz)/ √ 2 at x = 1, (Ex−Hz)/ √ 2 at y = 0, and (Ex+Hz)/ √ 2 at y = 1, we obtain a strongly stable approximation according to Proposition 3.3. For example, at the boundary x =−1 we have (by choosing σ = −1)
Σ =−A+=−1 2 ⎡ ⎣ 10 00 10 1 0 1 ⎤ ⎦ and A + Σ + ΣT = ⎡ ⎣ −1 00 0 00 0 0 −1 ⎤ ⎦ , which yields C0< 0. An error bounded scheme can be expected; see Proposition 3.4.
As an example we show what happens when the wave hits the discontinuity with an angle of incidence Θi = π6. A significant reflection at the interface between the
two mediums occurs; see Figure 4.2. In [16] it was shown that the accuracy in space is of design order.
Next we study the error development in time. If we specify the tangential electrical fields at the boundaries and measure the error as a function of time, we obtain the classical linear error growth shown in Figure 4.3. By specifying the characteristic variables we obtain the error bounded result shown in Figure 4.4. The results in Figures 4.3 and 4.4 are computed with 2, 4, and 6 orders of accuracy. The error levels using the characteristic boundary conditions are significantly lower than those obtained for noncharacteristic boundary conditions.
Remark. The result in Figure 4.4 is remarkable. By mesh refinement, we can obtain an arbitrarily high accuracy at any future time.
Remark. The error boundedness in this paper is due to the boundary treatment only; the internal discretization is not important. The same result was obtained using the spectral element code USEMe; see [10], [11].
5. Conclusions. We have shown that an error bound can be obtained for a
general time-dependent hyperbolic problem if a sufficiently dissipative (or sufficiently stable) boundary procedure can be constructed on a bounded domain.
Two examples, the one-dimensional one-way wave equation and the two-dimen-sional Maxwell’s equations in transverse mode, were studied and numerical calcula-tions were performed.
The results of the calculations using the error bounded schemes are remarkable. With accurate boundary data and mesh refinement, we can obtain an arbitrarily high accuracy at any future time.
The error boundedness in this paper is due to the boundary treatment only; the internal discretization (finite differences were used) is not important. The same result was obtained using a spectral element code.
REFERENCES
[1] S. Abarbanel, A. Ditkowski, and B. Gustafsson, On error bounds of finite difference
ap-proximations to partial differential equations—temporal behavior and rate of convergence,
J. Sci. Comput., 15 (2000), pp. 79–116.
[2] S. Abarbanel and D. Gottlieb, Optimal time splitting for two and three dimensional Navier–
Stokes equations with mixed derrivatives, J. Comput. Phys., 41 (1981), pp. 1–33.
[3] S. Abrabanel and A. Chertock, Strict stability of high-order compact implicit
finite-difference schemes: The role of boundary conditions for hyperbolic PDEs, I, J. Comput.
Phys., 160 (2000), pp. 42–66.
[4] S. Abrabanel, A. Chertock, and A. Yefet, Strict stability of high-order compact implicit
finite-difference schemes: The role of boundary conditions for hyperbolic PDEs, II, J.
[5] M. H. Carpenter, D. Gottlieb, and S. Abarbanel, Time-stable boundary conditions for
finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes, J. Comput. Phys., 111 (1994), pp. 220–236.
[6] M. H. Carpenter, J. Nordstr¨om, and D. Gottlieb, A stable and conservative interface
treatment of arbitrary spatial accuracy, J. Comput. Phys., 148 (1999), pp. 341–365.
[7] A. Ditkowski, K. Dridi, and J. S. Hesthaven, Convergent Cartesian grid methods for
Maxwell’s equations in complex geometries, J. Comput. Phys., 170 (2001), pp. 39–80.
[8] B. Gustafsson, H.-O. Kreiss, and J. Oliger, Time Dependent Problems and Difference
Methods, John Wiley & Sons, New York, 1995.
[9] J. S. Hestahaven and D. Gottlieb, A stable penalty method for the compressible Navier–
Stokes equations: I. Open boundary conditions, SIAM J. Sci. Comput., 17 (1996),
pp. 579–612.
[10] J. S. Hesthaven and T. Warburton, High-order nodal methods on unstructured grids. I.
Time-domain solution of Maxwell’s equations, J. Comput. Phys., 181 (2002), pp. 186–221.
[11] J. S. Hesthaven and T. Warburton, High-order nodal discontinuous Galerkin methods for
the Maxwell eigenvalue problem, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng.
Sci., 362 (2004), pp. 493–524.
[12] H.-O. Kreiss and G. Scherer, Finite element and finite difference methods for hyperbolic
partial differential equations, in Mathematical Aspects of Finite Elements in Partial
Dif-ferential Equation, C. de Boor, ed., Academic Press, New York, 1974, pp. 195–212. [13] C. Van Loan, Computational Frameworks for the Fast Fourier Transform, Frontiers Appl.
Math. 10, SIAM, Philadelphia, 1992.
[14] J. Nordstr¨om and M. H. Carpenter, High-order finite difference methods, multidimensional
linear problems and curvilinear coordinates, J. Comput. Phys., 173 (2001), pp. 149–174.
[15] J. Nordstr¨om, K. Forsberg, C. Adamsson, and P. Eliasson, Finite volume methods,
un-structured meshes, and strict stability, Appl. Numer. Math., 45 (2003), pp. 453–473.
[16] J. Nordstr¨om and R. Gustafsson, High order finite difference approximations of
electro-magnetic wave propagation close to material discontinuities, J. Sci. Comput., 18 (2003),
pp. 215–234.
[17] S. Tsynkov, V. Ryaben’kii, and V. Turchaninov, Long-time numerical computation of
