An efficient method to compute solitary wave solutions of fractional Korteweg-de Vries equations

TL;DR

This paper introduces an efficient numerical method combining fixed-point iteration and extrapolation to accurately compute solitary wave solutions of fractional Korteweg-de Vries equations, enabling detailed analysis of wave properties.

Contribution

The paper presents a novel, accelerated fixed-point iterative algorithm for computing solitary wave profiles of fractional KdV equations with improved convergence and accuracy.

Findings

Enhanced convergence speed and accuracy in computing wave profiles.

Validated numerical method through experiments and analysis.

Enabled detailed investigation of wave properties like speed-amplitude relation.

Abstract

Considered here is an efficient technique to compute approximate profiles of solitary wave solutions of fractional Korteweg-de Vries equations. The numerical method is based on a fixed-point iterative algorithm along with extrapolation techniques of acceleration. This combination improves the performance in both the velocity of convergence and the computation of profiles for limiting values of the fractional parameter. The algorithm is described and numerical experiments of validation are presented. The accuracy attained by the procedure can be used to investigate additional properties of the waves. This approach is illustrated here by analyzing the speed-amplitude relation.

Equations

$$u_t+u^p{u}_x-(D^{\alpha }u{)}_x=0.$$

$$\widehat{\left(D^{\alpha }g\right)}(\xi )=\beta (\xi )\hat{g}(\xi ),\beta (\xi )=|\xi {|}^{\alpha },$$

$$\hat{g}(\xi )={\int }_{-\infty }^{\infty }g(x)e^{-i\xi x}dx,\xi \in \mathbb{R}$$

$$L^{\infty }(\mathbb{R},H^{\alpha /2}(\mathbb{R}))=\{v:\mathbb{R}\to H^{\alpha /2}(\mathbb{R})/\underset{t\in \mathbb{R}}{\max }||v(t)|{|}_{H^{\alpha /2}(\mathbb{R})}<\infty \},$$

$$C(u)={\int }_{-\infty }^{\infty }u(x,t)dx,$$

$$M(u)={\int }_{-\infty }^{\infty }{u}^2(x,t)dx,$$

$$M(u)={\int }_{-\infty }^{\infty }\left(\frac{1}{2}|D^{\alpha /2}u(x,t){|}^2-\frac{1}{(p+1)(p+2)}{u}^{p+2}(x,t)\right)dx.$$

$$D^{\alpha }{\varphi }_c+c{\varphi }_c-\frac{{\varphi }_c^{p+1}}{p+1}=0.$$

$$u_t+(f(u){)}_x-(\mathcal{M}u{)}_x=0,$$

$$\beta (\xi )={\left(1+\gamma |\xi {|}^2\right)}^{1/2}{\left(\frac{\tanh \xi }{\xi }\right)}^{1/2},f(u)=\frac{u^2}{2}.$$

$$u_t+u_x+(f(u){)}_x+(\mathcal{M}u{)}_t=0,$$

$$\mathcal{L}\varphi =\mathcal{N}(\varphi ),\mathcal{L}=D^{\alpha }+c,\mathcal{N}(\varphi )=\frac{{\varphi }^{p+1}}{p+1},c>0.$$

$$\mathcal{N}(\lambda u)={\lambda }^{p+1}\mathcal{N}(u),\lambda ,u\in \mathbb{R}.$$

$${\mathcal{N}}^{\prime }(u)u=(p+1)\mathcal{N}(u),u\in \mathbb{R}.$$

$${\mathcal{L}}^{-1}{\mathcal{N}}^{\prime }({\varphi }_c){\varphi }_c=(p+1){\mathcal{L}}^{-1}\mathcal{N}({\varphi }_c)=(p+1){\varphi }_c.$$

$$\mathcal{L}{\varphi }_{n+1}=\mathcal{N}({\varphi }_n),n=0,1,\dots$$

$$m({\varphi }_n)=\frac{(\mathcal{L}{\varphi }_n,{\varphi }_n)}{(\mathcal{N}({\varphi }_n),{\varphi }_n)},$$

