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

The present work proposes a second-order time splitting scheme for a linear dispersive equation with a variable advection coefficient subject to transparent boundary conditions. For its spatial discretization, a dual Petrov–Galerkin method is considered which gives spectral accuracy. The main difficulty in constructing a second-order splitting scheme in such a situation lies in the compatibility condition at the boundaries of the sub-problems. In particular, the presence of an inflow boundary condition in the advection part results in order reduction. To overcome this issue a modified Strang splitting scheme is introduced that retains second-order accuracy. For this numerical scheme a stability analysis is conducted. In addition, numerical results are shown to support the theoretical derivations.


Introduction
The aim of this paper is to develop time splitting schemes in combination with transparent boundary conditions that have spectral accuracy in space. Splitting schemes are based on the divide and conquer idea; i.e. to divide the original problem into smaller sub-problems which are, hopefully, easier to solve. However, obtaining an approximation of the solution of the original problem from the solutions of the sub-problems is not always straightforward: order reductions or strong CFL conditions that destroy the convergence of the numerical scheme are known to arise, see e.g. [8,10,18]. Furthermore, transparent boundary conditions are non-local in time and depend on the B M. Residori mirko.residori@tu-dresden.de 1 Department of Mathematics, University of Innsbruck, Innsbruck, Austria solution. Imposing them with splitting methods poses a challenge in the derivation of stable numerical schemes of order higher than one.
In this paper, we show that it is possible to construct a second-order splitting scheme that performs well, in the context outlined above, and can be implemented efficiently. In particular, we show that the proposed numerical method is stable independent of the space grid spacing (i.e. no CFL type condition is needed). We focus our attention on a linearised version of the Korteweg-de Vries equation where T > 0. The same ideas, however, can be applied to a more general set of linear partial differential equations with variable coefficients. Note that the partial differential Eq. (1), despite being linear, finds many applications in a physical context. For example, it is used to model long waves in shallow water over an uneven bottom, see e.g. [17,21]. The goal of this work is to design a splitting scheme that is second order in time with spectral accuracy in space. This paper can be seen as an extension to [12], where a splitting scheme of order one in time and spectral accuracy in space is presented. When solving (1) one of the main difficulties one has to face is the unbounded domain R. Numerical simulations typically consider a finite domain that leads to boundary conditions. Our goal is to design a numerical scheme that retains the same dynamics as the original problem (1), but on a finite domain. This can be achieved by imposing transparent boundary conditions. The advantage of such boundary conditions is the zero-reflection property of the solution at the boundaries. Further, the solution can leave the finite domain and re-enter at a later time without any loss of information. On the downside, transparent boundary conditions are non-local in time (and space for two and three-dimensional problems), therefore, they become expensive for long time simulations. In particular, memory requirements grow proportionally with the number of time steps. While it is still possible to employ them in 1D, the multidimensional cases become impracticable. A remedy is to approximate transparent boundary conditions and obtain so-called absorbing boundary conditions. In this way, information at the boundaries is lost, but memory requirements remain constant. A lot of work has been done for the Schrödinger equation in recent years, see [1,3,5] and references therein. For third-order problems, we refer the reader to [6,7,12,23] and references therein.
In the present case, the third derivative in space renders any explicit integrator extremely expensive. Therefore, an implicit scheme should be implemented. While coupling an implicit time discretization with a spectral space discretization yields banded matrices for constant advection, they lead to full matrices if g varies in space. We therefore employ a time-splitting approach in order to separate the advection problem from the dispersive problem. Operator splitting methods for dispersive problems have been employed and studied before, we refer the reader to [9,11,12,15]. For splitting method with absorbing boundary conditions we cite the work [5]. Splitting methods allow us to design specific solvers for the variable coefficient problem. For example, [20] uses a technique based on preconditioning. However, a direct splitting of (1) is not advisable. The problem of separating the advection equation is the potential requirement of inflow conditions at the boundaries. The actual inflow, however is unknown and should be estimated for example by extrapolation methods. This leads to instabilities when spectral methods are applied unless a very restrictive CFL condition is satisfied. The idea to overcome this problem is to perform a modified splitting that allows us to treat the advection problem without prescribing any inflow condition. The boundary conditions are transferred to the dispersive problem only. In this case, we can compute the values we need with the help of the Z-transform, as has been done for a constant coefficient dispersive problem in [6,7]. Another popular technique to avoid reflections at the boundaries is the perfectly matched layer method (PML). This method has been introduced in [4] for Maxwell's equations. Subsequently, it has been adapted to the Schrödinger equation [22] and very recently a general PML approach in combination with pseudo-spectral methods has been proposed in [2]. To the best of our knowledge a PML method for a linearised Korteweg-de Vries equation is currently not available.
The paper is organized as follows. In Sect. 2 we derive the semi-discrete scheme, discrete in time and continuous in space, by applying the Strang splitting method. In Sect. 3 we impose transparent boundary conditions for the scheme derived in Sect. 2. In particular, we determine the proper values of the numerical solution at the boundaries with the help of the Z-transform. The stability of the resulting numerical method is then analyzed in Sect. 4. In Sect. 5 we describe a pseudo-spectral method for the spatial discretization which takes the transparent boundary conditions into account. Finally, in Sect. 6 we present numerical results that illustrate the theory.

