### 10.1 Outer boundary conditions

#### 10.1.1 Injection

Injecting boundary conditions is presumably the simplest way to numerically impose them. It implies simply overwriting, at every or some number of timesteps, the numerical solution for each incoming characteristic variable or its time derivative with the conditions that they should satisfy.

Stability of the injection approach can be analyzed through GKS theory [229], since energy estimates are, in general, not available for it (the reason for this should become more clear when discussing the projection and penalty methods). Stability analyses, in general, not only depend on the spatial approximation (and time integration in the fully-discrete case) but are in general also equation-dependent. Therefore, stability is, in general, difficult to establish, especially for high-order schemes and nontrivial systems. For this reason a semi-discrete eigenvalue analysis is many times performed. Even though this only provides necessary conditions (namely, the von Neumann condition (7.80)) for stability, it serves as a rule of thumb and to discard obviously unstable schemes.

As an example, we discuss the advection equation with “freezing” boundary condition,

where the boundary data is chosen to be constant, . As space approximation we use a Chebyshev collocation method at Gauss–Lobatto nodes. The approximation then takes the form
where is the corresponding Chebyshev differentiation matrix; see Section 9.6. The boundary condition is then imposed by replacing the equation by
For this problem and approximation an energy estimate can be derived [209], but we present an eigenvalue analysis as a typical example of those done for more complicated systems and injection boundary conditions. Figure 4 shows the semi-discrete spectrum using a different number of collocation points . As needed by the strong version of the von Neumann condition, no positive real component is present [cf. Eq. (7.82)]. We also note that, as discussed at the beginning of Section 9.8, the spectral radius scales as .

There are other difficulties with the injection method, besides the fact that stability results are usually partial and incomplete for realistic systems and/or high-order or spectral methods. One of them is that it sometimes happens that a full GKS analysis can actually be carried out for a simple problem and scheme, and the result turns out to be stable but not time-stable (see Section 7.4 for a discussion on time-stability). Or the scheme is not time-stable when applied to a more complicated system (see, for example, [116, 1]).

Seeking stable numerical boundary conditions for realistic systems, which preserve the accuracy of high-order finite-difference and spectral methods has been a recurring theme in numerical methods for time-dependent partial differential equations for a long time, especially for nontrivial domains, with substantial progress over the last decade – especially with the penalty method discussed below. Before doing so we review another method, which improves on the injection one in that stability can be shown for rather general systems and arbitrary high-order FD schemes.

#### 10.1.2 Projections

Assume that a given IBVP is well posed and admits an energy estimate, as discussed in Section 5.2. Furthermore, assume that, up to the control of boundary terms, a semi-discrete approximation to it also admits an energy estimate. The key idea of the projection method [317, 319, 318] is to impose the boundary conditions by projecting at each time the numerical solution to the space of gridfunctions satisfying those conditions, the central aspect being that the projection is chosen to be orthogonal with respect to the scalar product under which a semi-discrete energy estimate in the absence of boundaries can be shown. The orthogonality of the projection then guarantees that the estimate including the control of the boundary term holds.

In more detail, if the spatial approximation prior to the imposition of the boundary conditions is written as

and [for example, as is many times the case when SBP holds] there is a semi-discrete energy,
with respect to some scalar product, for which an estimate holds up to boundary terms:
with
for some constant independent of the initial data and resolution. Due to the SBP property, the boundary terms are exactly those present in the continuum energy estimate, except that so far they cannot be bounded because the numerical solution is not yet required to satisfy the boundary conditions.

The latter are then imposed by changing the semi-discrete equation (10.6) to

where projects to the space of gridfunctions satisfying the desired boundary conditions, is time-independent (which in particular requires the assumption that the boundary conditions are time-independent) and symmetric
(and, being a projection, ). Since the projection is assumed to be time-independent, projecting both sides of Eq. (10.6) implies that
and for any solution of (10.6) satisfying the boundary conditions, ,
Then
Since the solution now satisfies the boundary conditions, the boundary terms can be bounded as in the continuum and stability of the semi-discrete IBVP follows. In principle, this method only guarantees stability for time-independent, homogeneous boundary conditions. Homogeneity can be assumed by a redefinition of the variables; however, the restriction of a time-independent boundary condition is more severe.

Details on how to explicitly construct the projection can be found in [227]. The orthogonal projection method guarantees stability for a large class of problems admitting a continuum energy estimate. However, its implementation is somewhat involved.

#### 10.1.3 Penalty conditions

A simple and robust method for imposing numerical boundary conditions, either at outer or interpatch boundaries, such as those appearing in domain decomposition approaches, is through penalty terms. The boundary conditions are not imposed strongly but weakly, preserving the convergence order of the spatial approximation and leading to numerical stability for a large class of problems. It can be applied both to the case of FDs and spectral approximations. In fact, the spirit of the method can be traced back to Finite Elements-discontinuous Galerkin methods (see [147] and [28] for more recent results). Terms are added to the evolution equations at the boundaries to consistently penalize the mismatch between the numerical solution and the boundary conditions the exact solution is subject to.

##### Finite differences.
In the FD context the method is known as the Simultaneous Approximation Technique (SAT) [117]. For semi-discrete approximations of IBVPs of arbitrary high order with an energy estimate, both the order of accuracy and the presence of energy estimates are preserved when imposing the boundary conditions through it.

As an example, consider the half-space IBVP for the advection equation,

where is the characteristic speed in the direction.

We first consider a semi-discrete approximation using some FD operator satisfying SBP with respect to a scalar product , which we assume to be either diagonal or restricted full (see Section 8.3). As usual, denotes the spacing between gridpoints,

