Mathematical Physics and Quantum Field Theory, Electronic Journal of Differential Equations, Conf. 04, 2000, pp. 245–263 http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)

Localization of dependence for solutions of hyperbolic differential equations ^{∗}

Henry A. Warchall

Dedicated to Eyvind H. Wichmann

Abstract

We survey several results that localize the dependence of solutions to hyperbolic equations. These observations address questions that are central to numerical simulation of solutions on unbounded spatial domains. One result shows that in principle it is possible to numerically compute (the restriction of) a solution to a wave equation on an unbounded domain using only a bounded computational domain. Other results provide implementations of this fact in particular situations. In addition, we introduce a new diagrammatic way to generate explicit solutions to multiple-time initial-value problems for the wave equation in one space dimension.

1 Introduction

One of the many things I learned from Eyvind Wichmann is that it can be profitable to re-examine topics you know well; it is possible to see aspects of a subject that you missed the first time. In that spirit, I would like to revisit an equation we all know and love: the wave equation. This discussion reviews results on wave propagation in two articles of mine, and it introduces two new observations which I would like to present to Eyvind on this happy occasion.

In this first section, we will pose a question about the domain of dependence of solutions to the wave equation that is central to numerical simulation of solutions on unbounded domains. In the second section we will answer the question with a summary of the general results in [1] of wave propagation at computational domain boundaries, applied for simplicity to partial differential equations whose principal part is the d’Alembertian. In the third section we will discuss two schemes that implement the intent of the abstract result of the second section but not its precision. The first scheme has not been previously published, and the second scheme is found in [2]. The fourth section presents a new toy (also not previously published), a diagrammatic method to generate formulas for solutions of the (1 + 1)-dimensional wave equation with given multiple-time initial data.

∗Mathematics Subject Classifications: 35B30, 35L05, 35L15, 35L70, 35C10, 35A18.

Key words: Localization of dependence, wave equation, computational domain boundary,

exact nonreflecting boundary conditions.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University and University of North Texas.

Published November 3, 2000.

We begin with a gedanken experiment. Consider an infinite lake, with no shore or islands. We observe a fixed finite region of the surface. Suppose that the entire lake is initially quiet, except for a disturbance in the region under observation. We look away, then observe the same region at some later time. Does the data at the later time in the region alone suffice to predict the water’s future behavior in the region?

Intuitively speaking, we expect to be able to determine the future behavior of the surface in the region, given only data in the region at the intermediate time, because “all waves are outgoing.” On the other hand, we know that the domain of dependence for the solution at points near the edge of the region under observation extends outside the region, so it seems we may need data from outside the region to determine the solution inside the region in the future.

To make this question more precise, consider the computation of solutions to a (possibly nonlinear) wave equation on an unbounded spatial domain (say ℝ^{n}). Assume that the equation is autonomous, that the speed of propagation is finite, and that the initial data is supported inside a bounded domain on which values of the solution are to be computed. We desire to compute (only) the restriction to the bounded computational domain of the solution to the initial value problem on the unbounded domain. (In particular, the computational domain boundaries are not boundaries for the initial value problem, and do not affect the solution in any way.) Suppose we have an evolution scheme (for instance, numerical) to advance the solution stepwise in time on the spatial domain ℝ^{n}, starting from the given compactly-supported initial data. We apply the scheme to advance the solution in time until its spatial support reaches the boundary of the computational domain. It is not clear that we can continue computing the solution after that time using only data inside the computational domain, since the solution’s a priori domain of dependence for points near the computational domain boundary at later times extends outside the computational domain, as illustrated in Figure 1.

Usual (Dirichlet, Neumann, periodic) boundary conditions applied at the computational domain boundary generate spurious reflections that are not part of the solution on the unbounded domain. One can do better: approximatelyreflectionless boundary conditions have been developed (see [4] and references therein) but have varying degrees of accuracy. The application of such conditions at computational domain boundaries, yielding only approximations to solutions on the unbounded domain, is philosophically dissatisfying. In the situation outlined, data inside the computational domain at an intermediate time should suffice to determine the solution inside the domain at later times, because the solution is “outgoing” and its behavior outside the domain should not influence its later behavior inside. The remainder of our discussion is motivated by the two main questions:

To what extent can this intuition be made precise?

To what extent can it be made into a calculational tool?

Figure 1: Domain of dependence extends outside computational domain

Here we advocate an alternative to differential boundary conditions: propagation of waves through artificial domain boundaries.

An explicit example is in order. Consider the wave equation

u_{tt} − u_{xx} = 0

in one spatial dimension. The solution to the initial-value problem with data at time t = t_{0} is given explicitly by d’Alembert’s formula

2u(x, t) = u(x − Δt, t_{0}) + u(x + Δt, t_{0}) + ∫_{x−Δt}^{x+Δt} u_{t}(y, t_{0}) dy

where Δt ≡ t − t_{0}.

Suppose it is known that the support of the (t = t_{0}) initial data is to the left of x = b: supp {u(·, t_{0}), u_{t}(·, t_{0})} ⊂ (−∞, b). Suppose u and u_{t} are known to the left of x = b at time t_{1} > t_{0}. Given x ≤ b, is the solution at x and time t_{2} > t_{1} determined by u and u_{t} to the left of b at time t_{1} ?

The answer is yes:

Theorem 1.1 Suppose u(x, t) is a classical solution of u_{tt} − u_{xx} = 0 with support of initial (t = t_{0}) data to the left of x = b. Let t_{2} > t_{1} > t_{0} and set Δt ≡ t_{2} − t_{1}. If x ∈ (b − Δt, b + Δt) then

2u(x, t_{2}) = u(x − Δt, t_{1}) + u(b, t_{1}) + ∫_{x−Δt}^{b} u_{t}(y, t_{1}) dy.