Time discretization: modified splitting approach
In this section we derive a semi-discrete scheme by applying the Strang splitting method to problem (1) restricted to a finite interval [a, b], where a < b. Inspired by the ideas in [12], we perform a time splitting in order to separate the advection problem ∂ t u(t, x) + g(x)∂ x u(t, x) = 0 from the dispersive problem ∂ t u(t, x) + ∂ 3 x u(t, x) = 0. In the following, for brevity, time and space dependence for the unknown u = u(t, x) are omitted.
In Sect. 2.1 we present the canonical splitting of (1). This approach illustrates the difficulty to prescribe the inflow condition to the advection equation. In Sect. 2.2 we then propose the modified splitting and show how this problem can be avoided.

Canonical splitting
Before applying any splitting, a preliminary analysis shows us that the inflow conditions to the advection problem depend on the sign of g(x) at x = a and x = b. We summarise in Table 1 the four possible outcomes. For this presentation, we restrict our attention to g(x) > 0 for x ∈ [a, b]. This setting requires an inflow condition at x = a. Let M ∈ N, M > 0 be the number of time steps, τ = T /M the step size and t m = mτ , k = 0, . . . , M. We apply the Strang Table 1 This table summarises at which boundary points {a, b} the inflow condition needs to be prescribed for the advection problem, depending on the sign of g(x) at the boundaries splitting method to (1), which results in the two sub-problems Let ϕ [1] t be the flow of (2a) and let ϕ [2] t be the flow of (2b). Let u(t, x) be the solution of (1) at time t. Then, the solution to (1) at time t + τ is approximated by the Strang splitting In order to get a numerical scheme, we apply the Peaceman-Rachford scheme to (3). This consists in computing the first flow ϕ [1] τ 2 by the explicit Euler method, the middle flow ϕ [2] τ by the Crank-Nicolson method and the last flow by the implicit Euler method. Let u m (x) = u(t m , x). Then, we get The latter numerical scheme is known to be second-order in time due to its symmetry. Notice that Eq. (5) is a time approximation of The function f (t) encodes the inflow condition at x = a. For t ∈ (0, τ ] the inflow condition is unknown. It can be approximated by extrapolation methods which typically leads to instabilities. The idea to overcome this problem is to formulate the advection problem without any inflow condition. For this purpose we introduce next a modified splitting.

Modified splitting
Based on the observations in Sect. 2.1, we rewrite the governing equation in (1) as follows where p g (x) is the line connecting the points (a, g(a)) and (b, g(b)). We now apply a splitting method that results in the two sub-problems Let ϕ [1] t be the flow of (8a) and let ϕ [2] t be the flow of (8b). Let u(t, x) be the solution of (1) at time t. The solution to (1) at time t + τ is then approximated by the Strang splitting u(t + τ, ·) ≈ ϕ [1] τ 2 • ϕ [2] τ • ϕ [1] τ 2 (u(t, ·)) .
By applying the Peaceman-Rachford scheme to (9), we get where g * (x) = g(x) − p g (x). Notice that g * (a) = g * (b) = 0. This means that no inflow or outflow condition needs to be prescribed to Eq. (11). The modified splitting allows us to solve the advection equation only for the interior points, i.e. x ∈ (a, b).