In the SAT approach the boundary conditions are imposed through a penalty term, which is applied at the boundary point ,
where is the Kronecker delta, a real parameter to be restricted below, and is the -component of the SBP scalar product . Introducing the semi-discrete SBP energy
its time derivative is
For homogeneous boundary conditions () both numerical and time-stability follow if :
Strong stability follows if , including the case of inhomogeneous boundary conditions: for any such that
Eq. (10.20) implies
where
Integrating Eq. (10.23) in time,
thereby proving strong stability.

In the case of diagonal SBP norms it is straightforward to derive similar energy estimates for general linear symmetric hyperbolic systems of equations in several dimensions, simply by working with each characteristic variable at a time, at each boundary. A penalty term is applied to the evolution term of each incoming characteristic variable at a time, as in Eq. (10.18), where is replaced by the corresponding characteristic speed. In particular, edges and corners are dealt with by simply imposing the boundary conditions with respect to the normal to each boundary, and an energy estimate follows.

The global semi-discrete convergence rate can be estimated as follows. Define the error gridfunction as the difference between the numerical solution and the exact one evaluated at the gridpoints,

It satisfies the equation
where denotes the truncation error, here solely depending on the differentiation approximation:
with in general. For example, in the diagonal case , and in the restricted full one [see the discussion below Eq. (8.29)].

Using Eq. (10.26) and the SBP property, the norm of the error satisfies

where we have used and the Cauchy–Schwarz inequality in the second step. Dividing both sides of the inequality by and integrating we obtain
Since by Eq. (10.29) the norm of the truncation error is of the order , it follows that the -norm of the error converges to zero with rate . This also implies that the error converges pointwise to zero with rate,
In particular, this proves that the error at the boundary point converges to zero at least as fast as , and therefore, the penalty term in Eq. (10.18) is bounded. Two remarks are in order:
• In special cases it is possible to improve the error estimate exploiting the strong stability of the problem. Consider, for instance, the case of the operator defined in Example 51 with the associated diagonal scalar product . Then, Eq. (10.30) gives
where we have chosen . This implies that the error converges to zero in the -norm with second order and pointwise with order .
• At any fixed resolution, the error typically decreases with larger values of the penalty parameter , but the spectral radius of the discretization grows quickly in the process by effectively introducing dissipative eigenvalues (on the left half plane) in the spectrum, leading to demanding CFL limits (see, for example, [281]). Because the method is usually applied along with high-order methods, decreasing the error for a fixed resolution at the expense of increasing the CFL limit does not seem worthwhile. In practice values of in the range give reasonable CFL limits.

##### Spectral methods.
The penalty method for imposing boundary conditions was actually introduced for spectral methods prior to FDs in [198, 199]. In fact, as we will see below, the FD and spectral cases follow very similar paths. Here we only discuss its application to the collocation method. Furthermore, as discussed in Section 9.4, when solving an IBVP, Gauss–Lobatto collocation points are natural among Gauss-type nodes, because they include the end points of the interval. We restrict our review to them, but the penalty method applies equally well to the other nodes. We refer to [236] for a thorough analysis of spectral penalty methods.

Like the FD case, we summarize the method through the example of the advection problem (10.14, 10.15, 10.16), except that now we consider the bounded domain and apply the boundary condition at . Furthermore, we first consider a truncated expansion in Legendre polynomials. A penalty term with strength is added to the evolution equation at the last collocation point:

where are now the Gauss–Lobatto nodes; in particular, and .

Using the discrete energy given by Eq. (9.95) and SBP property associated with Gauss quadratures discussed in Section 9.4, a discrete energy estimate follows exactly as in the FD case if

where is the Legendre–Gauss–Lobatto quadrature weight at the right end point [cf. Eqs. (9.63, 9.65)]. Notice that this scaling of the penalty with the weight is the analog of the SBP weight term for the FD case [Eq. (10.18)]. The similarity is not accidental: both weights arrive from the discrete integration formulae for the energy. Notice:
• The approach for linear symmetric hyperbolic systems is the same as that we discussed for FDs: the evolution equation for each characteristic variable is penalized as in Eq. (10.33), where is replaced by the corresponding characteristic speed.

Devising a Chebyshev penalty scheme that guarantees stability is more convoluted. In particular, the advection equation is already not in the Chebyshev norm, for example (see [209]). Stability in the norm can be established using the Chebyshev–Legendre method [144], where the Chebyshev–Gauss–Lobatto nodes are used for the approximation, but the Legendre ones for satisfying the equation. In this approach, the penalty method is global, because it adds terms to the right-hand side of the equations not only at the endpoint, but at all other collocation points as well.

A simpler approach, where the penalty term is only applied to the boundary collocation point, as in the Legendre and FDs case, is to show stability in a different norm. For example, in [137] it is shown that a penalty term as in Eq. (10.33) is stable for the Chebyshev–Gauss–Lobatto case in the norm defined by the weight, cf. Eq. (9.17),

if
In Figure 5 we show the spectrum of an approximation of the advection equation (with speed and ) using the Chebyshev penalty method, collocation points and . There are no eigenvalues with positive real component. Even though this is only a necessary condition for stability, it suggests that it might be possible to lower the bound (10.36) on the minimum penalty needed for stability. On the other hand, for (and lower values) real positive values do appear in the spectrum. Compared to the injection case (Figure 4), the spectral radius for the cases shown in Figure 5 are about a factor of two larger. However, in most applications using the method of lines, the timestep is dictated by the condition of keeping the time integration error smaller than the spatial one, so as to maintain spectral convergence, and not the CFL limit.