\documentclass[twoside]{article} \usepackage{amsfonts, epsf} \pagestyle{myheadings} \markboth{Localization of dependence} { Henry A. Warchall } \begin{document} \setcounter{page}{245} \title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent Mathematical Physics and Quantum Field Theory, \newline Electronic Journal of Differential Equations, Conf. 04, 2000, pp. 245--263\newline http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)} \vspace{\bigskipamount} \\ % Localization of dependence for solutions of hyperbolic differential equations % \thanks{ {\em Mathematics Subject Classifications:} 35B30, 35L05, 35L15, 35L70, 35C10, 35A18. \hfil\break\indent {\em Key words:} Localization of dependence, wave equation, computational domain boundary, \hfil\break\indent exact nonreflecting boundary conditions. \hfil\break\indent \copyright 2000 Southwest Texas State University and University of North Texas. \hfil\break\indent Published November 3, 2000. } } \date{} \author{ Henry A. Warchall \\[12pt] {\em Dedicated to Eyvind H. Wichmann}} \maketitle \begin{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. \end{abstract} \renewcommand\theequation{\thesection.\arabic{equation}} \newtheorem{theorem}{Theorem}[section] \newtheorem{corollary}[theorem]{Corollary} \newcommand{\supp}{\mathop{\rm supp}\nolimits} \newcommand{\R}{{\mathbb R}} \newcommand{\eq}[1]{(\ref{#1})} \newcommand{\be}{\begin{equation}} \newcommand{\ee}{\end{equation}} \newcommand{\proof}{\smallskip\noindent{\bf Proof } } \section{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 \cite{WPCDB} 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 \cite{NW}. 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. 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 $\R^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 $\R^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. \begin{figure}[hbt] \begin{center} \epsfxsize=9cm \epsffile{fig1.eps} \end{center} \caption{Domain of dependence extends outside computational domain} \end{figure} 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: approximately-reflectionless boundary conditions have been developed (see \cite{T} 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? \noindent 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-\Delta t,t_0)+u(x+\Delta t,t_0)+\int_{x-\Delta t}^{x+\Delta t} {\,\,u_t(y,t_0)\,\,dy} $$ where $\Delta t\equiv 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\left\{ {u(\cdot ,t_0),\,\,u_t(\cdot ,t_0)} \right\}\subset (-\infty ,b)$. Suppose $u$ and $u_t$ are known to the left of $x=b$ at time $t_1 > t_0$. Given $x\le 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: \begin{theorem} 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 $\Delta t \equiv t_2 - t_1$. If $x\in (b-\Delta t,b+\Delta t)$ then $$ 2u(x,t_2)=u(x-\Delta t,t_1)+u(b,t_1)+\int_{x-\Delta t}^b {\,\,u_t(y,t_1)\,\,dy}. $$ \end{theorem} \begin{proof} There are two parts to the proof. \noindent First part: We claim that $u_t(x,t_1)=-u_x(x,t_1)$ for $x\ge b$. Since both $u$ and $u_t$ vanish at $t_0$ for $x\ge 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)+\int_{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\ge b$, as claimed. \noindent 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-\Delta t,t_1)+u(x+\Delta t,t_1)+\int_{x-\Delta t}^b {u_t(y,t_1)\,\,dy+} \int_b^{x+\Delta 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+\Delta t,t_1)$, and the assertion of the theorem follows immediately. QED \end{proof} 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-\Delta t,t_1)+u(b,t_1)+\int_{b-\Delta 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. \section{Abstract Uniqueness Results} Quite general results on wave propagation at computational domain boundaries are given in \cite{WPCDB} 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)\subset \R^n$ be the closed ball with center $x$ and radius $r$. Define the closed forward cone $V\equiv \left\{ {(x,t)\in \R^{n+1}\,\, \left| \quad{t\ge \,\,|x|} \right.\,\,} \right\}$. We denote by $\Omega^c$ the complement of the set $\Omega$. \begin{theorem} (Warchall) Let $P(x, D)$ be a linear differential operator on $\R^{n+1}$ with analytic coefficients and principal part $\partial _t^2-\Delta _x$. Let $\Omega \subset \R^n$ be an open convex set. Let $t_00$. To compute the restriction to ${\bar\Omega}$ of the solution to the initial-value problem \be \matrix{{P(x,t;D)u+G(x,t;u)=F(x,t)}\cr {\left. u \right|_{t=t_0}=u_0}\cr} \label{FullIVP} \ee we write the solution $u$ as the sum of the solutions $v^{(0)}$, $v^{(1)}$, and $w^{(1)}$ to the three initial-value problems $$ \matrix{P(x,t;D)v^{(0)}=F(x,t)\cr \left. {v^{(0)}} \right|_{t=t_0}=v_0^{(0)}\cr} $$ and $$ \matrix{P(x,t;D)v^{(1)}=0\cr \left. {v^{(1)}} \right|_{t=t_0}=v_0^{(1)}\cr} $$ and $$ \matrix{P(x,t;D)w^{(1)}+G(x,t;v^{(0)}+v^{(1)}+w^{(1)})=0\cr \left. {w^{(1)}} \right|_{t=t_0}=w_0^{(1)}\cr} $$ where $\supp w_0^{(1)}\subset B$, $\supp v_0^{(1)}\subset{\bar\Omega}$, 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\ge t_0$. The solution $w^{(1)}$ is to be advanced in time with a numerical algorithm. Let $t_1=t_0+\tau$ be the earliest time greater than $t_0$ at which the spatial support of $w^{(1)}$ can intersect $\partial\Omega$. Because $G(x,t;u)$ vanishes for $x$ outside $B$, and because $\supp w_0^{(1)} \subset 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 ${\bar\Omega}$. 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 \eq{FullIVP} is then given by $u=v^{(0)}+v^{(1)}+w^{(1)}$ in ${\bar\Omega}\times [t_0,t_1]$. We next continue the time evolution of $w^{(1)}$ by considering the two initial-value problems $$ \matrix{P(x,t;D)v^{(2)}=0\cr \left. {v^{(2)}} \right|_{t=t_1}=\left. {w^{(1)}} \right|_{t=t_1}\cr} . $$ and $$ \matrix{P(x,t;D)w^{(2)}+G(x,t;v^{(0)}+v^{(1)}+v^{(2)}+w^{(2)})=0\cr \left. {w^{(2)}} \right|_{t=t_1}=0\cr} $$ in terms of which $w^{(1)}=v^{(2)}+w^{(2)}$. The solution $v^{(2)}$ is by hypothesis an explicitly known function of the data $\left. {w^{(1)}} \right|_{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+\tau$, when its spatial support can reach $\partial\Omega$. The solution $u$ to the original initial-value problem \eq{FullIVP} is then given by $u=v^{(0)}+v^{(1)}+v^{(2)}+w^{(2)}$ in ${\bar\Omega}\times [t_1,t_2]$. This latter procedure can be continued arbitrarily many times to advance the solution by time $\tau$ at each step. In general, at the beginning of the $k^{\rm th}$ step, the solution $u$ is known in ${\bar\Omega}\times [t_0,t_{k-1}]$, being given in ${\bar\Omega}\times [t_{k-2},t_{k-1}]$ by $u=v^{(0)}+\cdots+v^{(k-1)}+w^{(k-1)}$. Here each function $v^{(j)}$, $j=1,2,\dots$, is known explicitly for all $t\ge t_{j-1}$, and the function $w^{(k-1)}$ has spatial support in ${\bar\Omega}$ at time $t_{k-1}=t_0+(k-1)\tau$. The $k^{\rm th}$ step of the computation proceeds with the explicit determination of $v^{(k)}$ for $t\ge t_{k-1}$ by the initial-value problem $$ \matrix{P(x,t;D)v^{(k)}=0\cr \left. {v^{(k)}} \right|_{t=t_{k-1}}=\left. {w^{(k-1)}} \right|_{t=t_{k-1}}\cr} $$ followed by numerical solution of $$ \matrix{P(x,t;D)w^{(k)}+G(x,t;v^{(0)}+v^{(1)}+\cdots+v^{(k)}+w^{(k)})=0\cr \left. {w^{(k)}} \right|_{t=t_{k-1}}=0\cr} $$ to determine $w^{(k)}$ in ${\bar\Omega}\times [t_{k-1},t_k]$. The solution $u$ to the original initial-value problem \eq{FullIVP} is given in ${\bar\Omega}\times [t_{k-1},t_k]$ by $u=v^{(0)}+v^{(1)}+\cdots+v^{(k)}+w^{(k)}$. In this way we can integrate this nonlinear equation using numerical computation only in $\bar\Omega$. \subsubsection*{Algorithm based on the wave equation} Morawetz' original result \cite{M} concerns the case where $P$ is the d'Alembertian $\Box\equiv\partial_t^2-\Delta$ in three spatial dimensions. In that case it is possible to do away with the second hypothesis concerning explicit representation of solutions to initial-value 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 $\R^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 \be \matrix{{\Box u+G(x,t;u)=0}\cr {\left. u \right|_{t=0}=u_0}\cr} \label{ivpG} \ee on $\R^{3+1}$ by computing only inside the concentric open ball $\Omega$ 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 $$ \matrix{\Box w^{(1)}+G(x,t;w^{(1)})=0\cr \left. {w^{(1)}} \right|_{t=0}=u_0\cr} $$ in $\Omega\times[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 $\Omega$ for all $t\in[0,t_1]$, by virtue of the unit speed propagation outside $B$. Thus we have $u=w^{(1)}$ in $\Omega\times[t_0,t_1]$. (In terms of the earlier notation, $F\equiv 0$, $t_0=0$, $\tau=2b$, $v_0^{(0)}=v_0^{(1)}=0$, $w_0^{(1)}=u_0$, and $v^{(0)}=v^{(1)}\equiv 0$.) We continue the time evolution of $w^{(1)}$ beyond time $t_1$ by solving the two initial-value problems \be \matrix{\Box v^{(2)}=0\cr \left. {v^{(2)}} \right|_{t=t_1}=\left. {w^{(1)}} \right|_{t=t_1}\cr} \label{waveivp} \ee and \be \matrix{\Box w^{(2)}+G(x,t;v^{(2)}+w^{(2)})=0\cr \left. {w^{(2)}} \right|_{t=t_1}=0\cr} \label{pertivp} \ee in terms of which $w^{(1)}=v^{(2)}+w^{(2)}$ for $t\ge t_1$, as before. Consider first the initial-value problem \eq{waveivp} for $v^{(2)}$. Because the support of the initial data is contained in $\Omega$, 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\equiv \{(0,\,5b)\}+V= \left\{ {(x,t)\in \R^{3+1}\,\,\left| \quad{t\ge \,\,|x|+5b} \right.\,\,} \right\}$. In particular, $v^{(2)}$ vanishes in $B$ for $t\ge 6b$. See Figure 3. \begin{figure}[hbt] \begin{center} \epsfxsize=9cm \epsffile{fig3.eps} \end{center} \caption{Spacetime regions in the piecewise-in-time algorithm based on the wave equation} \end{figure} But we can say more. We claim that $v^{(2)}$ vanishes in the cone $V_3\equiv \{(0,\,3b)\}+V$. To see this, let $\tilde{G}$ be $G$ cut off at time $t_1$: $$ \tilde{G}(x,t;u)=\left\{\matrix{G(x,t;u)&{\rm if}\quad t\le t_1\cr 0&{\rm if}\quad t>t_1}\right. $$ Thus $\tilde{G}$ has spacetime support in $B\times(-\infty,t_1]$. Then the solution to the initial-value problem \eq{ivpG} and the solution to the initial-value problem \be \matrix{{\Box \tilde{u}+\tilde{G}(x,t;\tilde{u})=0}\cr {\left. \tilde{u} \right|_{t=0}=u_0}\cr} \label{ivpGt} \ee are identical in $\Omega\times[0,t_1]$, and $\left. {w^{(1)}} \right|_{t=t_1}$ is the data at time $t_1$ for both problems. Further, the problem \eq{waveivp} gives the solution of \eq{ivpGt} for $t>t_1$, that is, $\tilde{u}=v^{(2)}$ for $t>t_1$. Let $(\xi,\tau)$ be a point in $V_3$. We can make use of Beltrami's formula \cite{CH} or Duhamel's principle \cite{F} to express the solution $\tilde{u}(\xi,\tau)$ in terms of the initial data $u_0$ on the sphere $\left| x-\xi\right|=\tau$ and an integral of $\tilde{G}(x,t;\tilde{u})$ over the truncated backward cone $A\equiv\{(\xi,\tau)\} -\left\{ {(x,t)\in \R^{3+1}\,\,\left| \quad{|x|\le t\le \tau} \right.\,\,} \right\}$. Because $u_0$ is supported in $B$ and thus vanishes on the sphere $\left| x-\xi\right|=\tau$, and because $\supp\tilde{G}$ does not intersect $A$, it follows that $v^{(2)}(\xi,\tau)=\tilde{u}(\xi,\tau)=0$ for all $(\xi,\tau)\in V_3$, as claimed. Thus in particular $v^{(2)}$ vanishes in $B$ for $t\ge t_2\equiv 4b$. Since we are interested only in the values of $v^{(2)}$ in $B$, we thus need only compute $v^{(2)}$ in $B\times[t_1,t_2]$. For this purpose, it more than suffices to apply a numerical algorithm to solve \eq{waveivp} in the truncated backward cone $A_5\equiv\{(0,5b)\}-\left\{ {(x,t)\in \R^{3+1}\,\,\left| \quad{|x|\le t\le 3b} \right.\,\,} \right\}$, which can be done in a straightforward fashion since for time evolution in that region the domain of dependence does not extend outside $\Omega$. 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\times[t_1,\infty)$ using a numerical algorithm and computing only inside $\Omega$. 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 $\Omega$ for all $t\in[t_1,t_2]$, by virtue of the unit speed propagation outside $B$. Then $u=v^{(2)}+w^{(2)}$ in $B\times[t_1,t_2]$. Succeeding steps of the algorithm proceed in the same way. At the beginning of the $k^{\rm th}$ step, the solutions $v^{(1)}$, $v^{(2)}$, $\dots$, $v^{(k-1)}$ have been computed in $B$, and all vanish in $B\times[t_{k-1},\infty]$, where $t_k=2kb$. The solution to \eq{ivpG} is given in $B\times[t_{j-1},t_j]$ by $u=v^{(j)}+w^{(j)}$ for $j=1,\dots,(k-1)$. The solution $w^{(k-1)}$ is supported in $\Omega$ 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 \be \matrix{\Box v^{(k)}=0\cr \left. {v^{(k)}} \right|_{t=t_{k-1}}=\left. {w^{(k-1)}} \right|_{t=t_{k-1}}\cr} \label{waveivpk} \ee and \be \matrix{\Box w^{(k)}+G(x,t;v^{(k)}+w^{(k)})=0\cr \left. {w^{(k)}} \right|_{t=t_{k-1}}=0\cr} \label{pertivpk} \ee in terms of which $w^{(k-1)}=v^{(k)}+w^{(k)}$ for $t\ge t_{k-1}$, in exactly the same fashion as in the second step. Then $u=v^{(k)}+w^{(k)}$ in $B\times[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 $\Omega$. \subsection*{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 {\it radius}. We consider the wave equation $\frac{\partial^2 u}{\partial t^2}-\Delta 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 $\Delta t \equiv t_2 - t_1$. Let $z\in\R^3$ be a point outside $B$, at distance $\Delta t$ or greater from $B$, and set $a\equiv \left| z \right|$. Thus $a>b+\Delta t$. The uniqueness result of Theorem 2.1 implies that $u(z,t_2)$ and ${{\partial u} \over {\partial 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 $\Delta t$ centered at $z$. This intersection is represented by the shaded region in Figure 4. \begin{figure}[hbt] \begin{center} \epsfxsize=9cm \epsffile{fig4.eps} \end{center} \caption{Spacetime regions in the spherical harmonic expansion algorithm} \end{figure} The spherical-harmonic expansion scheme does not quite make this dependence explicit, but instead furnishes an $\ell$-dependent one-sided propagation formula for the $\ell^{\rm 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-\Delta t$ at time $t_1$. This sphere is represented by the heavy circle in Figure 4. In \cite{NW} we build on the idea of Joseph Keller and Marcus Grote \cite{GK} 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. \subsubsection*{Outgoing wave condition} Let $(r,\theta,\phi)$ be the usual polar coordinates for $\R^3$. For $x\in \R^3\backslash B$, we expand the solution $u$ of the wave equation in terms of spherical harmonics: $$ u(x,t)=\sum\limits_{\ell=0}^\infty {\,\,\sum\limits_{m=-\ell}^\ell{\,\,\,u_{\ell m}(x,t)}} $$ where $u_{\ell m}(x,t)\equiv v_{\ell m}(r,t)Y_{\ell m}(\theta ,\phi )$, where $Y_{\ell m}$ is the usual spherical harmonic function. The coefficient function $v_{\ell m}$ of the partial wave solution $u_{\ell m}$ is given by $$ v_{\ell m}(r,t)\equiv \int_0^{2\pi } {\int_0^\pi {\overline {Y_{\ell m}(\theta ,\phi ) \,\,}u(r,\theta ,\phi,t)\,\,\sin \theta \,\,d\theta \,\,d\phi \,\,}} $$ and satisfies the reduced partial differential equation $$ \frac{\partial ^2v_{\ell m}}{\partial t^2}=v''_{\ell m}+{2 \over r}v'_{\ell m}- {{\ell(\ell+1)} \over {r^2}}v_{\ell m} $$ where the prime denotes differentiation with respect to $r$. For notational convenience we define the differential operator $$ L_\ell\equiv r^\ell\left( {-{\partial \over {\partial r}}{1 \over r}} \right)^\ell $$ and its formal adjoint $$ L_l^*\equiv \left( {{1 \over r}{\partial \over {\partial r}}} \right)^\ell \left[{r^\ell\cdot } \right]. $$ Set $$ w_{\ell m}(r,t)\equiv L_\ell ^*\left[ {r\,\,v_{\ell m}} \right](r,t); $$ then it is not difficult to show (\cite{NW}) that $w_{\ell m}$ satisfies the $(1+1)$-dimensional wave equation $$ {{\partial ^2w_{\ell m}} \over {\partial t^2}}={{\partial ^2w_{\ell m}} \over {\partial 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_{\ell m}$ satisfies the outgoing wave condition $$ {{\partial w_{\ell m}} \over {\partial \,t}}(r,t)+{{\partial w_{\ell m}} \over {\partial \,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_{\ell m}$. Knowledge of $w_{\ell m}$ alone, however, does not furnish an effective numerical algorithm, because recovery of $u$ from the $w_{\ell m}$ involves multiple integrations in the radial variable. In \cite{NW} we instead determine a single-point one-sided propagation algorithm that yields the coefficients $v_{\ell m}$ explicitly, by incorporating the outgoing wave condition in the propagation formula for the wave equation. \subsubsection*{One-sided propagation formula} To derive the single-point formula, we begin by noting that in terms of $v_{\ell m}$ the outgoing wave condition can be written as $$ \frac{\partial}{\partial \,r} L_\ell ^*\left[ {r\,\,v_{\ell m}} \right]+ L_\ell ^*\left[ {r\frac{\partial v_{\ell m}}{\partial t}} \right]=0 $$ for for $r>b$ and $t>t_0$. Now, because $u_{\ell m}(x,t)\equiv v_{\ell m}(r,t)Y_{\ell m}(\theta ,\phi )$ satisfies the $(3+1)$-dimensional wave equation for $|x|>b$, we may apply the standard integral propagation formula to propagate values of $u_{\ell m}(z,t)$ from time $t_1$ to time $t_2$: \begin{eqnarray*} 4\pi u_{\ell m}(z,t_2)&=&\oint \Big\{ u_{\ell m}(z+(\Delta t)\omega ,t_1) +(\Delta t)\big[ \dot u_{\ell m}(z+(\Delta t)\omega ,t_1)\\ &&+(\omega \cdot \nabla u_{\ell m})(z+(\Delta t)\omega ,t_1) \big] \Big\}\,d^2\omega \,. \end{eqnarray*} (Here the integration is over the unit sphere.) Without loss of generality (\cite{NW}) we assume that the direction of $z$ is along the north pole of the coordinate system. Then $u_{\ell m}(z,t)=0$ for $m\ne 0$ because $Y_{\ell m}(\theta =0,\phi )=0$ for $m\ne 0$. Thus $u(z,t)=\sum\limits_{\ell =0}^\infty {u_{\ell 0}(z,t)}$, and so we may consider only the time development of $u_{\ell 0}(x,t)=v_{\ell 0}(r,t)P_\ell (\cos \theta )$, where $P_\ell$ is the $\ell^{\rm th}$-order Legendre polynomial. Substituting this expression into the propagation formula and integrating by parts, we obtain \begin{eqnarray*} 2u_{\ell 0}(z,t_2)&=&{{s(s-\mu )} \over \tau }P_\ell (\mu )\,v(s) \left| {_{s=1-\tau }^{s=1+\tau }} \right.\\ &&+\int_{1-\tau }^{1+\tau } {\left\{ {sP_\ell (\mu )\dot v(s)-\tau \,P'_\ell (\mu )v(s)} \right\}\,ds}\\ \end{eqnarray*} where $v(s)\equiv v_{\ell 0}(as,\,\,t_1)$; $\dot v(s)\equiv a\,{{\partial v_{\ell 0}} \over {\partial \,t}}(as,\,\,t_1)$; $\mu \equiv \mu (s)\equiv {{s^2+1-\tau ^2} \over {2s}}$; $s\equiv {r \over a}$; $\tau \equiv {{\Delta t} \over a}$. It remains to apply the outgoing wave condition to this expression and simplify the result. To do so, we perform the following steps: \noindent (1) Manufacture the expression $L_\ell ^*\left[ {s\,\dot v} \right]$ from the term $s\,P_\ell (\mu )\dot v$ in the integrand above by determining a function $Q_l(s)$ such that $P_\ell (\mu )=L_\ell \left[ {Q_\ell } \right]$, and integrating by parts to convert the integrand term $sP_\ell (\mu )\dot v=L_\ell \left[ {Q_\ell } \right]\,\,s\dot v$ to $Q_\ell \,L_\ell ^*\left[ {s\dot v} \right]$. \noindent (2) Apply the outgoing wave condition: $Q_\ell \,L_\ell ^*\left[ {s\dot v} \right]=-\,Q_\ell \,{\partial \over {\partial s}}L_\ell ^*\left[ {sv} \right]$. \noindent (3) Integrate by parts again to convert this expression to $s\,v\,L_\ell \left[ {{\partial \over {\partial s}}Q_\ell } \right]$. \noindent (4) Note that $s\,v\,L_\ell \left[ {{\partial \over {\partial s}}Q_\ell } \right]$ is equal to the opposite of the only other integrand term, $-\,\tau \,P'_\ell (\mu )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+\tau$. The result is a single-point one-sided propagation formula that expresses $u_{\ell 0}(z,t_2)$ in terms of derivatives of $v_{\ell 0}$ at the single point $(r,t)=(a-\Delta t,\,\,t_1)$. Specifically, the radial derivatives of $v_{\ell 0}(r,\,\,t_1)$ to order $\ell$, and of $\dot v_{\ell 0}(r,\,\,t_1)$ to order $(\ell-1)$, are required only at $r=a-\Delta t$. Similar considerations hold (\cite{NW}) for general $z=(a,\theta,\phi)$, and we have finally that $$ u(z,t_2)=\sum\limits_{\ell =0}^\infty{\,\,\sum\limits_{m=-\ell }^\ell{\,\,\,v_{\ell m}(a,t_2)Y_{\ell m}(\theta ,\phi )}}, $$ where \be \begin{array}{rl} 2v_{\ell m}(a,t_2)=&(1-\tau )\,v(1-\tau )+\\ &+\left. {\left( {Q_\ell (s)L_\ell ^*\left[ {s\,v} \right]+\Gamma _\ell \left[ {{\textstyle{\partial \over {\partial s}}}Q_\ell ,\,sv} \right]- \Gamma _\ell \left[ {Q_\ell ,\,s\dot v} \right]\,} \right)} \right|_{s=1-\tau }\\ \end{array} \label{MainPropForm} \ee where $Q_\ell (s)\equiv {{(-1)^\ell } \over {2^\ell \ell !}}\left( {(s-\tau )^2-1} \right)^\ell $ and where the boundary-terms operator $\Gamma$ is given by $$ \Gamma _\ell \left[ {f,g} \right](s)\equiv -\sum\limits_{j=1}^\ell {\left\{ {\left( {-{\partial \over {\partial s}}{1 \over s}} \right)^{\ell -j}f(s)} \right\}\,\,\, {1 \over s}\left( {{1 \over s}{\partial \over {\partial s}}} \right)^{j-1}\left[ {s^\ell g(s)} \right]}. $$ In this formula, $v(s)\equiv v_{\ell m}(as,\,\,t_1)$ and $\dot v(s)\equiv a\, {{\partial v_{\ell m}} \over {\partial \,t}}(as,\,\,t_1)$. To better appreciate formula \eq{MainPropForm}, we list below the explicit expressions for the coefficients $v_{\ell m}$ for small values of $\ell$: $$ v_{00}(a,t_2)=(1-\tau )\,\,v_{00}(a-\Delta t,\,\,t_1) $$ \begin{eqnarray*}v_{1m}(a,t_2)&=&(1-\tau )\,\,\left\{ (1+\tau )v_{1m} (a-\Delta t,\,\,t_1)+ \right. \\ && \left. +(1-\tau )(\Delta t)\left[ \dot v_{1m}(a-\Delta t,\,\,t_1)+ v'_{1m}(a-\Delta t,\,\,t_1) \right] \right\} \end{eqnarray*} \begin{eqnarray*} v_{2m}(a,t_2)&=&(1-\tau )^2\,\, \left\{ (1+2\tau )v_{2m}(a-\Delta t,\,\,t_1)+ \right. \\ &&+(\Delta t)\left[ (1+2\tau )\dot v_{2m}(a-\Delta t,\,\,t_1)+ (1+3\tau )v'_{2m}(a-\Delta t,\,\,t_1) \right]+ \\ && \left. +(1-\tau )(\Delta t)^2 \left[ \dot v'_{2m}(a-\Delta t,\,\,t_1)+ v''_{2m}(a-\Delta t,\,\,t_1) \right] \right\} \end{eqnarray*} \begin{eqnarray*} v_{3m}(x_2,t_2)&=&(1-\tau )\left\{ (1+\tau )(1-5\tau ^2)v_{3m}\right. \\ &&+(1-\tau )(\Delta t)\left[ (1+\tau )^2\dot v_{3m}+ (1+3\tau -\tau ^2)v'_{3m} \right] \\ &&+ (1-\tau )^2(\Delta t)^2 \left[ {\textstyle{1 \over 3}}(3+10\tau )\dot v'_{3m}+(1+4\tau )v''_{3m} \right] \\ && \left. +{\textstyle{2 \over 3}}(1-\tau )^3 (\Delta t)^3 \left[ \dot v''_{3m}+v'''_{3m} \right] \right\} \end{eqnarray*} \begin{eqnarray*} \lefteqn{v_{4m}(a,t_2)}\\ &=&(1-\tau ) \left\{ (1+\tau -9\tau ^2-9\tau ^3+6\tau ^4)v_{4m}\right.\\ &&+(1-\tau )(\Delta t)\left[ {{\textstyle{1 \over 3}}(1-\tau ) (3+9\tau +8\tau ^2)\dot v_{4m}+(1+3\tau -5\tau^2-13\tau ^3)v'_{4m}} \right] \\ &&+{\textstyle{1 \over 3}}(1-\tau )^2(\Delta t)^2 \left[ (3+10\tau +11\tau ^2)\dot v'_{4m}+(3+12\tau +7\tau ^2)v''_{4m} \right] \\ &&+{\textstyle{1 \over 3}}(1-\tau )^3(\Delta t)^3\left[ {(2+9\tau )\dot v''_{4m}+ 2(1+5\tau )v'''_{4m}} \right] \\ &&\left. +{\textstyle{1 \over 3}} (1-\tau )^4(\Delta t)^4 \left[ \dot v'''_{4m}+v_{4m}^{(4)} \right] \right\} \end{eqnarray*} In the last two formulas, we have omitted the arguments for the functions $v_{lm}$ and $\dot v_{lm};$ they are, as always, $(a-\Delta t,\,\,t_1).$ We note that, because ${{\partial u} \over {\partial t}}(x,t)$ also satisfies the wave equation, we may obtain analogous propagation formulas for ${{\partial v_{lm}} \over {\partial 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 $\ell$ summation. Since the determination of $v_{lm}(a,t_2)$ is based on (spatial) derivatives of $v_{lm}(r,t_1)$ and $\dot v_{lm}(r,t_1)$ on the sphere $r=a-\Delta 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. \section{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 $\Omega^c$ at time $t_0$ and in $\Omega$ at time $t_1$. For the wave equation $$u_{tt}-u_{xx}=0$$ in one spatial dimension, we can obtain explicit solutions to certain multiple-time 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-\Delta t,t_0)+u(x+\Delta t,t_0)+\int_{x-\Delta t}^{x+\Delta t} {\,\,u_t(y,t_0)\,\,dy} $$ where now $\Delta t \equiv 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 $\int_{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\pm 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. \begin{figure}[hbt] \begin{center} \epsfxsize=4cm \epsffile{fig5.eps} \end{center} \caption{Representation of zero via d'Alembert's formula} \end{figure} By adding and subtracting 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, \begin{figure}[hbt] \begin{center} \epsfxsize=7cm \epsffile{fig6.eps} \end{center} \caption{Representation of zero via d'Alembert's formula employed twice} \end{figure} which relates values of $u$ at the four points shown and 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, \begin{figure}[hbt] \begin{center} \epsfxsize=7cm \epsffile{fig7.eps} \end{center} \caption{Multiple-time initial-value formula} \end{figure} from which we can read off the formula $$ \begin{array}{rl} 2u(b,t_2)=&u(b-(t_2-t_1),t_1)+u(b,t_1)+\int_{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)\\ &+\int_{b+(t_1-t_0)}^{b+(t_2-t_0)} {\,\,u_t(y,t_0)\,\,dy}\\ \end{array} $$ 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^2u=0$ with $m\ne 0$.) \paragraph{Acknowledgement} This work was partially supported by a University of North Texas Faculty Research Grant. \begin{thebibliography}{10} \bibitem{WPCDB} Henry A. Warchall, {\it Wave propagation at computational domain boundaries}, Commun. Partial Diff. Eq. {\bf 16} (1991) 31-41. \bibitem{NW} Jaime Navarro and Henry A. Warchall, {\it 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]. \bibitem{H} L. H\" ormander, {\it Uniqueness theorems and wave front sets for solutions of linear differential equations with analytic coefficients}, Comm. Pure Appl. Math. {\bf 24} (1971) 671-704. \bibitem{T} Semyon V. Tsynkov, {\it Numerical solution of problems on unbounded domains. A review}, Applied Numerical Mathematics {\bf 27} (1998) 465-532. \bibitem{M} Cathleen S. Morawetz, {\it Using Huygens' principle to avoid computation} (1995), unpublished. \bibitem{CH} Richard Courant and David Hilbert, Methods of Mathematical Physics, New York, Interscience Publishers (1953), Volume 2, Chapter VI, page 576. \bibitem{F} Gerald B. Folland, Introduction to Partial Differential Equations, Princeton University Press (1995), Chapter 5, page 174. \bibitem{GK} Marcus Grote and Joseph Keller, {\it Exact nonreflecting boundary conditions for the time dependent wave equation}, SIAM J. Appl. Math. {\bf 55} (1995) 280-297. \end{thebibliography} \medskip \noindent{\sc 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 \end{document}