Remark
Notice that both problems (8a), (8b) have a variable coefficient advection. However, as shown in Sect. 5 the matrix associated to the space discretization of Problem (8a), despite the space dependent coefficient p g , is still banded. This is a property of the spectral space discretization that we employ.

Discrete transparent boundary conditions
When it comes to numerical simulations, a finite spatial domain is typically considered. Problem (1) is then transformed into the following boundary value problem Due to the third order dispersion term, three boundary conditions are required. In particular, depending on the sign of the dispersion coefficient, we have either two boundary conditions at the right boundary and one at the left boundary or vice-versa. In this work we consider a positive dispersion coefficient. We assume g(x) constant for x ∈ R \ [a, b] and that u 0 (x) is a smooth initial value with compact support in whereas in the interval [b, ∞) we consider the problem The initial value u(0, x) is set to 0 because u 0 (x) has compact support in [a, b]. The boundary conditions at x → ±∞ are set to 0 because we ask for u ∈ L 2 (R). Therefore, the solution u must decay for x → ±∞. We focus our attention on (14) and impose discrete transparent boundary conditions at x = a. A similar procedure can be applied to (15).
The mathematical tool we employ in order to impose discrete transparent boundary conditions to (13) is the Z-transform. We recall the definition and the main properties of the Z-transform, which are used extensively in this section. For more details we refer the reader to [3]. The Z-transform requires an equidistant time discretization.
Given a sequence u = {u l } l , its Z-transform is defined bŷ where ρ is the radius of convergence of the series. The following properties hold where * d denotes the discrete convolution Remark The Peaceman-Rachford scheme given in (10) Discretizing (14) by the Crank-Nicolson method, gives Let u(x) = {u k (x)} k be the time sequence (x plays the role of a parameter) associated to the Crank-Nicolson scheme (17). Then its Z-transform is given bŷ Taking the Z-transform of (17) gives where we used the time advance property of the Z-transform and u 0 (x) = 0. In particular, (18) is an ODE in the variable x. It can be solved by using the ansatẑ where λ i , i = 1, 2, 3 are the roots of the characteristic polynomial associated to (18): The roots λ i can be ordered such that Re λ 1 < 0 and Re λ 2,3 > 0, see [6]. By the decay conditionû(x, z) → 0 for x → −∞, we obtain c 1 (z) = 0 and Since c 2 and c 3 are unknown, the way to compute the discrete transparent boundary conditions is to make use of the derivatives ofû to derive an implicit formulation.
Computing the first and second derivatives ofû gives We then have In the latter equation the z dependence is omitted. The roots λ i , i = 1, 2, 3 satisfy This allows us to rewrite Eq. (20) in terms of the root λ 1 to obtain We can finally determine the value of u m+1 (a) by evaluating (21) at x = a and taking the inverse Z-transform. Let where we used the convolution property of the Z-transform. We remark that to compute u m+1 (a) we need to know ∂ x u m+1 (a) and ∂ 2 x u m+1 (a). Similarly, for problem (15), we obtain where with σ 1 root of The time discrete numerical scheme to problem (13) becomes (for 0 ≤ m ≤ M − 1) where h m+1 Equations (25a)-(25c) are the Peaceman-Rachford scheme. Equation (25d) is the initial data and Equations (25e)-(25g) are the discrete transparent boundary conditions.
where S r is a circle with center 0 and radius r > ρ, where ρ is the radius of convergence in (16). An exact evaluation of the contour integrals might be too complicated or infeasible. Therefore, we employ a numerical procedure in order to approximate these quantities. In this work we use the algorithm described in [13,Sec. 2.3], which results in stable and accurate results.