Proof There are two parts to the proof. First part: We claim that u_{t}(x, t_{1}) = −u_{x}(x, t_{1}) for x ≥ b. Since both u and u_{t} vanish at t_{0} for x ≥ b, using the d’Alembert formula to advance from t_{0} to t_{1} gives

2u(x, t_{1}) = u(x − t_{1} + t_{0}, t_{0}) + ∫_{x−t_{1}+t_{0}}^{b} u_{t}(y, t_{0}) dy.

Note that the right-hand side depends only on (x − t_{1}), not (x + t_{1}). Thus u_{t}(x, t_{1}) = −u_{x}(x, t_{1}) for x ≥ b, as claimed. Second part: We use d’Alembert’s formula again to advance the solution from t_{1} to t_{2}, obtaining

2u(x, t_{2}) = u(x−Δt, t_{1})+u(x+Δt, t_{1})+∫_{x−Δt}^{b} u_{t}(y, t_{1}) dy+ ∫_{b}^{x+Δt} u_{t}(y, t_{1}) dy.

But, by virtue of the first part, the last integrand is just −u_{x}(y, t_{1}), so the last integral yields u(b, t_{1}) − u(x + Δt, t_{1}), and the assertion of the theorem follows immediately. QED

In particular, if we think of points x to the left of b as being inside a computational domain, and points to the right as being outside, the resultant formula

2u(b, t_{2}) = u(b − Δt, t_{1}) + u(b, t_{1}) + ∫_{b−Δt}^{b} u_{t}(y, t_{1}) dy

gives an explicit method for advancing the solution at the boundary, using only data inside the computational domain. Numerical implementation of this integral formula (which is equivalent to the boundary condition u_{t} = −u_{x}) works superbly to propagate waves through the computational domain boundary as if it were not there.

It is natural to ask how far this approach can be taken in treating higher numbers of spatial dimensions and other equations. In the next section we will see that in principle the same situation holds in a wide variety of cases. Note that for many computational purposes, it suffices to develop such “one-sided” propagation formulas for linear constant-coefficient equations. The reason is that in studies of nonlinear wave equations, an effective numerical simulation will have the large-amplitude part of the solution well inside the computational domain for times of interest, and the solution near the boundary will have small amplitude. Then propagation near the boundary is well-approximated by the linearized equation, and results for linear wave equations can be employed to approximate the solution near the boundary.

2 Abstract Uniqueness Results

Quite general results on wave propagation at computational domain boundaries are given in [1] for linear hyperbolic equations with analytic coefficients. Here we will specialize these results to equations that are close to the wave equation.

Let B(x, r) ⊂ ℝ^{n} be the closed ball with center x and radius r. Define the closed forward cone V ≡ {(x, t) ∈ ℝ^{n+1} | t ≥ |x| }. We denote by Ω^{c} the complement of the set Ω.

Theorem 2.1 (Warchall) Let P (x, D) be a linear differential operator on ℝ^{n+1} with analytic coefficients and principal part ∂_{t}^{2} − Δ_{x}. Let Ω ⊂ ℝ^{n} be an open convex set. Let t_{0} < t_{1} < t_{2}, and let z ∈ Ω̄.

Let f be a continuous function of t ≡ x_{n+1} with values in D′(ℝ^{n}) that is analytic in Ω^{c} × (−∞, ∞). Suppose f vanishes on an open set that contains {(x, t) ∈ {(z, t_{2})} − V | x ∈ Ω, t ≥ t_{1} }.

Suppose u ∈ D′(ℝ^{n+1}) is a solution to P (x, D)u = f such that:

(i) the data u(t_{0}) and u_{t}(t_{0}) have support in Ω, and (ii) the data u(t_{1}) and u_{t}(t_{1}) vanish on an open set that contains Ω∩B(z, t_{2}−t_{1}).

Then u vanishes in an open neighborhood of (z, t_{2}).

This result establishes that the solution (of the linear equation P u = f) in the domain at time t_{2} is uniquely determined by the restriction of the solution at earlier time t_{1} to the intersection of the computational domain with the a priori domain of dependence, as illustrated in Figure 2.

Figure 2: Solution at (z, t_{2}) is determined by data at time t_{1} in shaded region

Thus, if it is known that at some past time the data for the solution was supported in Ω, then the data in Ω at time t_{1} completely determines the solution in Ω at all later times. The drawback to this result is that the dependence is not made explicit, and in fact is known only for the wave equation in one space dimension.

About the proof of Theorem 2.1

The proof involves an argument about propagation of regularity. Let W F_{A}(u) be the analytic wave front set of the distribution u ∈ D′(ℝ^{N}). We make use of the following result. (Thanks to James Ralston for pointing it out.)

Theorem 2.2 (Hörmander [3]) Let P (x, D) be a differential operator of order m with analytic coefficients and real principal part P_{m}. Let L be a segment of a bicharacteristic strip of P, with \frac{∂P_{m}}{∂ξ} (x, ξ) ≠ 0 on L. Let f ∈ D′(ℝ^{N}) be such that Λ ∩ W F_{A}(f) is empty. Suppose that u ∈ D′(ℝ^{N}) solves P (x, D)u = f. Then either Λ ⊂ W F_{A}(u) or Λ ∩ W F_{A}(u) is empty.

This propagation of singularities along bicharacteristic strips also results in propagation of regularity:

Corollary 2.3 Let P (x, D) be a differential operator with analytic coefficients and real principal part P_{m}. Let u ∈ D′(ℝ^{N}), and set f = P (x, D)u. Suppose x ∈ ℝ^{N} is such that:

(1) f is real analytic at x;

(2) dP_{m}(x, ξ) ≠ 0 if (x, ξ) ∈ Char(P); and

