A pseudo-spectral splitting method for linear dispersive problems with transparent boundary conditions

TL;DR

This paper introduces a pseudo-spectral splitting method for linear dispersive equations with variable coefficients on unbounded domains, effectively implementing transparent boundary conditions to improve computational efficiency and accuracy.

Contribution

It develops a novel spectral splitting scheme that handles non-homogeneous transparent boundary conditions for dispersive problems, reducing computational cost while maintaining accuracy.

Findings

Efficient spectral discretization reduces grid points needed.

Implicit treatment of dispersive terms enhances stability.

Numerical simulations validate the scheme's effectiveness.

Abstract

The goal of the present work is to solve a linear dispersive equation with variable coefficient advection on an unbounded domain. In this setting, transparent boundary conditions are vital to allow waves to leave (or even re-enter) the, necessarily finite, computational domain. To obtain an efficient numerical scheme we discretize space using a spectral method. This allows us to drastically reduce the number of grid points required for a given accuracy. Applying a fully implicit time integrator, however, would require us to invert full matrices. This is addressed by performing an operator splitting scheme and only treating the third order differential operator, stemming from the dispersive part, implicitly; this approach can also be interpreted as an implicit-explicit scheme. However, the fact that the transparent boundary conditions are non-homogeneous and depend implicitly on the numerical solution presents a significant obstacle for the splitting/pseudo-spectral approach investigated here. We show how to overcome these difficulties and demonstrate the proposed numerical scheme by performing a number of numerical simulations.

Equations

$$u_t+g(x)u_x+u_{xxx}=h(t,x),t\in [0,T],x\in \mathbb{R}.$$