Stability of the semi-discrete scheme
For this section it is convenient to adopt a more compact notation. Thus, we write Then, the Peaceman-Rachford scheme (25a)-(25c) becomes The scheme can be rewritten separating the first step, i.e. when m = 0, as follows: Using the commutativity between I + τ 2 D 3 and I − τ 2 D 3 leads to We now show that the semi-discrete numerical scheme (26) is stable. The proof follows a similar approach as in [13].
Proof Let (·, ·) be the usual inner product on L 2 (a, b) and · the induced norm. We Applying the inner product with y m+1 + w gives or equivalently Integrating the right-hand side by parts gives Notice that ∂ x p g is constant since p g is a polynomial of degree 1. In order to complete the proof, a bound for w 2 is needed. By definition of w, we have Taking the inner product with w + y m gives Integrating by parts and using the fact that g * (a) = g * (b) = 0 gives Using the hypothesis τ < 4/ ∂ x g * ∞ leads to Combining (27) with (29) gives the bound where In the definition of B m we used p g (a) = g(a), p g (b) = g(b) and Multiplying both sides of (30) by 1 − τ 4 ∂ x g * ∞ and taking the sum over m gives By Lemma 4.1 the quantity B m is negative. Therefore, and stability follows by Gronwall's inequality since c 2 = O(τ ).
Inserting the discrete transparent boundary conditions (22) Let us extend the sequences B M a , B M b to infinity sequences by zero and apply Parseval's identity We obtain Notice that B M a and B M b are real values, therefore the imaginary parts of the right-hand sides in (31) must integrate to 0. Therefore,

Spatial discretization: pseudo-spectral approach
The spatial discretization of problem (25a)-(25g) is carried out by a dual Petrov-Galerkin method. In particular, we follow the approach given in [19] for the dispersive part and the approach given in [20] for the variable coefficient advection. It is very well known that pseudo-spectral methods achieve high accuracy even for a modest number of collocation points N , provided the solution is smooth. However, these methods have to be carefully designed in order to obtain sparse mass and stiffness matrices in frequency space. Then, the associated linear system can be solved in O(N ) operations.
In the following description we assume, without loss of generality, a = −1 and b = 1. The idea is to choose the dual basis functions of the dual Petrov-Galerkin formulation in such a manner that boundary terms from integration by parts vanish. Let us introduce a variational formulation for so that the discrete transparent boundary conditions are satisfied. To this goal, let P N be the space of polynomials up to degree N . For (33) and (35) we introduce the dispersive space The conditions in V d N collect the left-hand side of (25e)-(25g). Let (u, v) = b a u(x)v(x) dx be the usual L 2 inner product. The dual space V d, * N is defined in the usual way, i.e. for every φ d ∈ V d N and ψ d ∈ V d, * N it holds .
Proof Integrating ( p g ∂ x φ d , ψ d ) by parts and integrating (∂ 3 x φ d , ψ d ) by parts three times gives We want the boundary terms to vanish. For x = b we have The last equality is obtained by substituting p g (b) = g b and using the relations given by the space which leads to the boundary relations of the dual space V * N .
We proceed by introducing the advection space for (34): Notice that due to the variable coefficient g * the space V a N is free from inflow or outflow conditions. The dual space V a, * N is defined so that for every φ a ∈ V a N it holds for every ψ a ∈ V a, * N . The dual space V a, * N is P N . Let L j be the jth Legendre polynomial. We define where the coefficients α j , β j , γ j , α * j , β * j , γ * j are chosen in such a way that φ d j , ψ d j belong to V d N , V d, * N , respectively, see "Appendix A". The sequences {φ d j } N −3 j=0 and {ψ d j } N −3 j=0 are a basis of V d N and V d, * N respectively. We are now ready to consider the variational formulation.