(3) each bicharacteristic strip of P_{m} that contains (x, ξ) for some ξ ≠ 0 also contains some point q not in W F_{A}(u), with all points (y, η) of the strip between and including (x, ξ) and q satisfying (y, η) ∉ W F_{A}(f) and \frac{∂P_{m}}{∂η} (y, η) ≠ 0.

Then u is real analytic at x.

To establish this corollary, we note that by hypothesis (1), W F_{A}(u) ⊂ Char(P). Hypotheses (2) and (3) allow us to connect each point (x, ξ) ∈ Char(P) with a bicharacteristic strip to a point not in W F_{A}(u). Now apply Hörmander’s theorem to the strip to conclude (x, ξ) ∉ W F_{A}(u). Since this is true for all (x, ξ) ∈ Char(P), we have that u is real analytic at x, as claimed.

We may apply this Corollary to prove Theorem 2.1 as follows. Given z ∈ Ω̄, consider the timelike segment C ⊂ ℝ^{n+1} given by C ≡ {z} × [t_{1}, t_{2}]. Fix a point (z, t) ∈ C. Because Ω is convex, each backward characteristic (lightlike line) through (z, t) passes either through Ω^{c} × {t_{0}} or through (Ω ∩ B(z, t_{2} − t_{1})) × {t_{1}}. Because of the hypotheses on the data and on f, the solution u vanishes on an open neighborhood of each of these sets, and u is thus analytic at the (lower) endpoint of the segment between (z, t) and the corresponding set. Because the hypotheses insure that f is analytic along each such segment, and because the principal part P_{m} of the differential operator is the d’Alembertian, the hypotheses of the Corollary are satisfied, and we conclude that u is real analytic at (z, t).

The same argument holds at each point (z, t) ∈ C, and thus u is analytic on C. The hypotheses on the support of the data imply that u vanishes in an open neighborhood of the lower endpoint (z, t_{1}) of C, and thus by analytic continuation u vanishes on C. In particular, u vanishes in an open neighborhood of (z, t_{2}), as asserted. QED

3 Explicit but Nonlocal Schemes

The abstract result of Theorem 2.1 answers our first main question, but it does not provide an explicit algorithm by which to determine the future solution in the entire computational domain from data in the domain alone. In this section we discuss two explicit schemes for advancing the solution in the computational domain. Unfortunately, neither scheme implements the precise local dependence established by Theorem 2.1. The implementation of that dependence is an open question for all cases other than the wave equation in one spatial dimension.

Piecewise-in-time method

There is an explicit algorithm for using bounded computational domains to simulate certain problems on unbounded domains that is in principle exact and valid in arbitrary spatial dimension. (It does not, however, implement the local dependence established in Section 2.) This algorithm is inspired by an observation [5] by Cathleen Morawetz.

Algorithm based on explicit solutions

We consider a partial differential equation

P (x, t; D)u + G(x, t; u) = F (x, t)

where the linear hyperbolic differential operator P has uniformly bounded propagation speed, F is independent of u, and G is a function that can be nonlinear in u. To simplify notation, we will write initial conditions for the equation in the schematic first-order form u|_{t=t_{0}} = u_{0}.

Two hypotheses are crucial for the algorithm. First, we assume that G has spatial support in a fixed ball B ⊂ ℝ^{n}, that is, G(x, t; u) vanishes identically in the variables t and u when x is outside B. This condition is a strict realization of the idea that the equation under study is linear sufficiently far from the origin.

The second hypothesis is that the initial-value problem

P (x, t; D)v = F (x, t)
v|_{t=t_{0}} = v_{0}

for the linear system without the G term is explicitly solvable in the sense that the solution v is assumed to be readily computable everywhere from the initial data v_{0}, not requiring local numerical integration to advance v in time. For example, P might be a constant-coefficient operator such as the d’Alembertian or the Klein-Gordon operator, for which the solution v can be expressed in terms of F and the initial data v_{0} by an integral formula. The crux of this hypothesis is that v can be determined without advancing it in time with a local algorithm whose domain of dependence includes points outside the computational domain where the data is nonvanishing.

Let Ω ⊂ ℝ^{n} be an open, bounded (computational) domain containing B, with dist(B, ∂Ω) > 0. To compute the restriction to Ω̄ of the solution to the initial-value problem

\begin{matrix}P (x, t; D)u + G(x, t; u) = F (x, t) \\ u|_{t=t_{0}} = u_{0} \end{matrix}  (3.1)

we write the solution u as the sum of the solutions v^{(0)}, v^{(1)}, and w^{(1)} to the three initial-value problems

\begin{matrix}P (x, t; D)v^{(0)} = F (x, t) \\ v^{(0)}|_{t=t_{0}} = v_{0}^{(0)} \end{matrix}

and

\begin{matrix}P (x, t; D)v^{(1)} = 0 \\ v^{(1)}|_{t=t_{0}} = v_{0}^{(1)} \end{matrix}

and

\begin{matrix}P (x, t; D)w^{(1)} + G(x, t; v^{(0)} + v^{(1)} + w^{(1)}) = 0 \\ w^{(1)}|_{t=t_{0}} = w_{0}^{(1)} \end{matrix}

where supp w_{0}^{(1)} ⊂ B, supp v_{0}^{(1)} ⊂ Ω̄, and v_{0}^{(0)} + v_{0}^{(1)} + w_{0}^{(1)} = u_{0}. (It would be sufficient to take v_{0}^{(1)} = w_{0}^{(1)} = 0, but it may be convenient to make other choices.)

The solutions v^{(0)} and v^{(1)} are by assumption explicitly known functions of F and the initial data, for all time t ≥ t_{0}. The solution w^{(1)} is to be advanced in time with a numerical algorithm.

