A path integral method for solution of the wave equation with continuously-varying coefficients

TL;DR

This paper introduces a novel path integral method for solving the one-dimensional wave equation with smoothly varying coefficients, expressing solutions as convergent series of integrals involving initial data and PDE coefficients.

Contribution

It develops a new series-based approach for wave equations with variable coefficients, linking classical reflection/transmission concepts to Green's coefficients in smooth regions.

Findings

Series converges with proven bounds on truncation error

Method effectively approximates solutions with examples

Establishes combinatorial connection between coefficients

Abstract

A new method of solution is proposed for solution of the wave equation in one space dimension with continuously-varying coefficients. By considering all paths along which information arrives at a given point, the solution is expressed as an infinite series of integrals, where the integrand involves only the initial data and the PDE coefficients. Each term in the series represents the influence of paths with a fixed number of turning points. We prove that the series converges and provide bounds for the truncation error. The effectiveness of the approximation is illustrated with examples. We illustrate an interesting combinatorial connection between the traditional reflection and transmission coefficients for a sharp interface, and Green's coefficient for transmission through a smoothly-varying region.

Equations

$$u_{tt}$$

$$\begin{array}{rl}{p}_t(x,t)+K(x)u_x(x,t)& =0\\ u_t(x,t)+\frac{1}{\rho (x)}{p}_x(x,t)& =0.\end{array}$$

$$V(x)$$

$$w_t+\Lambda (x)w_x=0.$$

$$p(x,t)=w_1(x+ct,0)+w_2(x-ct,0).$$

$$p_0{v}_2^{+}$$

$$C_T(Z_{-},Z_{+})$$

$$C_R(Z_{-},Z_{+})$$

$$\begin{array}{rl}{w}_t+\Lambda (x)w_x& =(V^{-1}{)}^{\prime }q\\ & =(V^{-1}{)}^{\prime }Vw.\end{array}$$

$$(V^{-1}(x){)}^{\prime }V(x)$$

$$r(x)=\frac{Z^{\prime }(x)}{2Z(x)}$$

$$\frac{R}{\Delta x}\approx \frac{1}{\Delta x}\frac{Z(x+\Delta x)-Z(x)}{Z(x+\Delta x)+Z(x)}\approx \frac{1}{2}\frac{Z^{\prime }(x)}{Z(x)}.$$

$$(K(x),\rho (x))$$

$$\left[\begin{array}{c}p(x,0)\\ u(x,0)\end{array}\right]$$

$$X^{\prime }(t)$$

$$p^{\prime }(X(t))$$

$$p(X)$$

$$\frac{p_{+}}{p_0}={\left(\frac{Z_{+}}{Z_{-}}\right)}^{1/2}=C_G.$$

$$\begin{array}{rl}{p}_t(x,t)+{\mu }_x(x,t)& =0\\ {\mu }_t(x,t)+gh(x)p_x(x,t)& =0,\end{array}$$

$${\int }_0^{X(t/2)}r(x_1)d{x}_1=\frac{1}{2}\log \left(\frac{Z(X(t/2))}{Z_{-}}\right),$$

$$p(0,t)-p(-c_{-}t,0)$$

$$p(x_{+},t)$$

$$T_0(t)$$

$$R_1(t)$$

$$X_j^{\prime }(t)$$

$$x_j$$

$${\mathcal{P}}_n^{[\alpha ,\beta ]}:=\{\mathbf{x}\in [\alpha ,\beta {]}^n:\mathbf{x}\text{ is an alternating sequence.}\}.$$

$$R_{2m+1}(t)$$

$$T_{2m}(t)$$

$$\tau (x_1,x_2,\dots ,x_n)=\tau (\mathbf{x})=2\sum _{j=1}^n(-1{)}^{j+1}{\tau }_j,$$

$${\xi }_R(\mathbf{x},x,t)$$

$${\xi }_T(\mathbf{x},x,t)=-c_{-}(t-\tau (\mathbf{x})-\tau (x)).$$

Full text

1 Introduction and physical setting

We consider the Cauchy problem for the linear one-dimensional wave equation

[TABLE]

which can also be written in first-order form as

[TABLE]

Here we have used the notation of acoustics: $p$ is pressure, $u$ is velocity, $K$ is the bulk modulus, and $\rho$ is the density. Linear wave equations with the same mathematical structure arise in many other applications, with different interpretations of the material parameters, such as elasticity, electromagnetics, and linearized fluid dynamics or water waves. If the coefficients $(\rho(x),K(x))$ are constant or piecewise-constant, the problem may be solved exactly by the method of characteristics. On the other hand, for more general functions $\rho(x)$ and/or $K(x)$ the method of characteristics does not substantially simplify the problem since the solution varies at every point along a characteristic in a way that is coupled to other characteristics.

In this work we propose and demonstrate a method for approximately solving the general Cauchy problem for (2) in the presence of arbitrary variation in $\rho$ and $K$, by grouping characteristic paths according to the number of reflections. Our interest originated in a study of the shoaling of water waves over a continental shelf and our complementary works [7, 8] contain more discussion of this application and several illustrative examples using the linearized shallow water equations, a special case that is also discussed in Remark 1 below. Code to reproduce the numerical experiments in this paper is available online.111https://github.com/ketch/characteristics_rr

We focus on scattering of a localized pulse or front in a medium including a region of continuous variation. In this setting the standard analytical tool, the method of characteristics, leads to a system of ordinary differential equations [15]. In practice, to compute the solution at a single point from this expression, one must discretize the infinite system of ODEs and then solve them numerically. Herein we derive an expression for the solution at any point in terms of just a multidimensional integral. In practice this gives a method that is much more analytical than the classical method of characteristics, since the resulting expression only requires the evaluation of integrals (rather than the solution of an infinite system of ODEs). Our goal is not to provide an efficient computational method for (2), but rather a semi-analytical tool that may provide insight into solutions.

