Numerical Analysis and Approximate Travelling Wave Solutions for a Higher Order Internal Wave System

ABSTRACT

In this work we focus on the numerical solution of a higher order bidirectional nonlinear model of Boussinesq type involving a nonlocal operator. Based on a von Neumann stability analysis for the linearized problem, an efficient and stable scheme for the nonlinear system is proposed. Our method is based on a numerical scheme known from the literature that solves satisfactorily a lower order linear system. Additionally, approximate periodic travelling wave solutions profiles for the higher order nonlinear system are presented. Such approximate travelling wave solutions are obtained from a solitary wave family of solutions for the Intermediate Long Wave (ILW) equation and the regularized Intermediate Long Wave (rILW) equation.

Keywords:
spectral method; dispersive models; stability analysis; travelling waves

1 INTRODUCTION

Internal ocean waves are gravity waves that appear in stratified fluids. Stratification is a consequence of variations in density due mainly to differences in temperature and salinity. Abrupt variations in density justify the use of a layered model, the simplest of which is considered here and consists of a two-layer fluid configuration limited by a rigid lid at the top and a flat bottom. Each layer contains an inviscid, incompressible, irrotational fluid of constant density. The two fluids are immiscible and of different densities, the denser one is located below. An internal wave propagates at the interface between layers and the whole system evolves according to the Euler equations together with appropriate boundary conditions as detailed in ^{8), (15}. This simplified model retains the main features that enable the study of travelling wave solutions which reproduce the behaviour of well-identified disturbances that move with effectively constant speed for long periods of time as described in ^3),(10. Internal ocean waves are relevant for marine and submarine human activity. For example, the dead water phenomenon experienced by vessels when they travel through a relative thin layer of fresh water that does not mix with a denser layer below affects nautical operations, as originally reported by Nansen in ¹³. Internal waves are essential in ocean dynamics because of the amount of mass and energy transported by them, ultimately producing wave breaking and mixing and the distribution of nutrients among other kind of matter. For a broad discussion about internal waves, not limited to their presence in oceans, see ¹⁷.

Asymptotic analysis of the Euler equations is a successful method for the study of internal ocean waves. Several internal wave models described in the literature are derived by this process, for example, the ILW equation in the intermediate regime and the Benjamin-Ono (BO) equation in the deep water regime ^{5), (11}. For the case of intermediate depth for the lower layer and shallow upper layer, a strongly nonlinear model for internal waves was obtained in ^{15), (16}. It describes the evolution of the interface $\eta (x,t)$ between the fluids and the upper layer averaged horizontal velocity $u(x,t)$, where x and t represent the spatial and temporal variables, respectively. The coordinate system is set at the undisturbed interface between layers. Subscripts x and t stand for partial derivatives. In the flat bottom case and in nondimensional variables that system reads:

$\left\{\begin{array}{l}{\eta }_t={\left[(1-\eta )u\right]}_x,\\ u_t+u{u}_x+\left(1-\frac{{\rho }_2}{{\rho }_1}\right){\eta }_x=\sqrt{\beta }\frac{{\rho }_2}{{\rho }_1}{\mathcal{T}}_{\delta }{\left[(1-\eta )u\right]}_{xt}+\\ +\frac{\beta }{3(1-\eta )}{\left((1-\eta {)}^3\left(u_{xt}+u{u}_{xx}-u_x{u}_x\right)\right)}_x+\\ +\beta \frac{{\rho }_2}{{\rho }_1}{\left(\eta \left(\right(1-\eta )u{)}_{xt}+\frac{1}{2}((1-\eta )u{)}_x^2\right)}_x+\\ +\beta \frac{{\rho }_2}{{\rho }_1}{\mathcal{T}}_{\delta }{\left[\eta {\mathcal{T}}_{\delta }\left[\right(1-\eta )u{]}_x\right]}_{xt}+\frac{\beta }{2}\frac{{\rho }_2}{{\rho }_1}{\left({\mathcal{T}}_{\delta }{\left[\left(\right(1-\eta )u{)}_x\right]}^2\right)}_x.\end{array}\right.$ (1.1)

The thickness $h_1>0$ of the upper layer is much smaller than the characteristic wavelength $L>0$ at the interface. The thickness $h_2>0$ of the lower layer is comparable with L. The densities of the upper and lower fluid are denoted by ρ ₁ and ρ ₂, respectively. For a stable stratification, let ${\rho }_2>{\rho }_1>0$. Figure 1 illustrates this configuration. The nondimensional dispersion parameter $\beta =(h_1/L{)}^2$ is small. For numerical purposes we focus on the case where solutions are periodic functions in space with period 2l. In this case the nonlocal operator 𝒯_δ is defined by the symbol $\widehat{{\mathcal{T}}_{\delta }\left[f\right]}\left(k\right)=\mathit{i}\coth \left(\delta k\right)\hat{f}\left(k\right),k\in {\mathrm{\mathbb{Z}}}^{*}$, where δ = h ₂ /L and i denotes the imaginary unit. For more details about the operator 𝒯_δ see ¹⁵ or ¹.

Figure 1:
Two fluids configuration.

The second equation in system (1.1) is an approximation of order β ^3/2 while the first equation is exact, this means that it is a direct consequence of the Euler equations and no approximation from the asymptotic expansion was made ^{15), (16}. System (1.1) is a generalization of the model derived in ⁸ using a higher order asymptotic expansion that includes terms of order β which are not present in the model considered in ⁸.

Moreover, we can compare the dispersion relations of these models. In dimensional forms, being g the gravitational acceleration, the dispersion relation ω _h of system (1.1), the dispersion relation ω _l of the system proposed in ⁸ and the full dispersion relation ω _f that comes from the Euler equations are given, respectively, by

$$\begin{gathered} {\omega }_h^2=\frac{g({\rho }_2-{\rho }_1)k^2}{\frac{{\rho }_1}{h_1}+\frac{1}{3}{h}_1{\rho }_1{k}^2+{\rho }_{2}k\coth \left(k{h}_2\right)}, \\ {\omega }_l^2=\frac{g({\rho }_2-{\rho }_1)k^2}{\frac{{\rho }_1}{h_1}+{\rho }_{2}k\coth \left(k{h}_2\right)}, \\ {\omega }_f^2=\frac{g({\rho }_2-{\rho }_1)k^2}{{\rho }_{1}k\coth \left(k{h}_1\right)+{\rho }_{2}k\coth \left(k{h}_2\right)}. \end{gathered}$$

Noting that ${\rho }_{1}k\coth \left(k{h}_1\right)=1+\frac{(h_1{\rho }_1{)}^2}{3}+\mathcal{O}\left((h_1{\rho }_1{)}^4\right)$, we can see that ω _h is a better approximation to ω _f than ω _l . That is, system (1.1) approximates better the dispersive effect of the problem than the system proposed in ⁸.

Also, we can obtain the conservation of mass law for system (1.1) as follows: integrating the first equation of (1.1) on x and using that u is periodic on the x-variable, there exists an arbitrary constant d such that

${\int }_{-l}^l\eta (x,t)dx=d,$ (1.2)

which is valid for any value of t in a time interval where the solution is defined. As η comes from a perturbation of the interface at rest and the coordinate system is set at the undisturbed interface we must set $d=0$.

Considering a weakly nonlinear wave propagation regime, the nonlinearity parameter α, which is the ratio between the typical absolute wave amplitude value and the thickness h ₁, is introduced by scaling $\eta =\alpha {\eta }^{*},u=\alpha u^{*}$ where it is imposed that $\alpha =\mathcal{O}\left(\beta \right)$. Thus, omitting the asterisks and gathering the terms with order $\mathcal{O}\left(\alpha \sqrt{\beta }\right),\mathcal{O}\left(\alpha \beta \right)$ and $\mathcal{O}\left({\beta }^{3/2}\right)$, the following weakly nonlinear system with normalized shallow water velocity is obtained:

$\left\{\begin{array}{l}{\eta }_t-\left[\right(1-\alpha \eta )u{]}_x=0,\\ u_t+\alpha u{u}_x-{\eta }_x=\sqrt{\beta }\frac{{\rho }_2}{{\rho }_1}{\mathcal{T}}_{\delta }{\left[u\right]}_{xt}+\frac{\beta }{3}{u}_{xxt}+\mathcal{O}\left(\alpha \sqrt{\beta },\alpha \beta ,{\beta }^{3/2}\right).\end{array}\right.$ (1.3)

Besides, its linearized version around the zero equilibrium is

$\left\{\begin{array}{l}{\eta }_t=u_x,\\ u_t-\frac{{\rho }_2}{{\rho }_1}\sqrt{\beta }{\mathcal{T}}_{\delta }[u{]}_{xt}-\frac{\beta }{3}{u}_{xxt}={\eta }_x.\end{array}\right.$ (1.4)

A study about the existence and uniqueness of solution for both systems in periodic Sobolev spaces was presented in ⁶. Note that the conservation law (1.2) is also valid for systems (1.3) and (1.4).

As we see in ^{6), (15} the phase velocity $\omega /k=v\left(k\right)$ of system (1.4) is given by

$v\left(k\right)={\left(1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }\varphi \left(k\delta \right)+k^2\frac{\beta }{3}\right)}^{-1/2},$

where ϕ is defined as

$\varphi \left(k\right)=\left\{\begin{array}{ll}1,& k=0,\\ k\coth k,& k\ne 0.\end{array}\right.$ (1.5)

Thus, system (1.4) in the frequency domain can be written as

$\left[\begin{array}{c}{\hat{\eta }}_t\left(k\right)\\ {\hat{u}}_t\left(k\right)\end{array}\right]=\left[\begin{array}{cc}0& \mathit{i}k\\ \mathit{i}k{v}^2\left(k\right)& 0\end{array}\right]\left[\begin{array}{c}\hat{\eta }\left(k\right)\\ \hat{u}\left(k\right)\end{array}\right]=A\left(k\right)\left[\begin{array}{c}\hat{\eta }\left(k\right)\\ \hat{u}\left(k\right)\end{array}\right],k\in \mathrm{\mathbb{Z}},$

where matrix A(k) has pure imaginary eigenvalues ${\lambda }_k=\mathit{i}kv\left(k\right)$. Therefore, any numerical method must preserve this property of the spectrum. Because of this, in order to implement the method of lines, a fourth order finite difference scheme for spatial derivatives and a spectral approach for the dispersive terms are considered in the semi-discretization. The classical fourth order Runge-Kutta (RK4) scheme is used for time advancing. This combination proved to be the best since along the imaginary axis (where the eigenvalues of the spatial discretization operators must lie) the RK4 method has the largest stability interval if compared with the fifth order, four steps AdamsMoulton scheme and the fourth order, three steps Adams-Moulton scheme, see for instance ⁴.

The numerical method presented here for the linear system (1.4) is based in the one proposed in ¹ for the system

$$\begin{gathered} \left\{{\eta }_t=u_x, \\ {u}_t-\frac{{\rho }_2}{{\rho }_1}\sqrt{\beta }{\mathcal{T}}_{\delta }[u{]}_{xt}={\eta }_x,\right. \end{gathered}$$ (1.6)

which is the linearization around the zero equilibrium of the nondimensional form of the system proposed in ⁸ and differs from system (1.4) in the term involving u _xxt . The results of the von Neumann analysis for both systems are analogous, however the stability conditions obtained in this work are less restrictive than those presented in ¹. The improvement of the stability conditions is due to the term involving u _xxt which is also responsible for the better approximation of the Euler dispersion relation. Thus, the higher order linearized system has better physical and numerical properties.

The von Neumann analysis for the chosen method is performed together with a comparison with other schemes for spatial derivatives in semi-discretization, namely spectral differentiation and piecewise B-splines. The stability conditions are validated in numerical tests and extended to the scheme for the nonlinear system (1.3) which includes the discretization of the nonlinear terms $\alpha (\eta u{)}_x$ and αuu _x as presented in ¹². This allows us to study numerically the existence of approximate travelling waves.

This work is organized as follows. In Section 2 discretizations for both linear and nonlinear systems are presented including three possibilities for spatial derivatives and the von Neumann analysis for each method is performed. In addition, we present a comparison between stability conditions of the numerical methods for systems (1.4) and (1.6). In Section 3 the stability conditions obtained in the previous section are exemplified in numerical tests and extended for the nonlinear case. A convergence analysis for both methods is made in Subsection 3.1 and the solitary wave profiles of the ILW and rILW equations are tested as approximate periodic travelling wave solutions for our nonlinear system in Subsection 3.2. Conclusions appear in Section 4.

2 DISCRETIZATION AND STABILITY ANALYSIS

For the discretization of the linearized system (1.4) let us define the auxiliary function

$\psi =u-\frac{{\rho }_2}{{\rho }_1}\sqrt{\beta }{\mathcal{T}}_{\delta }[u{]}_x-\frac{\beta }{3}{u}_{xx},$ (2.1)

and rewrite system (1.4) as

$\left\{\begin{array}{l}{\eta }_t=u_x,\\ {\psi }_t={\eta }_x.\end{array}\right.$

Since u and η are 2l-periodic with respect to the variable x, we define a uniform grid on the interval [0, 2l], that is, $x_j=j\Delta x,j=1,\dots ,N$ where $N\in \mathrm{\mathbb{N}}$ is even and $\Delta x=2l/N$, the last element x _N is also identified with $x=0$ Define $\mathit{u}\left(t\right)=[u_1,\dots ,u_N{]}^T$ and $\mathit{\eta }\left(t\right)=[{\eta }_1,\dots ,{\eta }_N{]}^T$ where $u_j\approx u(x_j,t)$ and ${\eta }_j\approx \eta (x_j,t),j=1,\dots ,N$. The spatial discretization results in a system of Ordinary Differential Equations in the matrix form below

$\left\{\begin{array}{r}{\mathit{\eta }}_t=C\mathit{u},\\ P{\mathit{u}}_t=C\mathit{\eta }.\end{array}\right.$ (2.2)

Matrix C comes from the discretization of the first order x-derivative and matrix P is obtained using the Discrete Fourier Transform (DFT) in order to get $\mathit{\psi }=P\mathit{u}$ with $\mathit{\psi }\left(t\right)=[{\psi }_1,\dots ,{\psi }_N{]}^T$ where ψ _j approximates the expression (2.1) evaluated at $(x_j,t)$.

Defining the Fourier matrix F componentwise by $F_{m,j}=e^{(-2\pi \mathit{i}(m-N/2)j/N)},1\le m,j\le N$, the DFT of a vector $\mathit{w}=[w_1,\dots ,w_N{]}^T\in {\mathrm{ℝ}}^N$ is denoted by the vector $\hat{\mathit{w}}=[{\hat{w}}_{-N/2+1},\dots ,{\hat{w}}_{N/2}{]}^T\in {\mathrm{ℂ}}^N$ which is defined by the relation $\hat{\mathit{w}}=\Delta xF\mathit{w}$ and its inverse is given by $\mathit{w}=\frac{1}{2l}{\overline{F}}^T\hat{\mathit{w}}$.

Now, using that the composition of one spatial derivative with the operator 𝒯_δ has the symbol defined by equation (1.5) and some properties of the Fourier series, we can calculate $\hat{\psi }\left(k\right)=v(k\pi /l{)}^{-2}\hat{u}(k),k\in \mathrm{\mathbb{Z}}$. Defining the diagonal matrix $\hat{P}=\mathrm{diag}\left({\lambda }_{-N/2+1}\left(P\right),\dots ,{\lambda }_{N/2}\left(P\right)\right)$ whose entries are given by ${\lambda }_k\left(\hat{P}\right)=v(k\pi /l{)}^{-2},k=-N/2+1,\dots ,N/2$ we get $\hat{\mathit{\psi }}=\hat{P}\hat{\mathit{u}}$, thus

$\mathit{\psi }=\frac{1}{2l}\overline{F}{}^T\hat{\mathit{\psi }}=\frac{1}{2l}\overline{F}{}^T\hat{P}\hat{\mathit{u}}=\frac{1}{2l}\overline{F}{}^T\hat{P}\Delta xF\mathit{u}=\frac{1}{N}\overline{F}{}^T\hat{P}F\mathit{u}=P\mathit{u}.$ (2.3)

For the first order x-derivative, as in ¹, we choose numerical schemes such that $u_x\left(x_j\right)\approx (C\mathit{u}{)}_j$ where C is a real, skew-symmetric and Toeplitz circulant matrix. These assumptions are made in order to preserve the property of the spectrum described earlier and the periodicity of the spatial domain.

A Toeplitz circulant matrix C is defined by the elements $c_1,\dots ,c_N$ of its first row by

$C_{i,j}=\left\{\begin{array}{ll}{c}_{1+j-i},& i\le j,\\ c_{N+1+j-i},& i>j.\end{array}\right.$

It is shown in ¹ that the skew-symmetric Toeplitz circulant matrix C is diagonalized by matrix $(1/\sqrt{N})F$ and its eigenvalues are given by

${\lambda }_k\left(C\right)=\frac{\mathit{i}}{\Delta x}\gamma \left({\theta }_k\right),{\theta }_k=\frac{2\pi k}{N},\gamma \left(\theta \right)=2\Delta x\sum _{m=1}^{N/2-1}{c}_{1+m}\sin \left(m\theta \right).$

Therefore, from the diagonal matrix $\hat{C}=\mathrm{diag}\left({\lambda }_{-N/2+1}\left(C\right),\dots ,{\lambda }_{N/2}\left(C\right)\right)$ we recover $C=\frac{1}{\sqrt{N}}\overline{F}{}^T\hat{C}\frac{1}{\sqrt{N}}F$.

Three possibilities for matrix C are studied that correspond to each of the following discretizations for the first order spatial derivative: five points fourth order finite difference, piecewise linear B-splines and the spectral scheme. As shown in ¹, their respective functions γ are given by

${\gamma }^{\text{FD}}\left(\theta \right)=\frac{4}{3}\sin \theta -\frac{1}{6}\sin 2\theta ,{\gamma }^{\text{BS}}\left(\theta \right)=\frac{3\sin \left(\theta \right)}{2+\cos \left(\theta \right)/4},{\gamma }^{\text{SP}}\left(\theta \right)=\left\{\begin{array}{l}\theta ,\theta \in (-\pi ,\pi )\\ 0,\theta =\pm \pi .\end{array}\right.$

To integrate in time let us rewrite system (2.2) as

$\left[\begin{array}{c}{\mathit{\eta }}_t\\ {\mathit{u}}_t\end{array}\right]=D\left[\begin{array}{c}\mathit{\eta }\\ \mathit{u}\end{array}\right],D=\left[\begin{array}{cc}\mathbf{0}& C\\ P^{-1}C& \mathbf{0}\end{array}\right].$ (2.4)

Denoting the eigenvalues of P and D by ${\lambda }_k\left(P\right)$ and ${\lambda }_k\left(D\right),k=-N/2+1,\dots ,N/2$, respectively, it is not difficult to prove that ${\lambda }_k\left(D\right)=\pm {\lambda }_k\left(C\right)/\sqrt{{\lambda }_k\left(P\right)}$ since matrices C and P are diagonalized by the orthogonal matrix $(1/\sqrt{N})F$. Since ${\lambda }_k\left(D\right)$ is purely imaginary, the classic fourth order Runge-Kutta method is appropriate for this problem. Applying the RK4 method in system (2.4) we get

$\left[\begin{array}{c}{\mathit{\eta }}^{n+1}\\ {\mathit{u}}^{n+1}\end{array}\right]=\left(I+\Delta tD+\frac{\Delta t^2}{2}{D}^2+\frac{\Delta t^3}{3!}{D}^3+\frac{\Delta t^4}{4!}{D}^4\right)\left[\begin{array}{c}{\mathit{\eta }}^n\\ {\mathit{u}}^n\end{array}\right],$ (2.5)

where η ⁿ and u ⁿ are the approximations for $\mathit{\eta }\left(t_n\right)$ and $\mathit{u}\left(t_n\right)$, respectively.

Once the discretization of the linear system is complete, let us obtain the discretization of the nonlinear system. For that, let us rewrite system (1.3) as

$\left\{\begin{array}{r}{\eta }_t=E_1(\eta ,u),\\ {\psi }_t=E_2(\eta ,u),\end{array}\right.$

where $E_1(\eta ,u)=u_x-\alpha \eta u_x-\alpha u{\eta }_x$ and $E_2(\eta ,u)={\eta }_x-\alpha u{u}_x$.

Considering the spatial mesh defined before, the x-derivatives in E ₁ and E ₂ are approximated using one of the numerical schemes described previously, that is, five points finite difference, piecewise linear B-splines and the spectral scheme. The variable ψ is calculated by equation (2.3) as in the linear case. Integration in time is done with the fourth-order Runge-Kutta method. Note that in the problem that motivated this study α is of the same order of β, but setting $\alpha =0$ we obtain scheme (2.5).

2.1 Von Neumann Analysis

In order to prove that the numerical scheme (2.5) is stable, let us define the following norm

$\left|\right|\mathit{w}|{|}_{N,s}^2=2\pi \sum _{k=-N/2+1}^{N/2}{\left(1+(k\pi /l{)}^2\right)}^s|{\hat{w}}_k{|}^2,s\ge 0,$

where $\mathit{w}\in {\mathrm{ℝ}}^N$ and $\hat{\mathbf{w}}$ is its DFT. Our aim is to prove that there exists a positive constant C _s which does not depend on n and ∆t such that

${||[{\mathit{\eta }}^n,{\mathit{u}}^n{]}^T||}_{N,s,s+1}\le C_s{||[{\mathit{\eta }}^0,{\mathit{u}}^0{]}^T||}_{N,s,s+1},\forall n\Delta t\le T,$ (2.6)

where ${||[{\mathit{w}}_1,{\mathit{w}}_2{]}^T||}_{N,s,s+1}^2=\left|\right|{\mathit{w}}_1|{|}_{N,s}^2+|\left|{\mathit{w}}_2\right|{|}_{N,s+1}^2$, and both vectors w ₁ and w ₂ belong to ℝ^N .

Substituting matrix D in scheme (2.5) by its blocks described in (2.4) and using that $(1/\sqrt{N})F$ diagonalizes matrices C, P and P ⁻¹, we proceed as in ¹ to transform scheme (2.5) to the Fourier space:

$\left\{\begin{array}{l}{\hat{\mathit{\eta }}}^{n+1}=\left[I+\frac{\Delta t^2}{2}{\hat{P}}^{-1}{\hat{C}}^2+\frac{\Delta t^4}{4!}{\hat{P}}^{-2}{\hat{C}}^4\right]{\hat{\mathit{\eta }}}^n+\left[\Delta t\hat{C}+\frac{\Delta t^3}{3!}{\hat{P}}^{-1}{\hat{C}}^3\right]{\hat{\mathit{u}}}^n,\\ {\hat{\mathit{u}}}^{n+1}=\left[\Delta t{\hat{P}}^{-1}\hat{C}+\frac{\Delta t^3}{3!}{\hat{P}}^{-2}{\hat{C}}^3\right]{\hat{\mathit{\eta }}}^n+\left[I+\frac{\Delta t^2}{2}{\hat{P}}^{-1}{\hat{C}}^2+\frac{\Delta t^4}{4!}{\hat{P}}^{-2}{\hat{C}}^4\right]{\hat{\mathit{u}}}^n.\end{array}\right.$ (2.7)

For each frequency $k,k=-N/2+1,\dots ,N/2,$, system (2.7) returns

$\left[\begin{array}{c}{\hat{\eta }}_k^{n+1}\\ {\hat{u}}_k^{n+1}\end{array}\right]=\left[\begin{array}{cc}\text{c}({\theta }_k,\sigma ,\Delta x)& \mathit{i}{v}^{-1}\left(\frac{{\theta }_k}{\Delta x}\right)\text{s}({\theta }_k,\sigma ,\Delta x)\\ \mathit{i}v\left(\frac{{\theta }_k}{\Delta x}\right)\text{s}({\theta }_k,\sigma ,\Delta x)& \text{c}({\theta }_k,\sigma ,\Delta x)\end{array}\right]\left[\begin{array}{c}{\hat{\eta }}_k^n\\ {\hat{u}}_k^n\end{array}\right]=G_k\left[\begin{array}{c}{\hat{\eta }}_k^n\\ {\hat{u}}_k^n\end{array}\right],$ (2.8)

where $\sigma =\Delta t/\Delta x$ is the Courant number and

$$\begin{gathered} \text{c}(\theta ,\sigma ,\Delta x)=1-\frac{1}{2}{\sigma }^2{v}^2\left(\frac{\theta }{\Delta x}\right){\gamma }^2\left(\theta \right)+\frac{1}{4!}{\sigma }^4{v}^4\left(\frac{\theta }{\Delta x}\right){\gamma }^4\left(\theta \right), \\ \text{s}(\theta ,\sigma ,\Delta x)=\sigma v\left(\frac{\theta }{\Delta x}\right)\gamma \left(\theta \right)-\frac{1}{3!}{\sigma }^3{v}^3\left(\frac{\theta }{\Delta x}\right){\gamma }^3\left(\theta \right). \end{gathered}$$

It is easy to verify that we can write $G_k=V(k\pi /l){\tilde{G}}_{k}V(k\pi /l{)}^{-1}$, where

$V\left(\kappa \right)=\left[\begin{array}{cc}1& 0\\ 0& v\left(\kappa \right)\end{array}\right],{\tilde{G}}_k=\left[\begin{array}{cc}\text{c}({\theta }_k,\sigma ,\Delta x)& \mathit{i}\text{s}({\theta }_k,\sigma ,\Delta x)\\ \mathit{i}\text{s}({\theta }_k,\sigma ,\Delta x)& \text{c}({\theta }_k,\sigma ,\Delta x)\end{array}\right],$

and that the eigenvalues of G _k and ${\tilde{G}}_k$ are given by $g_k^{\pm }=g^{\pm }({\theta }_k,\sigma ,\Delta x)$, where

$g^{\pm }(\theta ,\sigma ,\Delta x)=\text{c}(\theta ,\sigma ,\Delta x)\mp \mathit{i}\text{s}(\theta ,\sigma ,\Delta x).$

Before enunciating the theorem which gives a sufficient condition for stability we need the following lemma which is proven in ⁶.

Lemma 2.1. There exist positive constants c ₁ and c ₂ such that$\forall y\in \mathrm{ℝ}$

$c_1\le v\left(y{)}^{-2}\right(1+y^2{)}^{-1}\le c_2.$ (2.9)

Theorem 2.1. Let$\beta >0$and$s\ge 0$. The numerical scheme (2.5) is stable, that is, inequality (2.6) holds, if$g^{\pm }(\theta ,\sigma ,\Delta x)$satisfy

$\left|g^{\pm }\right(\theta ,\sigma ,\Delta x\left)\right|\le 1,\forall \theta \in (-\pi ,\pi ].$ (2.10)

Proof. Let $[{\mathit{\eta }}^n,{\mathit{u}}^n{]}^T$ be the solution given by the method (2.5). For $s\ge 0$ we have

$\Vert [{\mathit{\eta }}^n,{\mathit{u}}^n{]}^T{\Vert }_{N,s,s+1}^2=\frac{1}{2l}\sum _{k=-N/2+1}^{N/2}{\left[1+(k\pi /l{)}^2\right]}^s{||B(k\pi /l)[{\hat{\eta }}_k^n,{\hat{u}}_k^n{]}^T||}_2^2,$

where ${||·||}_2$ denotes the usual Euclidean norm and $B\left(\kappa \right)=\mathrm{diag}(1,\sqrt{1+{\kappa }^2})$. Thus, it is enough to prove that there exists a positive constant C _s such that

${||B(k\pi /l)[{\hat{\eta }}_k^n,{\hat{u}}_k^n{]}^T||}_2^2\le C_s^2{||B(k\pi /l)[{\hat{\eta }}_k^0,{\hat{u}}_k^0{]}^T||}_2^2.$

Using equality (2.8) recursively we get $[{\hat{\eta }}_k^n,{\hat{u}}_k^n{]}^T=G_k^n[{\hat{\eta }}_k^0,{\hat{u}}_k^0{]}^T$, then

$$\begin{gathered} {||B(k\pi /l)[{\hat{\eta }}_k^n,{\hat{u}}_k^n{]}^T||}_2^2={||B(k\pi /l)\left(G_k{)}^n{B}^{-1}\right(k\pi /l\left)B\right(k\pi /l\left)\right[{\hat{\eta }}_k^0,{\hat{u}}_k^0{]}^T||}_2^2 \\ \le {||B(k\pi /l)\left(G_k{)}^n{B}^{-1}\right(k\pi /l)||}_2^2{||B(k\pi /l)[{\hat{\eta }}_k^0,{\hat{u}}_k^0{]}^T||}_2^2. \end{gathered}$$

Thereby, using that $G_k=V(k\pi /l){\tilde{G}}_{k}V(k\pi /l{)}^{-1}$ we get

${||B(k\pi /l)\left(G_k{)}^n{B}^{-1}\right(k\pi /l)||}_2^2\le {||B(k\pi /l)V(k\pi /l)||}_2^2\Vert {\tilde{G}}_k{\Vert }_2^{2n}{||{\left(B(k\pi /l)V(k\pi /l)\right)}^{-1}||}_2^2.$

Since matrix ${\tilde{G}}_k$ has eigenvalues $g^{\pm }({\theta }_k,\sigma ,\Delta x)$, by the hypothesis (2.10) we obtain $\Vert {\tilde{G}}_k{\Vert }_2^{2n}\le 1$.

On the other hand, using inequality (2.9) we get

$$\begin{gathered} \Vert B(k\pi /l)V(k\pi /l){\Vert }_2^2\Vert {\left(B(k\pi /l)V(k\pi /l)\right)}^{-1}{\Vert }_2^2= \\ \max \left\{v(k\pi /l{)}^2(1+(k\pi /l{)}^2),\frac{1}{v(k\pi /l{)}^2(1+(k\pi /l{)}^2)}\right\}\le \max \left\{1/c_1,c_2\right\}. \end{gathered}$$

Thus, we obtain ${||B\left(k\pi /l\right)[{\hat{\eta }}_k^n,{\hat{u}}_k^n{]}^T||}_2^2\le C_s^2{||B\left(k\pi /l\right)[{\hat{\eta }}_k^0,{\hat{u}}_k^0{]}^T||}_2^2$. Therefore we conclude that ${||[{\mathit{\eta }}^n,{\mathit{u}}^n{]}^T||}_{N,s,s+1}\le C_s{||[{\mathit{\eta }}^0,{\mathit{u}}^0{]}^T||}_{N,s,s+1}$. □

Condition (2.10) guarantees stability but it is not practical for numerical implementations. As shown in ¹, we can write the squared amplification factor as $|g_k^{\pm }{|}^2=1+p\left(\sigma v\left(\frac{{\theta }_k}{\Delta x}\right)\gamma \left({\theta }_k\right)\right)$, where $p\left(y\right)=y^6(y^2-8)/576$. Note that $p\left(y\right)\le 0$ if $\left|y\right|\le 2\sqrt{2}$, then the scheme is stable if

$\underset{\left|\theta \right|<\pi }{\max }\sigma \left|v\left(\frac{\theta }{\Delta x}\right)\gamma \left(\theta \right)\right|\le 2\sqrt{2}.$ (2.11)

In view of the previous computations, theorem 2.2 provides three practical conditions to guarantee stability.

Theorem 2.2. The numerical scheme (2.5) is stable if at least one of the following inequalities holds:

$\sigma =\frac{\Delta t}{\Delta x}\le {\gamma }_1\sqrt{1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }},{\gamma }_1=2\sqrt{2}{\left(\underset{\left|\theta \right|\le \pi }{\sup }\left|\gamma \right(\theta \left)\right|\right)}^{-1},$ (2.12)

$\mu =\frac{\Delta t}{\sqrt{\Delta x}}\le {\gamma }_2\sqrt{\sqrt{\beta }\left(1+\frac{{\rho }_2}{{\rho }_1}\right)},{\gamma }_2=2\sqrt{2}{\left(\underset{\left|\theta \right|\le \pi }{\sup }\left\{\frac{\left|\gamma \right(\theta \left)\right|}{\sqrt{\left|\theta \right|}}\right\}\right)}^{-1},$ (2.13)

$\Delta t\le {\gamma }_3\sqrt{\frac{\beta }{3}},{\gamma }_3=2\sqrt{2}{\left(\underset{\left|\theta \right|\le \pi }{\sup }\left\{\frac{\left|\gamma \right(\theta \left)\right|}{\left|\theta \right|}\right\}\right)}^{-1}.$ (2.14)

Proof. As previously stated, we just need to prove that each condition leads to inequality (2.11).

For the first condition we use the inequality $\varphi \left(y\right)\ge 1,\forall y\in \mathrm{ℝ}$ to obtain

$v\left(\kappa \right)={\left(1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }\varphi \left(\kappa \delta \right)+{\kappa }^2\frac{\beta }{3}\right)}^{-1/2}\le {\left(1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }\right)}^{-1/2}.$

Therefore,

$\left|v\left(\frac{\theta }{\Delta x}\right)\gamma \left(\theta \right)\right|\le \underset{\left|\theta \right|\le \pi }{\sup }v\left(\frac{\theta }{\Delta x}\right)\underset{\left|\theta \right|\le \pi }{\sup }\left|\gamma \left(\theta \right)\right|\le {\left[1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }\right]}^{-1/2}\underset{\left|\theta \right|\le \pi }{\sup }\left|\gamma \right(\theta \left)\right|.$ (2.15)

Multiplying inequality (2.15) by σ, applying condition (2.12) and making some simplifications we guarantee inequality (2.11).

For the second condition we use that $\varphi \left(y\right)\ge \left|y\right|,\forall y\in \mathrm{ℝ}$ and $1+\frac{\beta }{3}{y}^2>\sqrt{\beta }\left|y\right|,\forall y\in \mathrm{ℝ}$ for $\beta >0$ to get

$v\left(\kappa \right)={\left(1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }\varphi \left(\kappa \delta \right)+{\kappa }^2\frac{\beta }{3}\right)}^{-1/2}\le {\left(\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }\left|\kappa \delta \right|+\sqrt{\beta }\left|\kappa \right|\right)}^{-1/2},$

thus

$\left|v\left(\frac{\theta }{\Delta x}\right)\gamma \left(\theta \right)\right|\le \frac{\left|\gamma \right(\theta \left)\right|}{\sqrt{\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }|\theta \delta /\Delta x|+\sqrt{\beta }|\theta /\Delta x|}}\le \sqrt{\frac{\Delta x}{\sqrt{\beta }\left(1+\frac{{\rho }_2}{{\rho }_1}\right)}}\underset{\left|\theta \right|\le \pi }{\sup }\left\{\frac{\left|\gamma \right(\theta \left)\right|}{\sqrt{\left|\theta \right|}}\right\}.$

Multiplying inequality (2.13) $\left|v\left(\frac{\theta }{\Delta x}\right)\gamma \left(\theta \right)\right|/\sqrt{\Delta x}$ and applying the inequality above we obtain condition (2.11).

For the last condition we use that $\varphi \left(y\right)>0,\forall y\in \mathrm{ℝ}$, then

$\left|v\left(\frac{\theta }{\Delta x}\right)\gamma \left(\theta \right)\right|\le \underset{\left|\theta \right|\le \pi }{\sup }\left\{\frac{\left|\gamma \right(\theta \left)\right|}{\sqrt{\frac{\beta }{3}(\theta /\Delta x{)}^2}}\right\}=\left(\frac{\Delta x}{\sqrt{\frac{\beta }{3}}}\right)\underset{\left|\theta \right|\le \pi }{\sup }\left\{\frac{\left|\gamma \right(\theta \left)\right|}{\left|\theta \right|}\right\}.$

Using this inequality and condition (2.14) similarly to what was done previously we obtain inequality (2.11).

Therefore, if one of the conditions (2.12), (2.13) or (2.14) holds, the method is stable. □

The values of γ ₁, γ ₂ and γ ₃ depend on the spatial discretization used to obtain matrix C. Table 1 shows approximations for these values for each of the discretizations presented earlier in this section. The finite difference scheme provides larger values of γ ₁, γ ₂ and γ ₃ when compared to piecewise linear B-splines and the spectral scheme, thus its stability conditions are less restrictive than those provided by the other schemes.

Table 1:
Values of γ _i for each method.

	Finite Difference	B-Splines	Spectral Differentiation
γ 1	2.061	1.633	0.900
γ 2	2.651	2.300	1.595
γ 3	2.828	2.828	2.828

Now we can compare the scheme (2.5) for our dispersive system and the scheme presented in ¹ for system (1.6), which has the following stability conditions:

$$\begin{gathered} \sigma =\frac{\Delta t}{\Delta x}\le {\gamma }_1\sqrt{1+\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{\delta }},{\gamma }_1=2\sqrt{2}{\left({\sup }_{\left|\theta \right|\le \pi }\left|\gamma \right(\theta \left)\right|\right)}^{-1}, \\ \mu =\frac{\Delta t}{\sqrt{\Delta x}}\le {\gamma }_2\sqrt{\sqrt{\beta }\left(\frac{{\rho }_2}{{\rho }_1}\right)},{\gamma }_2=2\sqrt{2}{\left({\sup }_{\left|\theta \right|\le \pi }\left\{\frac{\left|\gamma \right(\theta \left)\right|}{\sqrt{\left|\theta \right|}}\right\}\right)}^{-1}. \end{gathered}$$

Due to the term $(\beta /3)u_{xxt}$ that leads to the term $(\beta /3)k^2$ in $v\left(k\right)$, the condition for $\Delta t/\sqrt{\Delta x}$ is less restrictive for system (1.4) than for system (1.6). Moreover, we can obtain a stability condition independent of ∆x for the scheme (2.5) which is not possible for the case of system (1.6) where there is no quadratic term like $(\beta /3)k^2$ in $v\left(k\right)$ and any inequality must rely on the linear growth of ϕ. The condition for the Courant number $\Delta t/\Delta x$ is the same in both cases.

3 NUMERICAL TESTS

In this section we will perform a sequence of numerical experiments in order to validate the methods and the stability conditions. The implementations were done in Octave. For the spatial derivative we choose the finite difference approximation because it has a low computational cost compared to the other schemes (piecewise linear B-splines and spectral scheme) and provides the largest values of γ ₁, γ ₂ and γ ₃, as shown in table 1.

In view of the conservation law (1.2) we choose as initial configuration the profile ${\eta }^0\left(x\right)=\overline{\eta }\left(x\right)-a$, where $\overline{\eta }\left(x\right)=\exp (-2x^2)$ and the constant a is calculated so that η ⁰ satisfies $0=\Delta x\sum _{j=1}^N{\eta }_j^0$ which is an approximation of (1.2) for $t=0$ by the trapezoidal rule.

We set u ⁰ so that η only propagates to the right in the linear system, that is, ${\hat{u}}^0\left(k\right)=v(k\pi /l){\hat{\eta }}^0\left(k\right)$. For the nonlinear system u ⁰ is set in the same way. For the numerical experiments we set ${\rho }_1=1,{\rho }_2=2,h_1=0.1$ and $h_2=3.505$ as in ⁸ and $l=10\pi$. The values of β, α, N, ∆x and ∆t are defined in each test. The other values are calculated by the relations $L=h_1/\sqrt{\beta }$ and $\delta =h_2/L$.

Figure 2 presents the graphics of the numerical solution of system (1.4) at different instants for $\beta =0.0001$ with $\Delta x=0.03068(N=2^{11})$ and $\Delta t=0.08043$ which satisfy the stability condition (2.13). We can see the dispersion acting on η and u because of the wave trains that form as time advances. In fact, as the phase velocity $v\left(k\right)$ decreases to zero when $\left|k\right|\to +\infty$ the high frequency components of the solution propagate more slowly and form the aforementioned wave trains.

Figure 2:
Solutions for the linear system (1.4) at times t = 0 (-), t = 1100∆t (− · −) and $t=2200∆t$ (···) with β = 0.0001, ∆x = 0.03068 and ∆t = 0.08043.

Although theorem 2.2 provides sufficient conditions for stability they are not necessary. In fact, if we choose $∆t=0.08847$ and keep the other parameters unaltered none of the conditions required for stability in theorem 2.2 is satisfied, but condition (2.10) still holds. However, for $∆t=0.09460$, condition (2.10) is not satisfied because $\left|g^{\pm }\right(\pi /2,\sigma ,\Delta x\left)\right|=1.0029$ and since no stability condition holds, high frequency components permeate the numerical solution as shown in figure 3 for $t=1800∆t$. Note that this phenomenon is different from dispersion in which the components propagate with different speeds. Here, besides the dispersion effect from system (1.4), the amplitude of high frequency components increases anomalously and the numerical solution is compromised. This example illustrates the relevance of a proper stability analysis even when the stability conditions are less restrictive.

Figure 3:
Solutions for the linear system (1.4) at times t = 0 (-), $t=900∆t$ (− · −) and $t=1800∆t$ (···) with β = 0.0001, ∆x = 0.03068 and ∆t = 0.09460.

Figure 4 presents the graphics of the numerical solution of system (1.4) for $\beta =0.001$ using $\Delta x=0.00383(N=2^{14})$ and $∆t=0.05163$ which satisfy the stability condition (2.14). We see that increasing the parameter β the phase velocity $v\left(k\right)$ decreases to zero faster, thus a wider wave train is formed. If $∆t=0.06196$, none of the stability conditions in theorem 2.2 is satisfied, but condition (2.10) still holds. For $∆t=0.07026$, condition (2.10) is not satisfied since we have $\left|g^{\pm }\right(1.12,\sigma ,\Delta x\left)\right|=1.0062$ and the solution obtained (keeping the rest of the parameters unaltered) is unstable as shown in figure 5.

Figure 4:
Solutions for the linear system (1.4) at times t = 0 (-), $t=1800∆t$ (− · −) and $t=3600∆t$ (···) with β = 0.001, ∆x = 0.00383 and ∆t = 0.05163.

Figure 5:
Solutions for the linear system (1.4) at times t = 0 (-), $t=1350∆t$ (− · −) and $t=2700∆t$ (···) with β = 0.001, ∆x = 0.00383 and ∆t = 0.07026.

The stability performance exhibited illustrates the fact that theorem 2.2 provides sufficient, not necessary restrictions for stability. Rather than the best possible condition, theorem 2.2 provides sufficient stability conditions which are useful for implementation codes.

Condition (2.10) only makes sense in the linear case where there exists a well defined amplification factor, but the stability conditions in theorem 2.2 remain valid for the scheme for the nonlinear system (1.3) as the following numerical experiments indicate. Figure 6 shows the graphics of the numerical solutions for system (1.3) for $\beta =\alpha =0.0001$, where $\Delta x=0.03068(N=2^{11})$ and $∆t=0.08043$ satisfy condition (2.13), also figure 7 presents the graphics of the numerical solutions for $\beta =\alpha =0.001$, with $\Delta x=0.00383(N=2^{14})$ and $∆t=0.05163$ satisfying the stability condition (2.14). We see that the conditions in theorem 2.2 still hold for these nonlinear cases. Figures 8 and 9 show that keeping all parameters used in figures 6 and 7, respectively, with the exception of ∆t, which is increased so that no stability condition in theorem 2.2 is satisfied, the numerical solution is unstable as in the linear case. The value used for ∆t is specified in the caption of each figure.

Figure 6:
Solutions for the nonlinear system (1.3) at times t = 0 (-), $t=1100∆t$ (− · −) and

Figure 7:
Solutions for the nonlinear system (1.3) at times t = 0 (-), $t=1800∆t$ (− · −) and

Figure 8:
Solutions for the nonlinear system (1.3) at times t = 0 (-), $t=900∆t$ (− · −) and

Figure 9:
Solutions for the nonlinear system (1.3) at times t = 0 (-), $t=1350∆t$ (− · −) and

Solutions for the nonlinear system (1.3) have a similar appearance to solutions for the linear system (1.4) with the same set of common parameters and starting from the same initial conditions. In order to compare them we compute the Euclidean norm of the vector that results from their difference at several instants. Table 2 shows that the Euclidean norm of the difference grows as time advances and it is larger for the larger value of the nonlinear parameter. In particular, the values at the intersection of the rows identified by $\beta =0.0001$ and the columns identified by $\overline{\eta }\left(x\right)=\exp (-2x^2)$ correspond to the two experiments that generate figures 2 and 6. Analogously, the values at the intersection of the rows identified by $\beta =0.001$ and the columns identified by $\overline{\eta }\left(x\right)=\exp (-2x^2)$ correspond to the two experiments that generate figures 4 and 7. Besides, table 2 shows that if $\overline{\eta }\left(x\right)$ changes from $0.1\exp (-2x^2)$ to $\exp (-2x^2)$ and therefore the initial values for η and u change by the same factor (keeping all the parameters equal), the Euclidean norm of the corresponding difference is multiplied approximately by one hundred, confirming that the nonlinear effects are present and are more noticeable for higher amplitude profiles.

Table 2:
Euclidean norm of the difference between the solutions of the linear system (1.4) and the nonlinear system (1.3).

		$\overline{\eta }\left(x\right)=0.1\exp (-2x^2)$		$\overline{\eta }\left(x\right)=\exp (-2x^2)$
β	time	norm (η)	norm (u)	norm (η)	norm (u)
0.0001	100∆t = 8.043	0.0000495	0.0000475	0.0049454	0.0047532
	500∆t = 40.213	0.0002119	0.0002041	0.0211847	0.0204099
	1100∆t = 88.468	0.0003597	0.0003474	0.0359573	0.0347279
	2200∆t = 176.937	0.0005106	0.0004941	0.0510226	0.0493752
0.001	150∆t = 7.745	0.0011563	0.0010590	0.1156249	0.1059029
	700∆t = 36.142	0.0025461	0.0023993	0.2542572	0.2395679
	1800∆t = 92.938	0.0031029	0.0029502	0.3092888	0.2940249
	3600∆t = 185.875	0.0033832	0.0032392	0.3365037	0.3221350

3.1 Convergence analysis

Now let us study the convergence of the method on the temporal and spatial variables. Since the classical Runge-Kutta method has fourth order of convergence, we just need to focus in the spatial variable.

Since the exact solution is not available for the nonlinear system, in order to compare the numerical approximations, we substitute it by an accurate numerical solution $(\overline{\mathit{\eta }},\overline{\mathit{u}})$ that is calculated with small values of $\Delta x=\Delta x^{*}$. For a fair comparison ∆t is kept constant and chosen to guarantee stability in all tests. For a given time $T=n\Delta t$ we define the time error $E_x(\eta ,\Delta x,\Delta x^{*},T)=\Vert {\mathit{\eta }}^n-{\overline{\mathit{\eta }}}^n{\Vert }_2,T=n\Delta t$ but considering only the points where all numerical solutions are calculated. We perform successive refinements of the spatial mesh dividing ∆x by 2 and estimating the spatial convergence rate by

$p\approx -{\log }_2\left(E_x\right(\eta ,\Delta x/2,\Delta x^{*},T)/E_x(\eta ,\Delta x,\Delta x^{*},T\left)\right).$

The error and rate for u are defined analogously.

A fourth order convergence rate is expected in the spatial variable since the finite difference formula has order 4 and the Fourier approximation has spectral convergence. The estimated convergence rate is confirmed in the experiments for the linear and nonlinear cases for different values of β and α as shown in tables 3-6.

Table 3:
Convergence analysis on spatial variable x for the linear system with β = 0.001, l = 10π,

∆x	error (η)	rate (η)	error (u)	rate (u)
0.122718	0.260102684392		0.231045716400
0.061359	0.019181498462	3.76129	0.016887753929	3.77413
0.030680	0.001209018675	3.98781	0.001064172067	3.98817
0.015340	0.000075683518	3.99771	0.000066613691	3.99777
0.007670	0.000004715024	4.00464	0.000004149941	4.00466
0.003835	0.000000277385	4.08730	0.000000244141	4.08730

Table 4:
Convergence analysis on spatial variable x for the linear system with β = 0.0001, l = 10π, ∆x ^* = 0.001917, ∆t = 0.020106 and T = 99.989358.

∆x	error (η)	rate (η)	error (u)	rate (u)
0.122718	0.275386408233		0.263332125055
0.061359	0.020799252938	3.72685	0.019828641788	3.73123
0.030680	0.001311657678	3.98707	0.001250336121	3.98720
0.015340	0.000082111042	3.99767	0.000078271354	3.99769
0.007670	0.000005115489	4.00463	0.000004876264	4.00464
0.003835	0.000000300945	4.08730	0.000000286871	4.08730

Table 5:
Convergence analysis on spatial variable x for the nonlinear system with $\beta =\alpha =0.001$, l = 10π, ∆x ^* = 0.001917, ∆t = 0.035755 and T = 99.970644.

∆x	error (η)	rate (η)	error (u)	rate (u)
0.122718	0.268262793175		0.238173684977
0.061359	0.019996662700	3.74582	0.017583749591	3.75970
0.030680	0.001261043561	3.98707	0.001108541009	3.98751
0.015340	0.000078944237	3.99764	0.000069394349	3.99770
0.007670	0.000004918219	4.00463	0.000004323218	4.00464
0.003835	0.000000289340	4.08730	0.000000254335	4.08730

Table 6:
Convergence analysis on spatial variable x for the nonlinear system with $\alpha =\beta =0.0001$, l = 10π, ∆x ^* = 0.001917, ∆t = 0.020106 and T = 99.989358.

∆x	error (η)	rate (η)	error (u)	rate (u)
0.122718	0.278408434936		0.266207168503
0.061359	0.021121230297	3.72044	0.020132847020	3.72493
0.030680	0.001332336906	3.98666	0.001269866417	3.98680
0.015340	0.000083407580	3.99764	0.000079495826	3.99765
0.007670	0.000005196285	4.00463	0.000004952569	4.00463
0.003835	0.000000305699	4.08730	0.000000291361	4.08730

3.2 Approximate travelling waves

Travelling wave solutions constitute a subject of major relevance in the study of nonlinear wave models. Important unidirectional models like Korteweg-de Vries (KdV), ILW and rILW equations and bidirectional models like Boussinesq systems exhibit travelling waves, see for example ^{2), (7), (11), (14), (18} and the references therein. These solutions do not change their shapes and propagate at constant speed by maintaining a balance between the nonlinear and the dispersive effects of the model. It is also important to know if the travelling wave solutions are stable to small perturbations, otherwise any physical or numerical disturbance will eventually destroy them, ⁹. In this section we will study the existence of travelling wave solutions for the nonlinear system (1.3) from a numerical point of view since theoretical results about their existence are not available.

As initial condition for the nonlinear system (1.3), we will use solitary wave profiles from the ILW and the rILW equations in order to measure if their shapes are well preserved while evolving according to the numerical scheme. Let us consider the rILW and the ILW equations given by equations (3.1) and (3.2), respectively:

${\eta }_t+{\eta }_x-\frac{3\alpha }{2}\eta {\eta }_x-\frac{\sqrt{\beta }}{2}\frac{{\rho }_2}{{\rho }_1}{\mathcal{T}}_{\delta }\left[{\eta }_{xt}\right]=0,$ (3.1)

${\eta }_t+{\eta }_x-\frac{3\alpha }{2}\eta {\eta }_x+\frac{\sqrt{\beta }}{2}\frac{{\rho }_2}{{\rho }_1}{\mathcal{T}}_{\delta }\left[{\eta }_{xx}\right]=0.$ (3.2)

Setting $c_1=-\frac{3}{2}\alpha ,c_2=\frac{{\rho }_2}{{\rho }_1}\frac{\sqrt{\beta }}{2}$, we can see in ¹⁵ that a travelling wave solution family for the rILW equation is

$\eta \left(y\right)=\frac{a{\cos }^2\left(\theta \right)}{{\cos }^2\left(\theta \right)+{\sinh }^2(y/\lambda )},y=x-ct,$ (3.3)

where

$a=\frac{4c{c}_2\theta \tan \theta }{\delta c_1},c=\frac{1}{1+\frac{2c_2}{\delta }\theta \cot \left(2\theta \right)},\lambda =\frac{\delta }{\theta },0<\theta <\pi /2.$

In ^{8), (11} we see that a travelling wave solution family for the ILW equation is also given by the expression (3.3) but with

$a=\frac{4c_2\theta \tan \theta }{\delta c_1},c=1-\frac{2c_2}{\delta }\theta \cot \left(2\theta \right),0<\theta <\pi /2.$

The numerical wave speed c _n is computed as follows: for each time $t_j=j\Delta t,j=1,2,\dots ,J$ big enough, we find the grid point $x_{l_j}$ that minimizes ${\eta }_l^j$ and make a linear regression on the points $\left\{\right(t_j,x_{l_j}\left)\right\}$ to estimate c _n . Using Fourier properties and the DFT, we define η * as the approximation of $\eta (x-c_{n}t,0)$ by ${\hat{\eta }}_k^{*}=\exp (-\mathit{i}k\pi c_{n}t/l){\hat{\eta }}_k^0$. We define the absolute error e _abs and the relative error e _rel at a chosen $t_n=n\Delta t$, respectively, by $e_{\text{abs}}=\Vert {\mathit{\eta }}^{*}-{\mathit{\eta }}^n{\Vert }_{\infty }$ and $e_{\text{rel}}=\Vert {\mathit{\eta }}^{*}-{\mathit{\eta }}^n{\Vert }_{\infty }/\Vert {\mathit{\eta }}^{*}{\Vert }_{\infty }$.

We will set the parameters, ${\rho }_1=1,{\rho }_2=2,h_1=0.1,h_2=3.505,l=10\pi$ and $\Delta x=0.03068(N=2^{11})$. The values of β, α, θ and ∆t are specified in each test. The other values are calculated by the relations $L=h_1/\sqrt{\beta }$ and $\delta =h_2/L$. We adjust the initial conditions using the approximation of the conservation law (1.2) as in section 3.

Figure 10 presents the graphics of η ⁿ , the evolution of the travelling wave initial condition η ⁰ and the corresponding approximate translation η * for $\alpha =\beta =0.0001,\theta =\pi /20$ and $n=1600$. Figure 11 presents the same type of results setting $\theta =\pi /40$. We can see that the ILW and rILW profiles preserve their shapes in both cases since their graphics coincide with their corresponding η *. This is confirmed by the small errors presented in table 7.

Figure 10:
Graphics of η * (− · −) and η ⁿ (-) for each initial condition at time $t=1600∆t$ with $\alpha =\beta =0.0001$, $\theta =\pi /20$, ∆x = 0.03068 and ∆t = 0.08043.

Figure 11:
Graphics of η * (− · −) and ηn (-) for each initial condition at time $t=1600∆t$ with $\alpha =\beta =0.0001$, $\theta =\pi /20$, ∆x = 0.03068 and ∆t = 0.08043.

Table 7:
Errors at time $t=1600∆t$ and wave velocities of η for each initial condition with $\alpha =\beta =0.0001$.

Profiles	θ	e abs	e rel	c n	c
rILW	$\pi /20$	0.0124475	0.0007270	0.97331	0.97315
ILW	$\pi /20$	0.0232071	0.0013190	0.97338	0.97241
rILW	$\pi /40$	0.0029185	0.0007436	0.97273	0.97248
ILW	$\pi /40$	0.0040683	0.0010080	0.97273	0.97170

Setting $\alpha =\beta =0.001$ the shape of the ILW and rILW profiles given by equation (3.3) are also preserved as the next experiments show. Figure 12 presents the graphics of the evolution η ⁿ for the travelling wave initial condition η ⁰ of each equation and the corresponding approximate translation η * for $\theta =\pi /20$ and $n=900$. Figure 13 presents the same type of results for $\theta =\pi /40$. Again, the graphics of η ⁿ coincide with their corresponding η * and the preservation of shape is confirmed by the small errors presented in table 8.

Figure 12:
Graphics of η * (− · −) and η ⁿ (-) for each initial condition at time $t=900∆t$ with

Figure 13:
Graphics of η * (− · −) and ηn (-) for each initial condition at time $t=900∆t$ with

Table 8:
Errors at time $t=900∆t$ and wave velocities of η for each initial condition with $\alpha =\beta =0.001$.

Profiles	θ	e abs	e rel	c n	c
rILW	$\pi /20$	0.0008796	0.0006144	0.97289	0.97315
ILW	$\pi /20$	0.0006704	0.0004557	0.97294	0.97241
rILW	$\pi /40$	0.0023305	0.0090672	0.97254	0.97248
ILW	$\pi /40$	0.0024069	0.0091069	0.97257	0.97170

Thus, the solitary wave family profiles of the ILW and rILW equations, for the corresponding values of α, β, ρ ₁, ρ ₂ and δ, perform satisfactorily as approximate travelling wave solutions for the nonlinear system (1.3), since their shapes are well preserved by the discrete scheme for the nonlinear system as time advances.

4 CONCLUSIONS

The numerical schemes proposed for both linear and nonlinear systems use a fourth order finite difference method for the first order x-derivatives, spectral schemes for the dispersive terms and the classical fourth order Runge-Kutta method for time evolution. Based on the results obtained in this work we conclude that the numerical methods for the linear and nonlinear systems are very robust. They present good stability conditions and convergence rates in temporal and spatial variables.

The shapes of the solitary waves of the ILW and rILW equations were well preserved when chosen as the initial condition for the discrete scheme of the nonlinear system. This may indicate that the nonlinear system (1.3) admits travelling wave solutions. A theoretical study of the existence of travelling waves is the subject of current research.

Acknowlegdments

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior Brasil (CAPES) Finance Code 001.

REFERENCES

• ¹
  D. Alfaro Vigo, S. Oliveira, A. Ruiz de Zárate & A. Nachbin. Fully discrete stability and dispersion analysis for a linear dispersive internal wave model. Computational and Applied Mathematics, 33 (2013), 203-221.
• ²
  J. Angulo, E. Cardoso Jr & F. Natali. Stability properties of periodic traveling waves for the Intermediate Long Wave Equation. Revista Matemática Iberoamericana, 33 (2015), 1-35.
• ³
  J. Apel, L. Ostrovsky, Y. Stepanyants & J. Lynch. Internal Solitons in the Ocean. Technical report, Woods Hole Oceanographic Institution (2006).
• ⁴
  U.M. Ascher & L.R. Petzold. “Computer methods for ordinary differential equations and differentialalgebraic equations”. Society for Industrial and Applied Mathematics (1998), 332 p.
• ⁵
  T.B. Benjamin. Internal waves of permanent form in fluids of great depth. Journal of Fluid Mechanics, 29(3) (1967), 559-592.
• ⁶
  J.S. Brodzinski. “Estudo de um modelo dispersivo na˜o linear para ondas internas”. Ph.D. thesis, UFPR (2016).
• ⁷
  M. Chen, N.V. Nguyen & S.M. Sun. Existence of traveling-wave solutions to Boussinesq systems. Differential Integral Equations, 24(9/10) (2011), 895-908.
• ⁸
  W. Choi & R. Camassa. Fully nonlinear internal waves in a two-fluid system. Journal of Fluid Mechanics , 396 (1999), 1-36.
• ⁹
  P.G. Drazin & R.S. Johnson. “Solitons: An Introduction”. Cambridge University Press (1989), 240 p.
• ¹⁰
  K.R. Helfrich & W.K. Melville. Long Nonlinear Internal Waves. Annual Review of Fluid Mechanics, 38 (2006), 395-425.
• ¹¹
  R.I. Joseph. Solitary waves in a finite depth fluid. Journal of Physics A Mathematical General, 10(12) (1977), L225-L227.
• ¹²
  W.C. Lesinhovski. “Análise de von Neumann de um modelo dispersivo de ondas internas”. Master’s thesis, UFPR (2017).
• ¹³
  F. Nansen. “Farthest North: The Epic Adventure of a Visionary Explorer”. Skyhorse Publishing Inc. (2008), 678 p.
• ¹⁴
  F. Oliveira. A note on the existence of traveling-wave solutions to a Boussinesq system. Differential Integral Equations , 29(1/2) (2016), 127-136.
• ¹⁵
  A. Ruiz de Zárate. “A Reduced Model for Internal Waves Interacting with Submarine Structures at Intermediate Depth”. Ph.D. thesis, IMPA (2007).
• ¹⁶
  A. Ruiz de Zárate , D. Alfaro Vigo , A. Nachbin & A. Choi. A Higher-order Internal Wave Model Accounting for Large Bathymetric Variations. Studies in Applied Mathematics, 122 (2009), 275-294.
• ¹⁷
  B.R. Sutherland. “Internal Gravity Waves”. Cambridge University Press (2010). doi:10.1017/CBO9780511780318.
  » https://doi.org/10.1017/CBO9780511780318
• ¹⁸
  G.B. Whitham. “Linear and Nonlinear Waves”. Wiley (2011), 660 p.

Publication Dates

• Publication in this collection
  04 Apr 2022
• Date of issue
  Jan-Mar 2022

History

• Received
  09 Nov 2020
• Accepted
  28 June 2021