Let t_{1} = t_{0} + τ be the earliest time greater than t_{0} at which the spatial support of w^{(1)} can intersect ∂Ω. Because G(x, t; u) vanishes for x outside B, and because supp w_{0}^{(1)} ⊂ B, the finite propagation speed of P implies that t_{1} is strictly larger than t_{0}. During the time interval [t_{0}, t_{1}] the spatial support of w^{(1)} is contained in the computational domain Ω̄.

The numerical integration of the initial-value problem for w^{(1)} is to be halted at time t = t_{1}. The solution u to the original initial-value problem (3.1) is then given by u = v^{(0)} + v^{(1)} + w^{(1)} in Ω̄ × [t_{0}, t_{1}].

We next continue the time evolution of w^{(1)} by considering the two initialvalue problems

\begin{matrix}P (x, t; D)v^{(2)} = 0 \\ v^{(2)}|_{t=t_{1}} = w^{(1)}|_{t=t_{1}} \end{matrix} .

and

\begin{matrix}P (x, t; D)w^{(2)} + G(x, t; v^{(0)} + v^{(1)} + v^{(2)} + w^{(2)}) = 0 \\ w^{(2)}|_{t=t_{1}} = 0 \end{matrix}

in terms of which w^{(1)} = v^{(2)}+w^{(2)}. The solution v^{(2)} is by hypothesis an explicitly known function of the data w^{(1)}|_{t=t_{1}}. Note that the evolution equation for w^{(2)} is driven by the known function v^{(0)} +v^{(1)} +v^{(2)}. As noted, the hypotheses on P and G insure finite propagation speed for w^{(2)} outside B. Again w^{(2)} is advanced in time with the numerical algorithm until time t_{2} = t_{1} + τ, when its spatial support can reach ∂Ω. The solution u to the original initial-value problem (3.1) is then given by u = v^{(0)} + v^{(1)} + v^{(2)} + w^{(2)} in Ω̄ × [t_{1}, t_{2}].

This latter procedure can be continued arbitrarily many times to advance the solution by time τ at each step. In general, at the beginning of the k^{th} step, the solution u is known in Ω̄ × [t_{0}, t_{k−1}], being given in Ω̄ × [t_{k−2}, t_{k−1}] by u = v^{(0)} + ⋯ + v^{(k−1)} + w^{(k−1)}. Here each function v^{(j)}, j = 1, 2, …, is known explicitly for all t ≥ t_{j−1}, and the function w^{(k−1)} has spatial support in Ω̄ at time t_{k−1} = t_{0} + (k − 1)τ.

The k^{th} step of the computation proceeds with the explicit determination of v^{(k)} for t ≥ t_{k−1} by the initial-value problem

\begin{matrix}P (x, t; D)v^{(k)} = 0 \\ v^{(k)}|_{t=t_{k−1}} = w^{(k−1)}|_{t=t_{k−1}} \end{matrix}

followed by numerical solution of

\begin{matrix}P (x, t; D)w^{(k)} + G(x, t; v^{(0)} + v^{(1)} + ⋯ + v^{(k)} + w^{(k)}) = 0 \\ w^{(k)}|_{t=t_{k−1}} = 0 \end{matrix}

to determine w^{(k)} in Ω̄ × [t_{k−1}, t_{k}]. The solution u to the original initial-value problem (3.1) is given in Ω̄ × [t_{k−1}, t_{k}] by u = v^{(0)} + v^{(1)} + ⋯ + v^{(k)} + w^{(k)}. In this way we can integrate this nonlinear equation using numerical computation only in Ω̄.

Algorithm based on the wave equation

Morawetz’ original result [5] concerns the case where P is the d’Alembertian □ ≡ ∂_{t}^{2}−Δ in three spatial dimensions. In that case it is possible to do away with the second hypothesis concerning explicit representation of solutions to initialvalue problems for the linear equation, by making use of (the strong) Huygens principle. (Thanks to Cathleen Morawetz for notifying me of this unpublished result.) In particular, let B be the open ball in ℝ^{3} of radius b, centered (without loss of generality) at the origin. We can compute the restriction to B of solutions to the initial-value problem

\begin{matrix}□u + G(x, t; u) = 0 \\ u|_{t=0} = u_{0} \end{matrix}  (3.2)

on ℝ^{3+1} by computing only inside the concentric open ball Ω with radius 3b. (Here we continue the notational abuse of writing initial conditions in schematic first-order form.) We proceed as follows.

We assume that the classical initial data u_{0} is supported in B, and, as before, that the continuous function G(x, t; u) has spatial support in B for all t and u. The first step consists of computing the solution to the initial-value problem

\begin{matrix}□w^{(1)} + G(x, t; w^{(1)}) = 0 \\ w^{(1)}|_{t=0} = u_{0} \end{matrix}

in Ω × [0, t_{1}], where t_{1} = 2b. This can be done with a straightforward numerical algorithm since the spatial support of the solution is contained in Ω for all t ∈ [0, t_{1}], by virtue of the unit speed propagation outside B. Thus we have u = w^{(1)} in Ω × [t_{0}, t_{1}]. (In terms of the earlier notation, F ≡ 0, t_{0} = 0, τ = 2b, v_{0}^{(0)} = v_{0}^{(1)} = 0, w_{0}^{(1)} = u_{0}, and v^{(0)} = v^{(1)} ≡ 0.)

We continue the time evolution of w^{(1)} beyond time t_{1} by solving the two initial-value problems

\begin{matrix}□v^{(2)} = 0 \\ v^{(2)}|_{t=t_{1}} = w^{(1)}|_{t=t_{1}} \end{matrix}  (3.3)

and

\begin{matrix}□w^{(2)} + G(x, t; v^{(2)} + w^{(2)}) = 0 \\ w^{(2)}|_{t=t_{1}} = 0 \end{matrix}  (3.4)

in terms of which w^{(1)} = v^{(2)} + w^{(2)} for t ≥ t_{1}, as before.