Among the vast literature on solutions of the wave equation in heterogeneous media there are other approaches that bear some relation to ours. In most cases, a piecewise-constant approximation is used; a combinatorial solution for scattering from arbitrary piecewise-constant media was developed in [9, 10]. For an approach similar to ours but in the setting of time-harmonic solutions, see [2, 14, 17].

The method developed in this paper results in an approximation series that has some similarity with the Born scattering series used in seismic imaging [20, 13]. Each term of the Born series arises as a perturbation expansion of the Green’s function solution depending on transmission or reflection at any point, while in our series each term arises solely from paths involving a particular number of turning points in a heterogeneous medium. The goal of forward seismic imaging is often to get a primary approximation of waves that have only suffered one reflection or to remove the effect of multiply reflected waves. So, methods developed to approximate primary waves such as in [12] and multiple removal algorithms like that explored in [13] bear some visual and physical similarities with the integrals developed in this paper.

In the remainder of this section we briefly review the mathematics of characteristics and reflection in one dimension.

1.1 The method of characteristics: homogeneous media

Defining $q=[p,u]^{T}$, the system (2) can be written as $q_{t}+A(x)q_{x}=0$, where $A$ has the eigenvalue decomposition $A=V(x)\Lambda V^{-1}(x)$ with

[TABLE]

Here $Z(x)=\sqrt{K(x)\rho(x)}$ is known as the impedance and $c(x)=\sqrt{K(x)/\rho(x)}$ is the sound speed. If $K(x)$ and $\rho(x)$ are constant (or more generally, if $Z(x)$ is constant) then $V(x)$ is also constant and, setting $w(x,t)=V^{-1}q(x,t)$, (2) can be rewritten as

[TABLE]

System (4) consists of two decoupled advection equations, indicating that one component of the solution ($w_{1}$) travels to the left (with velocity $-c$) while the other ($w_{2}$) travels to the right (with velocity $+c$). Lines of constant $x+ct$ and $x-ct$ are referred to as characteristics. The solution is simply the sum of the components transmitted along the two characteristic families:

[TABLE]

1.2 Piecewise-constant media: reflection and transmission

The method of characteristics can also be used to find the exact solution of (2) if $K(x)$ and $\rho(x)$ are piecewise-constant functions. Within each constant-coefficient domain the characteristic velocities are $\pm c(x)$. Consider a single interface where the impedance jumps from $Z_{-}$ on the left to $Z_{+}$ on the right. Let $v_{1}^{\pm},v_{2}^{\pm}$ denote the respective columns of $V(0^{\pm})$. For an incident right-going wave, the incident ($p_{0}$), transmitted ($p_{T}$), and reflected ($p_{R}$) wave pressures are related by

[TABLE]

Solving system (6) reveals that the transmitted and reflected waves are related to the incident wave by the transmission and reflection coefficients:

[TABLE]

1.3 Smoothly-varying media

Wherever the impedance $Z(x)$ is not constant, the system (2) cannot be decoupled as in (4) because the matrix $V(x)$ that relates $q$ and $w$ varies in space. If $Z(x)$ is differentiable, we have $w_{x}=(V(x)^{-1}q)_{x}=V^{-1}(x)q_{x}+(V^{-1}(x))^{\prime}q$ and we obtain instead of (4) the system

[TABLE]

Here

[TABLE]

Since the left-hand side of (8) is decoupled, we can write (8) as a system of ODEs using a simple coordinate transformation, and the solution of this system can be shown to be that of the original PDE (see e.g. [15, Section 4.1]). Using this approach, to find the solution $p(x,t)$ at a single point still requires solving an infinite number of ODEs. In this work we derive a semi-analytical method that only requires computing an iterated integral for each point value of the solution.

We see from (8) that information is still transmitted along characteristics, but the amplitude of each component is modified by the source terms that couple the characteristic variables through reflection. The coefficient

[TABLE]

gives the amplitude of these reflections and we refer to it as the infinitesimal reflection coefficient.

The infinitesimal reflection ceofficient $r(x)$ is related to the traditional reflection coefficient $R(Z_{-},Z_{+})$; if we take $Z(x)$ to be a continuous function with value $Z_{-}$ at $x$ and value $Z_{+}$ at $x+\Delta x$, the ratio $R/\Delta x$ approaches $r(x)$ as $\Delta x$ tends to zero:

[TABLE]

Characteristics for each of the three classes of media just discussed are illustrated in Figure 1. We see that in the presence of constant or piecewise-constant impedance, the number of characteristics that must be accounted for to compute the solution at a given point is finite. On the other hand, if $Z(x)$ varies continuously then there are in general infinitely many characteristics that influence a given point. The technique developed in the rest of this work is based on the hypothesis that the dominant contributions to the solution come from accounting for paths with relatively few reflections. Here a path is a continuous, piecewise smooth curve in the $x-t$ plane such that each smooth part follows a characteristic and $Z^{\prime}(x)\neq 0$ at each point of non-smoothness. This hypothesis is clearly reasonable when $|r(x)|<1$, since then each reflection must diminish the significance of the corresponding characteristic path. The motivation for this hypothesis more generally is given in Section 4.1.

2 Characteristics in continuously-varying media

In this section we develop an approximate solution to (2) in the form of an infinite series. We focus on the case of a finite region of variation in the spatial coefficients, as illustrated in Figure 2:

[TABLE]

Here $x_{+}$ is the width of the region of varying coefficients, and need not be small. For simplicity we consider the case of a right-going disturbance that is initially confined to $x<0$, and investigate the resulting reflected and transmitted disturbances. Thus

[TABLE]

We assume for simplicity that $Z(x)$ is continuous. Our method and results can be generalized in a natural way to arbitrary initial data and piecewise continuous media.

Outside of the region $[0,x_{+}]$, characteristics are straight lines in the $x$–$t$ plane. Let $X(t)$ denote the characteristic starting from $x=0$ at time zero; i.e., the solution of the initial value ODE

[TABLE]

Here $t_{+}$ is the crossing time so that $X(t_{+})=x_{+}$. It is convenient in what follows to extend $X(t)$ by defining $X(t)=0$ for $t<0$ and $X(t)=x_{+}$ for $t>t_{+}$.

2.1 Amplification or attenuation along characteristics: Green’s Law

In general the pressure is given by $p=w_{1}+w_{2}$; for the case of a pure right-going pulse (15), for which $w_{1}$ is zero, we have $p(x,0)=w_{2}(x,0)=p_{0}$. According to (8), along the path $X(t)$ the value of $w_{2}$ (and hence the value of $p$) satisfies the ODE

[TABLE]

with solution

[TABLE]

In particular, at $x=x_{+}$ we have

[TABLE]

Thus the amplitude of the unreflected part of the wave (for $x\geq x_{+}$) is $C_{G}p_{0}$ for any smoothly varying $Z(x)$, and depends only on the values $Z_{-}$ and $Z_{+}$; it is independent of how $Z$ varies over $[0,x_{+}]$. As we will see, (19) represents the first term in an infinite series that sums to the transmission coefficient $C_{T}$.

Remark 1**.**

We use $C_{G}$ for the quantity defined in (19) since this is the amplification factor given by Green’s law in the context of shoaling, as we discuss in more detail in [8]. The linearized shallow water equations used there can be put in the form (2) by introducing $p(x,t)$ as the depth perturbation of a small amplitude long wave on a background water depth $h(x)$, and $\mu(x,t)$ as the momentum perturbation. Then the linearized shallow water equations can be written in the nonconservative form

[TABLE]

where $g$ is the gravitational constant. This differs from the conservative form used in [8], and has the same form as (2) if we set $K(x)\equiv 1$ and $\rho(x)=1/(gh(x))$. Then the wave speed is $c(x)=\sqrt{gh(x)}$, the impedance is $Z(x)=1/\sqrt{gh(x)}$, and $C_{G}=(h_{-}/h_{+})^{1/4}$. This is the standard form of Green’s law used to estimate the amplification of a shoaling wave as it passes into shallower water, in which case $h_{-}>h_{+}$. Note that this particular application is a special case in that there is only a single variable coefficient $h(x)$, so it is not possible to vary the wave speed and impedance separately.

Both the amplification factor $C_{G}$ and the transmission coefficent $C_{T}$ defined in (7a) are related to the amplitude of transmitted waves. Their differing roles are illustrated in Figure 3, where we consider the propagation of a step function (taking $p_{0}(x)=1$, and with an impedance that grows linearly from $Z_{-}=1$ to $Z_{+}=3$ in the region $[0,1]$). Since $C_{G}$ governs the amplification along characteristics, the leading part of the transmitted wave (which is unaffected by paths with turning points, since they will emerge at later times) has amplitude $C_{G}$. Meanwhile, $C_{T}$ accounts for the cumulative effect of all paths (including those that have turning points) and so the amplitude of the transmitted wave at long times approaches $C_{T}$.

In Figure 4 we consider what happens as $x_{+}$ tends to zero, for fixed values of $Z_{-},Z_{+}$. We again take a step function as the initial condition (plotted as a dashed line). In this case the solution is invariant if $x_{+}$ and $t$ are scaled by the same factor, but in Figure 4 we have plotted solutions for different values of $x_{+}$ all at the same time $t$. We see that as the region $[0,x_{+}]$ shrinks, the width of the transmitted peak becomes increasingly narrow until in the limit $x_{+}=0$ (for which the impedance is discontinuous) the peak is gone and we have a single intermediate state dictated by the transmission coefficient. One way to think about this is that since the discontinuity has infinitesimally small width, the effects of all relevant paths must arrive in infinitesimally short time. The combinatorial relation between the transmission/reflection coefficients and $C_{G}$ is further explored in Section 4.

2.2 Approximating the reflected wave

Let us now turn our attention to the reflected wave in Figure 3. The main contribution to this wave comes from paths with one turning point, as illustrated in Figure 5. The figure on the left shows two paths that emerge at $x=0$ at the same time but started from different initial points and were reflected at different points. It is evident that at any time $t$ there will be such a path reaching $x=0$ that was reflected from point $x$ for each $x\in(0,X(t/2))$, since the path reflected from the rightmost point (the red path in the figure) must have traveled from $x=0$ to the point of reflection in time $t/2$. For the initial condition $p_{0}(x)=1$, the solution along each of these paths has the same initial amplitude. In this case the combined amplitude of these reflected waves is

[TABLE]

where $r(x)$ is defined in (12). Figure 5 shows this diagrammatically. Initially, the reflected wave only contains the contribution of paths reflected near $x=0$. After some time $2t_{+}$, reflections from the whole interval $[0,x_{+}]$ contribute, resulting in a constant asymptotic reflected amplitude.

3 General solution by integrating over all paths

In this section we construct a series for the solution at the boundaries of the variable region:

[TABLE]