$$\{\begin{array}{l}{u}_t+g(x)u_x+u_{xxx}=h(t,x),(t,x)\in [0,T]\times \mathbb{R},\\ u(0,x)=u_0(x).\end{array}$$

$$\{\begin{array}{l}{v}_t+g(x)v_x=0,\\ v(0,x)=v_0(x),\end{array}$$

$$\{\begin{array}{l}{w}_t+w_{xxx}=0,\\ w(0,x)=w_0(x).\end{array}$$

$$u^{m+1}(x)=\Phi^{\tau}\circ\Psi^{\tau}\big{(}u^{m}\big{)}(x).$$

$$v^{m+1/2}+\tau g{u}_x^m=u^m,u^{m+1}+\tau u_{xxx}^{m+1}=v^{m+1/2},m\ge 0.$$

$$u^{m+1}+\tau u_{xxx}^{m+1}=v^{m+1/2}=u^m-\tau g(x)u_x^m,u^0(x)=u(0,x).$$

$$u^{m+1}(x)=\Psi^{\tau}\circ\Phi^{\tau}\big{(}u^{m}\big{)}(x).$$

$$v^{m+1/2}+\tau v_{xxx}^{m+1/2}=u^m.$$

$$u^{m+1}=v^{m+1/2}-\tau g{v}_x^{m+1/2}.$$

$$u^m(x)\to 0,|x|\to \infty ,\text{for every }m\ge 0.$$

$$\hat{u}(z):=\mathcal{Z}(\mathbf{u})(z)=\sum _{l=0}^{\infty }{z}^{-l}{u}^l,z\in \mathbb{C},|z|>\rho \ge 1,$$

$$(\mathbf{u}{\ast }_d\mathbf{v}{)}^m:=\sum _{j=0}^m{u}^j{v}^{m-j},m\ge 0.$$

$$\mathbb{R}=(-\infty ,-1]\cup (-1,1)\cup [1,\infty ).$$

$$\mathcal{Z}(\{u^{m+1}(x){\}}_{m\ge 0})(z)=z\hat{u}(x,z)-z{u}^0(x)=z\hat{u}(x,z),x\in {\Omega }_{-}\cup {\Omega }_{+},$$

$$z(\hat{u}+\tau {\hat{u}}_{xxx})=\hat{u}-\tau g_{\pm }{\hat{u}}_x,x\in {\Omega }_{-}\cup {\Omega }_{+}.$$

$$\hat{u}(x,z)=\sum _{j=1}^3{c}_j(z)e^{r_j(z)x},$$

$$z\tau r^3+\tau g_{\pm }r+z-1=0.$$

$$\begin{split}&r_{j}(z)=\omega^{j-1}\,\zeta(z)-\frac{g_{\pm}}{3z\omega^{j-1}\,\zeta(z)},\quad\omega=e^{i2/3\pi},\quad j=1,2,3,\\ &\zeta(z)=-\bigg{(}\frac{G(z)}{2}\bigg{)}^{1/3},\\ &G(z)=\frac{z-1}{z\tau}+\sqrt{\bigg{(}\frac{z-1}{z\tau}\bigg{)}^{2}+\frac{4}{27}\bigg{(}\frac{g_{\pm}}{z}\bigg{)}^{3}}.\end{split}$$

$$\begin{array}{rl}& \hat{u}(x,z)=c_1(z)e^{r_1(z)x},x\in {\Omega }_{+},\\ & \hat{u}(x,z)=c_2(z)e^{r_2(z)x}+c_3(z)e^{r_3(z)x},x\in {\Omega }_{-}.\end{array}$$

$$\begin{array}{rl}& {\hat{u}}_x(x,z)=c_1(z)r_1(z)e^{r_1(z)x}=r_1(z)\hat{u}(x,z),x\in {\Omega }_{+},\\ & {\hat{u}}_x(x,z)=c_2(z)r_2(z)e^{r_2(z)x}+c_3(z)r_3(z)e^{r_3(z)x},x\in {\Omega }_{-}\end{array}$$

$$\begin{array}{rl}& {\hat{u}}_{xx}(x,z)=r_1^2(z)\hat{u}(x,z),x\in {\Omega }_{+},\\ & {\hat{u}}_{xx}(x,z)=c_2(z)r_2^2(z)e^{r_2(z)x}+c_3(z)r_3^2(z)e^{r_3(z)x},x\in {\Omega }_{-}.\end{array}$$

$$\hat{u}(x,z)=\frac{1}{r_1^2(z)}{\hat{u}}_{xx}(x,z),{\hat{u}}_x(x,z)=\frac{1}{r_1(z)}{\hat{u}}_{xx}(x,z),x\in {\Omega }_{+}.$$

$$\hat{u}(x,z)=\bigg{(}\frac{1}{r_{2}(z)}+\frac{1}{r_{3}(z)}\bigg{)}\,\hat{u}_{x}(x,z)-\frac{1}{r_{2}(z)r_{3}(z)}\,\hat{u}_{xx}(x,z),\quad x\in\Omega_{-}.$$

$$\begin{split}\hat{u}(-1,z)&=\bigg{(}\frac{1}{r_{2}(z)}+\frac{1}{r_{3}(z)}\bigg{)}\hat{u}_{x}(-1,z)-\frac{1}{r_{2}(z)r_{3}(z)}\hat{u}_{xx}(-1,z),\\ \hat{u}(1,z)&=\frac{1}{r_{1}^{2}(z)}\hat{u}_{xx}(1,z),\\ \hat{u}_{x}(1,z)&=\frac{1}{r_{1}(z)}\hat{u}_{xx}(1,z).\end{split}$$

$$\begin{split}\mathbf{u}(-1)&=\mathcal{Z}^{-1}\bigg{\{}z\mapsto\bigg{(}\frac{1}{r_{2}(z)}+\frac{1}{r_{3}(z)}\bigg{)}\hat{u}_{x}(-1,z)\bigg{\}}-\mathcal{Z}^{-1}\bigg{\{}z\mapsto\frac{1}{r_{2}(z)r_{3}(z)}\hat{u}_{xx}(-1,z)\bigg{\}},\\ \mathbf{u}(1)&=\mathcal{Z}^{-1}\bigg{\{}z\mapsto\frac{1}{r_{1}^{2}(z)}\hat{u}_{xx}(1,z)\bigg{\}},\\ \mathbf{u}_{x}(1)&=\mathcal{Z}^{-1}\bigg{\{}z\mapsto\frac{1}{r_{1}(z)}\hat{u}_{xx}(1,z)\bigg{\}}.\end{split}$$

$$\begin{split}\mathbf{Y}_{1}&=\mathcal{Z}^{-1}\bigg{\{}z\mapsto\frac{1}{r_{2}(z)}+\frac{1}{r_{3}(z)}\bigg{\}},\\ \mathbf{Y}_{2}&=\mathcal{Z}^{-1}\bigg{\{}z\mapsto-\frac{1}{r_{2}(z)r_{3}(z)}\bigg{\}},\\ \mathbf{Y}_{3}&=\mathcal{Z}^{-1}\bigg{\{}z\mapsto\frac{1}{r_{1}^{2}(z)}\bigg{\}},\quad\mathbf{Y}_{4}=\mathcal{Z}^{-1}\bigg{\{}z\mapsto\frac{1}{r_{1}(z)}\bigg{\}}.\end{split}$$

$$\begin{array}{rl}\mathbf{u}(-1)& ={\mathbf{Y}}_1{\ast }_d{\mathbf{u}}_x(-1)+{\mathbf{Y}}_2{\ast }_d{\mathbf{u}}_{xx}(-1),\\ \mathbf{u}(1)& ={\mathbf{Y}}_3{\ast }_d{\mathbf{u}}_{xx}(1),\\ {\mathbf{u}}_x(1)& ={\mathbf{Y}}_4{\ast }_d{\mathbf{u}}_{xx}(1).\end{array}$$

$$\begin{split}u^{m}(-1)-Y_{1}^{0}\,u^{m}_{x}(-1)-Y_{2}^{0}u^{m}_{xx}(-1)&=\sum_{j=1}^{m}\big{(}Y_{1}^{j}\,u^{m-j}_{x}(-1)+Y_{2}^{j}\,u^{m-j}_{xx}(-1)\big{)}=:h_{1}^{m},\\ u^{m}(1)-Y_{3}^{0}\,u_{xx}^{m}(1)&=\sum_{j=1}^{m}Y_{3}^{j}\,u^{m-j}_{xx}(1)=:h_{2}^{m},\\ u_{x}^{m}(1)-Y_{4}^{0}\,u_{xx}^{m}(1)&=\sum_{j=1}^{m}Y_{4}^{j}\,u_{xx}^{m-j}(1)=:h_{3}^{m}.\end{split}$$

$$\{\begin{array}{ll}& u^m+\tau u_{xxx}^m=u^{m-1}-\tau g(x)u_x^{m-1},x\in (-1,1),\\ & u^m(-1)-Y_1^0{u}_x^m(-1)-Y_2^0{u}_{xx}^m(-1)=h_1^m,\\ & u^m(1)-Y_3^0{u}_{xx}^m(1)=h_2^m,\\ & u_x^m(1)-Y_4^0{u}_{xx}^m(1)=h_3^m,\\ & u^0=u(0,x).\end{array}$$

Full text

Keywords— splitting methods, linear dispersive problems, pseudo-spectral methods, transparent boundary conditions

1 Introduction

In this paper, we consider a linear dispersive problem with a variable coefficient advection in one space dimension

[TABLE]

This partial differential equation (PDE) consists of an advection part given by $g(x)u_{x}$ and a dispersive part given by $u_{xxx}$. Applications can be found, for example, in modelling long waves in shallow water with an uneven bottom (see, for example, [13, 19]).

We face several difficulties in designing an efficient numerical method to solve this problem. The third order differential operator causes explicit methods to take excessively small time steps. This could be remedied by employing an implicit method, such as the Crank–Nicolson scheme that is used, for example, in [5]. In this paper, our goal is to use a spectral approach to discretize space. This drastically reduces the number of grid points needed. However, the variable coefficient advection results in a full matrix. Thus, applying an implicit method to the entire problem is very costly in terms of computational effort. We, therefore, choose to employ a splitting approach to separate the dispersive part, which will be treated implicitly, from the advection part, which will be treated explicitly. Operator splitting schemes have been employed for dispersive problems before and we refer the reader to the literature [9, 10, 11, 14, 17]. Finally, problem (1) is posed on an unbounded domain. However, to perform numerical simulations the domain has to be restricted to a compact subset of the real line. Our goal is to derive a numerical scheme that, despite this truncation, retains the dynamics of the original problem. This is accomplished using so-called transparent boundary conditions (TBCs). Transparent boundary conditions are derived in such a manner that a wave propagating to either the left or the right can leave the computational domain without producing reflections. A previous work which combines splitting scheme with absorbing boundary conditions for the Schrödinger equation can be found in [4]. We further note that combining TBCs with a pseudo-spectral approach is not entirely straightforward as TBCs can not be formulated as Dirichlet or Neumann boundary conditions and the value imposed depends on the history of the numerical solution at the boundary.

This work is inspired by the paper of Besse et al. [5]. In their work the $\mathcal{Z}\text{-transform}$ is used to compute transparent boundary conditions for a third-order linear problem with constant coefficients. The temporal discretization is carried out by the Crank–Nicolson scheme and the spatial discretization by finite differences. Differently from their approach we will follow a splitting strategy in order to divide the full equation into its dispersive part $u_{t}+u_{xxx}=0$ and its advection part $u_{t}+g(x)u_{x}=0$. The transparent boundary conditions are then derived in a semi-discrete setting using the backward Euler and the explicit Euler method for the dispersive and the advection equation, respectively. Let us also remark that, since we limit ourselves to the first order Lie–Trotter splitting in this work, the resulting numerical time integrator can also be written as an implicit-explicit (IMEX) scheme. The space discretization is performed by a pseudo-spectral approach similarly to what is described in [18, 20]. This results in super-polynomial convergence in space. Consequently very accurate numerical solutions are achieved by using only a small number of grid points.

Due to the non-locality in time (in 1-D) and space (2-D or higher dimensions), transparent boundary conditions are expensive to compute. While it is possible to employ them in 1-D, the 2-D case becomes impracticable and one needs an approximation of these conditions, the so called absorbing boundary conditions (ABCs). In recent years a lot of work has been done for the Schrödinger equation coupled with transparent boundary conditions, see [2, 3]. For third order problems less literature is available. We refer the reader to [5, 6, 20].

The paper is organized as follows. In Section 2 we derive a semi-discrete scheme, discrete in time and continuous in space, by applying the Lie–Trotter splitting. In Section 3 we impose transparent boundary conditions for the scheme derived in Section 2. In particular, we determine the proper values of the numerical solution at the boundaries with the help of the $\mathcal{Z}$-transform. In Section 4 we describe a pseudo-spectral method for the spatial discretization which takes the TBCs into account. Finally, in Section 5 we present some numerical results that illustrate the theoretical findings.

2 Time discretization: a splitting approach

In this section we derive a semi-discrete scheme by applying Lie–Trotter splitting to (1). Splitting methods give us the possibility to choose different numerical methods in order to solve the advection part and the dispersive part. For the dispersive part it is convenient to employ an implicit time integrator. Indeed, the third derivative in space makes the Courant–Friedrichs–Lewy (CFL) condition very strict for any explicit scheme. In this work we choose the implicit Euler scheme for the dispersive part. For the advection part we use a different time integration scheme. A possible choice is the explicit Euler method, since the CFL condition given by the advection term is generally not prohibitive. Both the implicit and the explicit Euler method converges with order one. Numerical methods of order one are a sufficient choice in our case because of the bottleneck given by the Lie–Trotter splitting, which converges with order one only.

We consider problem (1) supplemented with initial data $u_{0}$

[TABLE]

The initial data $u_{0}$ and the source term $h$ are assumed to be compactly supported. The coefficient $g(x)$ is constant outside a bounded interval. For simplicity of exposition we assume $h=0$. Let $u^{m}(x)$ be a numerical approximation to the exact solution $u(t,x)$ of (2) at time $t=t_{m}=m\tau$. We split problem (2) into two sub-problems:

[TABLE]

[TABLE]

Let us denote by $\Psi^{\tau}(v_{0})$ and $\Phi^{\tau}(w_{0})$ the flows of (3) and (4), respectively. We approximate the solution $u(t,x)$ at time $t=t_{m+1}$ starting from $u^{m}(x)\leavevmode\nobreak\ \approx\leavevmode\nobreak\ u(t_{m},x)$ by applying a Lie–Trotter splitting, i.e.

[TABLE]

The semi-discrete numerical schemes for (3) and (4), respectively, take the form

[TABLE]

Composing the two flows by Lie–Trotter splitting we obtain the numerical scheme

[TABLE]

Remark*.*

In this work the order of composition of the flows $\Psi^{\tau}$ and $\Phi^{\tau}$ is important. When we exchange the flow’s composition we obtain

[TABLE]

The numerical solution for \Phi^{\tau}\big{(}u^{m}\big{)} is given by implicit Euler, we have

[TABLE]

The numerical solution for \Psi^{\tau}\big{(}v^{m+1/2}\big{)} is given by explicit Euler, we have

[TABLE]

It is not possible from (7)-(8) to obtain a simple formula which relates $u^{m+1}$ to $u^{m}$ as in (6). Therefore, the numerical scheme in (6) is easier to handle.

3 Transparent boundary conditions

Transparent boundary conditions are non-local in time. Therefore, they cannot be computed separately for the two sub-problems of the Lie–Trotter splitting. The idea is first to derive the numerical scheme (6) and then to compute the boundary conditions for that particular scheme. The drawback of this approach is that different time discretizations require different boundary conditions. On the other hand, the computed boundary conditions produce (theoretically) no reflections since they are perfectly designed for the employed numerical scheme.

The mathematical tool we use to derive the transparent boundary conditions is the $\mathcal{Z}$-transform. The $\mathcal{Z}$-transform requires an equidistant time discretization. Further, we assume that

[TABLE]

We refer to (9) as the decay condition. Given a sequence $\mathbf{u}=\{u^{l}\}_{l\geq 0}$ its $\mathcal{Z}$–transform is defined by

[TABLE]

where $\rho$ is the radius of convergence of the series.We recall the main properties of the $\mathcal{Z}$-transform that are used in the following.

Linearity: for $\alpha,\beta\in\mathbb{R}$, $\mathcal{Z}(\alpha\mathbf{u}+\beta\mathbf{v})(z)=\alpha\hat{u}(z)+\beta\hat{v}(z)$;

• 2.

Time advance: for $k>0$, $\mathcal{Z}(\{u^{l+k}\}_{l\geq 0})(z)=z^{k}\hat{u}(z)-z^{k}\sum_{l=0}^{k-1}z^{-l}u^{l}$;

• 3.

Convolution: \mathcal{Z}\big{(}\mathbf{u}*_{d}\mathbf{v}\big{)}(z)=\hat{u}(z)\hat{v}(z);

where $*_{d}$ denotes the discrete convolution

[TABLE]

For more details and properties about the $\mathcal{Z}$-transform we refer the reader to [5].

Now we determine the transparent boundary conditions for the scheme (6). For this we assume that the coefficient $g(x)$ is constant for $x\leq\alpha$ and $x\geq\beta$ and that $u^{0}$ is continuous and its support is contained in $[\alpha,\beta]$. Without loss of generality, we further assume $\alpha=-1$ and $\beta=1$. We then split the domain into three parts

[TABLE]

Let $\Omega_{-}=(-\infty,-1]$ and $\Omega_{+}=\leavevmode\nobreak\ [1,\infty)$ be the two outer domains. Note that $u^{0}(x)\leavevmode\nobreak\ =0$ for $x\in\Omega_{-}\cup\Omega_{+}$. We now consider the problem in the two outer domains. We denote by $g_{\pm}$ the constant values that $g$ assumes for $x\in\Omega_{-}$ and $x\in\Omega_{+}$, respectively. Computing the $\mathcal{Z}$-transform of the sequence $\{u^{m}(x)\}_{m\geq 0}$, where $x\in\Omega_{-}\cup\Omega_{+}$ plays the role of a parameter, we get

[TABLE]

where we used the time advance property of the $\mathcal{Z}$-transform and the assumption of a supported $u^{0}(x)$ in $\Omega$. By applying the $\mathcal{Z}$-transform to (6), we obtain an ordinary differential equation in the variable $x$:

[TABLE]

The solution of this equation for $x\in\Omega$ is obtained by employing an exponential ansatz. It is given by

[TABLE]

where $c_{j}$ are coefficients and $r_{j}$, $j=1,2,3$ are the roots of the characteristic polynomial associated to (10), i.e.

[TABLE]

The roots $r_{j}$ can be computed analytically. They are given in [5] by

[TABLE]

In the same paper [5] the following property is shown:

Theorem 3.1**.**

The roots $r_{j}(z)$ for $j\in\{1,2,3\}$ are such that $\mathrm{Re}\,r_{1}(z)<0$ and $\mathrm{Re}\,r_{2,3}(z)>1$, $z=\rho\mathrm{e}^{i\theta}$, $\theta\in[0,2\pi)$ and for a fixed $\rho>1$.

The solution (11) must be coupled with the decay condition $u^{m}(x)\to 0,\,|x|\to\leavevmode\nobreak\ \infty$. Therefore, we choose the coefficients $c_{j}(z)$ such that the decay condition is satisfied. This requires $c_{1}(z)=0$ for $x\in\Omega_{-}$. Otherwise the solution would grow to infinity as $\mathrm{Re}\,r_{1}(z)<0$. Similarly for $x\in\Omega_{+}$ we impose $c_{2}(z)=c_{3}(z)=0$. Otherwise the solution would grow to infinity as $\mathrm{Re}\,r_{2,3}(z)>0$. Finally we obtain

[TABLE]

The coefficients $c_{j}$ are not yet determined. Notice that differentiating the equation (13) w.r.t. $x$ we obtain

[TABLE]

and differentiating once more we get

[TABLE]

Rewriting formulae (13)–(15), we can express the function $\hat{u}$ of (10) and $\hat{u}_{x}$ for $x\in\Omega_{+}$ in terms of its second derivative

[TABLE]

The function $\hat{u}$ for $x\in\Omega_{-}$ can be expressed in terms of its first and second derivative

[TABLE]

The transparent boundary conditions in the $\mathcal{Z}$-transformed space are simply computed by evaluating $\hat{u}$ at $x=\pm 1$ and $\hat{u}_{x}$ at $x=1$, namely

[TABLE]

The solution in the physical space is then obtained by applying the inverse $\mathcal{Z}\text{-transform}$ to (16). From (16) we get

[TABLE]

By the convolution property of the $\mathcal{Z}$-transform, we have $\mathbf{u}*_{d}\mathbf{v}=\mathcal{Z}^{-1}\{z\mapsto\hat{u}(z)\hat{v}(z)\}$. Let us denote by

[TABLE]

Equations (17) take the form

[TABLE]

We rewrite the boundary conditions in order to highlight the part that depends on $m$ (collected on the left hand-side) and the history that depends on previous time steps (collected on the right-hand side)

[TABLE]

The quantities $\mathbf{Y}_{k}$, $k=1,\dots,4$ can be computed numerically as we briefly describe in Section 5. For more details about the $\mathcal{Z}$-transform and its numerical implementation as well as the inverse $\mathcal{Z}$-transform, we refer the reader to [5, 21].

Summing up, we obtain the following boundary value problem to solve:

[TABLE]

The task of solving this problem will be carried out in the next sections.

4 Pseudo-spectral space discretization

The space discretization of problem (19) is carried out by a pseudo-spectral method. We implement a dual-Petrov–Galerkin method. This method has been employed for the pure dispersive equation in [20]. We remark that other spatial discretizations are also possible. For example, in [5] a finite difference method is employed to discretize in space. Finite differences are easy to implement, but an accurate numerical solution is achieved only when a fairly large amount of grid points are used. On the other hand, pseudo-spectral methods require a modest number of points and provide very accurate solutions. These methods have to be carefully designed in order to get banded spatial discretization matrices so that the resulting linear system is cheap to solve.

Let $a<b$ and let $\mathcal{P}_{N}$ be the space of polynomials of degree less or equal than $N$ on the interval $[a,b]$. We define

[TABLE]

where $Y^{0}_{k}$, $k=1,\dots,4$ are the coefficients described in Section 3.

Let $(u,v)=\int_{a}^{b}uv\,\text{d}x$ be the usual $L_{2}$ inner product. $V_{N}^{*}$ is defined as the dual space of $V_{N}$ so that for every $u\in V_{N}$ it holds $(u_{xxx},v)=-(u,v_{xxx})$ for every $v\in V_{N}^{*}$.

Lemma 4.1**.**

The dual space $V_{N}^{*}$ of $V_{N}$ is given by

[TABLE]

Proof.

Integrating the quantity $(u_{xxx},v)$ by parts three times we get

[TABLE]

The boundary terms in (21) must vanish. For $x=b$ we have

[TABLE]

The last equality is obtained substituting in place of $u_{x}(b)$ and $u(b)$ the relations with $u_{xx}(b)$ given by the space $V_{N}$. Similarly for $x=a$ we have

[TABLE]

From equations (22) and (23) we get the boundary relations for the dual space $V_{N}^{*}$. ∎

Let $w^{m}\in V_{N}$ be the solution of the homogeneous problem associated to (19), i.e. where we set $h_{i}^{m}=0$ for $i=1,2,3$. The solution $u^{m}$ of the inhomogeneous problem is then given by $u^{m}=w^{m}+p_{2}^{m}$, where $p_{2}^{m}$ is the unique polynomial of degree $2$ such that

[TABLE]

We remark that the polynomial $p_{2}^{m}$ depends on $h_{i}^{m}$, $i=1,2,3$. Therefore, it must be computed at each time step. For more details we refer the reader to [18, 20].

The variational formulation of (6) reads: find $u^{m}=w^{m}+p_{2}^{m}\in P_{N}$ with $w^{m}\in V_{N}$ such that for every $v\in V_{N}^{*}$ it holds

[TABLE]

In order to numerically evaluate the $L_{2}$ inner products $(\cdot,\cdot)$, we have to choose some interior collocation points. A typical choice are the Gauss–Lobatto points. These points are efficient for solving second order differential equations. However, due to the lack of symmetry of the considered third order problem, a better option is to choose as interior collocation points $\{x_{k}\}_{k=2}^{N-1}$ the roots of the Jacobi polynomial $P^{(2,1)}_{N-2}(x)$, $x_{1}=a$, $x_{N}=b$, see [1, 16].

The discrete inner product associated to the Gauss–Jacobi quadrature rule is given by

[TABLE]

where the associated weights are

[TABLE]

The employed quadrature rule integrates exactly polynomials up to degree $2N-\leavevmode\nobreak\ 2$. A detailed description of the generalized quadrature rule for third order problems used in this work can be found in [12].

We substitute $u^{m}=w^{m}+p_{2}^{m}$ and we use the fact that $w^{m}_{xxx}=u^{m}_{xxx}$, so we obtain

[TABLE]

For any $0\leq k\leq N-3$, set

[TABLE]

where $L_{k}$ is the $k$th Legendre polynomial and $\alpha_{k},\beta_{k},\gamma_{k},$ are uniquely determined so that $\phi_{k}$ and $\psi_{k}$ belong to $V_{N}$ and $V_{N}^{*}$, respectively. The sequences $\{\phi_{k}\}_{k=0}^{N-3}$ and $\{\psi_{k}\}_{k=0}^{N-3}$ constitute a basis of $V_{N}$ and $V_{N}^{*}$, respectively.

Lemma 4.2**.**

The basis functions $\phi_{k}$ and $\psi_{k}$ satisfy

[TABLE]

where $\delta_{ij}$ is the Kronecker delta and $(\cdot,\cdot)$ is the usual $L_{2}$ inner product.

Proof.

By definition, we have

[TABLE]

The Legendre polynomials $L_{i}$ satisfy

[TABLE]

For $j>i$, we have $(\partial_{x}^{3}\phi_{i},\psi_{j})=0$ since $\partial_{x}^{3}\phi_{i}$ is a linear combination of Legendre polynomials whose highest degree Legendre polynomial is $L_{i}(x)$. The conclusion follows by the orthogonality property of the Legendre polynomials. On the other hand, for $i<j$, we have

[TABLE]

The only case left is $i=j$, where we have

[TABLE]

By using the relations (26), we obtain

[TABLE]

The only non-zero term is given by $k=i+1,\,j=k-1=i$, so we have

[TABLE]

By using the orthogonality property of the Legendre polynomials, one can easily prove

[TABLE]

which is the desired result. ∎

We remark the fact that the $L_{2}$ inner product $(p,q)$ coincides with the discrete inner product $(p,q)_{N}$ for all polynomials $p$, $q$ such that $\deg p+\deg q\leq 2N-2$. This is a property of the quadrature rule.

Since $w^{m}\in V_{N}$ we can express it in terms of $V_{N}$ basis functions

[TABLE]

Therefore, the variational formulation (25) with $v=\psi_{k}$ becomes

[TABLE]

for $j,k=0,\dots,N-3$.

Let us define the mass matrix

[TABLE]

and the stiffness matrix

[TABLE]

for $j,k=0,\dots,N-3$.

The discrete inner product $(\phi_{j},\psi_{k})_{N}$ coincides with the $L_{2}$ inner product $(\phi_{j},\psi_{k})$ for $j+k+6\leq 2N-2$ since $\phi_{k}$ and $\psi_{k}$ are polynomials of degree $k+3$. We can compute the entries $m_{j,k}$ using the $L_{2}$ inner product except for

[TABLE]

where the discrete inner product must be used. The entries $s_{j,k}$ can be computed using the $L_{2}$ scalar product only, since $\partial_{x}^{3}\phi_{j}$ is a polynomial of degree at most $N-3$ and $(\partial_{x}^{3}\phi_{j})\psi_{k}$ does not exceed the degree $2N-2$. The matrices $M$ and $S$ are seven-diagonal and diagonal, respectively. Therefore, the linear system

[TABLE]

where

[TABLE]

has a small bandwidth and is solvable in $O(N)$ operations.

Remark*.*

For a constant $g$ the variational formulation (25) becomes

[TABLE]

The matrix $D=(d_{j,k})$ given by $(\partial_{x}\phi_{j},\psi_{k})_{N}$ is not a banded matrix, but a lower triangular matrix as a straightforward computation shows. The linear system (28) is, therefore, expensive to solve.

Remark*.*

For arbitrary $g$ the linear system (27) can be rewritten as

[TABLE]

where $T$ is $m$ independent. The second term on the right-hand side is expensive to compute, since it involves $O(N^{2})$ operations. This number can be reduced to $O(N\log N)$ if the Chebyshev–Legendre dual-Petrov–Galerkin method is implemented, see [8]. Namely, we introduce $I^{C}_{N}$, the interpolation operator based on the Chebyshev–Gauss–Lobatto points, and we replace $g(x)u_{x}^{m-1}$ with I^{C}_{N}\big{(}g(x)u_{x}^{m-1}\big{)}. Now, we can use the fast Chebyshev–Legendre transform to compute the coefficients of the second term in $O(N\log N)$ operations. For more details we refer the reader to [16].

5 Implementation and numerical results

In this section, we discuss some implementation aspects. We also present numerical results that show the properties of the numerical scheme derived in the previous sections. In particular, we illustrate the super-polynomial spatial convergence of the employed dual-Petrov–Galerkin method. We show first-order convergence in time and the effect of the transparent boundary conditions on the numerical scheme. We also perform numerical simulations for different choices of the function $g$ and show that the transparent boundary conditions do not cause any reflections at the boundaries.

5.1 Stability

The numerical scheme (6) is derived by discretizing in time the dispersive part using the implicit Euler scheme and the advection part by the explicit Euler scheme. In the simpler case where we assume periodic boundary conditions and $g(x)=g$ constant we can easily show stability via Fourier analysis. In particular, we write $u^{m}$ as its Fourier series and we indicate the coefficients of the series with $(\hat{u}^{m}_{k})_{k\in\mathbb{Z}}$. We obtain

[TABLE]

The squared amplification factor given by

[TABLE]

is always less or equal than $1$ for $|g|\leq 1$, and can be bounded by $1+c\tau$ for $|g|>1$ with a constant $c$ depending only on $g$. We thus obtain

[TABLE]

which shows that the scheme is stable. We remark that using explicit Euler for the semi-discrete advection problem would be unstable, but the overall scheme is stable due to the use of the implicit Euler method for the dispersive part. Further, we remark that for $|g|\leq 1$ the problem is unconditionally stable. This simple analysis gives an indication that also for the problem considered in this work, where $g(x)$ is no longer constant and transparent boundary conditions are employed, the same might hold. The numerical experiments support the above observations.

5.2 Inverse $\mathcal{Z}$-transform

A numerical procedure to compute the inverse $\mathcal{Z}$-transform is given in [5]. Here we recall the main aspects. Given a sequence $(u^{l})_{l\geq 0}$ its $\mathcal{Z}$-transform reads

[TABLE]

To recover the $l$th element of the sequence starting from $\hat{u}$ one can use Cauchy’s integral formula

[TABLE]

where $S_{r}$ is the circle with center [math] and radius $r$. Approximating the integral by trapezoidal rule with $N$ points, we get

[TABLE]

where $\mathcal{F}^{-1}$ denotes the inverse discrete Fourier transform and $\mathbf{U}=\{U_{k}\}_{0\leq k\leq N}$ is the $N$-periodic sequence with $U_{k}=\hat{u}(re^{2\pi ik/N})$.

The radius $r$ and the number of points $N$ should be chosen in order to guarantee good approximations of the coefficients $u^{l}$. The choice of the radius $r$ has been studied in [15, 21]. Notice that $r>1$ must be sufficiently close to $1$ so that the inverse $\mathcal{Z}$-transform does not become numerically unstable when $l$ grows. A possibility to suppress this instability problem is to compute the inverse $\mathcal{Z}$-transform with a quadruple precision algorithm as proposed in [5]. A more effective remedy to overcome the instability has been proposed in [6]. There the authors exploit the relation between the coefficients of the polynomial (12) and its roots. For simplicity, we only consider the former approach in the present work. We set $r=\mathrm{e}^{C\tau}$ with appropriately chosen $C>0$ and $N=m\lceil|\log(10^{-7})|\rceil$, where $m$ is the number of time steps. A complete analysis of the optimal choice of the radius is not the goal of this work. For a thorough discussion we refer the reader to [15, 21].

The number of points $N$ discretizing (30) should be large enough in order to achieve good approximations of the coefficients $u^{l}$ in (31). It is known that the discrete Fourier transform using $N$ points is a good approximation of the Fourier transform for all the $u^{l}$ such that $|l|\leq N/2$. Therefore, $N\geq 2l$ is a good choice.

In formulae (18) convolution products appear for computing the transparent boundary conditions. The computation of such quantities becomes expensive as $l$ grows. Therefore, for sufficiently long time simulations, the computational cost will be dominated by the number of time steps. In order to keep the computational cost of the convolution sums to a minimum, one can approximate them, for example with a sum of exponentials approach. For a thorough discussion about the discrete transparent boundary conditions and the sum of exponential approach we refer the reader to [3, 5].

5.3 Numerical results

For the numerical tests we consider the problem

[TABLE]

subject to transparent boundary conditions. The spatial domain is chosen sufficiently large so that the initial data $u(0,x)$ is close to [math] for $x=\pm 6$, in particular $|u(0,\pm 6)|<10^{-15}$. Therefore, even if the initial data has not compact support, the effects produced by $u(0,x)$ at the boundaries are negligible.

As first example we consider problem (32) with

[TABLE]

The same problem is considered in [5, 20]. In this case the exact solution is available and given by

[TABLE]

where $*$ is the convolution on the entire real line and $\mathrm{Ai}$ denotes the Airy function. Notice that the spatial domain is $[-6,6]$ whereas the pseudo-spectral approach derived in Section 4 works only for $x\in[-1,1]$. Therefore, we first scale problem (32) to $[-1,1]$. Then, we compute the numerical solution to the scaled problem and finally we scale back to the original domain. In Fig. 1 we show the result of the numerical simulation. No reflections can be seen at the boundaries.

In Fig. 2 we present the numerical solution for

[TABLE]

This example is also considered in [5], where a numerical solution was provided by employing a finite difference spatial discretization. Our approach achieves very accurate spatial result using a modest number of grid points. Therefore, it gives an improvement compared to the finite difference approach, where a fairly large amount of spatial gird points are used. We notice the effect of the advection term that shifts the solution to the right. The exact solution can be obtained via Fourier transform, see [5]. Again no reflections can be observed at the boundaries.

In Fig. 3 we present the numerical solution for

[TABLE]

This example has a negative velocity $g$ and the exact solution can be obtained via Fourier transform. We notice that the advection term shifts the solution to the left, so no dynamic is happening at the right boundary. The solution leaves the domain through the left boundary without reflections.

Finally, we consider a third example where $g$ is no longer constant and given by

[TABLE]

The function $g$ is chosen in such a way that its extension to the whole domain $\mathbb{R}$ is smooth. Therefore, since $g$ is assumed to be constant outside the finite computational domain $[-6,6]$, we ask for $g_{x}(\pm 6)=0$. In Fig. 4 we plot the function $g$ as well as the numerical solution $u(T,x)$. We notice an advection for $x\in[-6,0]$ and a decay of $g$ for $x\in[0,6]$. One expects that the solution will be shifted to the right due to the advection term for $x\in[-6,0]$ and that there will be almost no effects in the second part of the domain, i.e., for $x$ close to $6$. This behaviour can be clearly seen in the plot, where, at the right boundary, the solution stays close to [math].

In Fig. 5 (left picture) we show the behaviour of the spatial error of the numerical solution compared to a reference solution for a fixed time $T=0.5$ and time step $\tau=2^{-14}$ for the three different choices of $g(x)$. As a reference solution we consider the numerical approximation employing $N=64$ grid points. In this setting the spatial error dominates the time error. Let

[TABLE]

be the relative $\ell^{2}$ spatial error computed at time $t_{m}=\tau m$ for $N\in[16,48]$. The grid points $x_{i}$ are the collocation points as in (24). From $\text{err}^{m}$ we compute the $\ell^{2}$ error in time as

[TABLE]

We observe in a semilogy plot the typical supergeometric convergence $\exp(-cN^{2})$ for analytic functions by spectral methods, see [7]. The error plot confirms the theoretical derivations of the previous sections. In particular, with less then $50$ spatial points we can achieve numerical solutions with an error less than $10^{-7}$.

In Fig. 5 (right picture) we show the time convergence for $g(x)=0$, $g(x)=\pm 6$, g(x)\leavevmode\nobreak\ =\leavevmode\nobreak\ \pi\big{(}1+\cos(\frac{\pi(x+6)}{12})\big{)} and $N=64$ space grid points. For the case where $g$ is constant, we consider the exact solution. For $g=g(x)$ we compute a reference solution employing $m=2^{14}$ time steps and compare the $\ell^{2}$ norm evaluated at $T=0.5$ for different $m=2^{7},\dots,\leavevmode\nobreak\ 2^{10}$. The order is one, as expected.