Variational formulation
The dual Petrov-Galerkin formulation of (33) reads: find u * ∈ P N such that holds for every ψ d j ∈ V d, * N , j = 0, . . . , N − 3. In general the function u m does not belong to the space V d N . Indeed u m satisfies the discrete transparent boundary conditions However, we can write u m = u m h + p m 2 , where u m h ∈ V d N and p m 2 is the unique polynomial of degree two such that The function u * also does not belong to the space. Similarly, we can write u * = u * h + p * 2 . We assume that u * satisfies the same boundary conditions as u m , therefore p * 2 = p m 2 . We thus obtain We proceed with the dual Petrov-Galerkin formulation of (34). Find u m+1/2 ∈ P N such that holds for every ψ a j ∈ V a, * N , j = 0, . . . , N . Notice that u * , u m+1/2 ∈ V a N . So, differently from the dispersive case, we obtain the solution u m+1/2 without performing any shift.

Implementation in frequency space
This section is dedicated to compute the mass and stiffness matrices for (39)-(41). When it comes to numerical implementation, the L 2 inner product (u, v) needs to be approximated. We use two different discrete inner products for the spaces V d N and V a N . This choice is motivated by the fact that the spaces V d N and V a N satisfy different boundary conditions. Definition (Dispersive inner product) Let ·, · d N be the dispersive inner product defined as where y l are the roots of the Jacobi polynomial P (2,1) N −2 (y) and w l the associated weights. Definition (Advection inner product) Let ·, · a N be the advection inner product defined as where y l are the roots of the Jacobi polynomial P For more details about generalized quadrature rules, we refer the reader to [16].

Stiffness and mass matrices for
The first step is to obtain the frequency coefficientsũ m,d h,k in (44). We take the dispersive inner product on both sides The mass matrix is Using the orthogonality relation between φ d k and ψ d j gives φ d k , ψ d j d N = 0 if |k− j| > 3 and j +k ≤ 2N −8, see "Appendix B". Then, M d is a 7-diagonal matrix. Equation (45) in matrix form reads The left-hand side of (46) can also be written in matrix form: We obtain the frequency coefficients The second step is to compute the stiffness matrix and the frequency coefficients of the second term of the addition in (39). The stiffness matrix is

Lemma 5.2 S d is a 7-diagonal matrix.
Proof We now that φ d k is a polynomial of degree k + 3. Therefore, q(x) := is a polynomial of degree ≤ k + 3. We write q as a linear combination of Legendre polynomials up to degree k + 3: Let us consider the dispersive inner product q, ψ d j d N and k + 3 < j. Then, The last equation follows from the definition of ψ d j and the orthogonality property of the Legendre polynomials. Let j < k+3 with k+ j ≤ 2N −8, then (see "Appendix B") Similarly to q, we obtain φ d k ,q d N = 0 and the result follows.
The frequency coefficients of the second term on the right-hand side of (39) are given byp Notice that p g ∂ x p m 2 is a polynomial of degree 2. Therefore, it can be written as a linear combination of the Legendre polynomials L 0 , L 1 and L 2 . Using the orthogonality property of Legendre polynomials we obtain A similar procedure applies to (41), where we obtain Both linear systems (47) Similarly to the dispersive case, we need the frequency coefficients in (49), (50). We take the advection inner product in (49), (50) on both sides Using the orthogonality relation between φ a k and ψ a j gives φ a k , ψ a j a N = 0 if k = j. Then, the mass matrix with the stiffness matrix S a ∈ R (N +1)×(N +1) defined by The stiffness matrix S a is in general a full matrix. A direct inversion of (53) requires O(N 3 ) operations, thus is not advisable. Applying an iterative scheme is preferable, but multiplying the matrix S a with a vector costs O(N 2 ) operations. A more efficient way is to compute g * ∂ x u * (and g * ∂ x u m+1/2 ) in the physical space at the advection collocation points. The point-wise multiplication with g * costs only O(N ) operations. The result is then transformed back to the frequency space. Transforming back and forth to the frequency space can be done efficiently by employing the discrete Lagrange transform (DLT) and the inverse discrete Lagrange transform (IDLT) developed in [14], see "Appendix C".
Remark For the special case where g * is a polynomial of degree n, the stiffness matrix S a is banded with bandwidth less or equal to 2n. This implies that for a small n the linear system (53) is sparse and can be solved in O(N ) operations without switching from the frequency to the physical space.
Transition matrices In order to connect (53) to (47) and (48), it is necessary to transfer information from the dispersive space to the advection space and vice-versa. In particular, the aim is to translate the frequency coefficients from the dispersive space to the advection space in an efficient way. Let Then The coefficientsũ m+1/2,d can be directly obtained fromũ m+1/2,a in O(N ) operations.
Full discretization The implementation in frequency space results in The solution u m+1 (x) can be reconstructed by

