\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}