Here $R_{n}$ and $T_{n}$ denote contributions from paths involving $n$ reflections. We have effectively computed $T_{0}(t)$ and $R_{1}(t)$ already in the previous sections; from our derivations of (19) and (21) it is straightforward to obtain the more general expressions

[TABLE]

which give the part of the transmitted solution due to paths with no reflections and the part of the reflected solution due to paths with a single reflection, respectively. Here and below, $\tau_{j}$ denotes the time for a characteristic to reach $x_{j}$ from $x=0$.

The function $X(t)$ defined in (16) gives a characteristic; i.e., a path with no points of reflection. More generally, consider a path involving reflection at the sequence of points ${\mathbf{x}}=\{x_{1},x_{2},\dots,x_{n}\}\in[0,x_{+}]$, which we refer to as the reflection point sequence for this path. This path is a union of curves $X_{j}(t)$ ($j=0,1,2,\dots$), each of which is the solution of an initial value problem:

[TABLE]

Here $x_{0}=0$ and $x_{n+1}$ is either zero (for reflected paths) or $x_{+}$ (for transmitted paths). The value of $t_{j}$ is the time at which the path reaches $x_{j}$. For a given medium, a path is determined completely by the reflection points ${\mathbf{x}}$ and the initial time $t_{0}$. Some examples of such paths are given in Figures 6b and 6a. Notice that the shape of the curves $X_{j}$ depends on the variation of $c(x)$, but all can be obtained by applying a temporal offset to $X(t)$ and (for left-going segments) reflecting the curve $X(t)$ vertically in the $x$–$t$ plane.

In keeping with the method of characteristics, we would like to add up the contributions of all paths arriving at a given place and time $(x,t)$. One way to do this is to sum over all valid reflection point sequences. Notice that the reflection point sequence cannot be an arbitrary sequence of points in $[0,x_{+}]$. We need to sum over all paths with an alternating sequence of reflection points, as defined by:

Definition 1**.**

A sequence ${\mathbf{x}}=\{x_{1},x_{2},\cdots x_{n}\}$ is an alternating (down-up) sequence if

[TABLE]

Henceforth we use the term alternating to mean, specifically, down-up sequences. Let

[TABLE]

Then an integral over all paths with $n$ reflection points in $[0,x_{+}]$ is an integral over $\mathcal{P}_{n}^{[0,x_{+}]}$. Thus the terms in (22) are given by the following iterated integrals. Note that $x_{1}$ in the outermost integral can be anywhere in $[0,x+]$; then $x_{2}$ must be chosen in $[0,x_{1}]$, and $x_{3}$ in $[x_{2},x+]$, etc.

[TABLE]

Here $\xi_{R}({\mathbf{x}},x,t)$ is the starting point for the path with reflection points ${\mathbf{x}}$ arriving eventually at $(x,t)$, while $\xi_{T}({\mathbf{x}},x,t)$ is the starting point for the path with reflection points ${\mathbf{x}}$ arriving eventually at $(x,t)$. The factor $(-1)^{m}$ appears because the reflection coefficient for a characteristic initially going left is $-r(x)$, so every even-numbered reflection involves a factor of $-1$. Thus $p_{0}(\xi)$ gives the initial solution value corresponding to a given path, and the product of reflection coefficients gives the part of that value that eventually contributes to the reflected or transmitted wave. The limits of integration take into account that the reflection points must be an alternating sequence.

The factor $C_{G}$ appearing in (24b) is due to variation along a characteristic as described by the solution of (17). It is absent in (24a) because the value of the solution along a characteristic traveling left changes by exactly the reciprocal factor and so there is no net change in amplitude for a path that returns to $x=0$.

We can compute the full solution (to any desired accuracy) by considering the contributions from all paths involving $n=1,2,\dots,N$ reflection points. To complete this approach we only need to determine how $\xi$ depends on ${\mathbf{x}}$ and $t$, which we do in the next two subsections.

3.1 Reflection

Let us work out the initial location $\xi({\mathbf{x}},x,t)$ for a path passing through $x=0$ (going to the right) that is subsequently reflected at the points $x_{1},x_{2},\dots,x_{n}\in(0,x_{+})$ and eventually arrives going to the left at $(x,t)$, for some $x\in[0,x_{+}]$ (see fig. 6b). It is convenient to define

[TABLE]

where again $\tau_{j}$ is the travel time from $x=0$ to $x_{j}$, The travel time for the path is $\tau({\mathbf{x}})-\tau(x)$. Thus this path must have first passed through $x=0$ at time $t-\tau({\mathbf{x}})+\tau(x)$. Hence it must have originated at time zero from

[TABLE]

We can compute the contribution of all paths that are eventually reflected, for any initial condition $p_{0}$, using (24a) with $\xi_{R}({\mathbf{x}},0,t)$ given by (26).

3.2 Transmission

Consider a path starting at $x=0$ (going to the right) that is reflected at the points $x_{1},x_{2},\dots,x_{n}\in(0,x_{+})$ and arrives at $(x,t)$ going to the right, for some $x\in[0,x_{+}]$ (see fig. 6a). The total time to traverse this path is $\tau({\mathbf{x}})+\tau(x)$. Hence it must have originated at time zero from

[TABLE]

Each path contributing to the transmitted wave starts at $x\leq 0$ (with impedance $Z_{-}$) and ends at $x=x_{+}$ (with impedance $Z_{+}$), so the net change in the value of the solution along this path due to Green’s law is given by the factor $C_{G}$ defined in (19). Hence the contribution to the solution is given by $C_{G}r(x_{1})r(x_{2})\cdots r(x_{n})$, leading to the integral (24b) for the total contribution of all paths with $n$ reflections.