$$\mathcal{L}{\varphi }_{n+1}=m({\varphi }_n{)}^{ϵ}\mathcal{N}({\varphi }_n),n=0,1,\dots ,$$

$${\mathcal{L}}_h={\mathcal{N}}_h({\varphi }_h),$$

$${\mathcal{L}}_h=\gamma c{I}_N+D_N^{\alpha },{\mathcal{N}}_h({\varphi }_h)=\frac{{\varphi }_h{.}^{p+1}}{p+1},$$

$$D_N^{\alpha }=F_N^{-1}{\Lambda }_N^{\alpha }{F}_N,$$

$$m({\varphi }_n)=\frac{{\sum }_{k=0}^{N-1}\left(c+|k\pi /l{|}^{\alpha }\right)|\widehat{{\varphi }_n}(k){|}^2}{{\sum }_{k=0}^{N-1}\widehat{\mathcal{N}({\varphi }_n)}(k)\overline{\widehat{{\varphi }_n}(k)}},$$

$$\widehat{{\varphi }_{n+1}}(k)=m({\varphi }_n{)}^{ϵ}\frac{\widehat{\mathcal{N}({\varphi }_n)}(k)}{c+|k\pi /l{|}^{\alpha }},$$

$$|1-m_n|,n=0,1,\dots$$

$$ERRO{R}_c(n)=||{\varphi }_n-{\varphi }_{n-1}||,n=0,1,\dots$$

$$RES(n)=||{\mathcal{L}}_h{\varphi }_n-{\mathcal{N}}_h({\varphi }_n)||,n=0,1,\dots ,$$

$$Amp={\left(\frac{c}{2}(p+1)(p+2)\right)}^{\frac{1}{p}},$$

Full text

††articletype: ARTICLE TEMPLATE

1 Introduction

The paper is concerned with the computation of solitary wave solutions of the fractional Korteweg-de Vries (fKdV) equation

[TABLE]

In (1), $u=u(x,t)$ is a real-valued function of $x\in\mathbb{R},t\geq 0$, $p\in\mathbb{N},\alpha\in\mathbb{R}$ and $D^{\alpha}$ stands for the linear operator represented by the symbol

[TABLE]

where

[TABLE]

is the Fourier transform, defined on the space of squared integrable functions $g\in L^{2}(\mathbb{R})$. Equation (1) is relevant as a dispersive and nonlinbear perturbation of Burgers’ inviscid equation. The parameters $\alpha$ and $p$ govern, respectively, the dispersion and the nonlinear effects and in this sense (1) is suitable to investigate the relations between nonlinearity and dispersion that lead to different dynamics, such as existence and stability of solitary waves, blow-up phenomena, etc, [23]. The case $\alpha=2$ corresponds to the classical Korteweg-de Vries (KdV, when $p=1$) and generalized Korteweg-de Vries (gKdV, when $p\geq 2$) equations, while $\alpha=1$ leads to the Benjamin-Ono (BO) equation and generalized (gBO) versions.

The parameters $\alpha$ and $p$ also determine several mathematical properties of the initial value problem for (1), (2). The main results on well-posedness in the literature concern the case $p=1$, for which the Cauchy problem, when $\alpha\geq 1$, is proved to have global solutions in suitable functional spaces, [19]. When $-1<\alpha<0$, suitable smooth initial data lead the corresponding solution to blow up at finite time, [23]. As far as the case $0<\alpha<1$ is concerned, see [26] for local well-posedness results. Global weak solutions, without uniqueness, in

[TABLE]