Consider first the initial-value problem (3.3) for v^{(2)}. Because the support of the initial data is contained in Ω, and because the wave equation in odd spatial dimensions three or greater enjoys the strong form of Huygens’ principle, the solution vanishes in the forward cone V_{5} ≡ {(0, 5b)} + V = {(x, t) ∈ ℝ^{3+1} | t ≥ |x| + 5b }. In particular, v^{(2)} vanishes in B for t ≥ 6b. See Figure 3.

Figure 3: Spacetime regions in the piecewise-in-time algorithm based on the wave equation

But we can say more. We claim that v^{(2)} vanishes in the cone V_{3} ≡ {(0, 3b)} + V. To see this, let G̃ be G cut off at time t_{1}:

G̃(x, t; u) = \{ \begin{matrix}G(x, t; u) & if t ≤ t_{1} \\ 0 & if t > t_{1} \end{matrix}

Thus G̃ has spacetime support in B ×(−∞, t_{1}]. Then the solution to the initialvalue problem (3.2) and the solution to the initial-value problem

\begin{matrix}□ũ + G̃(x, t; ũ) = 0 \\ ũ|_{t=0} = u_{0} \end{matrix}  (3.5)

are identical in Ω×[0, t_{1}], and w^{(1)}|_{t=t_{1}} is the data at time t_{1} for both problems. Further, the problem (3.3) gives the solution of (3.5) for t > t_{1}, that is, ũ = v^{(2)} for t > t_{1}.

Let (ξ, τ) be a point in V_{3}. We can make use of Beltrami’s formula [6] or Duhamel’s principle [7] to express the solution ũ(ξ, τ) in terms of the initial data u_{0} on the sphere |x − ξ| = τ and an integral of G̃(x, t; ũ) over the truncated backward cone A ≡ {(ξ, τ)} − {(x, t) ∈ ℝ^{3+1} | |x| ≤ t ≤ τ }. Because u_{0} is supported in B and thus vanishes on the sphere |x − ξ| = τ, and because supp G̃ does not intersect A, it follows that v^{(2)}(ξ, τ) = ũ(ξ, τ) = 0 for all (ξ, τ) ∈ V_{3}, as claimed.

Thus in particular v^{(2)} vanishes in B for t ≥ t_{2} ≡ 4b. Since we are interested only in the values of v^{(2)} in B, we thus need only compute v^{(2)} in B × [t_{1}, t_{2}]. For this purpose, it more than suffices to apply a numerical algorithm to solve (3.3) in the truncated backward cone A_{5} ≡ {(0, 5b)} − {(x, t) ∈ ℝ^{3+1} | |x| ≤ t ≤ 3b }, which can be done in a straightforward fashion since for time evolution in that region the domain of dependence does not extend outside Ω. We have no need for, and do not compute, values of v^{(2)} in the complement of A_{5}.

We can thus determine the solution v^{(2)} in B × [t_{1}, ∞) using a numerical algorithm and computing only inside Ω. To complete the second step, we employ a numerical algorithm to advance w^{(2)} in time from t_{1} to t_{2}. As in the first step, this can be done in a straightforward manner since the spatial support of the solution is contained in Ω for all t ∈ [t_{1}, t_{2}], by virtue of the unit speed propagation outside B. Then u = v^{(2)} + w^{(2)} in B × [t_{1}, t_{2}].

Succeeding steps of the algorithm proceed in the same way. At the beginning of the k^{th} step, the solutions v^{(1)}, v^{(2)}, …, v^{(k−1)} have been computed in B, and all vanish in B × [t_{k−1}, ∞], where t_{k} = 2kb. The solution to (3.2) is given in B × [t_{j−1}, t_{j}] by u = v^{(j)} + w^{(j)} for j = 1, …, (k − 1). The solution w^{(k−1)} is supported in Ω at time t_{k−1}, and we advance w^{(k−1)} from time t_{k−1} to time t_{k} by solving the two initial-value problems

\begin{matrix}□v^{(k)} = 0 \\ v^{(k)}|_{t=t_{k−1}} = w^{(k−1)}|_{t=t_{k−1}} \end{matrix}  (3.6)

and

\begin{matrix}□w^{(k)} + G(x, t; v^{(k)} + w^{(k)}) = 0 \\ w^{(k)}|_{t=t_{k−1}} = 0 \end{matrix}  (3.7)

in terms of which w^{(k−1)} = v^{(k)} + w^{(k)} for t ≥ t_{k−1}, in exactly the same fashion as in the second step. Then u = v^{(k)} + w^{(k)} in B × [t_{k−1}, t_{k}].

Thus for this class of problems the procedure furnishes a numerical scheme that computes the restriction of the solution to B, using only computation in Ω.

Spherical harmonic expansion

The second scheme to determine the future solution in the entire computational domain from data in the domain alone is based on a (truncated) expansion in spherical harmonics. This scheme implements the evolution with a computation that is local in radius.

We consider the wave equation \frac{∂^{2}u}{∂t^{2}} −Δu = f with source f in 3+1 spacetime dimensions. We assume that at all times t the source f has spatial support in the ball B of radius b centered, say, at the origin. The initial data for u at time t_{0} is supported in B.

Let t_{1} and t_{2} be two later times with t_{2} > t_{1} > t_{0}, and set Δt ≡ t_{2} − t_{1}. Let z ∈ ℝ^{3} be a point outside B, at distance Δt or greater from B, and set a ≡ |z|. Thus a > b+Δt. The uniqueness result of Theorem 2.1 implies that u(z, t_{2}) and \frac{∂u}{∂t}(z, t_{2}) are completely determined by the data at time t_{1} in the intersection of the ball of radius a centered at the origin, and the ball of radius Δt centered at z. This intersection is represented by the shaded region in Figure 4.

Figure 4: Spacetime regions in the spherical harmonic expansion algorithm

The spherical-harmonic expansion scheme does not quite make this dependence explicit, but instead furnishes an ℓ-dependent one-sided propagation formula for the ℓ^{th} partial wave in the spherical harmonic decomposition of u outside of B. This gives a formula for u(z, t_{2}) in terms of (radial derivatives of) the data on a sphere centered at the origin of radius a − Δt at time t_{1}. This sphere is represented by the heavy circle in Figure 4.

In [2] we build on the idea of Joseph Keller and Marcus Grote [8] to expand u using spherical harmonics to get partial waves, and then to determine an operator that converts the partial waves into solutions of the (1+1)-dimensional wave equation. (Thanks to Joseph Keller for notifying me of this work prior to its publication.) Our contribution is to use a differential operator instead of an integral operator, and combine the resultant “outgoing wave condition” with the explicit propagation formula for the wave equation in 3 + 1 dimensions to obtain a single-point propagation formula.

Outgoing wave condition

Let (r, θ, ϕ) be the usual polar coordinates for ℝ^{3}. For x ∈ ℝ^{3}\B, we expand the solution u of the wave equation in terms of spherical harmonics:

u(x, t) = ∑_{ℓ=0}^{∞} ∑_{m=−ℓ}^{ℓ} u_{ℓm}(x, t)

where u_{ℓm}(x, t) ≡ v_{ℓm}(r, t)Y_{ℓm}(θ, ϕ), where Y_{ℓm} is the usual spherical harmonic function. The coefficient function v_{ℓm} of the partial wave solution u_{ℓm} is given by

v_{ℓm}(r, t) ≡ ∫_{0}^{2π} ∫_{0}^{π} \overline{Y_{ℓm}(θ, ϕ)} u(r, θ, ϕ, t) sin θ dθ dϕ

and satisfies the reduced partial differential equation

\frac{∂^{2}v_{ℓm}}{∂t^{2}} = v′′_{ℓm} + \frac{2}{r}v′_{ℓm} − \frac{ℓ(ℓ + 1)}{r^{2}}v_{ℓm}

where the prime denotes differentiation with respect to r.

For notational convenience we define the differential operator

L_{ℓ} ≡ r^{ℓ} (−\frac{∂}{∂r} \frac{1}{r})^{ℓ}

and its formal adjoint

L_{l}^{∗} ≡ (\frac{1}{r} \frac{∂}{∂r} )^{ℓ} [r^{ℓ}·] .

Set

w_{ℓm}(r, t) ≡ L_{ℓ}^{∗} [r v_{ℓm}] (r, t);

then it is not difficult to show ([2]) that w_{ℓm} satisfies the (1 + 1)-dimensional wave equation

\frac{∂^{2}w_{ℓm}}{∂t^{2}} = \frac{∂^{2}w_{ℓm}}{∂r^{2}}

for r > b. As in Section 1, the hypotheses on the supports of the data and f imply (via Duhamel’s principle) that w_{ℓm} satisfies the outgoing wave condition

\frac{∂w_{ℓm}}{∂ t} (r, t) + \frac{∂w_{ℓm}}{∂ r} (r, t) = 0

for r > b and t > t_{0}.

The initial-value problem for this simple equation can of course be solved explicitly for w_{ℓm}. Knowledge of w_{ℓm} alone, however, does not furnish an effective numerical algorithm, because recovery of u from the w_{ℓm} involves multiple integrations in the radial variable. In [2] we instead determine a single-point one-sided propagation algorithm that yields the coefficients v_{ℓm} explicitly, by incorporating the outgoing wave condition in the propagation formula for the wave equation.

One-sided propagation formula

To derive the single-point formula, we begin by noting that in terms of v_{ℓm} the outgoing wave condition can be written as

\frac{∂}{∂ r} L_{ℓ}^{∗} [r v_{ℓm}] + L_{ℓ}^{∗} [r \frac{∂v_{ℓm}}{∂t} ] = 0

for for r > b and t > t_{0}.

Now, because u_{ℓm}(x, t) ≡ v_{ℓm}(r, t)Y_{ℓm}(θ, ϕ) satisfies the (3 + 1)-dimensional wave equation for |x| > b, we may apply the standard integral propagation formula to propagate values of u_{ℓm}(z, t) from time t_{1} to time t_{2}:

4πu_{ℓm}(z, t_{2}) = ∮ {u_{ℓm}(z + (Δt)ω, t_{1}) + (Δt)[u̇_{ℓm}(z + (Δt)ω, t_{1})
+(ω · ∇u_{ℓm})(z + (Δt)ω, t_{1})]} d^{2}ω .

(Here the integration is over the unit sphere.) Without loss of generality ([2]) we assume that the direction of z is along the north pole of the coordinate system. Then u_{ℓm}(z, t) = 0 for m ≠ 0 because Y_{ℓm}(θ = 0, ϕ) = 0 for m ≠ 0.

Thus u(z, t) = ∑_{ℓ=0}^{∞} u_{ℓ0}(z, t), and so we may consider only the time development

of u_{ℓ0}(x, t) = v_{ℓ0}(r, t)P_{ℓ}(cos θ), where P_{ℓ} is the ℓ^{th}-order Legendre polynomial. Substituting this expression into the propagation formula and integrating by parts, we obtain

2u_{ℓ0}(z, t_{2}) = \frac{s(s − μ)}{τ}P_{ℓ}(μ) v(s) |{s=1+τ \atop s=1−τ}
+ ∫_{1−τ}^{1+τ} {sP_{ℓ}(μ)v̇(s) − τ P′_{ℓ}(μ)v(s)} ds

where v(s) ≡ v_{ℓ0}(as, t_{1}); v̇(s) ≡ a \frac{∂v_{ℓ0}}{∂ t} (as, t_{1}); μ ≡ μ(s) ≡ \frac{s^{2}+1−τ^{2}}{2s}; s ≡ \frac{r}{a}; τ ≡ \frac{Δt}{a}.

It remains to apply the outgoing wave condition to this expression and simplify the result. To do so, we perform the following steps:

(1) Manufacture the expression L_{ℓ}^{∗} [s v̇] from the term s P_{ℓ}(μ)v̇ in the integrand above by determining a function Q_{l}(s) such that P_{ℓ}(μ) = L_{ℓ} [Q_{ℓ}], and integrating by parts to convert the integrand term sP_{ℓ}(μ)v̇ = L_{ℓ} [Q_{ℓ}] sv̇ to Q_{ℓ} L_{ℓ}^{∗} [sv̇].

(2) Apply the outgoing wave condition: Q_{ℓ} L_{ℓ}^{∗} [sv̇] = − Q_{ℓ} \frac{∂}{∂s}L_{ℓ}^{∗} [sv].

(3) Integrate by parts again to convert this expression to s v L_{ℓ} [\frac{∂}{∂s}Q_{ℓ}].

(4) Note that s v L_{ℓ} [\frac{∂}{∂s}Q_{ℓ}] is equal to the opposite of the only other integrand term, − τ P′_{ℓ}(μ)v.

Remarkably, this procedure reduces the propagation formula to surface terms alone. We can furthermore choose the function Q_{l}(s) such that all surface terms vanish at s = 1 + τ. The result is a single-point one-sided propagation formula that expresses u_{ℓ0}(z, t_{2}) in terms of derivatives of v_{ℓ0} at the single point (r, t) = (a − Δt, t_{1}). Specifically, the radial derivatives of v_{ℓ0}(r, t_{1}) to order ℓ, and of v̇_{ℓ0}(r, t_{1}) to order (ℓ − 1), are required only at r = a − Δt.

Similar considerations hold ([2]) for general z = (a, θ, ϕ), and we have finally that

u(z, t_{2}) = ∑_{ℓ=0}^{∞} ∑_{m=−ℓ}^{ℓ} v_{ℓm}(a, t_{2})Y_{ℓm}(θ, ϕ),

where

\begin{aligned}2v_{ℓm}(a, t_{2}) = & (1 − τ) v(1 − τ)+ \\ & + (Q_{ℓ}(s)L_{ℓ}^{∗} [s v] + Γ_{ℓ} [\frac{∂}{∂s}Q_{ℓ}, sv] − Γ_{ℓ} [Q_{ℓ}, sv̇] )|_{s=1−τ} \end{aligned} (3.8)

where Q_{ℓ}(s) ≡ \frac{(−1)^{ℓ}}{2^{ℓ}ℓ!} ((s − τ)^{2} − 1)^{ℓ} and where the boundary-terms operator Γ is given by

Γ_{ℓ} [f, g] (s) ≡ − ∑_{j=1}^{ℓ} \{(−\frac{∂}{∂s} \frac{1}{s})^{ℓ−j} f(s)\} \frac{1}{s} (\frac{1}{s} \frac{∂}{∂s})^{j−1} [s^{ℓ}g(s)].

In this formula, v(s) ≡ v_{ℓm}(as, t_{1}) and v̇(s) ≡ a \frac{∂v_{ℓm}}{∂ t} (as, t_{1}).

To better appreciate formula (3.8), we list below the explicit expressions for the coefficients v_{ℓm} for small values of ℓ:

v_{00}(a, t_{2}) = (1 − τ) v_{00}(a − Δt, t_{1})
v_{1m}(a, t_{2}) = (1 − τ) {(1 + τ)v_{1m}(a − Δt, t_{1})+
+(1 − τ)(Δt) [v̇_{1m}(a − Δt, t_{1}) + v′_{1m}(a − Δt, t_{1})]}

v_{2m}(a, t_{2}) = (1 − τ)^{2} {(1 + 2τ)v_{2m}(a − Δt, t_{1})+

+(Δt) [(1 + 2τ)v̇_{2m}(a − Δt, t_{1}) + (1 + 3τ)v′_{2m}(a − Δt, t_{1})] +

+(1 − τ)(Δt)^{2} [v̇′_{2m}(a − Δt, t_{1}) + v′′_{2m}(a − Δt, t_{1})]}

v_{3m}(x_{2}, t_{2}) = (1 − τ) {(1 + τ)(1 − 5τ^{2})v_{3m}
+(1 − τ)(Δt) [(1 + τ)^{2}v̇_{3m} + (1 + 3τ − τ^{2})v′_{3m}]
+(1 − τ)^{2}(Δt)^{2} [\frac{1}{3}(3 + 10τ)v̇′_{3m} + (1 + 4τ)v′′_{3m}]
+\frac{2}{3}(1 − τ)^{3}(Δt)^{3} [v̇′′_{3m} + v′′′_{3m}]\}

v_{4m}(a, t_{2})

= (1 − τ) {(1 + τ − 9τ^{2} − 9τ^{3} + 6τ^{4})v_{4m}

+(1 − τ)(Δt) [\frac{1}{3}(1 − τ)(3 + 9τ + 8τ^{2})v̇_{4m} + (1 + 3τ − 5τ^{2} − 13τ^{3})v′_{4m}]

+\frac{1}{3}(1 − τ)^{2}(Δt)^{2} [(3 + 10τ + 11τ^{2})v̇′_{4m} + (3 + 12τ + 7τ^{2})v′′_{4m}]
+\frac{1}{3}(1 − τ)^{3}(Δt)^{3} [(2 + 9τ)v̇′′_{4m} + 2(1 + 5τ)v′′′_{4m}]
+\frac{1}{3}(1 − τ)^{4}(Δt)^{4} [v̇′′′_{4m} + v_{4m}^{(4)}]\}

In the last two formulas, we have omitted the arguments for the functions v_{lm} and v̇_{lm}; they are, as always, (a − Δt, t_{1}).

We note that, because \frac{∂u}{∂t}(x, t) also satisfies the wave equation, we may obtain analogous propagation formulas for \frac{∂v_{lm}}{∂t} (a, t_{2}). Thus an approximate numerical algorithm for propagation near a computational domain boundary could be based on this propagation formula with a truncated ℓ summation. Since the determination of v_{lm}(a, t_{2}) is based on (spatial) derivatives of v_{lm}(r, t_{1}) and v̇_{lm}(r, t_{1}) on the sphere r = a − Δt, which is inside the computational domain boundary, it is conceivable that a numerical routine for the interior time development could be devised to maintain sufficient accuracy to allow accurate approximation of these radial derivatives.

4 The Wave Equation in One Dimension, Revisited

The uniqueness result in Theorem 2.1 can be reinterpreted as a statement about uniqueness of solutions to a multiple-time “initial” value problem in which data is presented in Ω^{c} at time t_{0} and in Ω at time t_{1}. For the wave equation

u_{tt} − u_{xx} = 0

in one spatial dimension, we can obtain explicit solutions to certain multipletime initial value problems. To do so, it is convenient to introduce a pictorial representation of d’Alembert’s formula, which we rewrite as

0 = −2u(x, t) + u(x − Δt, t_{0}) + u(x + Δt, t_{0}) + ∫_{x−Δt}^{x+Δt} u_{t}(y, t_{0}) dy

where now Δt ≡ t−t_{0}. In the pictorial representation, a closed dot at the spacetime point (x, t) represents the value u(x, t), an open dot at the spacetime point (x, t) represents the value −u(x, t), and a solid horizontal line between spacetime points (x_{1}, t) and (x_{2}, t) represents the value of the integral ∫_{x_{1}}^{x_{2}} u_{t}(y, t) dy along the segment between the points. Additionally, a constant c adjacent to a closed (open) dot at (x, t) represents cu(x, t) (−cu(x, t)). In the pictures, a dashed line between spacetime points (x, t) and (x ± k, t + k) serves only to indicate lightlike separation between the points and does not represent a value.

With these conventions, we can regard d’Alembert’s formula as saying that the number 0 has the representation shown in Figure 5. By adding and sub-

Figure 5: Representation of zero via d’Alembert’s formula

tracting various (differently-sized and -located) such representations of 0, we can get many formulas involving a solution of the wave equation. Here we will establish one formula that furnishes the solution to a multiple-time initial value problem.

We start with Figure 5 and subtract a smaller version of zero in the lower left corner to obtain Figure 6, which relates values of u at the four points shown and

Figure 6: Representation of zero via d’Alembert’s formula employed twice

the indicated integral of u_{t}. We note that if the data vanishes along the bottom solid line segment then the solution at the apex depends only on the value of the solution at the single location at the intermediate time. This single-point influence is at the heart of the propagation formula in Section 3.2.

We can obtain the one-sided propagation formula of Section 1 by adding another small representation of zero to the “notch” in Figure 6, to obtain Figure 7, from which we can read off the formula

Figure 7: Multiple-time initial-value formula

2u(b, t_{2}) = u(b − (t_{2} − t_{1}), t_{1}) + u(b, t_{1}) + ∫_{b−(t_{2}−t_{1})}^{b} u_{t}(y, t_{1}) dy

−u(b + (t_{1} − t_{0}), t_{0}) + u(b + (t_{2} − t_{0}), t_{0}) + ∫_{b+(t_{1}−t_{0})}^{b+(t_{2}−t_{0})} u_{t}(y, t_{0}) dy

for the solution at time t_{2} in terms of data at times t_{0} and t_{1}. This is an explicit formula for the solution, in terms of data presented at two different times. It shows the precise dependence of the solution of a multiple-time initial-value problem. We invite the reader to play with this new toy to generate many fascinating formulas!

The wave equation in higher space dimensions does not behave this nicely. (Nor does even the one-space-dimension Klein-Gordon equation u_{tt} − u_{xx} + m^{2}u = 0 with m ≠ 0.)

Acknowledgement This work was partially supported by a University of North Texas Faculty Research Grant.

References

[1] Henry A. Warchall, Wave propagation at computational domain boundaries, Commun. Partial Diff. Eq. 16 (1991) 31-41.

[2] Jaime Navarro and Henry A. Warchall, Reflectionless boundary propagation formulas for partial wave solutions to the wave equation, Electronic Journal of Differential Equations, Vol. 1995 (1995) No. 17, 1-14, [http://ejde.math.swt.edu/Volumes/1995/17/abstr.html].

[3] L. Hörmander, Uniqueness theorems and wave front sets for solutions of linear differential equations with analytic coefficients, Comm. Pure Appl. Math. 24 (1971) 671-704.

[4] Semyon V. Tsynkov, Numerical solution of problems on unbounded domains. A review, Applied Numerical Mathematics 27 (1998) 465-532.

[5] Cathleen S. Morawetz, Using Huygens’ principle to avoid computation (1995), unpublished.

[6] Richard Courant and David Hilbert, Methods of Mathematical Physics, New York, Interscience Publishers (1953), Volume 2, Chapter VI, page 576.

[7] Gerald B. Folland, Introduction to Partial Differential Equations, Princeton University Press (1995), Chapter 5, page 174.

[8] Marcus Grote and Joseph Keller, Exact nonreflecting boundary conditions for the time dependent wave equation, SIAM J. Appl. Math. 55 (1995) 280-297.

Henry A. Warchall (e-mail: hankw@unt.edu) Department of Mathematics University of North Texas Denton, TX 76203-1430 and Division of Mathematical Sciences National Science Foundation 4201 Wilson Blvd., Suite 1025 Arlington, VA 22230