We can compute the contribution of all paths that are eventually transmitted, for any initial condition $p_{0}$, using (24b) with $\xi_{T}(x,x_{+},t)$ given by (27). In the next two subsections we specialize the formula above for the important cases of a step (Heaviside) function and a $\delta$-function.

3.3 Scattering of a step function

In this section we apply the approach just outlined to the propagation of an initial condition consisting of a step:

[TABLE]

Since the step function is the integral of a $\delta$-function, the resulting solution gives the integral of the Green’s function for the problem, and can be used as a basis to obtain solutions for arbitrary initial data.

Straightforward calculation shows that the required values of $\xi$ in this case are simply

[TABLE]

for reflected components and

[TABLE]

for transmitted components, where $\tau$ is defined in (25).

The integrals (24a) and (24b) for the reflected and transmitted components can thus be written

[TABLE]

Here the limits of integration $\hat{x}_{2m+1}(t;{\mathbf{x}})$ and $\hat{x}_{2m}(t;{\mathbf{x}})$ impose the condition that the path must reach $x=x_{+}$ by time $t$:

[TABLE]

3.4 Scattering of a delta function

Next we consider initial data consisting of a $\delta$-function:

[TABLE]

Since this is the distributional derivative of the step function (28), we can obtain the solution by differentiating (31b), which yields

[TABLE]

4 Relation between Green’s coefficient and the transmission/reflection coefficients

Let us consider what happens for long times; let $T_{2m}^{\infty}=\lim_{t\to\infty}T_{2m}(t)$. Then $\hat{x}_{2m}(t;{\mathbf{x}})=x_{+}$ (for all ${\mathbf{x}}$) and it is straightforward but tedious to evaluate the multiple integral (31b); the result depends only on $C_{G}$ and $m$. For each value of $m$, $T_{2m}^{\infty}=(-1)^{m}a_{2m}C_{G}(\log(C_{G}))^{2m}$, where the constants $a_{2m}$ for $m=1,2,3,\dots$ are

[TABLE]

We now explain where this sequence comes from.

For $t\geq(n+1)t_{+}$, the integral (31b) for $T_{2m}$ is over all of $\mathcal{P}_{n}$. It can be simplified using the substitution $y(x)=\log(Z(x))/2$. Also let $y_{+}=\log(Z_{+})/2,y_{-}=\log(Z_{-})/2$. For simplicity we assume that $Z(x)$ is monotone increasing. Then

[TABLE]

Let $n=2m$; then (ignoring the sign for the moment) this integral is the volume of some subset of the $n$-dimensional hypercube $[y_{-},y_{+}]^{n}$; namely, the volume of the set $\mathcal{P}^{[y_{-},y_{+}]}_{2m}$ (see (23)). It is does not include the full hypercube because the reflection points are required to be an alternating sequence (this requirement is enforced by the limits of integration). Notice that since $Z(x)$ is monotone increasing, this is equivalent to the condition that the sequence ${y_{1},y_{2},\dots,y_{n}}$ be alternating. The integral in (34) gives the volume of the subset of the hypercube that satisfies this alternating condition. The volume of the whole hypercube is of course $(y_{+}-y_{-})^{n}=(\log(C_{G}))^{n}$.

To determine the value of the integral (34), let us partition the hypercube into $n!$ equal parts, where each part is defined by a particular ordering of the $y_{j}$. For instance, with $n=4$ we would write

[TABLE]

where $(i,j,k,l)$ ranges over all permutations of $(1,2,3,4)$. Each of the sets $V_{ijkl}$ must have the same volume since there is nothing to distinguish a particular coordinate direction. Thus each has volume $(\log(C_{G}))^{n}/n!$. The value of the integral (34) is determined by how many of the $V_{ijkl}$ satisfy the alternating condition. With $n=4$ there are 5 alternating sequences:

[TABLE]

so the integral yields $(5/24)(\log(C_{G}))^{4}$. In general, the number of alternating sequences of length $n$ is known as the $n$th Euler zigzag number (or just zigzag number); for even $n$ these are also known as secant numbers or simply zig numbers [1, 16]. We have proved

Lemma 1**.**

Let $\mathcal{P}_{n}$ be defined by (23). Then

[TABLE]

where $A_{n}$ is the $n$th zigzag number; i.e., the number of alternating permutations of a sequence of length $n$.

An immediate consequence is

Theorem 1**.**

Let $Z(x)$ be monotone and define

[TABLE]

Then the asymptotic contributions for the step are given by

[TABLE]

where $i$ denotes the imaginary unit.

The name zigzag seems eminently appropriate for numbers that appear in the context of Figure 6. Nevertheless, it is worth noting that the original meaning of the name was a reference to zigzags in the discrete setting and had nothing to do with space or paths. There are many recursive formulas for the zigzag numbers; in the course of this work we rediscovered the following formula by evaluating the multiple integrals (31b) [16]. Let $a_{n}=A_{n}/n!$; then the $a_{n}$ are generated by setting $a_{0}=a_{1}=1$ and computing

[TABLE]

These formulas recover the values (33) and the corresponding sequence for the reflection terms. We recall the following combinatorial result due to André [1]:

Theorem 2** (André’s Theorem).**

Let $b_{n}(z)$ be defined by (36). Then

[TABLE]

Comparison of this result with our series (37a)-(37b) leads immediately to

[TABLE]

Further comparing with the expressions for the transmission and reflection coefficients yields

Corollary 1**.**

Let $e^{-\pi}<Z_{+}/Z_{-}<e^{\pi}$. Then

[TABLE]

Proof.

