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 , 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,, 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.
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
The latter are then imposed by changing the semi-discrete equation (10.6) to
Details on how to explicitly construct the projection can be found in . The orthogonal projection method guarantees stability for a large class of problems admitting a continuum energy estimate. However, its implementation is somewhat involved.
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  and  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.
As an example, consider the half-space IBVP for the advection equation,
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 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,
Using Eq. (10.26) and the SBP property, the norm of the error satisfiespointwise to zero with rate,
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:
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
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 ). Stability in the norm can be established using the Chebyshev–Legendre method , 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  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),
Living Rev. Relativity 15, (2012), 9
This work is licensed under a Creative Commons License.