Dimensional splitting of hyperbolic partial differential equations using the Radon transform

TL;DR

This paper presents a novel dimensional splitting method for hyperbolic PDEs using the Radon transform, enabling advanced multi-dimensional techniques like large time-step methods and boundary conditions.

Contribution

It introduces a new splitting approach based on the Radon transform's properties, expanding the toolkit for solving multi-dimensional hyperbolic PDEs.

Findings

Enables large time-step extensions for hyperbolic PDEs

Facilitates implementation of absorbing boundary conditions

Supports multi-dimensional transport reversal

Abstract

We introduce a dimensional splitting method based on the intertwining property of the Radon transform, with a particular focus on its applications related to hyperbolic partial differential equations (PDEs). This dimensional splitting has remarkable properties that makes it useful in a variety of contexts, including multi-dimensional extension of large time-step (LTS) methods, absorbing boundary conditions, displacement interpolation, and multi-dimensional generalization of transport reversal.

Figures

Equations

$$q_t+A{q}_x=0$$

$$\hat{\phi }(\omega ,s)=\mathcal{R}\phi (\omega ,s)={\int }_{x\cdot \omega =s}\phi (x)\text{d}m(x),$$

$$\check{\psi }(x)={\mathcal{R}}^{\#}\psi (x)={\int }_{S^{n-1}}\psi (\omega ,\omega \cdot x)\text{d}S(\omega ),$$

$${\mathcal{S}}_H({\mathbb{P}}^n)=\{\begin{array}{ll}& \text{for each }k\in {\mathbb{Z}}^{+},{\int }_{\mathbb{R}}F(\omega ,p)p^k\mathrm{d}p\text{ is}\\ F\in \mathcal{S}({\mathbb{P}}^n):& \text{a homogeneous polynomial in }{\omega }_1,...,{\omega }_n\\ & \text{of degree }k\end{array}\}.$$

$${\left(\frac{\partial }{\partial x_i}\phi (x)\right)}^{\wedge }={\omega }_i\frac{\partial }{\partial s}\hat{\phi }(\omega ,s).$$

$$q_t+\theta \cdot \nabla q=0\text{ where }\theta \in S^1.$$

$${\hat{q}}_t+(\theta \cdot \omega ){\hat{q}}_s=0,$$

$${\left[\begin{array}{c}p\\ u\\ v\end{array}\right]}_t+\left[\begin{array}{ccc}0& K_0& 0\\ 1/{\rho }_0& 0& 0\\ 0& 0& 0\end{array}\right]{\left[\begin{array}{c}p\\ u\\ v\end{array}\right]}_{x_1}+\left[\begin{array}{ccc}0& 0& K_0\\ 0& 0& 0\\ 1/{\rho }_0& 0& 0\end{array}\right]{\left[\begin{array}{c}p\\ u\\ v\end{array}\right]}_{x_2}=0.$$

$${\left[\begin{array}{c}\hat{p}\\ \hat{u}\\ \hat{v}\end{array}\right]}_t+\left[\begin{array}{ccc}0& {\omega }_1{K}_0& {\omega }_2{K}_0\\ {\omega }_1/{\rho }_0& 0& 0\\ {\omega }_2/{\rho }_0& 0& 0\end{array}\right]{\left[\begin{array}{c}\hat{p}\\ \hat{u}\\ \hat{v}\end{array}\right]}_s=0.$$

$${\left[\begin{array}{c}\hat{p}\\ \hat{\mu }\end{array}\right]}_t+\left[\begin{array}{cc}0& K_0\\ 1/{\rho }_0& 0\end{array}\right]{\left[\begin{array}{c}\hat{p}\\ \hat{\mu }\end{array}\right]}_s=0.$$

$$\hat{q}(T,\omega ,s)=\mathcal{K}(T;\omega ){\hat{q}}_0(\omega ,s).$$

$$c_{n}q(T,x)=\{\begin{array}{rlr}& {\mathcal{R}}^{\#}\frac{d^{n-1}}{d{s}^{n-1}}\hat{q}(T,\omega ,s)& \text{ if }n\text{ is odd,}\\ & {\mathcal{R}}^{\#}{H}_s\frac{d^{n-1}}{d{s}^{n-1}}\hat{q}(T,\omega ,s)& \text{ if }n\text{ even,}\end{array}$$

$$q_0(x)$$

$$p_0(x)$$

$$\hat{q}(t,\omega ,s)=\frac{1}{2}\left(r_1{\hat{p}}_0(\omega ,s-t)+r_2{\hat{p}}_0(\omega ,s+t)\right)\text{ where }{r}_1=\left[\begin{array}{c}1\\ {\omega }_1\\ {\omega }_2\end{array}\right]\text{ and }{r}_2=\left[\begin{array}{c}1\\ -{\omega }_1\\ -{\omega }_2\end{array}\right].$$

$${\left({\int }_0^{4\sqrt{2}}|p(\rho ,t)-p_{\text{ref}}(\rho ,t){|}^{\mathfrak{p}}\rho \text{d}\rho \right)}^{\frac{1}{\mathfrak{p}}}$$

$$\text{ where }\rho =\sqrt{x_1^2+x_2^2},p(\rho ,t)=p\left(\frac{\rho }{\sqrt{2}},\frac{\rho }{\sqrt{2}},t\right),\mathfrak{p}=1,2.$$

$$q_0(x)$$

$$p_1(x)$$

$$q_t+f(q{)}_{x_1}+g(q{)}_{x_2}=0,$$

$${\hat{q}}_t+{\left[{\omega }_{1}f(q)+{\omega }_{2}g(q)\right]}_s^{\wedge }=0.$$

$$h(q)={\omega }_{1}f(q)+{\omega }_{2}g(q).$$

$${\hat{q}}_t+h(q{)}_s^{\wedge }=0.$$

$${\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\end{array}\right]}_t+{\left[\begin{array}{c}\rho u\\ \rho u^2+\frac{1}{2}\bar{g}{\rho }^2\\ \rho uv\end{array}\right]}_{x_1}+{\left[\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+\frac{1}{2}\bar{g}{\rho }^2\end{array}\right]}_{x_2}=0.$$

$${\left[\begin{array}{c}\rho \\ \rho \mu \\ \rho \nu \end{array}\right]}_t^{\wedge }+{\left[\begin{array}{c}\rho \mu \\ \rho {\mu }^2+\frac{1}{2}\bar{g}{\rho }^2\\ \rho \mu \nu \end{array}\right]}_s^{\wedge }=0.$$

$${\hat{Q}}_{i,j}^n\approx {\int }_{{\xi }_{i,j}}q(t_n,x)\text{d}m(x).$$

$${\hat{Q}}_{i_0,j}^{n+1}={\hat{Q}}_{i_0,j}^n-\Delta t(F_{i_0,j+\frac{1}{2}}^n-F_{i_0,j-\frac{1}{2}}^n).$$

$$D_{2n}(h,2s)$$

$$D_{2n}(h,2s+1)$$

$${\left({\mathcal{R}}_N^{a}A\right)}_{h,s}$$

$${\left({\mathcal{R}}_N^{b}A\right)}_{h,s}$$

$${\left({\mathcal{R}}_N^{c}A\right)}_{h,s}$$

$${\left({\mathcal{R}}_N^{d}A\right)}_{h,s}$$

$${\mathcal{R}}_{N}A=({\mathcal{R}}_N^{a}A,{\mathcal{R}}_N^{b}A,{\mathcal{R}}_N^{c}A,{\mathcal{R}}_N^{d}A).$$

Full text

1 Introduction

Dimensional splitting provides the simplest approach to obtaining a multi-dimensional method from a one-dimensional method [28, 37, 14, 25]. Although extremely powerful, existing splitting methods do not preserve a special feature that is easily obtained for 1D methods. For 1D hyperbolic partial differential equations (PDEs) of the type

[TABLE]

where $A$ is a constant diagonalizable matrix with real and distinct eigenvalues, one can devise large time-step (LTS) methods that allow the solution to be solved up to any time without incurring excessive numerical diffusion [21, 23, 22]. Previous splitting methods do not lead to such LTS methods in multi-dimensions.

In this paper, we introduce a dimensional splitting method that allows multi-dimensional linear constant coefficient hyperbolic problems to be solved up to desired time. The method relies on the intertwining property of Radon transforms [18, 29], thereby transforming a multi-dimensional problem into a family of one-dimensional ones. Simply by applying an 1D LTS method on each of these one-dimensional problems, one obtains a multi-dimensional LTS method. While this intertwining property is well-known and is utilized to analyze PDEs in standard references [20], it has not been used for constructing multi-dimensional numerical methods, to the best of our knowledge.

The method also has implications for the problem of imposing absorbing boundary conditions, a problem that has received sustained interest over many decades [15, 3, 11, 4]. By using the Radon transform, the splitting decomposes multi-dimensional waves into planar ones, thereby allowing a separate treatment of each incident planar wave near the boundary. This yields the desired absorbing boundary conditions in odd dimensions, and in even dimensions one obtains an approximation up to $\mathcal{O}(1/t)$ that does not cause spurious reflections.

Another useful application of this dimensional splitting is in displacement interpolation, a concept that arises naturally in optimal transport [40]. Our interest in displacement interpolation is motivated by model reduction. To construct reduced order models for typical hyperbolic problems, one cannot rely solely on linear subspaces [1, 10], and it is necessary to interpolate over the Lagragian action [36, 35, 32, 34]. In a single spatial dimension this can be done in a relatively straightforward manner, owing to the LTS methods available for 1D [34]. The multi-dimensional LTS method is useful also for the multi-dimensional extension of displacement interpolation, and this in turn will yield a straightforward way for low-dimensional information to be extracted for multi-dimensional hyperbolic problems.

For the dimensional splitting to be computationally successful, one requires an algorithm for computing the Radon transform and its inverse efficiently. Throughout this paper we use the approximate discrete Radon transform (ADRT), also called simply the discrete Radon transform (DRT), devised in [8, 17]. We will refer to ADRT as DRT. It is a fast algorithm with the computational cost of $\mathcal{O}(N^{2}\log N)$ for an $N\times N$ image or grid111The term grid (cell) is a more appropriate term for our PDE applications, but DRT originally comes from imaging literature so we will sometimes also use the term image (pixel), interchangeably., and the efficiency is obtained through a geometric recursion of so-called digital lines. The inversion algorithm using the full multi-grid method appeared in [30], but here we adopt a simpler approach by making use of the conjugate gradient algorithm [16] for the inversion.

This paper is organized as follows. In Section 2, we give a review of the intertwining property of the Radon transform and introduce the dimensional splitting method. In Section 3, we give a brief introduction to the DRT algorithm and discuss its inversion. In Section 4, we discuss its applications in absorbing boundary and in displacement interpolation. In this paper we will fully implement only constant coefficient linear problems in spatial dimension two, although we will also discuss how the splitting can be extended to fully nonlinear problems and to higher spatial dimensions. Further investigations into these and other related topics will be mentioned in Section 5.

The Radon transform was introduced by Johann Radon [31] and has been a major subject of study, primarily due to its use in medical imaging but also as a general mathematical and computational tool.

2 Dimensional splitting using the Radon transform

In this section, we briefly review the intertwining property of the Radon transform [18, 29] then show that it can be used as a dimensional splitting tool that extends the large time-step (LTS) operator to multiple spatial dimensions. It preserves the ability to take large time-steps without loss of accuracy in the constant coefficient case. Moreover, this splitting can potentially be used for fully nonlinear problems as well, in a similar manner to the other splitting methods, with the usual CFL condition for the time-step.

2.1 Intertwining property of Radon transforms

The Radon transform $\hat{\varphi}:S^{n-1}\times\mathbb{R}\to\mathbb{R}$ of the function $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ is defined as

[TABLE]

in which $\mathrm{d}m$ is the Euclidean measure over the hyperplane. For any fixed pair $(\omega,s)\in S^{n-1}\times\mathbb{R}$, the set $\{x=(x_{1},x_{2},...,x_{n})\in\mathbb{R}^{n}:x\cdot\omega=s\}$ defines a hyperplane, so the transform is simply an integration of the function over this hyperplane. We will denote the space of hyperplanes parametrized above by $\mathbb{P}^{n}$. In effect, $\hat{\varphi}$ decomposes $\varphi$ into planar waves in the direction of $\omega$.

The back-projection is defined as the dual of $\mathcal{R}$ with respect to the obvious inner product over $S^{n-1}\times\mathbb{R}$. For $\psi:S^{n-1}\times\mathbb{R}\to\mathbb{R}$ the back-projection $\check{\psi}$ is

[TABLE]

where $\textrm{d}S$ is the measure on $S^{n-1}$.

The Radon transform $\mathcal{R}$ is a linear one-to-one map between $\mathcal{S}(\mathbb{R}^{n})$ and $\mathcal{S}_{H}(\mathbb{P}^{n})$ [18] in which $\mathcal{S}(\mathbb{R}^{n})$ denotes the Schwartz class and

[TABLE]

The correspondence can be naturally extended to distributions, and we refer the reader to standard references for further details.

The Radon transform has a remarkable property, that it intertwines a partial derivative with a univariate derivative. The $i$-th partial derivative $\partial/\partial x_{i}$ of $\varphi$ is now transformed to the derivative of $\hat{\varphi}$ with respect to $s$ multiplied by $\omega_{i}$,

[TABLE]

This is the key property that allows us to transform a multi-dimensional hyperbolic problem into a collection of one-dimensional problems. For example, let us apply the Radon transform to the transport equation in $\mathbb{R}^{2}$, in which the scalar state variable $q:\mathbb{R}^{+}\times\mathbb{R}^{2}\to\mathbb{R}$ satisfies,

[TABLE]

The transformation produces a family of 1D advection equations

[TABLE]

whose coefficient varies for each $\omega$. Similarly, consider the acoustic equations for $p,u,v:\mathbb{R}^{+}\times\mathbb{R}^{2}\to\mathbb{R}$, where the state variable $p$ denotes pressure, $u$ the velocity in $x_{1}$-direction, $v$ the velocity in $x_{2}$-direction,

[TABLE]

After the transform, we obtain

[TABLE]

This PDE has one spatial dimension in variable $s$. Letting $\mu=\omega_{1}u+\omega_{2}v$ and $\nu=-\omega_{2}u+\omega_{1}v$, (2.7) can be rewritten as three equations for new states $\hat{p},\hat{\mu}$ and $\hat{\nu}$. If one omits the trivial equation $\nu_{t}=0$, the equation (2.7) is reduced to the 1D acoustic equations,

[TABLE]

In this case, the equation depends on $\omega$ through the variable $\mu$. However, the equation itself is invariant over all $\omega$, owing to the fact that the problem (2.7) is isotropic. Moreover, note that this is exactly the same equation obtained in the physical space if you consider the case of a plane wave where the data varies only in the direction $\omega$ so that derivatives in the orthogonal direction vanish.

The Radon transform therefore transforms $n$-dimensional hyperbolic problems such as (2.5) and (2.7) into their $1$-dimensional counterparts (2.6) and (2.9), respectively.

2.2 Multi-dimensional extension of large time-step (LTS) methods

Previous dimensional splitting methods [28, 37, 14, 25] such as Strang splitting do not allow a natural extension of large time-step (LTS) methods [21, 23, 22] to multiple spatial dimensions. In order to take large time-steps for constant coefficient multi-dimensional hyperbolic problems, one can use the Fourier transform, for example. Upon taking the Fourier transform, one is left with a set of ordinary differential equations (ODEs) different from the original problem [38, 7]. On the other hand, using the Radon transform, one obtains a dimensional splitting that reduces the multi-dimensional problem into a family of one-dimensional counterparts of similar (if not identical) form. This allows 1D LTS methods to be applied for each of these problems, and the multi-dimensional solution is obtained by computing the inverse of the Radon transform. Moreover, the Radon transform provides an intuitive geometrical interpretation as a decomposition into planar waves and thus yields other useful applications. These applications will be illustrated in Section 4.

This multi-dimensional extension of the LTS method for the constant coefficient case is very straightfoward. Taking the Radon tranform of the problem as above, one obtains a set of $1$D problems such as (2.6) or (2.9). Then one applies the $1$D LTS solution operator $\mathcal{K}$ to evolve the initial data $\hat{u}_{0}(\omega,s)$ for each $\omega$ up to desired final time $T$. The operator $\mathcal{K}$ may depend on the direction $\omega$, so we denote the dependence as a parameter by writing $\mathcal{K}=\mathcal{K}(T;\omega)$. This yields the Radon transform of the solution at time $T$,

[TABLE]

Then, to compute the solution $q$ we can apply the inversion formula

[TABLE]

where the constant $c_{n}=(4\pi)^{(n-1)/2}\Gamma(n/2)/\Gamma(1/2)$ and $H_{s}$ denotes the Hilbert transform. Much is known about the inversion; see standard texts such as [18, 29] for more details.

This splitting can also be related to the Strang splitting, if one views it as a decomposition of the multi-dimensional problem into planar wave propagation. In Strang splitting one constructs the planar waves emanating in varying directions by dividing a single time-step into multiple successive planar wave propagations. The Radon transform decomposes the multi-dimensional directions by explicitly discretizing the sphere $S^{n-1}$.

Let us consider a concrete example, the 2D acoustic equation (2.7). Let us set $K_{0}=\rho_{0}=1$, so that we have the sound speed $c=1$, and impose the initial conditions

[TABLE]

The initial pressure profile is a cosine hump supported in a disk of radius $1$ centered at the origin, and the initial velocity profile is identically zero. We will also set absorbing boundary conditions in the manner to be described in Section 4.1.

On the transformed side (2.9), the evolution for any fixed direction $\omega\in S^{1}$ is given by the d’Alembert solution (2.13),

[TABLE]

This reduces to simple shifts at corresponding speeds, which can be computed easily up to any time $t$. This is precisely the 1D LTS solution for the constant coefficient case.

The solution to the acoustic equation computed on the domain $[-4,4]\times[-4,4]$ is shown in the left column of Figure 1. The Radon transform of the pressure term $\hat{p}$ is plotted in the right column of the same figure. Note that this problem is radially symmetric about the origin. A consequence of this is that the Radon transform is invariant with respect to the variable $\omega$, hence the Radon transform of the solutions at different times all appear as horizontal stripes. (There is a small amount of shift, following from the fact that for an image of even size $N$, the origin is chosen as the $(N/2,N/2)$-pixel, slightly off center.)

A key observation is that the evolution of the solution in the transformed variables is a sum of two shifting horizontal stripes, although the wave profile in the spatial domain propagates radially. For each fixed angle $\omega$, one only need solve the d’Alembert solution (2.13), which is easy to solve to any time $t$ by shifting the initial profile twice each according to two opposite speeds, and summing them. Intuitively, the shifts correspond to the propagation of decomposed planar waves for any fixed normal directions in $S^{1}$.

The actual computational did not make use of the continuous Radon transform (2.1), but rather a completely discrete approximation called the DRT, which will be introduced and discussed in further detail in Section 3. Here it will suffice to mention that a grid of size $128\times 128$ was used and prologation of $p=2$ was used for the DRT, and that the continuous transform can be obtained by an easy change of variables (3.5) which scale the domain and amplitude of the DRT, and that the change of variables do not affect the intertwining property. The 1D LTS method can still be used on the DRT, just as in the case of the continuous transform. The computational cost for this solution is conjectured to be $\mathcal{O}(N^{5/2}\log N)$: $\mathcal{O}(N^{2}\log N)$ for the forward DRT, $\mathcal{O}(N^{2})$ for the 1D LTS solution, and $\mathcal{O}(N^{5/2}\log N)$ for the inverse DRT (see Section 3). We note that the 1D LTS method can be applied for each angle in parallel.

Since the problem is radially symmetric, we can compare the solution to a 1D reference solution of high accuracy. We computed the 1D problem using Godunov flux with 4000 grid cells, as implemented in the Clawpack software package [12]. We compared the diagonal slice of our 2D solution at angle $\pi/4$ with the reference solution at time $t=0$ and $t=3$. To observe the accuracy of the solution with respect to the grid-size, solutions of sizes $N=8,16,32,64,128,256,512$ are also compared. The error was computed for the pressure variable $p$ using the weighted $L^{1}$ and $L^{2}$ norms,

[TABLE]

The comparison results are displayed in Figure 3. The error at the later time is at the level of the initial discretization error. We also observe that the convergence rate is between first and second-order with respect to the cell diameter $1/N$.

While this particular problem was radially symmetric, the splitting is by no means restricted to problems with radial symmetry. Let us modify the initial conditions above so that it is the sum of two cosine humps of different radii and heights,

[TABLE]

The splitting solution and its continuous Radon transform is plotted in Figure 2. In the first row of the figure, the initial condition and its Radon transform are shown. The two cosine humps in the initial condition each correspond to a sinusoidal signal on the transformed side. Recall the horizontal line centered at $s=0$ from the previous example (Figure 1). The sinusoidal shift away from $s=0$ is due to the fact that translation is an anisotropic operation. This can also be deduced from the transformed transport equation (2.6) in which the transport speed is $\theta\cdot\omega$, that is, $\cos\phi$ where $\phi$ is the angle between transport direction $\theta$ and the direction of the transform $\omega$. For example, when the cosine hump at the origin $p_{0}$ (2.12) is transported away from the origin by $r(\theta\cdot\omega)$, $\hat{p_{0}}$ is shifted by $\hat{p_{0}}(\omega,s-r\cos(\phi))$.

In any case, the solution is still given by the d’Alembert solution (2.13) and the acoustics equation can be solved exactly the same way as before. The DRT used in the actual computations are plotted in Figure 7. Each corresponds to a continuous transform in Figure 2.

2.3 Splitting for the nonlinear case

Here we will discuss how the splitting above to can be applied to a fully nonlinear system of hyperbolic equations. For a state vector $q$, such a PDE is given in the form

[TABLE]

where $f$ and $g$ are flux functions that can be nonlinear. Taking the Radon transform as before, we obtain

[TABLE]

Let us define the directional flux function as

[TABLE]

Then one obtains the nonlinear $1$D equations,

[TABLE]

As in the acoustics equation (2.9), the dependence on $\omega$ enters through the flux function $h(q)$, while the form of the equation is invariant with respect to $\omega$.

As an example, consider the shallow water equations in 2D, in which $\rho,u,v:\mathbb{R}^{+}\times\mathbb{R}^{2}\to\mathbb{R}$ denote water height, velocity in the $x_{1}$-direction and velocity in the $x_{2}$-direction, respectively,

[TABLE]

Here $\bar{g}$ denotes the gravitational constant. The Radon transform as above yields $1$D equation in the form (2.19), in the transformed velocity variables $\mu=\omega_{1}u+\omega_{2}v$ and $\nu=-\omega_{2}u+\omega_{1}v$,

[TABLE]

Note that the first two equations of (2.21) are just the shallow water equation in a single dimension in the normal direction of the hyperplane, whereas the third equation is the conservation of momentum in the transverse direction.

We observe that the transformed equations resemble a finite volume discretization. Let us say that $\xi_{i,j}$ is a discretization of the hyperplane $\{x\in\mathbb{R}^{n}:x\cdot\omega_{i}=s_{j}\}$. The specific discretization for the hyperplanes can take on many different forms, but here we will leave it in a general form. We denote the approximation to $\hat{q}(t_{n},\omega_{i},s_{j})$ at time-step $t_{n}$ by $\hat{Q}^{n}_{i,j}$,

[TABLE]

For each fixed direction $i=i_{0}$, the collection of hyperplanes $\{\xi_{i_{0},j}\}$ form a partition of the domain. We can consider these hyperplanes to be finite volume cells. In the equation (2.19) the flux function $h(q)$ assigns the flux between $\xi_{i_{0},j}$ and $\xi_{i_{0},j+1}$. If the cell boundary between $\xi_{i_{0},j}$ and $\xi_{i_{0},j+1}$ is denoted by $\xi_{i_{0},j+\frac{1}{2}}$, we define the numerical flux $F_{i_{0},j+\frac{1}{2}}^{n}$ to be an approximation to the flux at $\xi_{i_{0},j+\frac{1}{2}}$, valid from time-step $t_{n}$ to $t_{n+1}$. Then we have the finite volume update

[TABLE]

Once these updates are made for all $i$, the updated $\hat{Q}^{n+1}_{i,j}$ are combined through the inversion formula (2.11) to yield the numerical solution at time $t_{n+1}$.

The dimensional splitting strategy would be to compute the numerical flux $F^{n}_{i,j+\frac{1}{2}}$ by solving only the 1D Riemann problems in the $x$ and $y$ directions. Since the flux function $h(q)$ is a linear combination of normal fluxes $f(q)$ and $g(q)$ (2.18), we can compute $h$ once we have the approximation for these normal fluxes. In other words, we can solve the $1$D Riemann problems for piecewise constant jumps locally in $x$ and $y$ directions, then sum these fluxes across the cell boundary $\xi_{i,j+\frac{1}{2}}$ to obtain the flux between hyperplanes. One thereby decomposes the multi-dimensional Riemann problem into a set of single-dimensional ones, to be combined together by the inversion formula (2.11).

Unlike in the constant coefficient case, the flux function $h(q)$ must be updated at every time step, as is usually done for finite volume methods, although one may apply the nonlinear LTS method on the transformed problem regardless. This would be based on the 1D analogues studied in [21, 23, 22]. The fully nonlinear splitting will not be implemented here, but will be investigated in a future work.

3 Discrete Radon transform (DRT)

There are many different discretizations of the Radon transform and its inverse [27, 6, 9, 19, 2], arguably the most well-known being the filtered backprojection (FBP) algorithm [29]. However, its reliance on Fourier transforms and spherical harmonics lead to some filtering of high-frequency content, causing Gibbs phenomenon near the sharp edges in the solution. This is not suitable for use in hyperbolic PDEs, which are known to develop shock discontinuities.

Instead, we consider the use of a completely discrete analogue, namely the approximate discrete Radon transform (ADRT), which we refer to simply as the discrete Radon transform (DRT), introduced in [17, 8]. Rather than interpolating pixel values onto straight lines passing near it, DRT sums one entry for each row or column, along so-called digital lines or d-lines. The d-lines are defined recursively, allowing for a fast computation in $\mathcal{O}(N^{2}\log N)$ for an image of size $N\times N$. The back-projection is given by reversing the recursion, and is also fast with the same computational cost of $\mathcal{O}(N^{2}\log N)$. The precise definitions are given below.

3.1 Recursive definition of DRT

The d-lines of length $N$ are denoted by $D_{N}(h,s)$ with two parameters $h$ and $s$ (see Figure 4.) $h$ denotes the height ($x$-intercept) and $s$ the slope ($x$-displacement), and the pair corresponds to $s$ and $\boldsymbol{\omega}$ for the continuous transform (2.4), respectively. Although the same notation $s$ is used here again after having been used in the continuous setting (2.1) we will keep the notation in order to follow the intuitive notation of [30], and mark the continuous variable with a subscript $s_{c}$ whenever the two are used simultaneously. The definition uses the recursion in which d-lines of length $2n$ are split into left and right d-lines of half its length,

[TABLE]

The recursion (3.1) defines only a quarter of the possible d-lines, as the slope $s$ will range from 0 to $N$, corresponding to angles [math] to $\pi/4$ starting from the $x$-axis in the counter-clockwise direction. This is referred to by saying that the d-lines cover one quadrant, for the full transform one needs to cover the angles from [math] to $\pi$. The other d-lines can be computed by transposing and flipping the indices $h$ and $s$. We will denote the d-lines and DRT corresponding to the angular intervals $[0,\pi/4]$,$[\pi/4,\pi/2]$,$[\pi/2,3\pi/4]$, and $[3\pi/4,\pi]$ by $a$,$b$,$c$, and $d$.

The DRT of an array $A\in\mathbb{R}^{2}$ for the quadrants $a,b,c$ and $d$ are given by the summation of entries of $A$ over the d-lines,

[TABLE]

The full DRT is simply the ordered tuple of all quadrants, and we write

[TABLE]

See Fig. 5 for a visual illustration. Due to the recursive form of (3.1), the transform can be computed in $\mathcal{O}(N^{2}\log N)$. The parameters $h$ and $s$ belong to the range $[-s+1,N]$ and $[0,N]$ so $R^{a}_{N}A\in\mathbb{R}^{\left(\frac{3}{2}N+\frac{1}{2}\times N\right)}$. Therefore, $\mathcal{R}_{N}:\mathbb{R}^{N\times N}\to\mathbb{R}^{(6N+2)\times N}$. For example, the DRT of the solutions displayed in Figs. 1 and 2 are plotted in Figs. 6 and 7, respectively.

There is a simple relationship between the DRT and the continuous Radon transform. First let us say that $s_{c}\in[-1,1]$ (perhaps through proper scaling) and parameterize $\omega$ by $\omega=(\cos\theta,\sin\theta)$ where $\theta\in[0,\pi]$ . The relation to continuous variables $(s_{c},\omega)$ is given by

[TABLE]

Then the explicit relation between the DRT $\mathcal{R}_{N}$ and the continuous transform $\mathcal{R}$ are given after the density of the lines are also transformed depending on the angle by $\cos\theta$,

[TABLE]

We note that the DRT approximates the continuous transform with first order accuracy with respect to the grid cell width $1/N$ [8].

The back-projection is the dual of this transform with respect to the usual dot product in $\mathbb{R}^{N^{2}}$. We will denote the back-projection by $B_{N}$ or $\mathcal{R}_{N}^{T}$. If one explicitly forms the matrix for the linear transforms $\mathcal{R}_{N}$ and $\mathcal{R}_{N}^{T}$ they are indeed transposes of each other. $\mathcal{R}_{N}^{T}$ is the discrete analogue of $\mathcal{R}^{\#}$ in (2.2), a summation of all values assigned to d-lines passing through a point.

$\mathcal{R}^{T}_{N}$ is computed by reversing the sweep (3.1) above and computing a sequence of back-projections of decreasing size. Given a matrix $\hat{A}\in\mathbb{R}^{\left(\frac{3}{2}N+\frac{1}{2}\right)\times N}$, the reverse sweep for one quadrant is given by

[TABLE]

where the initial array $B_{N}(h,s)=\hat{A}_{h,s}$, and $n$ is set to $N/2$. This summation is repeated for the two half-images $B_{n}^{L}$ and $B_{n}^{R}$ on the LHS, until $n$ reaches $1$. Again, this summation is only for one quadrant, and we denote the end result as $\left(B_{N}^{a}\hat{A}\right)_{i,j}$.

The full back-projection is given by

[TABLE]

Its inversion algorithm using a full multi-grid method was demonstrated in [30] along with convergence analysis. In this paper, we use the conjugate gradient (CG) algorithm as will be discussed below in Section 3.2.

We end this section with the remark that the recursion (3.1) need not be in two and can be in any prime number, much like the fast Fourier transform [13].

3.2 Inversion of DRT with Conjugate Gradient Method

In order to use the dimensional splitting method to solve PDEs, a method for computing the inverse of a DRT (2.11) is needed. An inversion algorithm using a full multi-grid method appeared in [30]. Here we explore the application of the conjugate gradient method [16] to the least-squares problem

[TABLE]

The matrices for the transforms $\mathcal{R}^{T}_{N}$ and $\mathcal{R}_{N}$ are never explicitly formed, as we can use the fast algorithm. The computational cost of a DRT inversion is conjectured to be $\mathcal{O}(N^{5/2}\log N)$ for an $N\times N$ image [33]. Note that this is slightly more costly than $\mathcal{O}(N^{2}(\log N)^{3})$ that was conjectured for the full multi-grid method [30]. A more careful study of this inversion is of interest on its own right, and will appear elsewhere.

The inversion of the Radon transform, be it continuous or discrete, is mildly ill-posed [29]. Numerically speaking, this means the matrix $\mathcal{R}^{T}_{N}\mathcal{R}_{N}$ operator to be inverted will be ill-conditioned. One approach commonly used to deal with this issue is to use regularizations, for example in medical imaging applications. However, there is an important distinction to the tomography setting, namely that it is feasible to make additional measurements. In our setting, making more measurements from the original image $X$ in (3.8) would only incur additional computational effort, whereas in medical imaging it would require more physical measurements.

For example, in this work we perturb the range of the Radon transform and there is the possibility the perturbed function on the space of hyperplanes no longer lies in the range of the transform. The inversion (3.8) is exact only when $B$ lies in the range of $\mathcal{R}_{N}$, and this assumption cannot be satisfied in general once $B$ is evolved with respect to the dynamics of the transformed variables, as in (2.10). Therefore, changes in the transformed variables will cause $\hat{q}$ to depart from the range of $\mathcal{R}_{N}$. This becomes a source of error, incurring numerical artifacts in the computed inverse. The DRT employed in this paper also does incur these artifacts.

A brute force way to avoid this error is to make more measurements, i.e., oversample. We simply prolong the original image $q$ before manipulation, and restrict after the back-projection. This would correspond to making additional measurements in the tomography setting. This enlarges the range of the transform, and allows one to control the error. Therefore $\mathcal{R}_{N}$ will be replaced by $\mathcal{R}_{2pN}\mathcal{P}_{2p}$ where $\mathcal{P}_{p}$ is the [math]-th order prologation (where the value of each cell in the original grid is assigned to a $2p\times 2p$ cells in the enlarged grid) and $\mathcal{R}^{T}_{N}$ by $\mathcal{S}_{2p}\mathcal{R}^{T}_{2pN}$ where $\mathcal{S}$ is the restriction operator. The oversampling strategy does not affect the overall complexity.

3.3 DRT in dimension three

Just as the continuous Radon transform was defined in (2.1) for arbitrary number of dimensions $n$, the DRT can also be generalized to higher dimensions [26]. Here we treat the 3D case as an example. The recursive definitions (3.1) for the d-planes parametrized by three parameters $(h,s_{1},s_{2})$ can be derived for each hexadecant in a straightforward manner, as follows

[TABLE]

The DRT over one hexadecant (a quarter of a quadrant) is defined as the sum over these d-planes as in (3.2), and now the full transform in 3D is given by applying these to each of the hexadecant

[TABLE]

Hence, via transposing and flipping the indices as necessary, the full DRT is the ordered tuple

[TABLE]

The corresponding back-projection operation for a hexadecant is given by

[TABLE]

The full back-projection is then the average of back-projections $B_{N}$ over all hexadecants in $\mathcal{H}$,

[TABLE]

The computational cost of both operations would be $\mathcal{O}(N^{3}\log N)$.

4 Applications of the dimensional splitting

The dimensional splitting described in Section 2 above is a decomposition of hyperbolic solutions into evolution of planar waves. This decomposition can be useful in diverse settings. Here we discuss two applications: the absorbing boundary conditions and the displacement interpolation.

4.1 Absorbing boundary conditions

It is well-known that imposing absorbing boundary conditions to emulate infinite domains in multi-dimensional wave propagation is a challenging problem [15, 3, 11, 4]. On the other hand, the $1$D extrapolation boundary condition is much more tractable [25]. A major advantage of this splitting method is that the 1D extrapolation boundary conditions can be used on the transformed side at the computational boundary to avoid any reflections. This yields exactly the desired absorbing boundary conditions in the odd-dimensional case. Therefore, the dimensional splitting introduced in the previous section can be used directly to impose absorbing boundary conditions in 3D. On the contrary, there is an error caused by such an extrapolation in the even-dimensional case. This is due to the Huygens’ principle, evident in the presence of the Hilbert transform in the inversion formula (2.11). In this section, we discuss the type of error caused by imposing such extrapolation boundary conditions via the Radon transform in even dimensions.

In the true infinite domain, the non-zero values in transformed variables beyond the computational boundary of $\mathcal{S}^{n-1}\times\mathbb{R}$ affect the solution within the computational domain in the original variables $\mathbb{R}^{n}$. For example, the vertical translation of horizontal strips in Figure 1 should continue beyond the finite computational boundary, and by imposing a 1D extrapolation boundary condtion we would be neglecting this infinite propagation. To make this more precise, denote the computational (finite) transformed domain by $\Omega=\{(\omega,s)\in S^{n-1}\times(-b,b)\}$ for some $b>0$. Let $\chi_{\Omega}$ be the characteristic function of the finite domain and $\chi_{\mathbb{R}^{n}\setminus\Omega}=1-\chi_{\Omega}$. For $n$ even, the exact solution $q$ can be written as,

[TABLE]

Recall that $\mathcal{R}^{\#}$ is the back-projection (2.2). The first term in (4.2) is the approximate solution one would obtain if extrapolation boundary was set up at the boundaries $s=\pm b$. Let us call this approximate solution $q_{h}(x)$. Then the error is

[TABLE]

where $\mathrm{p.v}$ denotes the principal value integral. Note that in hyperbolic problems in free space, wave profiles will radiate outwards, that is, the support of $\hat{q}$ will be transported towards $r=\pm\infty$. This causes the RHS above to decay with time. Furthermore, the principal value integral is a smooth function of $s$ as long as $\partial^{n-1}\hat{q}/\partial r^{n-1}$ is integrable. Since we will also apply $\mathcal{R}^{\#}$, we expect the error to be smoother than $q(x)$.

Let us revisit the acoustic equations example (2.7) from Section 2.2, with initial conditions (2.12). Since $\hat{q}_{0}(x_{1},x_{2})$ is supported in $\{(\omega,s)\in S^{1}\times\mathbb{R}:\left\lvert s\right\rvert\leq 1\}$, we have a simple estimate for the case when $t$ is sufficiently large so that $(-1+t,1+t)\subset(b,\infty)$,

[TABLE]

where $\lVert\cdot\rVert_{1}$ is the $\ell^{1}$-norm for $\mathbb{R}^{3}$. Therefore we see that the effect of the extrapolation boundary decays relatively slowly, at the rate of $\mathcal{O}(1/t)$.

The solution is now compared with a fully 2D reference finite volume solution computed on a larger domain $[-8,8]\times[-8,8]$, using the wave propagation algorithm [24] with Lax-Wendroff flux and Van Leer limiter [39, 25], implemented in Clawpack. The reference solution was computed on a $1024\times 1024$ finer grid-cells of uniform size, then corresponding cells have been summed and compared with coarser cells of the DRT solution. The $L^{1}$ and $L^{2}$-norms of the difference over time is displayed below in Figure 8.

The error is in the order of $10^{-3}$ up to time $t=2.5$, before the profile starts approaching the boundary. The error from the truncation (4.4) begins to appear around time $t=2.5$ and peaks around time $t=5$, then decays to zero with time. The solution at time $t=3$, as it has begun to interact with the boundary, is shown in Figure 9. Note that there are no reflections from this boundary condition. On the other hand, a thin layer appears at the computational boundary. The layer is clearly non-physical, but is localized and has limited affect on the solution further in the interior. The DRT of the solution is also shown to its right, and we can see that the two pulses from the d’Alembert solution are hitting the 1D extrapolation boundaries (the top and bottom boundaries of the polygonal region in dark blue). The pulses first arrive at the DRT boundary at the angles $0,\pi/2$ and $\pi$, and this agrees with the solution plot to the left.

Comparing this solution to the reference solution, one discovers that the bulk of this error is concentrated near the thin layer which appears near the boundary. In Figure 10, we have plotted the difference between our solution and the reference solution on the computational domain $[-4,4]\times[-4,4]$ to the left. When we restrict the contour plot to the interior portion of the domain $[-3,3]\times[-3,3]$ as is shown to the right, we see that the error is significantly smaller as we move away from the boundary. The estimate (4.4) helps us understand this behavior. Note the decay of the principal value integral in (4.3): the further away the interior point is from the boundary, the smaller is the effect of the trunctation by $\chi_{\Omega}$. If we denote this distance by $d$, then the decay will be $\mathcal{O}(1/(d+t))$. As $t$ increases and the waves leave the domain, the error also decays, at the rate $\mathcal{O}(1/t)$ estimated by (4.4).

The observations above suggest several potential approaches in further improvements to this approach. For example, one can exploit the fact that the Hilbert transform in (4.3) commutes with translation, to emulate the effect of the infinite domain. One may also exploit the decay with respect to the distance to the boundary by placing finer mesh along the boundary in an adaptive fashion.

4.2 Displacement interpolation

In projection-based model reduction, the solution to a parametrized PDE is projected into a low-dimensional subspace, yielding a fast solver with significantly lower computational cost without compromising accuracy. To discover this low-dimensional subspace, the popular approach is to use proper orthogonal decomposition (POD) [5]. For hyperbolic PDEs, however, the solutions do not lie in a low-rank linear subspace, even for the simplest problems [36, 1, 10, 32, 34]. For instance, the d’Alembert solution (2.13) is a linearly independent function of $s$ for each $t>0$. It is easy to see that a linear projection of this solution to a low-dimensional basis would not yield a good approximation of the solution. Naturally, methods to remove translational symmetry [36, 32, 34] are being actively explored. This is also intimately related to the concept of displacement interpolation, a term we borrow from the optimal transport literature [40], in which one aims to minimize the Wasserstein distance, although we will not make the connection more explicit here.

Let us first illustrate how displacement interpolation arises naturally, with a simple 1D example. Suppose $\phi_{0}$ is a hat function, given by

[TABLE]

for some $h>0$. Let $\phi_{1}$ and $\phi_{2}$ be translation and scaling of $\phi_{0}$,

[TABLE]

For $h=0.1$, the two functions are shown in the first row of Figure 11. The linear interpolation $\psi$ of the two functions with weights $(1-\tau)$ and $\tau$ is given by

[TABLE]

whereas a displacement interpolation between the two functions under a simple transport map (1D translation) would be given by

[TABLE]

The two interpolants for $\tau=0.25$ are plotted in the bottom row in Figure 11 . Since $\phi_{1}$ and $\phi_{2}$ are both translates of a scalar multiple of $\phi_{0}$, the displacement interpolation reveals the low-rank nature of the two functions, whereas the linear interpolant remains rank two for $\tau\in(0,1)$.

In practice, one must be able to deduce that $\phi_{1}$ and $\phi_{2}$ above lie in the translates of $\text{span}\{\phi_{0}\}$ without a priori knowledge. To achieve this, one may apply the template-fitting procedure [36] which solves the minimization problem

[TABLE]

where $\mathcal{K}$ is the translation operator, $\mathcal{K}(\tau)[\phi_{1}(x)]=\phi_{1}(x-\tau)$, then perform a singular value decomposition (SVD) on $\{\phi_{2},\mathcal{K}(\tau_{*})\phi_{1}\}$ [36, 32]. However, this simple formulation does not take into account multiple traveling speeds or heavily deforming profiles, which limits its applicability. Transport reversal was introduced in [34] to overcome these limitations. The algorithm is a greedy iteration over a generalized form of (4.9), which decomposes the 1D function $\phi_{2}(x)$ into multiple traveling structures. To be more precise, given two functions $\phi_{1}$ and $\phi_{2}$ as in (4.9), transport reversal yields the decomposition

[TABLE]

where $\eta_{k}$ is a scaling function and $\rho_{k}$ a cut-off function. For more detailed treatment of this decomposition in 1D, we refer the reader to [34]. Let us suppose we have computed this decomposition. The displacement interpolant $\psi_{D}$ resulting from this decomposition is set to satisfy,

[TABLE]

As an example, let us assume we are given two functions $\phi_{1}$ and $\phi_{2}$ as shown in Figure 12(a) and 12(b). These are taken from the 1D slice located at $s=\tan(\frac{5}{8}\pi)(N-1)$ of the DRT from Figure 7. The transport reversal would decompose $\phi_{2}$ into a superposition of two traveling profiles,

[TABLE]

each plotted in Figure 13. In exact arithmetic, the two iterations of transport reversal would pick up exactly the d’Alembert solution (although in practice numerical error would require further iterations to pick off the residuals). That is, we would obtain $h_{1}=h_{2}=1/2$ and $\rho_{1}=\rho_{2}=1$ with $\nu_{1}=-\nu_{2}=1$ and $K=2$ in (4.10). Now, the displacement interpolation for $\tau=1/2$ can be computed, yielding $\psi_{D}$ shown in Figure 12(c). The exact evolution of the two iterates (4.12) are shown in Figure 13.

Now, this displacement interpolation was done for a single slice of the fixed $\omega$ in the transformed variables. Suppose we are given a function $\varphi$ in 2D. Then by performing the same transport reversal on its Radon transform $\hat{\varphi}$ for all $\omega$ as functions of the variable $s$, we obtain an extension of the 1D displacement interpolant (4.8) to higher spatial dimensions. For each fixed angle $\omega_{i}\in S^{1}$ we obtain the transport reversal in terms of the traveling structures,

[TABLE]

These can be used for displacement interpolation as above for each $\omega$. The inverse transform can be taken to obtain the displacement interpolant $\psi_{D}$.

Let us clarify the implication. For the acoustic equation example with the initial condition (2.15), we were given a snapshot of the solution $q$ at time $t_{1}=0.5$ and $t_{2}=1.5$. From the two snapshots, we were able to accurately approximate the solution for all time, without additional information about the dynamics, without even knowing the PDE. Thus this interpolant can be more useful than linear interpolation: the linear subspace spanned by $\{q(t_{1},x),q(t_{2},x)\}$ does not a contain a good approximation for representing the evolving solution.

This ability to compute the displacement interpolation by exploiting the simple dynamic on the transformed side will be useful in the future development of transport reversal as a model reduction tool in multi-dimensional settings.

5 Conclusion and future work

We have introduced a dimensional splitting method using the intertwining property of the Radon transform. Its applications in solving hyperbolic PDEs, imposing absorbing boundary conditions, and computing displacement interpolations were discussed. For the inversion of DRT the conjugate gradient method was used.

As noted in Section 4.1, the dimensional splitting proposed here used with DRT in 3D (Section 3.3) allows one to impose absorbing boundary conditions for 3D problems without incurring any error of the type (4.3) that appears in 2D. This will be verified in future work. The application of this splitting to fully nonlinear hyperbolic PDEs as discussed in Section 2.3 will be studied as well. The utility of the Radon transform for displacement interpolation (Section 4.2) will be much more compelling when used in conjunction with the fully multi-dimensional transport reversal [34] as a model reduction tool for general hyperbolic PDEs, and work is underway for such an extension.

The number of CG iterations for the inversion (3.8) can be estimated to justify the conjectured $\mathcal{O}(N^{5/2}\log N)$ cost for inversion: this and other inversion results will appear elsewhere. While the prologation used in the inversion (Section 3.2) causes expense only of a constant factor, it can be of significant computational cost. Other approaches to reduce the amount of computational effort will be explored. The DRT is essentially a structured matrix multiplication and may be amenable to parallelization, and its performance on graphical processing units (GPUs) will be a future topic of research.

Acknowledgments

The author expresses gratitude to Randall J. LeVeque for carefully reviewing this manuscript. He also thanks Gunther Uhlmann for helpful discussions.