We prove the transmission coefficient part; the proof for the reflection coefficient is similar. From (37a) we have

[TABLE]

This is the Maclaurin series for $\operatorname{sech}(z)$ with $z=\log(C_{G})$; the sequence is convergent for $|z|<\pi/2$, which is equivalent to the condition $e^{-\pi}<Z_{+}/Z_{-}<e^{\pi}$. Meanwhile, we can express the transmission coefficient in terms of the Green’s coefficient as follows:

[TABLE]

Substituting $z=-i\log(C_{G})$ (so $C_{G}=e^{iz}$) we find

[TABLE]

∎

Corollary 1 gives simple expressions for the transmission and reflection coefficients in terms of the Green’s coefficient. It also says that if we add up all the long-time asymptotic contributions from paths with any even number of reflections, we obtain the same value given by the transmission coefficient. Similarly, if we add up all contributions from paths with any odd number of reflections, we obtain the same value as the reflection coefficient. Thus the asymptotic state near $x=0$ for the reflection of the step is just the middle state resulting from the Riemann problem. In fact, corollary 1 could instead be proven directly, using PDE-based arguments to show that the net effect of all terms asymptotically depends only on $Z_{+},Z_{-}$ and so must sum to the traditional transmission and reflection coefficients.

Remark 2**.**

The connection between paths and the zigzag numbers can be seen also in the following way. The path with reflection points $\{x_{1},\dots,x_{n}\}$ can be associated with a tree, essentially as is done in [9, Section 2.3], where each node corresponds to a reflection point. For each admissible ordering of the reflection points, the resulting tree (with nodes labeled in the order in which they are visited) is an increasing 0-1-2 tree, and it can be shown that all increasing 0-1-2 trees of $n$ nodes correspond to some ordering of the reflection points. The number of increasing 0-1-2 trees with $n$ nodes is known to be equal to $A_{n}$ [3].

4.1 Convergence

In this section we consider the convergence of the series (22). The analysis above and the examples in Section 4.2 provide evidence that the series (22) derived in the previous sections approximates the solution of the initial boundary value problem. On the other hand, in the long-time limit, the series may diverge for large impedance ratios. Here we show that the series (22) always converges for any finite time. The main result is Theorem 3; a more specialized result is given in Theorem 4 that has much stronger assumptions but also much smaller error bounds.

Let us extend formulas (24), which approximate solution values at $x=0$ and $x=x_{+}$ respectively, to approximate the solution also for $x\in[0,x_{+}]$. The same reasoning used in the previous sections leads to the more general formulas

[TABLE]

where

[TABLE]

Note that

[TABLE]

Lemma 2**.**

Let $p_{0}(x),Z(x)\in C^{1}$. Let $w^{1}_{2m+1},w^{2}_{2m}$ be defined by (40) for $m\geq 0$ and define $w^{1}_{-1}(x,t)=0$. Then their partial derivatives exist and satisfy

[TABLE]

Proof.

The proof is by direct computation. We illustrate by taking $m=0$. First, note that $C_{G}^{\prime}(x)=r(x)C_{G}(x)$ and $\tau^{\prime}(x)=c(x)$. Thus we have

[TABLE]

so that

[TABLE]

Next we have

[TABLE]

so that

[TABLE]

∎

We remark that this result extends in a natural way to more general initial data by using the theory of distributions.

Let us define formally

[TABLE]

Using Lemma 2 and assuming for the moment that the sums converge, we can formally write

[TABLE]

showing that the limiting functions $\overline{w}^{1}(x,t),\overline{w}^{2}(x,t)$ satisfy (8) and therefore give the solution to the wave equation (2). It remains to show that $\overline{w}^{1},\overline{w}^{2}$ exist. To this end, we require the following Lemma whose proof is deferred to the appendix.

Lemma 3**.**

Consider the set

[TABLE]

Then

[TABLE]

Lemma 4**.**

Consider problem (2) with $C^{1}$ coefficients (14) and initial data (15). Let $w^{1,2}_{n}$ be defined as in (40) and let $t<\infty$ be fixed. Then the sums (42) are uniformly convergent, as are the sums

[TABLE]

Furthermore, we have

[TABLE]

for some $t_{*}\in[0,t]$, where

[TABLE]

Proof.

We bound the magnitude of each term of each series by the product of the volume of integration (using Lemma 3) and the maximum magnitude of the integrand. Notice that in the limits of integration for (40a) we can replace $\mathcal{P}_{n}$ by $\mathcal{P}_{n}(t+\tau(x))$ since for paths outside the latter set we have $p_{0}(\xi_{R}({\mathbf{x}},x,t))=0$. Thus

[TABLE]

Therefore

[TABLE]

Similarly, we obtain

[TABLE]

where $D=\max_{x}|p_{0}^{\prime}(x)|$. The error bounds in the theorem then follow from Taylor’s theorem. ∎

Remark 3**.**

Using (1) without replacing $x_{+}$ in the arguments above leads to estimates that are independent of $t$ but blow up when $\max Z(x)/\min Z(x)$ is too large.

Finally, we obtain

Theorem 3**.**

Consider problem (2) with $C^{1}$ coefficients (14) and initial data (15). Let $w^{1,2}_{n}$ be defined as in (40) and let $t<\infty$ be fixed. Let

[TABLE]

Then $\overline{p}(x,t)$ is the solution of the initial value problem.

Proof.

Lemma 4 shows that $\overline{p}(x,t)$ is well-defined and also (with Lemma 2) that $\overline{w}^{1,2}$ satisfy (8). ∎