Numerical results
In this section, we present numerical results that illustrate the theoretical investigations of the previous chapters. For that purpose, we consider with final time T = 1 and initial value u 0 (x) = e −x 2 . We restrict (61) to the interval (−6, 6) and impose transparent boundary conditions at x = ±6. The initial data is chosen such that |u 0 (±6)| ≤ 10 −15 .
For the numerical simulations, we employ a time discretization with constant step size and a space discretization given by the dual-Petrov-Galerkin variational formulation with N collocation points. We consider the error 2 of the full discretization defined as

Example 1 (Constant advection)
We consider (61) with constant advection g(x) = 6. This is the same problem which is considered in [6]. The setting reduces the advection equation in the modified splitting (8b) to the identity map. Even if the time-splitting is trivial, for this particular problem the exact solution can be computed via Fourier transform, see [6]. Consequently the constant advection problem offers a good benchmark for testing the convergence of the proposed numerical method in that context.
In Fig. 1 snapshots of the numerical solution u m N for t = 1 4 , 1 2 , 3 4 , 1 and τ = 2 −12 are shown. Notice that the numerical solution "leaves" the domain at the boundary x = −6 without any reflection. As time increases the solution moves to the right and re-enters the computational domain. Finally, the solution matches the boundary at x = 6 without any reflection.  Table 2 the full discretization error between the numerical solution and the exact solution varying N and M is reported. In particular, in Table 2 (left) the number of time steps M is fixed to 2 12 and the number of collocation points N is varying from 24 to 40. In this way the time discretization error is small enough to be negligible with respect to the spatial error. The value α denotes the slope of the line obtained by connecting two subsequent error values and varying N in a semi-logarithmic plot. More specifically, let N 1 and N 2 with N 1 < N 2 be two subsequent values of N and err 1 2 , err 2 2 the associated error values. Then Notice that α remains constant when N is varying, which confirms the spectral accuracy of the numerical scheme.
In Table 2 (right) the number of collocation points is fixed to 2 6 and the number of time steps M is varying from 2 5 to 2 8 . In this way the space error is small enough to be negligible with respect to the time error. The value β denotes the slope of the line obtained connecting two subsequent error values and varying M in a doublelogarithmic plot. More specifically, let M 1 and M 2 with M 1 < M 2 be two subsequent On the left side M is fixed to 2 12 so that the time error is negligible w.r.t. the spatial error. On the right side N is fixed to 2 6 so that the spatial error is negligible w.r.t. the time error. In both tables, errors are obtained testing u m N against the exact solution computed via Fourier transform as in [6]. The fact that α remains constant confirms the spectral accuracy of the proposed method, while the fact that β ≈ 2 confirms the second order in time values of M and err 1 2 , err 2 2 the associated error values. Then We can clearly see β ≈ 2, which confirms second order accuracy in time.

