A scattering-based algorithm for wave propagation in one dimension

TL;DR

This paper introduces a new explicit numerical scheme for solving the one-dimensional variable coefficient wave equation, offering minimal restrictions on coefficients and initial data for wave propagation simulations.

Contribution

The paper proposes a novel scattering-based explicit algorithm that improves flexibility and applicability in 1D wave propagation modeling.

Findings

Efficient numerical solution for variable coefficient wave equations

Minimal restrictions on initial data and coefficients

Potential for improved accuracy in wave simulations

Abstract

We present an explicit numerical scheme to solve the variable coefficient wave equation in one space dimension with minimal restrictions on the coefficient and initial data.

Equations

$$u_{tt}-\frac{1}{\zeta }(\zeta u_x{)}_x=0$$

$$u(x,0)=f(x),u_t(x,0)=g(x)$$

$$\alpha (x)=\frac{1}{2}\left(f(x)-\frac{1}{\zeta (x)}{\int }_{-\infty }^x\zeta (s)g(s)ds\right),\beta (x)=\frac{1}{2}\left(f(x)+\frac{1}{\zeta (x)}{\int }_{-\infty }^x\zeta (s)g(s)ds\right).$$

$$\zeta^{P}(x)=\left\{\begin{array}[]{ll}\zeta_{-}&\mbox{ if }x<x_{1}^{\ast}\\ \zeta(x_{j})&\mbox{ if }x_{j-1}^{\ast}\leq x<x_{j}^{\ast}\quad(2\leq j\leq n)\\ \zeta_{+}&\mbox{ if }x_{n}^{\ast}\leq x\end{array}\right..$$

$$u_{tt}-\frac{1}{{\zeta }^P}({\zeta }^P{u}_x{)}_x=0$$

$$u(x,0)=f(x),u_t(x,0)=g(x),$$

$$u(x_j^{\ast }-,t)=u(x_j^{\ast }+,t)\text{mbox}and{\zeta }^P(x_j)u_x(x_j^{\ast }-,t)={\zeta }^P(x_{j+1})u_x(x_j^{\ast }+,t).$$

$$v^{k+1}=v^k{A}_P+(v_1^{k+1},0,\dots ,0,v_{2n+2}^{k+1}),$$

$$A_{P}=\left(\begin{array}[]{cccccccccc}0&r_{1}&r_{1}+1&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&r_{2}&r_{2}+1&0&0&0&0&0\\ 0&1-r_{1}&-r_{1}&0&0&0&0&0&0&0\\ 0&0&0&0&0&r_{3}&r_{3}+1&0&0&0\\ 0&0&0&1-r_{2}&-r_{2}&0&0&0&0&0\\ 0&0&0&0&0&0&0&r_{4}&r_{4}+1&0\\ 0&0&0&0&0&1-r_{3}&-r_{3}&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&1-r_{4}&-r_{4}&0\end{array}\right).$$

$$\begin{split}&U_{P}=\\ &\left(\begin{array}[]{cccccccccc}0&r_{1}&\sqrt{1-r_{1}^{2}}&&&&&&&\\ 1&0&0&0&&&&&\mbox{\Large$0$}&\\ 0&0&0&r_{2}&\sqrt{1-r_{2}^{2}}&&&&&\\ &\sqrt{1-r_{1}^{2}}&-r_{1}&0&0&0&&&&\\ &&0&0&0&\ddots&\ddots&&&\\ &&&\sqrt{1-r_{2}^{2}}&-r_{2}&\ddots&0&0&&\\ &&&&\ddots&\ddots&0&r_{n}&\sqrt{1-r_{n}^{2}}&\\ &&&&&\sqrt{1-r_{n-1}^{2}}&-r_{n-1}&0&0&0\\ &\mbox{\Large$0$}&&&&&0&0&0&1\\ &&&&&&&\sqrt{1-r_{n}^{2}}&-r_{n}&0\end{array}\right).\end{split}$$

$${\sigma }_j=\prod _{k=j}^n\sqrt{\frac{1+r_k}{1-r_k}}=\sqrt{\frac{\zeta (x_j)}{\zeta (x_{n+1})}}(1\le j\le n),$$

$$D_P=\mathrm{diag}({\sigma }_1,{\sigma }_1,{\sigma }_2,{\sigma }_2,\dots ,{\sigma }_n,{\sigma }_n,1,1).$$

$$J=\mathrm{diag}(0,1,1,\dots ,1,0).$$

$$A_P=D_P(U_{P}J)D_P^{-1}.$$

$$|\lambda |\Vert v\Vert =\Vert U_{P}Jv\Vert =\Vert Jv\Vert \le \Vert v\Vert .$$

$$\begin{array}{rl}{\alpha }_k(x)& =\frac{1}{2}\left(u^P(x,t_k)-\frac{1}{\zeta (x)}{\int }_{-\infty }^x\zeta (s)u_t^P(s,t_k)ds\right)\text{mbox}and\\ {\beta }_k(x)& =\frac{1}{2}\left(u^P(x,t_k)+\frac{1}{\zeta (x)}{\int }_{-\infty }^x\zeta (s)u_t^P(s,t_k)ds\right).\end{array}$$

$$\begin{array}{rl}{\alpha }_{k+1}(x)& =(1+r_{j-1}){\alpha }_k(x-\Delta )-r_{j-1}{\beta }_k(x_j^{\ast }-(x-x_{j-1}^{\ast }))\text{mbox}and\\ {\beta }_{k+1}(x)& =r_j{\alpha }_k(x_j^{\ast }-(x-x_{j-1}^{\ast }))+(1-r_j){\beta }_k(x+\Delta ).\end{array}$$

$${\alpha }_{k+1}(x_j)$$

$${\beta }_{k+1}(x_j)$$

$$\zeta(x)=\left\{\begin{array}[]{cc}1&\mbox{ if }x\leq 1\\ \rule[-10.0pt]{0.0pt}{28.0pt}2+\tanh\left(\frac{8(x-5)}{16-(x-5)^{2}}\right)&\mbox{ if }1<x<9\\ 3&\mbox{ if }9\leq x\end{array}\right..$$

$$\alpha (x)=2e^{-.05(x+10{)}^2}$$

Full text

MSC 35L05, 65M; Keywords: one dimensional wave equation, numerical methods, scattering

1 The algorithm

A range of physical and biological applications involve the one dimensional wave equation222Note that the equation $u_{tt}-c^{2}u_{zz}=0$ (with $\zeta=1/c$) and the coupled system $\rho u_{t}+p_{z}=0$ and $\frac{1}{K}p_{t}+u_{z}=0$ (with $c=\sqrt{K/\rho}$, $\zeta=\sqrt{K\rho}$) transform to (1.1a) under the change of variables $x=\int_{0}^{z}\frac{1}{c(s)}\,ds$.

[TABLE]

where: there exist $x_{-}<x_{+}$ such that $\zeta$ takes constant values $\zeta_{-}$ and $\zeta_{+}$ on the respective intervals $(-\infty,x_{-})$ and $(x_{+},\infty)$; and $0<c<\zeta<C$ for some $c<C\in\mathbb{R}$. Applications include imaging of layered media such as seismic imaging and the acoustic imaging of laminated structures, microwave imaging of skin tissue, and modelling the human vocal tract and cochlea; see [8],[3], [17], [16], [18]. Certain of these applications have a long history, yet there seems to be no clear consensus on how best to compute solutions to (1.1) in the general setting where $\zeta$ has complicated structure on $[x_{-},x_{+}]$ such as rapid oscillations or discontinuities, or where the initial data is singular (e.g., a Dirac comb), or both (for background see [14, Ch. 10], [1, Ch. 5], [12, Ch. 9]). Recent theoretical insights into the case where $\zeta$ is piecewise constant offer a fresh perspective from which to consider the numerical analysis of (1.1) for general $\zeta$ ([9], [10]).

The purpose of this note is to present a simple, explicit numerical scheme to solve (1.1) on a region $(x,t)\in[a,b]\times[0,T]\subset\mathbb{R}^{2}$ where $[x_{-},x_{+}]\subset[a,b]$. The scheme is based on a natural idea from scattering theory, but with crucial differences from established methods in both the details of the implementation and the interpretation. These differences yield an order of magnitude reduction in computational expense coupled with remarkable versatility.

Algorithm 1**.**

Preliminary step. Convert the initial data $f(x)$ and $g(x)$ from (1.1b) to a pair of functions

[TABLE]

[TABLE]

2 Discussion

Rather than discretize derivatives, the idea behind Algorithm 1 is to replace $\zeta$ with a step function and solve the resulting equation exactly. This echoes the piecewise constant approximations of Goupillaud and Godunov and their many subsequent manifestations (cf. [2], [11], [5], [7], [8, Ch. 3]); however, Algorithm 1 differs markedly from previously established methods. In detail, let $x_{j}^{\ast}=(x_{j}+x_{j+1})/2$ $(1\leq j\leq n)$ denote the midpoints of consecutive spatial grid points, write $P=\{x_{1},x_{1}^{\ast},x_{2},x_{2}^{\ast},\ldots,x_{n},x_{n}^{\ast},x_{n+1}\},$ and set

[TABLE]

Then $\zeta^{P}$ is a step function having equally spaced jump points $x_{1}^{\ast},\ldots,x_{n}^{\ast}$, and it coincides with $\zeta$ on the spatial grid. Replacing $\zeta$ by $\zeta^{P}$ in (1.1) yields

[TABLE]

the first equation of which reduces to $u_{tt}-u_{xx}=0$ on intervals of the form $(-\infty,x_{1}^{\ast})$, $(x_{j-1}^{\ast},x_{j}^{\ast})$ or $(x_{n}^{\ast},\infty)$, where $\zeta^{P}$ is constant. Its solution within each such interval is therefore a pair of left and right travelling waves. One-sided limits of $u$ and $\zeta^{P}u_{x}$ are required to agree at each jump point $x_{j}^{\ast}$

[TABLE]

By a standard computation (see [8, Ch. 3]), this yields reflection and transmission factors: a right-travelling wave in $(x_{j-1}^{\ast},x_{j}^{\ast})$ gets partially reflected back off $x_{j}^{\ast}$, with reflection factor $r_{j}$ as defined in Algorithm 1, and partially transmitted into $(x_{j}^{\ast},x_{j+1}^{\ast})$, with transmission factor $1+r_{j}$. In particular, characteristic lines for the equation (2.2a) bifurcate at jump points of $\zeta^{P}$. The corresponding domains of dependence and influence for a grid point $(x_{j},t_{k})$ are depicted in Figure 1. In Algorithm 1 the pair $v^{k}_{2j-1},v^{k}_{2j}$ encodes the respective amplitudes of the right and left moving parts of the solution to (2.2) at the grid point $(x_{j},t_{k})$, and the iterative step accounts for the reweighted reflected and transmitted waves comprising the domain of dependence of $(x_{j},t_{k+1})$ at time $t_{k}$. In this way the algorithm computes the exact solution to (2.2) at the grid points.

The matrix formulation of the dependence of (row vector) $v^{k+1}$ on $v^{k}$ as expressed in the scattering algorithm is

[TABLE]

where $A_{P}$ is a $(2n+2)\times(2n+2)$ pentadiagonal matrix having $4n$ non-zero entries. For example, if $n=4$,

[TABLE]

Solutions to (2.2) outlined in, e.g., [5] and [8] involve upper or lower triangular systems requiring approximately $Tn^{3}/(b-a)$ flops, while the sparse structure of $A_{P}$ reduces this to $7Tn^{2}/(b-a)$ flops. Thus Algorithm 1 is qualitatively faster than previously established methods, allowing computation on a vastly finer spatial grid.

The following facts justify Algorithm 1 theoretically.

Theorem 1**.**

The eigenvalues of $A_{P}$ lie properly inside the unit circle, implying stability of the scattering algorithm.

Proof.   Let $U_{P}$ denote the $(2n+2)\times(2n+2)$ pentadiagonal matrix

[TABLE]

$U_{P}$ is evidently unitary since its columns (and rows) are pairwise orthogonal. Set

[TABLE]

and let $D_{P}$ denote the $(2n+2)\times(2n+2)$ diagonal matrix

[TABLE]

Let $J$ denote the $(2n+2)\times(2n+2)$ diagonal projection of rank $2n$ obtained by zeroing out the first and last diagonal entries of the identity matrix,

[TABLE]

Then by direct computation,

[TABLE]

Thus the spectrum of $A_{P}$ is the same as that of $U_{P}J$. Note that for any column vector $v$, $\lVert Jv\rVert\leq\lVert v\rVert$, with equality only if $Jv=v$. If $U_{P}Jv=\lambda v$ for some scalar $\lambda$, then, since $U_{P}$ is unitary,

[TABLE]

Hence the spectrum of $U_{P}J$ belongs to the closed unit disk. If $|\lambda|=1$, then $Jv=v$, as noted above, and hence $U_{P}v=\lambda v$. Given that $|\lambda|=1$, this forces $v=0$, as follows. Note that $v_{1}=v_{2n+2}=0$ since $Jv=v$. Since $U_{P}v=\lambda v$, the structure of the second row of $U_{P}$ in (2.5) implies that $v_{1}=\lambda v_{2}$, so $v_{2}=0$. Looking at the first row of (2.5), $r_{1}v_{2}+\sqrt{1-r_{1}^{2}}v_{3}=\lambda v_{1}$, forcing $v_{3}=0$. In general, if $v_{2k}=v_{2k+1}=0$ for some $1\leq k\leq n-1$, then inspection of row $2k+2$ of $U_{P}$, followed by $2k+1$ yields that $v_{2k+2}=v_{2k+3}=0$. Since $v_{2n+2}=0$ this proves $v=0$ as claimed. Thus if $|\lambda|=1$ then $\lambda$ is not an eigenvalue, so the spectrum of $U_{P}J$ and hence that of $A_{P}$, is contained in the interior of the unit disk.

Theorem 2**.**

Let $u^{P}$ denote the solution to (2.2), $(x_{j},t_{k})$ an arbitrary grid point as defined in Algorithm 1, and $\tilde{u}$ the output of the algorithm.

If $\alpha$ and $\beta$ are regular functions (as opposed to distributions), then $u^{P}(x_{j},t_{k})=\tilde{u}^{k}_{j}$. 2. 2.

If $\alpha(x)=\sum_{j\in\mathbb{Z}}c_{j}\delta(x-a-j\Delta)$ and $\beta(x)=\sum_{j\in\mathbb{Z}}d_{j}\delta(x-a-j\Delta)$ are Dirac combs supported on an extension of the spatial grid, and values of $\alpha$ and $\beta$ are replaced in Algorithm 1 by the corresponding coefficients $c_{j}$ and $d_{j}$, then $u^{P}(x_{j},t_{k})=\tilde{u}^{k}_{j}\;\delta(x-a-j\Delta)$.

Proof.   Let $\alpha_{k}$ and $\beta_{k}$ denote the respective right and left-moving components of $u^{P}$ at time $t_{k}$ analogous to (1.2),

[TABLE]

The standard laws of reflection and transmission described above imply that for $x\in(x_{j-1}^{\ast},x_{j}^{\ast})$,

[TABLE]

See Figure 1. These formulas can be given a valid interpretation for distributions (evaluated on appropriately chosen test functions) as well as ordinary functions. Specializing to the point $x=x_{j}\in(x_{j-1}^{\ast},x_{j}^{\ast})$ yields the formulas

[TABLE]

The latter formulas make sense for ordinary functions $\alpha_{k},\beta_{k}$. But distributions do not in general have meaningful restrictions to a single point—unless they are supported at a single point, such as is the case for Dirac functions. Thus formulas (2.7), which are precisely the iterative step of the scattering algorithm, apply both to ordinary functions and Dirac distributions, proving that the algorithm coincides with the exact solution to (1.1) in these cases.

Theorem 3**.**

If there is a sequence of partitions $\{P\}$ such that $\zeta^{P}\rightarrow\zeta$ uniformly, then $u^{P}\rightarrow u$ on $[a,b]\times[0,T]$ in an appropriate sense depending on the smoothness of the initial data.

The convergence $u^{P}\rightarrow u$ is uniform if initial data are smooth; for certain distributional initial data the convergence may be in the Sobolev space $H^{-1}$. Proofs and full technical details will be presented elsewhere. In essence, Theorem 3 follows from abstract functional analysis such as [15, Ch. 3.8] or [4, Ch. 10]. See [13, Theorem 3.5] for a particular version, proved by way of the Hille-Yosida theorem. The heuristic upshot of Theorem 3 is that uniform convergence of the coefficient $u^{P}$ drives convergence of the solution.

Combining Theorems 2 and 3 one can exploit a separation of scales phenomenon to compute both the regular and singular parts of the solution $u(x,t)$ to (1.1) in the case of Dirac function initial data. The solutions to the approximations $u^{P}$ are Dirac combs whose components have the form $\gamma\delta(x-x_{j}-t_{k})$. If $(x_{j},t_{k})$ is a common grid point for a succession of partitions $P$, then the rescaled coefficient $\gamma/\Delta$ diverges if $(x_{j},t_{k})$ belongs to the singular support of $u$; otherwise $\gamma/\Delta$ converges to the regular value of $u(x_{j},t_{k})$.333For technical reasons, one in fact computes the solution at two successive spatial grid points and divides by $2\Delta$. On the other hand, the unscaled coefficient $\gamma$ goes to zero at regular points and stabilizes at singular points. See Figure 2.

Algorithm 1 also handles the contrasting scenario in which initial data is smooth and $\zeta$ is a highly oscillatory step function. If the jump points of $\zeta$ are equally spaced, one can choose $P$ such that $\zeta^{P}=\zeta$, whereby the algorithm is exact up to roundoff error (by Theorem 2). Comparison with the (much slower) exact methods of [9] shows that a good approximation is also obtained if the jump points of $\zeta$ are not equally spaced, provided the underlying grid is sufficiently fine. Indeed, Algorithm 1 reveals a remarkable qualitative phenomenon whereby smooth initial data can lead to discontinuities, as illustrated in Figure 3.

The gist of Theorems 1, 2 and 3 is that Algorithm 1 is guaranteed to approximate the solution to the original equation (1.1) in a wide range of cases, including with distributional initial data. Up to roundoff error, the approximation is exact for piecewise constant $\zeta$, provided its jump points are included in those of $\zeta^{P}$, irrespective of whether the initial data are smooth or singular. The general requirement that $\zeta$ be a uniform limit of step functions (i.e., a regulated function) is very weak, allowing $\zeta$ to have essentially arbitrary structure on $[x_{-},x_{+}]$ (see [6, Ch. VII]). Thus Algorithm 1 differs essentially from existing methods in its interpretation—initial data may be non-smooth or purely distributional, and the coefficient $\zeta$ may have discontinuities.

It does not appear that Algorithm 1 can be generalized to higher dimensions in a way that retains the exactness results of Theorem 2—in this sense the method is purely one dimensional. Within the context of one-dimensional evolutionary systems governed by (1.1) however, the algorithm marks a qualitative improvement over existing methods in terms of both speed and generality. In summary, Algorithm 1 comprises an easily-implemented, versatile and accurate computational tool applicable to various imaging modalities and to modelling of various biological or built structures.

Appendix A Computational data

The present section gives further details on Figures 2 and 3 to facilitate replication of the results. The smooth ramp in Figure 2 is given by the formula

[TABLE]

The source is a right moving unit Dirac function initially centered at $x=-5$.

In the case of Figure 3, $\zeta$ is a step function with 40 randomly generated jump points lying in the interval $[1,10]$. $\zeta(x)=1$ for $x<1$ and $\zeta(x)=2/3$ for $x>10$. The source wave is a right-travelling gaussian initially of the form

[TABLE]

in the first example; the source is shifted 15 units to the right in the second example. The qualitative nature of the propagating wavefield does not depend on the detailed structure of the step function.