The error estimate given in Lemma 4 is typically too large to be useful. As we will see in the examples of Section 4.2, the series often converges much faster. The next theorem gives an example of conditions under which more rapid convergence can be guaranteed.

Theorem 4**.**

Consider problem (2) with coefficients (14) and unit step function initial data (28). Let $R_{n}(t)$ and $T_{n}(t)$ be defined as in (24). Let $Z(x)$ be monotone with $e^{-2\sqrt{2}}<Z_{+}/Z_{-}<e^{2\sqrt{2}}$. Then for any time $0\leq t<\infty$ the following limits exist:

[TABLE]

Furthermore, the terms $|R_{2m+1}(t)|$ and $|T_{2m}(t)|$ decrease monotonically with $m$ and the approximation error can be bounded as follows:

[TABLE]

Note that our assumption on the impedances gives $|\log(Z_{+}/Z_{-})|<{2\sqrt{2}},$ and hence we have $\frac{1}{2}C_{G}^{2}<1$ (and approaching $1/2$ as $Z_{+}/Z_{-}\rightarrow 1$), giving exponentially fast convergence.

Proof.

From (24) we see that if $Z(x)$ is monotone then the series $R_{2m+1}$ and $T_{2m}$ are alternating series (i.e., successive terms in each series have opposite sign). It is sufficient to prove that the terms $|R_{2m+1}(t)|$ and $|T_{2m}(t)|$ decrease monotonically with $m$; then the rest of the theorem follows from standard results for alternating series. We prove convergence of the transmission series $T_{2m}(t)$. The proof for the reflection series is similar. For simplicity, we consider the case in which $Z(x)$ is increasing.

Let $m$ and $t$ be fixed and let $Z(x)$ be as stated in the Theorem. As discussed already $T_{2m}(t)$ is given by integrating over $\mathcal{P}_{2m}(t+t_{+})$. For clarity, in the remainder of the proof we write $\mathcal{P}_{2m}$ with no argument; it is implicitly $t+t_{+}$.

[TABLE]

The second equality holds because, since $Z(x)$ is monotone, the integrand has the same sign for all paths. This also means that if $\mathcal{P}_{2m}$ is replaced by a larger set of paths, the resulting integral provides an upper bound on $|T_{2m}(t)|$.

Notice that every path in $\mathcal{P}_{2m+2}$ can be obtained in exactly one way by taking a particular path in $\mathcal{P}_{2m}$ and appending two (admissible) reflection points $x_{2m+1},x_{2m+2}$. Admissibility of the resulting path involves a restriction in the total path length (travel time $\tau({\mathbf{x}})\leq t$) and the condition that $x_{2m+1}\geq\max(x_{2m},x_{2m+2})$. Let us consider the larger set $\widehat{\mathcal{P}}_{2m+2}$ obtained by omitting the path length restriction and requiring only that $x_{2m+1}\geq x_{2m+2}$. In other words, $\widehat{\mathcal{P}}_{2m+2}$ is obtained by appending, for each path in $\mathcal{P}_{2m}$, all pairs $(x_{2m+1},x_{2m+2})$ such that $0\leq x_{2m+1}\leq x_{2m+2}\leq x_{+}$. Clearly $\mathcal{P}_{2m+2}\subset\widehat{\mathcal{P}}_{2m+2}$, so we have

[TABLE]

Since $|\log(Z_{+}/Z_{-})|<{2\sqrt{2}},$ we have $\frac{1}{2}C_{G}^{2}<1$, so $|T_{2m+2}(t)|<|T_{2m}(t)|$, so the alternating series is convergent. ∎

4.2 Examples

In this section we illustrate, through numerical examples, the method just proposed. For comparison, we compute reference solutions using the finite volume solver Clawpack [4]. Code for reproducing these results is available online.222https://github.com/ketch/characteristics_rr In each case the reference solution is computed with a discretization sufficiently fine so that further refinement produces no visible change in the solution.

We take $x_{+}=1$ in all examples. In the first three examples we take the functions $c(x),Z(x)$ to be linear in the interval $(0,1)$:

[TABLE]

Let $s=c_{+}-c_{-}$. Then a right-going characteristic starting from $x=0$ at $t=0$ satisfies the ODE

[TABLE]

with solution

[TABLE]

The total time to cross from $x=0$ to $x=x_{+}$ is thus

[TABLE]

For each example, we show the solution corresponding to an initial step function ($p_{0}(x)=1$ for all $x<0$) and a Dirac $\delta$-function ($p_{0}(x)=\delta(x)$). For the $\delta$-function examples, the $\delta$-function part of the transmitted wave is represented by a larger red circle that also indicates the mass of the transmitted $\delta$-function.

A first example, with very mild variation in $Z$, is shown in fig. 7. The solution involving only terms up to $T_{2}$ is already highly accurate. In the second example, shown in fig. 8, $Z$ varies by a factor of 8. In this case it can be seen that the approximation using terms up to $T_{4}$ gives a significant improvement.

Both of the previous examples satisfy the conditions given in Theorem 4. The next two examples do not. In the third example, we take $Z_{-}=1$ and $Z_{+}=20$. It can be seen that in this case the convergence for large times is much slower and the series including terms up to $T_{4}$ is a good approximation only for short times.

In the final example, $Z(x)$ is non-monotone:

[TABLE]

The solution given by including terms up to $T_{4}$ captures the oscillating solution well. This example also illustrates that when $Z(x)$ is a non-monotone function, the transmitted wave amplitude can exceed $C_{G}$ at some points.

Remark 4**.**