Example 2
We consider (61) with g(x) = −x 3 /54 + x + 3. As mentioned in Sect. 5.2, for g being a low degree polynomial, the stiffness matrix S a results in a banded matrix. Therefore, the linear system associated to the advection equation can be solved in O(N ) operations. The exact solution for this problem is not known, so we test the numerical solution u m N against a reference solution u m ref computed using a significantly greater number of points (both in time and space).
In Fig. 2 snapshots of the numerical solution u m N for t = 1 4 , 1 2 , 3 4 , 1 are shown. The solution is dragged to the right with an increasing speed. No appreciable reflections can be seen at the boundaries. Similarly to example 1, we report in Table 3  Example 3 We consider (61) with g(x) = e −(x+6) 2 + e −x 2 + e −(x−6) 2 − 1 2 . This example is interesting because of g is not polynomial and its sign alternates. The produced effects are a concentration of mass at the pointsx such that g(x) = 0, ∂ x g(x) < 0 and a thinning out where g(x) = 0, ∂ x g(x) > 0. Snapshots of the numerical solution that illustrate this phenomena are shown in Fig. 3. No reflections are detected at the boundaries, as expected.
Similarly to example 2, we report in Table 4 Table 4 (left) we observe a smaller value α with respect to Tables 2 and  3. Therefore, spatial convergence is slower with respect to examples 1 and 2, but still spectral accuracy is achieved. The slower convergence rate is related to the variations of the function g * , which are greater in magnitude than in examples 1 and 2.  In Fig. 4 we collect error plots for examples 1, 2 and 3. For all numerical tests we observe second order in time and the typical exponential convergence exp(−α N 2 ), α > 0 in space.
The numerical experiments confirm that the proposed approach performs well in the one dimensional case. However, the extension to higher dimensions is not straightforward. Transparent boundary conditions together with the pseudo-spectral discretization become more involved to compute. This poses a real challenge and is object of future studies.

A Finding the coefficients in (36)
The Legendre polynomials L j (x) satisfy the following orthogonality relation (L j , L k ) = δ jk 2 2 j + 1 .
Further, at x = ±1 we have where ( j) k = j( j + 1), . . . ( j + k − 1). Inserting the dispersive basis function φ d j given in (36) in the boundary relations of the space V d N leads to the following linear system for (α j , β j , γ j ) T : This means that all entries of M d k j except for ( j, k) = {(N − 4, N − 3), (N − 3, N − 4), (N − 3, N − 3)} can be analytically pre-computed. For the last three entries the discrete inner product defined in (42) must be used. This implies that for the last three entries the orthogonality relation between φ d k and ψ d j might not hold for the dispersive inner product. However, the bandwidth of the matrix will not change. A similar analysis applies for the stiffness matrix S d . For M ad we have Since 0 ≤ k ≤ N − 3 and 0 ≤ j ≤ N all entries can be computed analytically. A similar analysis applies for the advection mass matrix M a .

C DLT and IDLT
We recall briefly the definitions of DLT and IDLT. For more details we refer the reader to [14] . Given N + 1 valuesũ 0 ,ũ 1 , . . .ũ N the discrete Legendre transform is defined by where y k are the roots of the Legendre polynomial L N +1 (y). The inverse discrete Legendre transform computesũ 0 ,ũ 1 , . . . ,ũ N for given u 0 , u 1 , . . . , u N . It takes the formũ where w k , k = 0, . . . N are the Gauss-Legendre quadrature weights. (Notice that y k and w k are the same collocation points and weights as defined in the advection inner product (43)). Both the DLT and the IDLT can be computed in O(N (log N ) 2 / log log N ) operations, see [14]. For our application, let s * := IDLT(g * ∂ x u * ) ands m+1/2 := IDLT(g * ∂ x u m+1/2 ).
Clearly, to computes * (ands m+1/2 ) we need to reconstruct ∂ x u * (and ∂ x u m+1/2 ) starting from the frequency coefficientsũ * ,a (andũ m+1/2,a ). This can be done as follows. Note that we have The first relation is obtained by simply taking the derivative with respect to x in (49). The second relation comes from the fact that ∂ x u * is a polynomial of degree up to N and thus it belongs to V a N . Therefore, it can be written as a The choice of the test functions is motivated by the fact that the resulting matrices are both banded matrices with bandwidth two and three, respectively. In matrix form, (63) reads F Tũ * ,a = G T ∂ x u * a , from which we obtain ∂ x u * a in O(N ) operations. A similar procedure applies to the frequency coefficients of ∂ x u m+1/2 . Finally ∂ x u * and ∂ x u m+1/2 are obtained by applying the DLT to the corresponding frequency coefficients. Summarizing, we havẽ s * = IDLT g * · DLT G −T F Tũ * ,a , s m+1/2 = IDLT g * · DLT G −T F Tũm+1/2,a .
Thus, (40) can be solved in O(N (log N ) 2 / log log N ) operations for a general function g * .