Simulation of water wave equations has been a challenging problem due to the bad wellposedness of the equation. To understand this difficulty, we will first review previous remarkable theoretic work on water wave equations. In both Wu's work [28] and Beale and Hou's work [5], the authors used Riemann mappings to find the right variables and rewrote the water wave equation. Both versions of the equation in [28] and [5] have the common leading order structure: a non-local hyperbolic equation, say (1.6). By analyzing this simplified model, we derive an optimal discretization strategy in the form of a CFL condition. This condition is rigorously proved for system (1.6) and numerically verified for water wave equations. Following this road map, we will first derive the simplified model.
The unsteady system of incompressible free surface flow in two-dimension has attracted much theoretic and numerical attention [6,10,11,28]. Governed by the irrotational Euler equation, this free surface flow problem is also referred to as the water wave problem which dates back to the early 20th century [23,25]. By observing the equation in both Eulerian coordinate and Lagrangian coordinate, insightful analytical results were derived since then. One can refer to [11] for a review of recent related results. In [5], Beale, Hou, and Lowengrub formulated the water wave equation in Lagrangian coordinates and considered the linearization of it which is a nonlocal system. Later in [28], by reducing the system to a nonlocal hyperbolic equation in Eulerian coordinate, Wu derived impressive results on the well-posedness of the water wave equation. These two works directed us to focus on a simplified system that inherits the common dominating structure shared by both works: a nonlocal hyperbolic model, which played a critical role in the proof of well-posedness in [28].
The water wave equation is formulated as follows (see [5]) in Lagrangian coordinates. Consider a 2π-periodic two-dimensional fluid with infinite depth whose surface is described by z : R×[0,∞) → C: z(α,t) = x(α,t)+iy(α,t).
(1.1)
Here α∈R is a material coordinate that parametrizes the undisturbed surface. Periodicity of the fluid wave implies that s(α,t) := z(α,t)-α is a 2π-periodic function in α. Because the fluid is inviscid and irrotational, the velocity can be written as ∇Φ where Φ(x,y,t) is the velocity potential. Let φ : R×[0,∞) → R, (α,t) → φ(α,t) := Φ(x(α,t),y(α,t),t)
be the evaluation of the velocity potential at the surface so that φ(α,t) is a periodic function in α with period 2π. In [5], starting from the irrotational Euler equation, Beale, Hou, and Lowengrub derived the equation for the interface of fluid with infinite depth which is parametrized by z(α,t) in Lagrangian coordinate, i.e., z(α,t) and φ(α,t) satisfy the following system of equations:
2z α (α) =: w(α,t),
|w| 2 -gy,
where z means the complex conjugate of z, γ(α,t) quantifies the derivative of the dipole strength and g is the gravity. This derivation could also be found in [3,6] (Equations ( 1)-( 3) in [6]). Authors of [5] proved that the linearization of (1.2) leads to (1.3) (Equation (2.8) in [5]), namely for (α,t) ∈ R×(0,∞):
(
Here σ and c are positive, which depend on the solution of the water wave equation, but independent of η and ζ. η is the normal component of the perturbation of the position of the interface. δ is a combination of the tangential and normal components of the perturbation of the position. ζ describes the perturbation of the potential. g 1 and g 2 are some extra terms in the linearization which contain linear terms in η,δ. The operator Λ = (-∆) 1/2 = H∂ α (1.4) is the 1/2-fractional Laplacian with Fourier symbol |k|, where H is the Hilbert transform whose Fourier symbol is -isgn(k). On R, the Hilbert transform H is given by
x-y dy.
As we shall see, system (1.3) is L 2 stable and dispersive.
Another system with the same structure of (1.3) is derived in [28] in the Eulerian coordinates. Wu achieved remarkable results in [28] in proving the well-posedness of the water wave problem in Sobolev spaces. Wu used a conformal mapping formulation and reduced the water wave system to a quasi-linear hyperbolic system (see (4.6) and (5.8ǫ) in [28] and let w = -v) for (β,t) ∈ R×(0,∞) :
where σ > 0, c > 0 and β is the Eulerian coordinate instead of a material coordinate. Let X be the x-coordinate of the interface, then u = X tt and v = -X t in (1.5). Due to the extra time derivative, this equation can also be viewed as a linearization of the water wave equations in [5]. The usage of conformal mappings in this work was successful. It transformed all the nonlocal terms on a time-dependent domain into the half Laplacian, i.e. Λ = (-∆) 1/2 = H∂ β on R which is time-invariant. Wu referred to this system as the 'hyperbolic system' in [28]. We will also preserve this description in the present paper.
The slight difference between system (1.5) and (1.3) is that transport terms appear in the former but disappear in the latter. This is because: in (1.5), variable β is not the material coordinate but a variable associated with the conformal mapping. Except for this difference, both equation (1.5) and (1.3) share the same dominating structure, i.e. the system
for (x,t)∈R×(0,∞). This system is intrinsic to the water wave equation since it is detected in both Eulerian and Lagrangian coordinates. Motivated by the preceding discussion, we will focus on this nonlocal hyperbolic system in the rest of the paper.
As in [6], we impose periodic boundary condition and study the system on the torus, i.e.
with θ ∈T=R/2πZ and t∈(0,∞). The Hilbert transform H still has the symbol -isgn(k) but the formula now is given by
We will focus on equation (1.7) in the following sections. Consider the special case of (1.7) where σ,c are constant and
We derived much insight into the numerical simulation of the water wave from this case. In this constant coefficient case, the system is reduced to the following second -order (in time) nonlocal hyperbolic equation
where µ=σc. For heuristic purposes, we carry out some preliminary analysis and present the basic properties of (1.8) in Section 2.1. A more careful analysis for the constantcoefficient case is conducted in Section 3. Numerical studies of water waves have been performed in many papers [6,8,12,13,19,20,27]. The numerical methods can roughly be divided into two classes. In the first class, conformal mapping was not used. In [6], water wave problems were solved by the discretization of an integral formulation (see Section 5 for more information). However, the convergence was merely proved with time being kept continuous. The discussion of the fully discretized system seems challenging. In the second class [8,12,13,27], conformal mappings are used for numerical simulations but no rigorous numerical analysis for conformal mapping formulation has been performed. Meanwhile, although the analytical properties of the nonlocal hyperbolic system (1.6) are relatively well understood in [5,28], the numerical studies of such equations have not been thoroughly investigated.
Thus to shed light on the distinct properties of such hyperbolic systems and water wave simulations, we intend to focus on numerical analysis of the simplified model (1.6). Because (1.6) has nonlocal terms and the nonlocal terms have a simple Fourier symbol (say -isgn(k)), we select the pseudo-spectral approximation in the spatial discretization, which is often favored by wave equations (see [4,26]).
The primary goal of this paper is to analyze Runge-Kutta methods when applied to such nonlocal wave equations. In particular, we emphasize two insightful properties and implications that our analysis of Runge-Kutta methods provides:
First, we explore the optimal time step sizes in terms of a CFL type condition, when certain explicit Runge-Kutta methods are used. As we shall show in Section 2.1, the hyperbolic system (1.6) is also dispersive and may exhibit multi-scale behavior. Therefore, the time step constraint is more severe in the high-frequency regime. Consequently, finding optimal time steps with respect to the wave numbers is naturally desired [4,15,16,21].
In detail, we have systematically analyzed stability conditions of general Runge-Kutta methods for the hyperbolic system (1.3) with constant coefficients, including the highfrequency regime. In constant coefficient case, we have shown that, naive time discretization of system (1.7) results in the familiar hyperbolic CFL constraint ∆t = O(∆x). If we use Runge-Kutta methods whose absolute stable region contains a part of the imaginary axis, this CFL constraint can be relaxed to ∆t = O( √ ∆x) which is huge progress. See Theorem 3.1 for detail.
In high-frequency regime, this relaxation on CFL condition provides an optimal time discretization strategy which reduces computational load when simulating (1.7). In this regime, we consider the equation of u ′ (x,t) = u(ǫx,ǫt),v ′ (x,t) = v(ǫx,ǫt) which is the rescaling of (u,v). Still in the constant coefficient case, the equation of (u ′ ,v ′ ) is rewritten as
Due to the factor ǫ before v t , careful treatment of CFL conditions is necessary. In theorem 3.2, we conclude that Runge-Kutta schemes whose stability regions cover part of the imaginary axis are stable as long as both time step and spatial grid size resolve the wave oscillation, i.e.
This result is sharp in the view that one cannot capture the accurate wave function without resolving its oscillations. Second, as a corollary of the stability analysis, we investigated the energy stable property of Runge-Kutta methods. The key point is that the nonlocal hyperbolic system (1.8) with constant coefficients is energy preserving, here the energy is given by E = T |v t | 2 + 1 2 µvΛv dx. When applying pseudo-spectral approximation in spatial discretization, this Hamiltonian is still conserved by Parseval's equality. We will prove that, when applying Runge-Kutta methods in time discretization, the Hamiltonian is nonincreasing as long as the applied Runge-Kutta method has an absolute stable region that covers part of the imaginary axis. Available numerical discretization include explicit-k Runge-Kutta methods for k ≥ 3. See Corollary 3.1 for detail.
In the variable coefficient case, we discuss the extension of stability analysis from the constant-coefficient case and to the full water wave simulations. The proposed CFL conditions (say ∆t = O( √ ∆x)) are verified in numerical experiments. Notice that all preceding discussion is established for water waves of infinite depth, we will also consider an analog in finite depth case.
The rest of the paper is organized as follows: in Section 2.2, we introduce basic notations and the setup for the numerical analysis. In Section 3, we discuss thoroughly the discretization of the nonlocal system with Runge-Kutta (both explicit and implicit) methods in time and the Fourier spectral method in space. We then study the discretization of the system with variable coefficients in Section 4. We prove the convergence for the semi-discrete schemes using the Fourier spectral method or the filtered Fourier spectral method and then discuss the time discretization using Runge-Kutta methods. We then connect the nonlocal hyperbolic system to water wave equations in Section 5. Analog in water waves of finite depth is also discussed. Lastly, in Section 6, we perform numerical experiments. The stability conditions for the nonlocal hyperbolic system with variable coefficients and water wave equations are confirmed numerically. Numerical experiments suggest possible caustics for the system in the high-frequency regime.
In this section, we discuss the special case (1.8) and then introduce necessary notations for numerical analysis later.
In this section, we present a concise review of basic properties of (1.8), a special case of the hyperbolic system (1.3). Multiplying by u t on both sides of (1.8), and integrating over x, we derive
This means that the energy
is conserved in time. Because v also satisfies equation (1.8), we know that
is also conserved. To derive the dispersion relation, suppose that the plane wave u(x,t)= Ae i(2πkx-wt) is the solution to the equation (1.8). On the Fourier side, (1.8) reads as ûtt = -µ|ξ| û.
(2.2)
Notice that the Fourier transform of plane wave u(x,t) = Ae i(2πkx-wt) is û(ξ,t) = Aδ(kξ)e -iwt , substituting it into the equation yields the dispersive relation:
-
Because ω is real and ω = const×ξ, the system is dispersive. Next, we derive the Green function of (1.8). The explicit expression of the Green function also suggest the dispersion relationship. In fact, the solution u(x,t) to the initial value problem
(2.4) can be written as u(x,t) = f (x) * F(x,t)+g(x) * G(x,t).
(2.5)
Here * represents convolution in space, i.e. f (x) * F(x,t) = R f (x-y)F(y,t)dy. Function G(x,t) is the Green's function and F(x,t) is the time derivative of it. In fact, G(x,t) and F(x,t) satisfy following Cauchy problems respectively:
Remember that the nonlocal operator Λ has Fourier symbol |ξ|, thus G and F are respectively given by
Here, F -1 is the inverse Fourier transform, i.e. (F
For the sake of completeness, we provide the details and more discussions of the Green's function in Section A. Consider the special case where f (x) = cos(kx),g(x) = 0. Then the solution of (2.4) is u(x,t) = cos(kx)cos( µ|k|t). Thus, the angle velocity in space and time has the squareroot relation in k: k and |µ|k. This can be explained by the dispersion relation (2.3). In fact, this also suggest a relaxed CFL condition as we shall explain later.
Remark 2.1. Eq. (1.8) is reminiscent of the surface quasi-geostrophic equations (SQG) studied in [7,17]. However, the surface SQG equation is dissipative while (1.8) is dispersive.
In this work, we consider the one-dimensional nonlocal hyperbolic system (1.7) on T = R/2πZ. The spatial discretization is selected as the Fourier pseudo-spectral method or the filtered Fourier pseudo-spectral method.
We discretize the spatial domain with grid size h = 2π/N, and we denote grid points by
where N ∈N is even. We denote the time step size by τ, and denote t n = nτ. The notation u n j represents the numerical value of u(θ,t) at (θ j ,t n ), and u n represents the vector u n = (u n j ). Given any N-vector f = ( f j ), we expand each component as a sum of discrete Fourier modes via
2 N} and the discrete Fourier transform f = ( fk ) are given by fk = 1
Note that the Hilbert transform H and differentiation operators become certain multipliers when the Fourier transform is applied. When projected onto a uniform grid, those transforms between two functions reduce to corresponding relations between the discrete Fourier transforms of two functions confined on the grid. We define the projected differential operator and the projected Hilbert transform H in the following. For two N-vector f and g, we write
)
We introduce the notation L = DH as the projected Λ = ∂H, so that
Recall the discrete inner product between two N-vectors is defined as
where ḡ means the complex conjugate. The discrete ℓ 2 and ℓ ∞ norms are defined by
For discrete case, we still have Parserval's equality:
Lemma 2.1. The discrete Parserval's equality holds
Consider the constant coefficient case of the simplified hyperbolic model, i.e.
where (θ,t) ∈ T×(0,∞). On the Fourier side, this equation reads as
In this section, we will thoroughly analyze the stability of Runge-Kutta methods applied to (3.2). First, we will conduct the Von Neumann analysis [18] on Runge-Kutta methods. Both explicit and implicit ones will be systematically investigated. As a result, we will derive stability conditions in terms of CFL conditions. These CFL conditions provide necessary guidance for the simulation of water wave equations which is conducted in Section 6. Eq. (3.2) in the high-frequency regime is also considered. In the interest of avoiding aliasing error [24], an optimal discretization strategy is developed as a consequence of the analysis of Runge-Kutta methods. The strategy is optimal in the sense that it resolves oscillations at the lowest requirement. In what follows, we denote the Butcher tableau of a certain n-step Runge-Kutta method by
where G is the Runge-Kutta matrix, w are the weights and p are the nodes.
Direct calculation shows that A has 2 complex eigenvalues λ 1,2 =± i cσ|k|. Thus, matrix A is similar to the diagonal matrix D := diag{λ 1 ,λ 2 }:
where
This similarity relation reduces the coupled system (3.2) into the following decoupled system:
Hence, the Von Neumann analysis on the original coupled system 3.2 reduces to that on each scalar equation in (3.4). Notice that P 2 has entry σ|k|/c which contains Fourier variable k, so L 2 -stability of (3.4) does not directly imply L 2 -stability of (3.2). However, this fact does not affect stability of Runge-Kutta methods on (3.2). As we will explain later, one can still ensure stability of certain Runge-Kutta methods of (3.2) under a weighted norm, which actually implies a weaker stability in L 2 -norm. Therefore, we first consider stability analysis of (3.4). Notice that the 2 eigenvalues of matrix D have the same absolute value, so we can focus on the stability analysis for only the first equation which is exactly the ODE system y ′ = λ 1 y.
In the following lemma, an explicit formula of the growth index is derived when a certain Runge-Kutta method is employed. This facilitates the Von Neumann analysis.
Here τ is the time step size, G the Runge-Kutta matrix, w the weights in the Runge-Kutta method (see (3.3)) and ν = cσ|k|.
The proof merely requires standard techniques in numerical analysis textbooks, hence attached to Appendix B for the sake of completeness.
An important case of this lemma is for explicit Runge-Kutta methods. For explicit Runge-Kutta methods, matrix G is a lower triangular matrix, thus det(I +τ 2 νG 2 ) = 1. Then formula (3.5) reduces to
This implies that ψ 2 (τ,ν) is a constant coefficient polynomial of τ 2 ν which depends on G,w. This fact can also be derived from properties of explicit methods. If one uses explicit RK-p method to discretize the first component of system (3.4), the growth index in Lemma 3.1 has the following expression [18]:
Moreover, the absolute stable region of RK-p method contains a part of imaginary axis, i.e. [-iC 1 (p),iC 1 (p)] if p ≥ 3 [18]. Thus requiring
is sufficient to ensure the stability of explicit RK-p methods when solving system (3.4). Denote the growth matrix for a certain Runge Kutta method by B(τ). According to [18], we recall the following definitions. Definition 3.1. A scheme is called weakly stable if at least for τ sufficiently small, there exists α >0 such that |B(τ)| ≤1+ατ holds. A scheme is called strongly stable if |B(τ)| ≤1 holds at least for τ sufficiently small. Here |B(τ)| is the norm of matrix induced by the norm |•| defined on R n .
One can see that even with weak stability, the numerical scheme is convergent as τ → 0 for any fixed time T > 0. Thus, we use (3.7) and Lemma 3.1 to prove the following theorem which provides the stability condition: Theorem 3.1. There are 2 cases in this theorem:
1. For any Runge-Kutta scheme, if there exists a positive real number C such that τ Ch holds, then the scheme is stable for system (3.4) or (3.2). More specifically, (a) if tr(G 2 ) > tr (G-ew T ) 2 , then the scheme is strongly stable;
(b) if tr(G 2 ) tr (G-ew T ) 2 , then the scheme is weakly stable.
Here e is the vector of ones and G,w is defined in (3.3).
then the scheme is strongly stable.
Proof. We first prove statement (a) and (b) in 1. Recall that ν = cσ|k| and |k|h π, thus we have
Therefore, τν=O(1) and τ 2 ν=O(τ) as τ→0. Using Lemma 3.1 and the fact det(I +τX)= 1+tr(X)τ +O(τ 2 ), we have
Then we consider 2 cases respectively:
1. If tr(G 2 ) > tr (G-ew T ) 2 , then ψ(τ,ν) < 1 when τ goes to 0, which means strong stability. This proves statement (a).
) tr (G-ew T ) 2 , by the estimate (3.10), we have ψ(τ,ν) 1+ Lτ +O(τ 2 ) when τ goes to 0, where
This implies weak stability which proves statement (b).
For the third statement, by estimate (3.10) again, we have τ
then τ √ ν stays in the absolutely stable region of the numerical scheme which indicates strong stability. Remark 3.1. In the proof of Theorem 3.1, we proved stability in the weighted norm P 2 • , which actually does not directly imply stability in the L 2 -norm. In fact, in this case, we still have stability in the L 2 -norm by Theorem 3.1: notice that |k| ≤ √ π/ √ h, so there exists constants C 1 and C 2 such that for all k ≥ 1,
for the mode k = 0, any numerical scheme is always stable since (3.1) reduces to ∂ t u = 0,
= -cu(0)t, which is linear. So the mode k = 0 is always stable, no matter which scheme is utilized. Therefore, if the total error of a numerical scheme is O(h n ), n ≥ 1 in norm P 2 • , then the total error will be O(h n-1/2 ) in the L 2 -norm.
We have to emphasize that statement (a) and (b) in the preceding theorem are not sharp. This is because that they merely ensure zero stability of a scheme while some of them can be even unconditionally absolutely stable. Despite of this theoretic unsharpness, these two statements are still general and practical because they are conclusive for any Runge-Kutta method. Moreover, they provide necessary guidance for numerical simulations of water wave equations in Section 6.
To illustrate this theorem, we analyze two examples of Runge-Kutta methods: a weighted Euler method which is semi-implicit, and the explicit RK-4 method. Both methods are stable under the sufficient condition τ Ch as stated in the first part of the theorem. Nevertheless, this CFL condition is unnecessary for the former method under certain weights. For the second method, the condition τ C √ h claimed in the second statement of Theorem 3.1 is observed.
(a) Semi-implicit RK2: weighted Euler method. The tableau of a weighted Euler method is:
Suppose that τ ≤Ch for some positive constant C, one can verify that ψ(τ,ν)≤1+C 1 τ holds as τ → 0 for some constant C 1 . Thus by Theorem 3.1, this scheme is stable.
However, if δ 1 2 , then ψ(τ,ν) 1 holds without any requirement on τ,h. This implies unconditionally strong stability. This example illustrates that condition τ Ch is sufficient but not necessary.
(b) Explicit RK4. The tableau of explicit RK4 is:
, where
Remember that the absolute stable region of explicit RK-4 contains a part of the imaginary axis, thus by Theorem 3.1, we expect to acquire a CFL condition of form τ C √ h. Substitute G,w T into ψ(τ,ν) and denote z = τ 2 ν, we have
Here c,σ is the constant in the system (3.2). Thus one can take
πσc in Theorem 3.1 to ensure stability.
A corollary of Theorem 3.1 is that the discretization of energy E (see (2.1)) is stable for any explicit Runge-Kutta method whose absolute stable region contains part of the imaginary axis. Remember that in proving Theorem 3.1, we focused on vector ( û′ , ŵ) in (3.4). Also notice that û′ 2 2 + ŵ 2 2 = P 2 (u,v) T 2 and
Lv,v , so the stability analysis in Theorem 3.1 is actually established on the L 2 -norm of (u ′ ,w) (or P 2 (u,v) T ), which is a combination of Ḣ1/2 -norm of v and L 2 -norm of u. This L 2 -norm is exactly the discretization of energy
with a constant factor on the Fourier side by substituting v t = -cu. Thus strong stability implies that the energy is non-increasing which is exactly the following corollary:
Corollary 3.1. For any RK method whose absolute stable region contains a part of the imaginary axis [-C 1 i,C 1 i], there exists a positive constant C 2 such that when
the discretization of the energy E in (2.1) is non-increasing, i.e.
is non-increasing.
In fact, one has
So by the second statement in Theorem 3.1, we know that the scheme is strongly stable and E 1 is non-increasing. This proves the corollary.
To investigate the behavior of high frequency waves in (3.1), we consider the rescaling (x,t) → (ǫx,ǫt) which yields the following system:
where T 1 = R/Z. Here ǫ = ℓ/K is a small number, ℓ is the considered length scale and K ≫ 1 is the typical wave number. Eq. (3.20) is equivalent to the following equation of second order in time:
Again, this rescaled system is Hamiltonian: the energy T 1 ǫ|u t | 2 + 1 2 µuΛu dθ is conserved. To see this, one can multiply by u t on both sides of (3.21) and integrate over θ to derive
On the Fourier side, the first order system reads as
Analysis of RK methods for this system can be conducted in the same way as in Section 3.1 if we replace c by c/ǫ. This concludes the following theorem of stability: Theorem 3.2. There are 2 cases in this theorem:
1. For any Runge-Kutta scheme, if there exists a positive real number C such that τ ǫCh holds, the scheme for system (3.22) is stable. More specifically, (a) if tr(G 2 ) > tr (G-ew T ) 2 , then the scheme is strongly stable;
(b) if tr(G 2 ) tr (G-ew T ) 2 , then the scheme is weakly stable.
Here e is the vector of ones and ψ(τ,ν) is defined in (3.3).
[-C 1 i,C 1 i], there exists a positive constant C 2 such that when
the scheme is strongly stable.
Energy is also non-increasing if the scheme is strongly stable, i.e., Corollary 3.2. For any RK method whose absolute stable region contains a part of the imaginary axis [-C 1 i,C 1 i], there exists a positive constant C 2 such that when
the discretization of the energy E in (2.1) is non-increasing, i.e.
An implication of Theorem 3.2 is that formula (3.23) theoretically provides guidance on selecting step size when simulating the behavior system (3.22). A notorious difficulty in numerical simulation of high-frequency waves is the spatial aliasing error [24] and its consequential strict requirement of stability. To avoid aliasing error, one needs to ensure h = O(1/K) [24]. Recall that ǫ = ℓ/K where ℓ is the considered length scale, thus h = O(1/K) = O(ǫ). As in Theorem 3.2, although τ Cǫh can ensure stability, if this CFL condition is rigorously obeyed, the numerical simulation would suffer from a heavy load of computation since the time step size τ satisfies
This small time step size indicates that O(K 2 ) times of computation will be conducted which could be unacceptable. However, still in Theorem 3.2, we have proved that τ Cǫh can be improved to τ C 2 √ ǫh for typical Runge-Kutta schemes whose stable regions cover a part of the imaginary axis. In this case, the requirement on τ is reduced from
Hence, both time steps and spatial grid sized only need to resolve the wave oscillation to ensure stability, i.e.:
This result is sharp in the view that one cannot capture the accurate wave function without resolving its oscillations. Therefore, Theorem 3.2 and Eq. (3.23) direct design of numerical schemes of (3.20) by validating an optimal stability condition. Moreover, we expect that (3.23) can be employed in simulations of water waves, which is in fact realized in Section 6.
In this section, we consider system (1.6) with variable coefficients. We will prove the convergence of the semi-discretization and provide some evidence for the validity of the CFL condition, namely τ C √ h in the variable-coefficient case. However, we have to admit that rigorous proof of convergence for fully discrete approximation is not given, this would be considered in our further work. Nevertheless, careful numerical experiments are conducted to verify the convergence and stability of τ C √ h for both system (1.6) and the water wave equation. Thus, our analysis is still meaningful in the sense that it enables us to shed light on the convergence and stability conditions of the numerical simulation of the water wave equation.
We assume that all the coefficients are smooth on T, σ ≥ σ 0 > 0 and c(θ,t) ≥ c 0 > 0. Consider the energy functional
Taking the derivative and using (1.7), we find
Compared with the energy in constant-coefficient case (2.1), the additional v,v term assists on estimating the energy: it is needed to control the linear terms like
hold, thus these terms can be controlled by
By Young's inequality, we have
The Gr önwall inequality yields
This implies that E is bounded in finite time. This is consistent with the hyperbolicity.
In this subsection, we discretize spatial variables first and then show that the semidiscretization is convergent. As mentioned in Section 2.2, we use the Fourier pseudospectral method for spatial discretization. However, for the variable-coefficient case, we desire to damp high-frequency modes to gain sufficient smoothness so that some desired properties hold (see Lemma 4.4 and Section 4.3 below). Therefore, we introduce the filter function ρ on the Fourier side [6]. We use N-vectors U h and V h to approximate u(τ,t) and v(τ,t) respectively. Let σ,g 1 ,g 2 ,c be restricted to the grid points. Given a filter function ρ: (-π,π]→R, we denote the operator with symbol ρ h (k) = ρ(hk) by ρh , i.e.,
Then, we have the filtered version of action of operators on N-vectors:
The hugest advantage of this filter function is its generalization and flexibility. First, it includes some typical finite difference schemes. For instance, the centered difference on torus (D c u) j = 1 2h (u j+1 -u j-1 ) can be regarded as a filtered Fourier differentiation with filter ρ(ξ) = sin(ξ) ξ . Second, if a filter is unnecessary, one can simply set ρ = 1. We will assume the following conditions for the filter function: • ρ ≥ 0, even and ρ ∈ C 2 (-π,π] (Note that ρ may not be C 2 on torus).
• There exists r ∈ N + such that sup ξ∈(0,π)
Because ρ is non-negative, we can then define the natural discrete Sobolev norms associated with ρ to be
From Lemma 2.1, we have following properties:
Lemma 4.1. Suppose f ,g are two N-vectors. Then integration by parts formulas hold:
These formulas are guaranteed by Parserval's equality. In fact:
Other equalities can be similarly checked and we omit the details.
With the filter function, we discretize the system in space with the filtered pseudo Fourier spectral method, while keeping time continuous:
Here σ,c,λ 1 ,λ 2 , f 1 , f 2 are also discretized in space, but continuous in time, i.e., they are evaluated at the grid points.
To prove convergence, we first check the consistency of this discretization. By Fourier analysis and the aliasing formula, we have the following lemma: Lemma 4.2. Let ϕ ∈ C ∞ (T) and N ∈N. Then, the restriction f = ( f j ) = (ϕ(θ j )) of ϕ to the grid points satisfies
where R i : T → R (i = 1,2) are functions with |∂ α θ R i (θ,h,r)| bounded uniformly in θ and h, for any α ∈ N.
Proof of this lemma can be found in [26]. As a corollary of Lemma 4.2, we have the following consistency result which is direct:
) and the filter satisfies (4.4). Let U e = (u(θ j ,t)), V e = (v(θ j ,t)), i.e. the restriction of exact solutions on grids. Then we have
where R i (•,•;h)h r (i = 3,4) are the local truncations errors and R i (•,•;h) (i = 3,4) are two smooth functions on T×[0,T] with W α,∞ norms uniformly bounded in h for any α ∈ N.
Now we show the convergence of the semi-discretized equations (4.6).
Proposition 4.1. Consider (1.6) with σ ≥ σ 0 > 0 and c ≥ c 0 > 0. Assume that all the coefficients are smooth. Let the exact solution to (1.6) be (u,v) ∈ C ∞ (T×[0,T]). Let r be the constant in (4.4). Let (U e ,V e ) be the restriction of the exact solution to grid and (U h ,V h ) be the numerical solution given by the pseudo-spectral method (4.6) with the same initial values. Then there exists a constant M(T) > 0, such that ∀t ∈ [0,T]:
Proof. Define the error vectors
Taking the difference of Eqs. (4.6) and (4.8), we find the error functions satisfy the following equations
Consider the energy functional for this ODE system (analogy to (4.1))
). Therefore,
According to Eq. (4.11),
In the first estimate, we used the fact that L ρ R 4 is uniformly bounded by the smoothness of the error (by Lemma 4.3). The last term d e v ,e v dt is straightforward:
By Gr önwall's inequality, we finally obtain
which leads to our estimate for the error directly.
In this section, we aim to study the spatial operators on the right-hand side of (4.6). In particular, we are interested in the eigenvalues of the operators, which will shed light on the time discretization of the ODE system (4.6). The strategy is similar as that in Section 3, i.e. we introduce a similar transformation given by P(t) so that the stiff part of the operator scaling as 1 √ h becomes anti-symmetric. For the convenience of further discussion, we introduce the notion of smoothing operators which is an analogy to the big-O notation introduced in [6]: Theorem 4.1. Suppose that the filter satisfies Condition 4.1 and ρ(π) = 0. Then, we can decompose the operator A(t) (Eq. (4.16)) as
(recall h = 2π/N), where the linear operators A 1 (t
(ii) The eigenvalues of P(t) -1 A 1 (t)P(t) are purely imaginary whose norms are bounded by a constant C > 0 independent of N. The eigenvalues of P(t) -1 A 2 (t)P(t) are bounded by a constant C > 0s independent of N.
Proof. We consider the operator B(t) whose domain is Q, defined by
Then, it is given by:
We then define A 1 (t) as
and
We can directly verify that the ranges of A i are in E N ×Q. Moreover, A 1 is bounded and anti-symmetric. Now we focus on A 2 . Note that
To see this, let us consider the continuous version of the equations, and consider the same energy functional (4.1). One can estimate Ė similarly as before except for term Λv,b∂ θ v . To estimate this term, we find
where [H,b] = Hb-bH is the commutator. Note that b is smooth. Thus the commutator [H,b]∂ θ v gives a convolution between a smooth function and
Therefore, Λv,b∂ θ v can also be bounded by the
Unfortunately, the same estimate does not work again for the discrete case. A filter function is necessary for a stable energy. In fact, for the discretized Hilbert transform, [H,b] is not smooth in general. In [6], the authors found that [H,b] may not even be A -1 . By Lemma 4.4, one needs a filter ρ so that the commutator [H,b] ρ = A -2 , which has the smoothing effect to ensure that
(ii) On the other side, conventional numerical treatments on the transport terms require the CFL condition of the form
This condition is much more restrictive compared with (3.9), which could be resolved using semi-Lagrangian method [16,21].
In Section 1, we explained that the nonlocal hyperbolic systems (1.6) and (4.19) are closely related to water wave equations with infinite depth. In this section, we provide sufficient insight into how our results of the nonlocal system imply the stability conditions for the simulation of water wave equations.
Recall the linearization of (1.2), which leads to (1.3). Indeed, this linearization happens for the numerical schemes as well. In [6], the authors proposed a filtered pseudospectral differentiation method to discretize the spatial variables. On the basis of condition 1, they also assumed that the filter satisfies (i) ρ(π) = 0 and ρ ′ (π) = 0; (ii) there exists r ≥ 4 such that |ρ(ξ)-1| ≤ C|ξ| r for ξ ∈ (-π,π]. Let j ∈ [N] and (z j ,φ j ,γ j ) be the numerical solutions at the grid points. Then, the discretization to (1.2) is given by (Eqs. ( 7)-( 9) in [6])
(5.1)
To control the numerical error, the authors in [6] introduced the following functions:
ζ j = (φ j (t)-φ(α j ,t))-Re(w j (z j (t)-z(α j ,t)).
(5.2)
Moreover, they assumed the strong Taylor sign condition (Eq. ( 88) in [6]) c(α,t) := -∂ n p ≥ c 0 > 0 and that the true solutions are smooth with
Based on these assumptions, the authors in [6] found that these variables for errors satisfy the following semi-linear non-local hyperbolic system (Eqs. ( 89)-( 91)):
(5.3) See Definition 4.1 for A 0 . The leading order behavior of the semi-linear system (5.3) is exactly (1.3), the nonlocal hyperbolic system. Making use of this hyperbolic structure, the authors showed that the semi-discrete system (5.1) converges to the original water wave problem (1.2). However, there was no discussion about time discretization. Therefore, given the same linearization effect on both continuous equations and numerical schemes, we expect that similar stability conditions hold for both (1.6) and the water wave equation. Because the errors satisfy the leading order equation of the nonlocal hyperbolic system (1.3) or (1.6), previous discussion in (4.2) convinced us that the In the variable-coefficient case, because Λ 1 can be regarded as the composition of the Hilbert transform and the pseudo-derivative with symbol ktanh(H 0 |k|), the commutator estimates in Theorem 4.1 will not change either. This new commutator is smooth at k = 0 but it does not help us to damp the high frequency for controlling aliasing errors. Hence, all the conclusions will be similar to the infinite depth case. We expect that when using the energy stable Runge-Kutta methods to solve them, the CFL condition is again
This will be left for future numerical study.
In this section, we present some simulations to verify our conclusion and carry out careful numerical experiments. We will introduce our results in the following way: in Section 6.1 and 6.4, we utilize the explicit RK-4 method for the temporal discretization to discretize both the simplified nonlocal hyperbolic system (1.6) and the water wave equation. We verify that the stability condition agrees with (3.9) for both systems which are compatible with our previous analysis. Convergence of the discretization of the nonlocal hyperbolic system is demonstrated in Section 6.2 and exploration of the nonlocal system in the highfrequency regime is performed in Section 6.3. Moreover, a turn-over wave example in [6] is recovered to verify the correctness of our simulation program run in Section 6.4. For all simulations in this section, periodic boundary conditions are selected. The (filtered) Fourier spectral method is conducted for spatial discretization.
In this example, we test stability conditions for the nonlocal hyperbolic system (1.6) with g 1 = 0, g 2 = 0. We consider both the constant-coefficient case with c = 3, σ = 1 and the variable-coefficient case with c(θ,t) = exp(cos(θ +t)), σ(θ,t) = 2+sin(θ +t).
We perform simulations for various step sizes using Fourier spectral method in space and forward Euler (FE) and Runge-Kutta 4 (RK4) for temporal discretization. The numerical solutions are calculated up to T = 10. Results are presented in Fig. 1 where the blue part represents the unstable region while the yellow part represents the stable region. The stability in both Section 6.1 and Section 6.4 was determined by the amplitude at a certain breaking point T 0 , before the terminal time T. For a certain group of temporal step size and grid size, If the L ∞ norm of the solution at T 0 is larger than a threshold (or even diverges, say 'NaN'), then it is determined as instable. Otherwise, it is stable. In the top half of Fig. 1, we plot stability in the √ h-τ plane. Notice that the borders for RK4 are more flat and similar to lines. This plot indicates that the stability condition for RK4 is really (3.9). Meanwhile, the borders for FE are some convex curves similar to a parabola. Thus the stability condition for FE should be the linear CFL condition τ ≤ Ch.
In order to illustrate this point clearly, we zoom in on the figures and re-plot the stable region in h-τ plane. To check the condition for FE, we re-plot the variable-coefficient case as shown at the bottom of Fig. 1. The new figure shows that the stability condition for FE is τ ≤ Ch. In all cases, the stable region of RK4 is much larger than that of FE, which indicates that RK4 is more stable than FE.
These results verify our analysis in Sections 3 and 4. Recall that the stability region for FE only intersects the imaginary axis at z = 0, while RK4 contains some part of the imaginary axis. Therefore, by Theorem 3.1, the CFL condition for RK4 should be (3.9); the linear CFL condition sufficiently ensures stability of FE.
In this subsection, we verify the convergence numerically for the nonlocal hyperbolic system (1.6). For transport terms, we take g 1 = 0, g 2 = 0. The constant coefficients are given by c = 3, σ = 1.
The initial conditions are given by u 0 (θ) = e sin(θ) +cos(θ), v 0 (θ) = cos 2 (θ).
For the temporal discretization, forward Euler (FE), backward Euler (BE), Crank-Nicolson (CN) and Runge-Kutta 4 (RK4) are used.
All simulations are computed up to T = 2. The reference solution (or 'accurate solution') is computed using Runge-Kutta 4 with h = 2π/2 7 and τ = 10 -5 . The error plots are shown in Fig. 2. In Fig. 2(a), spectral accuracy is clearly observed in spatial discretization: when h ≈ 0.2, the errors have already been dominated by the temporal error. Fig. 2(b) indicates that the temporal errors are of the order as expected. Therefore, our discretization schemes indeed converge.
In this section, we investigate the system in the high-frequency regime. In particular, we aim at verifying whether there is a caustic phenomenon that is similar to [4].
In research of high-frequency waves, a WKB kind initial value is typically considered. Meanwhile, we rescale the system as in Section 3.3. Therefore, we consider (3.21) with a selected initial value:
e ilog(20cosh(5x-2.5))/ǫ , u t (x,0) = e -100(x-0.5) 2 e ilog(20cosh(5x-2.5))/ǫ . (6.1)
The initial value consists of a Gaussian function and a high-frequency term. The former is used to control the support and the latter is a WKB type function. In numerical experiments, we select µ = 1 in (6.1) and different rescaling factors: ǫ = 2 -i , i = 4,5,••• ,12. We plot the snapshot of the amplitude at t = 0.0625 for ǫ = 2 -4 ,2 -6 ,2 -8 ,2 -10 in Fig. 3. The figure indicates that the amplitude increases when ǫ decreases. Remember that the initial value u(x,0) is of an amplitude no greater than 1 for all different ǫs', thus this growing trend of the amplitude provides evidence for caustic phenomenon. To obtain a more careful observation, we check the maximum amplitude before time t = 0.0625 in the whole domain [0,1], and plot the following log-log figure in Fig. 4. From Fig. 4, we observe that when ǫ is sufficiently small, the curve is almost a line whose slope is 1, which indicates that the maximum amplitude is approximately proportional to 1/ǫ. This observation also supports the existence of the caustic phenomenon.
In this section, we simulate the water wave equation (Eq. (1.2) with α∈T). As in previous simulations, the spatial discretization is implemented using the filtered Fourier spectral method. The filter function we use in this section is the same as in [6, Section 6]: ρ(ξ) = exp(-10(|ξ|/π) 25 ), ξ ∈ (-π,π]. (6.2) This filter function numerically satisfies the condition ρ(π) = 0 and ρ ′ (π) = 0: |ρ(π)| and |ρ ′ (π)| are sufficiently small. To verify that our simulation programs, we first recover the same example in [6, Section 6] with turn-over phenomenon. The initial data are given by:
x(α,0) = α, y(α,0) = 0.6cos(α), γ(α,0) = 1+0.6sin(α). In the simulation, the spatial grid size is h = 1/512 and the time step size is τ = 1/4000. we use RK4 for time discretization here. The snapshots of the waves at different times are shown in Fig. 5. As one can see in the figure below, the same numerical results in [6] are recovered. Now, we numerically investigate the stability condition for the discretization for the water wave equation. We use the following initial data:
x(α,0) = α, y(α,0) = 0.3cos(α), γ(α,0) = 1+0.3sin(α). (6.4) The spatial discretization is performed using the filtered Fourier spectral method with the same filter in Eq. (6.2). For temporal discretization, we select the FE method and the RK-4 method. The simulations are performed up to time T = 4.
The results are presented in Fig. 6. Same as in Section 6.1, the blue part represents the unstable region while the yellow part represents the stable region. Also, partially zoomed-in versions of the figure are presented in the case of FE to clearly check the stability condition. From this figure, one can tell that the stability condition for RK4 is close to τ √ h, while the stability condition for FE is more similar to τ h. This observation is in accordance with our conclusions in Section 5 and previous analysis.