Although we have focused on media with continuous coefficients $Z(x),c(x)$, it is possible to extend this approach to piece-wise continuous media by incorporating the effect of reflection and transmission at points of discontinuity. fig. 11 shows the approximation $R_{1}(t)$ for an example with $x_{+}=1$ and a single discontinuity at $x=1/2$. Let $Z^{\pm}_{1/2}$ denote the impedance just to the left and right of $x=1/2$, and assume that $Z(x)$ varies continuously over each interval $[0,1/2]$ and $[1/2,1]$. Then the transmission and reflection coefficients are given by $C_{T}(Z^{-}_{1/2},Z^{+}_{1/2})$, $C_{R}(Z^{-}_{1/2},Z^{+}_{1/2})$ for waves incident from the left and by $C_{T}(Z^{+}_{1/2},Z^{-}_{1/2})$, $C_{R}(Z^{+}_{1/2},Z^{-}_{1/2})$ for waves incident from the right (see (7)). We can write $R_{1}(t)=R_{1}^{\textup{cont}}(t)+R_{1}^{\textup{discont}}(t)$ where the contribution from the discontinuity is given by the reflection coefficient. The continuous contribution, $R_{1}^{\textup{cont}}$, is obtained by integrating over the two continuous portions. For a step function initial condition, this is

[TABLE]

which can also be written as

[TABLE]

The term with the two transmission coefficient factors accounts for the paths that pass through the discontinuity (once in each direction). The number of separate integrals that must be evaluated increases for higher-order terms and for media with more discontinuities.

5 Conclusions

We have developed a new approximation to the solution of the wave equation in one space dimension in the presence of a region of continuously-varying coefficients, by accounting for all paths along which information can reach a given point. This can be extended in a straightforward way to other linear hyperbolic systems in one dimension. This approximation gives an explicit expression for the solution as an infinite sum of integrals depending only on the initial data and the material properties. We have shown that the series converges to the solution of the wave equation. We have also demonstrated that it can give theoretical insight, by studying the propagation of a Heaviside function and elucidating the relation between transmission and reflection coefficients expected in the limiting case of a sharp interface, and the Green’s law behavior expected for sufficiently smooth transitions in material properties.

It is natural to ask how the approach described in this paper compares, as a computational tool, to traditional numerical PDE discretizations like finite difference or finite volume methods. It is difficult to give a meaningful answer to this question, because the approaches are fundamentally different:

• •

Numerical methods begin with discretization: the medium is approximated, generally in a piecewise-constant manner and the solution is approximated by a representation in some finite basis. The fundamental approximation in our approach is instead truncation of the series (22); discretization is eventually required for numerical evaluation of integrals, but this can be done to machine precision if desired.

• •

Numerical discretizations require computation of the solution at a large number of points (in $(x,t)$, and this number must be increased in order to obtain higher accuracy. In our approach the solution can be computed at a single point to any desired accuracy without computing the solution at other points.

• •

The importance of various factors influencing the size of the error are very different in numerical discretizations versus our approach. For instance, numerical discretizations have difficulty in accurately capturing narrow peaks such as the leading part of the transmitted or reflected waves above. In order to capture these, we had to use especially fine grids in Clawpack. But the path integral method is most accurate at these points; in fact, the first term in the infinite series already gives the exact solution.

Because of these differences, it is easy to construct situations in which one approach or the other is vastly more efficient. For instance, the method described here can be more efficient if the solution is needed only at one or a few points and if the initial data is not smooth. On the other hand, if the ratio $\max Z(x)/\min Z(x)$ is large and/or solution values are needed at very many points, the approach described here may be much more costly than traditional numerical discretizations. We have not investigated techniques for reducing the computational cost or made any detailed comparisons.

It is natural to expect that the series (22) may converge because paths involving many reflections contribute in successively smaller amounts to the solution. Examining (24), this viewpoint makes sense only if $|r(x)|<1$. However, our examples and analysis show that (22) converges quite independently of any such condition. Theorem 4 indicates that in general (22) converges for a completely different reason: the number of contributing paths (more precisely, the volume they occupy in an appropriate space) becomes vanishingly small as $n\to\infty$.

Acknowledgments

We are grateful to Ernst Hairer for a comment that led us to the connection with zigzag numbers, and to Lajos Lóczi for reviewing an early draft of this work. We also thank an anonymous referee for very helpful comments and suggestions.

Appendix A Proof of Lemma 3

Proof.

First, for simplicity take $c(x)=1$ so that the travel time between two points is just the distance between them. Because $c(x)=1$, in this case the set $\mathcal{P}_{n}(t)$ (defined in (43)) is just

[TABLE]

the set of alternating sequences with path length at most $t$. Define the mapping $f:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ by

[TABLE]

This mapping can be represented by a lower-triangular matrix whose diagonal entries are $\pm 1$, so it preserves volume. Note also that

[TABLE]

and for any alternating sequence ${\mathbf{x}}\geq 0$ we have $f({\mathbf{x}})\geq 0$. Let ${\mathcal{B}}_{n+}^{1}(t)$ denote the intersection of the $n$-dimensional $L_{1}$ ball of radius $t$ with the positive orthant:

[TABLE]

For any ${\mathbf{x}}\in{\mathcal{X}}_{n}(t)$, we have $f({\mathbf{x}})\in{\mathcal{B}}_{n+}^{1}(t)$, so

[TABLE]

The value of the last integral is a classical result due to Dirichlet [6, p. 168].

To extend the proof to arbitrary $c(x)$, let $C=\max_{x}|c(x)|$. Then the length of a path emerging at time $t$ is no greater than $Ct$, so $\mathcal{P}_{n}(t)\subset{\mathcal{X}}_{n}(Ct)$. Thus

[TABLE]

∎