for initial data in the Sobolev space $H^{\alpha/2}(\mathbb{R})$ (with the usual norm $||\cdot||_{H^{\alpha/2}(\mathbb{R})}$, are proved to exists in [36] (see also [26]) for $\alpha>1/2$. Also, Klein and Saut, [23] conjecture several cases of blow-up with different structure: first, no hyperbolic blow-up (blow-up of the spatial gradient with bounded sup-norm) exists; when $1/2<\alpha<1$, the solution is global; when $1/3<\alpha<1/2$, there is a sort of nonlinear dispersive blow-up, [28, 29, 30], and, finally, blow-up of different type occurs when $0<\alpha<1/3$.

On the other hand, at least formally, the following quantities are preserved by smooth enough, decaying solutions of (1)

[TABLE]

[TABLE]

[TABLE]

The quantity (3) is well defined when $\alpha\geq 1/3$ and provides a Hamiltonian structure to (1), see [23] for the case $p=1$ and [5] for $p>1$.

An additional, relevant point on the dynamics of (1) concerns the existence and stability of solitary wave solutions. They are solutions of the form $u(x,t)=\phi(x-ct)$ for some speed $c>0$ and profile $\phi_{c}=\phi_{c}(X)$ with $\phi_{c}(X)\rightarrow 0$ as $|X|\rightarrow\infty$. Substituting into (1) and integrating once, the profile $\phi_{c}$ must satisfy

[TABLE]

Some known results on existence and stability of solutions of (4) are now summarized. For the case $p=1$, solitary waves are proved to exist when $\alpha>1/3$, [13, 14, 15], with an asymptotic decay as $1/x^{1+\alpha},|x|\rightarrow\infty$; orbital stability, [4, 17], for $\alpha>1/2$ is proved in [26], while some spectral instability analysis can be obtained from [22, 32, 5]. In the case $p>1$, as mentioned in [5] (see also references therein), existence of solutions of (4) holds for $\alpha>1$ and any $p$, with orbital stability when $p<2\alpha$. For $0<\alpha<1$, existence can be derived for $1<p<2\alpha/(1-\alpha)$, [42, 41, 6]. Results on orbital stability when $1/2<\alpha<1,p<2\alpha$, as well as a linear instability criterium, can be seen in [5].

Since no explicit formulas for solitary wave solutions of (1) are, except in the classical cases $\alpha=1,2$, unknown, then some numerical method for the generation of approximate profiles is required. In this sense, Klein and Saut, [23], solve numerically (4) for $p=1$ to construct approximate solitary wave profiles. The numerical method to this end is based on transforming (4) into the corresponding algebraic equations for the Fourier transform of the profile, which are iteratively solved by Newton iteration. The procedure is implemented by approximating each profile with a trigonometric interpolant polynomial on a long enough interval, in such a way that the Newton method is applied to the system for the corresponding Fourier components. As mentioned in [23], the algebraic decrease of the modulus of the Fourier coefficients, due to the loss of smoothness of the periodic approximations of the profiles at the boundary of the computational domain (the solitary wave profile decreases slowly at infinity) leads to a slow convergence of the iteration. This is addressed by taking a large computational domain and a high number of Fourier modes. The resolution is also performed by using GMRES, [35], to compute iteratively the inverse of the Jacobian matrix.

In this paper, an alternative to construct numerically solitary wave profiles of (4) is proposed. The technique was successfully applied to the numerical generation of periodic traveling wave solutions of the fKdV equation in [3]. The main points of this approach here are the following:

• •

The method is also based on the implementation of (4) in Fourier space for the periodic approximation on a long enough interval.

• •

The algebraic system for the discrete Fourier coefficients of the trigonometric interpolant is however iteratively solved by the Petviashvili method, [34], a fixed point type algorithm which may overcome some of the limitations of the Newton iteration.

• •

In order to improve the slow convergence due to the periodic approximation, the Petviashvili method is complemented with the use of acceleration techniques based on extrapolation, [38, 39], which have shown a relevant performance in the numerical generation of solitary waves, [2].

The main contributions are the following:

• •

The combination of the Petviashvili method with extrapolation improves the computation of the solitary wave profiles. The improvement is observed in mainly two points: the first one is that the method is able to generate numerical profiles for limiting values of $\alpha$. A second point of improvement is found in the efficiency, since with a high number of Fourier modes and on a long interval, the extrapolation technique accelerates the convergence in a relevant way.

• •

These new advantages can be used to study computationally additional properties of the waves, such as the speed-amplitude relation.

• •

The method can also be applied to compute approximate solitary wave solutions of other generalizations of (1) of the form, [5]

[TABLE]

where ${\cal M}$ is a linear, pseudo-differential operator associated to a continuous, even, real-valued Fourier symbol $\beta(\xi)$ and $f$ is a smooth, nonlinear, real-valued function. A relevant example of (5) is the extended Whitham equation, [24, 25, 11], for which

[TABLE]

In (6), the parameter $\gamma\geq 0$ controls the surface tension effects in the model. The case $\gamma=0$ leads to the classical Whitham equation, [43, 31]. The computation of traveling-wave solutions of equations of the form (5), kncluding the extended Whitham equation (6), has been made in the literature with different techniques, see e. g. the references in [3] and, more recently, the method introduced and performed in [21], based on continuation with spectral projection.

The technique can also be applied to study the solitary wave solutions of fractional BBM type equations

[TABLE]

with ${\cal M}$ and $f$ defined as in (5), see [5].

The structure of the paper is as follows. In Section 2 the numerical method to compute approximate solitary profiles, based on the Petviashvili method and accelerating techniques, is described, along with some implementations details. The purpose of Section 3 is two-fold: a first group of experiments validates the efficiency of the method and studies its performance. This is used, in a second part, to analyze computationally additional properties of the solitary waves. The illustration is focused on the speed-amplitude relation and its dependence on the parameters $\alpha$ and $p$.

2 Numerical generation of solitary waves

2.1 The Petviashvili method

In order to describe the numerical method to compute solitary waves of (1), observe that (4) can be written in the form

[TABLE]

Note now that the nonlinear term $\mathcal{N}$ is homogeneous of degree $p+1$, in the sense that

[TABLE]

Thus, differentiating (8) with respect to $\lambda$ and evaluating at $\lambda=1$, we have

[TABLE]

On the other hand, the operator $\mathcal{L}$ is invertible for $c>0$ and if $\phi=\phi_{c}$ satisfies (7) then, using (9), we have

[TABLE]

This means that $\phi_{c}$ is an eigenfunction of the iteration operator $\mathcal{L}^{-1}\mathcal{N}^{\prime}(\phi_{c})$ with eigenvalue $p+1>1$. Therefore, for a given initial $\phi_{0}$, the classical fixed point iteration

[TABLE]

will not be, in general, convergent. An alternative iterative method of fixed-point type is the Petviashvili method, [34], which is formulated as

[TABLE]

for some parameter $\epsilon$. The term (10a) is called the stabilizing factor. The convergence of (10) for equations of the form (1) was studied in [33]. Pelinovsky and Stepanyants prove the convergence for $1<\epsilon<(p+2)/p$, under some hypotheses on the spectrum of the linearized operator $\mathcal{L}-\mathcal{N}^{\prime}(\phi_{c})$ at the profile $\phi_{c}$, with the fastest rate of convergence given by $\epsilon^{*}=(p+1)/p$. The inclusion of the stabilizing factor modifies the spectrum of $\mathcal{L}^{-1}\mathcal{N}^{\prime}(\phi_{c})$ in such a way that the eigenvalue $\lambda=p+1$ becomes, for the values of $\epsilon$ considered, an eigenvalue of the iteration operator of (10) at $\phi_{c}$ with magnitude below one (and which is equals zero in the case of choosing $\epsilon=\epsilon^{*}$). The rest of the spectrum is preserved, see [1] and references therein.

The implementation of (10) is typically carried out ´by Fourier pseudospectral approximation, [1, 2]. Let us consider the periodic problem of (7) on a sufficiently long interval $(-l,l)$. This is discretized by a uniform grid $x_{j}=-l+jh,j=0,\ldots,N-1,h=2l/N$, in such a way that (7) can be approximated by the discrete system

[TABLE]

where $\phi_{h}$ is a $N$-vector approximation to the profile $\phi$ at the collocation points $x_{j}$, $I_{N}$ is the $N\times N$ identity matrix and $D_{N}^{\alpha}$ defined as the $N\times N$ matrix

[TABLE]

with $F_{N}$ the discrete Fourier transform matrix on $\mathbb{C}^{N}$ and $\Lambda_{N}^{\alpha}$ the $N\times N$ diagonal matrix with diagonal entries of the form $|k\pi/l|^{\alpha},k=0,\ldots,N-1$, [10, 8]. The dot in (11b) stands for the Hadamard product. Finally, system (11a) is implemented for the discrete Fourier coefficients, $\widehat{\phi}_{h}=(1/N)F_{N}\phi_{h}$, of $\phi_{h}$ and the resulting algebraic system is iteratively solved with the discrete version of (10)

[TABLE]

for $k=0,\ldots,N-1,n=0,1,\ldots$ and where $\widehat{\phi_{n}}=(\widehat{\phi_{n}}(0),\ldots,\widehat{\phi_{n}}(N-1))^{T}$.

A final point concerns the way how the iteration is controlled. This is done by using three strategies:

• •

Since in the case of convergence, the sequence of the stabilizing factors $m_{n}:=m(\phi_{n})$, must go to one, see (12a), then a first control is given by the differences

[TABLE]

• •

A second group of control parameters is given by the sequence of the Euclidean errors between two consecutive iterations

[TABLE]

• •

Finally, the sequence of the residual errors (also in Euclidean norm)

[TABLE]

is also considered.

Thus, the iteration is run up to one of (13)-(15) is below a fixed tolerance parameter $tol$ which, for the experiments below, has been taken as $10^{-10}$.

2.2 Acceleration techniques

As mentioned in the Introduction, the loss of smoothness at the boundary due to the periodic approximation to the equations for the profiles implies an algebraic decrease of the modulus of the Fourier coefficients and, consequently, a slow iteration. In order to improve the velocity of convergence, our proposal here is including some acceleration method in the iterative process, [9, 39]. In this sense, the so-called Vector Extrapolation Methods (VEM), [37, 40, 20], introduce a final stage of extrapolation at the end of each iteration of (12). This is usually carried out in a cycling way: from the last iterate $\psi_{0}=\phi_{n}$ at stage $n$, a number $mw$ (called width of extrapolation) of iterations $\psi_{1},\ldots,\psi_{mw}$ of (12) is computed, and the next iteration $\phi_{n+1}$ is derived as a suitable extrapolation formula from $\psi_{0},\ldots,\psi_{mw}$, see [39] for details. The coefficients of the extrapolation are functions of previous steps of the iteration and their derivation makes use of different criteria, leading to different methods. The one considered for the experiments in this paper is the minimal polynomial extrapolation (MPE), a polynomial method which computes the coefficients by setting orthogonality conditions on the generalized residual, see [37, 20].

The procedure (12), accelerated with MPE is now illustrated by the following numerical results.

3 Numerical experiments

3.1 Some experiments of validation

Some experiments on the performance of the iterative scheme are first presented. Figure 1 shows the form of approximate solitary wave profiles of (4) with $p=1,c=1,mw=6$ and several values of $\alpha$. The initial iteration is a squared hyperbolic secant and for the implementation, an interval with $l=2048$ and $N=2^{18}$ Fourier modes were taken. As it is known (at least for $p=1$, [23]), the more peaked the profile the smaller $\alpha$ is. (This behaviour is independent of the nonlinearity parameter, although as $p$ increases, the amplitude of the corresponding profile has been observed to decrease, see figures 7 and 8 below.) In all the computations, the minimum value of the profile is below $10^{-5}$.

From Figure 1, note that the limiting case $\alpha_{l}=1/3$, for $p=1$, of $\alpha$ looks to be computable. This also holds for any $p>1$, for which $\alpha_{l}=p/(p+2)$, see Figure 2.

The algebraic decay at infinity of some of the profiles can be observed in Figure 3, which displays the corresponding phase portraits and where the derivative has been computed by using pseudospetral differentiation, see [10, 8].

In order to check the accuracy of the computed profiles, several experiments are made. Figure 4 displays the behaviour of the residual error and the stabilizing factor as function of the number of iterations and for each of the waves computed in Figure 1. The results confirm the convergence of the sequences (13), (15) and consequently of the iteration. (The second control sequence (14) behaves even better and the corresponding results are not shown.) As far as the performance is concerned, note that in all the cases, the tolerance $tol=10^{-10}$ is attained in less than $50$ iterations.

A second experiment to check the accuracy is as follows. The computed profiles were taken as initial condition of a numerical method to approximate the periodic initial value problem associated to (1). The numerical scheme is based on a pseudospectral Fourier discretization in space and a fourth order, diagonally implicit, Runge-Kutta composition method, described in [16] (see also references therein) as time integrator. The scheme has relevant geometric properties, [18], and has been shown to be efficient in nonlinear wave problems, [21].

The evolution of the corresponding numerical solution was monitored and, in the case of $\alpha=0.7$, is represented at several times in Figures 5(a),(b). They show that the initial approximate profile evolves in a solitary way without relevant disturbances, suggesting that the computed profile represents a solitary wave of (1) with a high degree of accuracy. This is also confirmed by the evolution of the amplitude and speed of the numerical approximation during the integration, displayed in Figures 5(c) and (d). (The computation of amplitude and speed was made in the standard way, see e. g. [12].)

The following experiments complement the illustration of the performance of the method. Figure 6 shows the number of iterations required to get a residual error below the tolerance as function of the fractional parameter $\alpha$ and for several values of the width of extrapolation $mw$, in the approximation to a solitary wave profile of (4) with $p=1$ and $c=1$. Note that, for moderate values of $mw$, as $mw$ increases, the number of iterations decreases and in this sense the performance improves. However, observe that the parameter $mw$ cannot be fixed a priori in this sort of nonlinear problems, [39, 40]. Also, from some value of $mw$, the number of iterations and the computational time will not improve anymore. For one of the values of $mw$ considered, the number of iterations does not increase as $\alpha$ decreases; it is maximum when $\alpha$ is close to $0.6$.

3.2 Application to study additional properties of the solitary waves

The final group of experiments is concerned with the model (1). Once the accuracy and performance of the iterative method have been checked, the procedure can be used to obtain some additional information about the solitary waves. This idea is focused here on the search for the speed-amplitude relation and its dependence on the parameters $\alpha$ and $p$. Figure 7(a) displays this relation for $\alpha=0.8$ fixed and several values of $p$. The amplitude is an increasing function of the speed but the increment depends on $p$. When $\alpha=2$ (gKdV case) the solitary waves are known explicitly and the amplitude has the exact formula, [7]

[TABLE]

see Figure 7(b).

A comparison of both figures suggests to consider the experiment of studying the speed-amplitude relation for a fixed value of $p$ and as function of $\alpha$. This is illustrated by Figure 8. Note that for a fixed value of $c$ and $\alpha$, increasing $p$ typically leads to a solitary wave of smaller amplitude. When $c$ and $p$ are fixed, the taller the wave the smaller value of $\alpha$ is (cf. Figures 1 and 7).

Due to the peaked form of the profiles, the maximum is, for small values of $\alpha$, not easy to compute with a minimum of accuracy. However, the most reliable numerical results suggest that a similar relation to (16) holds for any $\alpha$. This is observed in Table 3.2 where, for several values of $p$ and $\alpha$, the numerical data speed-amplitude were fitted to a power function $f(x)=ax^{b}$. The accuracy of the results was guaranteed by a goodness of fit where some statistical parameters are around some fixed tolerance. Explicitly, the following statistics were used: the sum of squares due to error SSE (with a tolerance threshold of $10^{-6}$), the R-squared (which in all the cases is equals $1$) and the root mean squared error RSME (around $10^{-4}$).

Note that while the coefficients $b$ suggest an amplitude as a power $1/p$ of the speed as in (16), the dispersion parameter $\alpha$ looks to affect only the coefficient $a$ of the fit.

Acknowledgement(s)

This work was supported by Spanish Ministerio de Economía y Competitividad under the Research Grant MTM2014-54710-P.