Communication in Numerical Analysis

Volume 2012 (2012), Article ID cna-00059, 28 Pages

doi: 10.5899/2012/cna-00059


Research Article


Numerical simulation of GEW equation using RBF collocation method


Hamid Panahipour *


Safahan Institute of Higher Education, Isfahan 81747-43196, Iran


* Corresponding author. h.panahipour@gmail.com Tel: +989131706401


Received: 03 August 2011; Accepted: 19 June 2012


Copyright © 2012 Hamid Panahipour. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The generalized equal width (GEW) equation is solved numerically by a meshless method based on a global collocation with standard types of radial basis functions (RBFs). Test problems including propagation of single solitons, interaction of two and three solitons, development of the Maxwellian initial condition pulses, wave undulation and wave generation are used to indicate the efficiency and accuracy of the method. Comparisons are made between the results of the proposed method and some other published numerical methods.


Keywords: Equal width equation, modified equal width equation, RBF collocation method

1 Introduction

Various methods [1]-[2] have been devised to find the exact and approximate solutions of nonlinear partial differential equations in order to provide more information for understanding many physical phenomena arising in numerous scientific and engineering fields. GEW equation is an important nonlinear wave equation of the form


\begin{aligned} u_t+\epsilon {u^p}{u_x}-\mu u_{xxt}=0, \end{aligned}

where p is a positive integer, \epsilon and \mu are positive constants which require the boundary conditions u\rightarrow 0 as x\rightarrow \pm\infty. It is related to the generalized regularized long wave (GRLW) equation [22, 23] and the generalized Korteweg-de Vries (GKdV) equation [17], based on the equal width (EW) equation [31, 29, 35, 36], and has solitary solutions in pulse-like form. GEW equation was formulated by Peregrine [26, 27], and Benjamin et al. [4] as an alternative to GRLW equation and GKdV equation [22, 23, 5, 6] for studying soliton phenomena and as a model for small-amplitude long waves on the surface of water in a channel. The study of GEW equation provides the opportunity of investigating the creation of secondary solitary waves and/or radiation to get insight into the corresponding processes of particle physics [12, 24]. This equation has many applications in physical situations such as unidirectional waves propagating in a water channel, long waves in near-shore zones, and many others. When $p=1,$ we get a special case, EW equation [31, 29, 35, 36] and when $p=2,$ we get a modified EW (MEW) equation [15]. EW equation was solved numerically by explicit finite difference methods [28], Galerkin method using cubic B-spline finite elements [15], Galerkin and Petrov-Galerkin methods using quadratic B-spline finite elements [16, 32], spectral method based on Chebyshev polynomials [3], a least square technique using linear space-time finite elements and Petrov-Galerkin finite element scheme with shape functions taken as quadratic B-spline functions [35, 36] and the collocation method with quintic B-spline, quartic B-spline and cubic, quadratic, quintic B-spline in [29, 30, 31], respectively. Zaki in [37] solved the MEW equation numerically by a Petrov-Galerkin method using quintic B-spline finite elements and Raslan solved it in [13] using the collocation method by quadratic B-spline at mid points as element shape functions and recently, Saka has applied quintic B-spline collocation algorithms for numerical solution of it in [33].


Unlike the traditional numerical methods in solving partial differential equations, the meshfree or meshless methods need no mesh generation. The collocation methods using RBFs which proposed by Kansa [19, 20], are without any connectivity requirement, hence they are truly meshless and simple to implement as are interesting methods for modelling rather high dimensional problems because of spatial dimension independency [34, 8, 9, 10].


This paper is organized as follows: In the next section, RBF collocation method for solving GEW equation is presented. In section 3, the stability of the method is examined by using a linearized stability analysis. Section 4 is devoted to numerical results. A brief conclusion is presented in the last section of this paper.


2 Construction of the method

We consider GEW equation as


u_t+\epsilon {u^p}{u_x}-\mu u_{xxt}=0,\qquad x\in\Omega=(a,b)\subset\mathbb{R},\qquad t>0,

(2.1)


with the initial condition


u(x,0)=f(x),\qquad x\in\overline{\Omega}=[a,b],

(2.2)


and the boundary conditions


u(a,t)=g_a(t),\qquad u(b,t)=g_b(t),\qquad t\geq0,

(2.3)


where \epsilon and \mu are positive parameters and $f(x),~g_a(t)$ and $g_b(t)$ are known functions. We discretize time derivative of GEW equation using a classic finite difference formula and space derivatives by the \theta-weighted $(0 \leq \theta \leq 1)$ scheme between successive two time levels n and n+1 as


\frac{u^{n+1}-u^n}{\delta t}+\theta (\epsilon u^pu_x)^{n+1}+(1-\theta)(\epsilon u^pu_x)^n-\frac{\mu}{\delta t} \bigg{(}u_{xx}^{n+1}-u_{xx}^{n}\bigg{)}=0,

(2.4)


where $u^n=u(x,t^n),~t^n=t^{n-1}+\delta t$ and \delta t is a time step size.


The nonlinear term $(u^pu_x)^{n+1}$ in Eq. (2.4) can be approximated by using the following formulas which obtained by applying the Taylor expansion,


\begin{aligned} (u^p)^{n+1}&\approx(u^p)^{n}+\delta t\Big{(}(u^p)_t\Big{)}^n\approx(u^p)^{n}+\delta t p\Big{(}u^{p-1}\Big{)}^n\Big{(}\frac{u^{n+1}-u^n}{\delta t}\Big{)}+O(\delta t^2),\\ (u_x)^{n+1}&\approx(u_x)^{n}+\delta t\frac{u_x^{n+1}-u_x^n}{\delta t}+O(\delta t^2), \end{aligned}

thus


\begin{aligned} (u^pu_x)^{n+1}&\approx(u^pu_x)^n+\delta t\bigg{(}(u^p_t)^nu_x^n+(u^p)^nu_{xt}^n\bigg{)}+O(\delta t^2)\\ &=(u^pu_x)^n+\delta t\bigg{(}\frac{(u^p)^{n+1}-(u^p)^{n}}{\delta t}u_x^n+(u^p)^{n}\frac{(u_x)^{n+1}-(u_x)^{n}}{\delta t}\bigg{)}\\ &+O(\delta t^2)\\ &=(u^p)^{n+1}u_x^n+(u^p)^nu_x^{n+1}-(u^p)^nu_x^{n}+O(\delta t^2)\\ &=(u^{p})^nu_x^{n+1}+p(u^{p-1})^{n}u_x^{n}u^{n+1}-p(u^{p})^nu_x^{n}+O(\delta t^2). \end{aligned}

(2.5)


So Eq. (2.4) can be rewritten as


\begin{aligned} u^{n+1}+\delta t\theta\epsilon\bigg{(}(u^n)^pu_x^{n+1}+p(u^{n})^{p-1}u_x^{n}u^{n+1}\bigg{)}-\mu u_{xx}^{n+1}\\ =u^{n}+\delta t\epsilon\bigg{(}((p+1)\theta-1)(u^n)^pu_x^{n}\bigg{)}-\mu u_{xx}^{n} \end{aligned}

(2.6)


Now, let us choose the collocation points $x_i,~i=1,\ldots,N$ over $\overline{\Omega}$ such that $x_i,~i=2,\ldots,N-1$ are interior points and $x_i,~i=1,N$ are boundary points and then apply the following approximation


u(x,t^n)=u^n(x)\simeq\sum_{j=1}^N\lambda_j^n\phi_j(x),

(2.7)


where \lambda_j^n are the coefficients to be determined later and $\phi_j(x)=\varphi(\parallel x-x_j\parallel).~\varphi(x)$ is a form of the RBFs presented in Table 1, where c is known as the shape parameter of the RBFs which is found experimentally for each RBF and for each problem separately in order to obtain accurate solution and m=2 in our case in TPS function. Moreover $\parallel\cdot\parallel$ denotes the Euclidean norm.


Table 1: some well-known radial basis functions


Therefore for each collocation point $x_i,$ Eq. (2.7) can be written as


\begin{aligned} u^n(x_i)=\sum_{j=1}^N{\lambda_j^{n}\phi_j(x_i)},\qquad i=1,\ldots,N, \end{aligned}

(2.8)


with the following matrix form


\begin{aligned} u^n={\bf A}\lambda^n, \end{aligned}

(2.9)


where


\begin{aligned} {\bf A}=[\phi_j(x_i):i=1,\ldots,N,\ j=1,\ldots,N]_{N\times N}, \qquad \lambda^{n}=[\lambda_1^{n},\cdots,\lambda^{n}_N]^T. \end{aligned}

Now we can spilt the matrix \bf{A} into two matrices \bf{A}_{d} and \bf{A}_{b} corresponding to N-2 interior points and two boundary points as follows


\begin{aligned} {\bf A}={\bf A}_{d}+{\bf A}_{b}, \end{aligned}

where


\begin{aligned} {\bf A}_{d}&=[\phi_j(x_i):i=2,\ldots,N-1,\ j=1,\ldots N\mbox{\and 0 elsewhere}]_{N\times N},\\ {\bf A}_{b}&=[\phi_j(x_i):i=1,N,\ j=1,N\mbox{\ and 0 elsewhere}]_{N\times N}. \end{aligned}

Now by substituting Eq. (2.8) into Eqs. (2.6) and (2.3) and using the collocation points $x_i,~i=1,\ldots,N,$ we obtain the following equations


\begin{aligned} &\sum_{j=1}^N{\lambda_j^{n+1}\phi_j(x_i)}-\mu{\sum_{j=1}^N{{\lambda_j^{n+1}}{\phi''_j(x_i)}}}+\delta t\theta\epsilon\Bigg{(}\Big{(}{\sum_{j=1}^N{\lambda_j^{n}\phi_j(x_i)}}\Big{)}^p\Big{(}\sum_{j=1}^N{\lambda_j^{n+1}\phi'_j(x_i)}\Big{)}\\ &+p\Big{(}\sum_{j=1}^N{\lambda_j^{n}\phi_j(x_i)}\Big{)}^{p-1}\Big{(}\sum_{j=1}^N{\lambda_j^{n}\phi'_j(x_i)}\Big{)}\Big{(}\sum_{j=1}^N{\lambda_j^{n+1}\phi_j(x_i)}\Big{)}\Bigg{)}\\ &=\sum_{j=1}^N{\lambda_j^{n}\phi_j(x_i)}-\mu{\sum_{j=1}^N{{\lambda_j^{n}}{\phi''_j(x_i)}}}\\ &+\delta t\theta\epsilon\Bigg{(}((p+1)\theta-1)\Big{(}{\sum_{j=1}^N{\lambda_j^{n}\phi_j(x_i)}}\Big{)}^p\Big{(}\sum_{j=1}^N{\lambda_j^{n}\phi'_j(x_i)}\Big{)}\Bigg{)}, \qquad i=2,\ldots,N-1 \end{aligned}

and


\begin{aligned} &\sum_{j=1}^N{{\lambda_j^{n+1}}{\phi_j(x_1)}}=g_a(t^{n+1}),\\ &\sum_{j=1}^N{{\lambda_j^{n+1}}{\phi_j(x_N)}}=g_b(t^{n+1}), \end{aligned}

where $\phi'_j(x_i)=\frac{d}{dx}\phi_j(x)|_{x=x_i}$ and $\phi''_j(x_i)=\frac{d^2}{dx^2}\phi_j(x)|_{x=x_i}.$


These equations generate a system of N linear equations in N unknown parameters \lambda_j^{n+1} which can be expressed in a matrix form as follows


\begin{aligned} {\bf M}\lambda^{n+1}={\bf R}, \end{aligned}

(2.10)


where


\begin{aligned} {\bf M}&=[{\bf A}_d+{\bf A}_b-\mu{\bf C}+\theta\delta t\epsilon\{{\bf E}+{\bf D}\}],\\ {\bf R}&=[{\bf A}_d-\mu{\bf C}+\delta t\{\epsilon((p+1)\theta-1){\bf E}\}]\lambda^n+{\bf F}^{n+1}, \end{aligned}

and


\begin{aligned} &{\bf B}=[\phi'_j(x_i):i=2,\ldots,N-1,\ j=1,\ldots N\mbox{\ and 0 elsewhere}]_{N\times N},\\ &{\bf C}=[\phi''_j(x_i):i=2,\ldots,N-1,\ j=1,\ldots N\mbox{\ and 0 elsewhere}]_{N\times N},\\ &{u_x^n}={\bf B}\lambda^n,{\bf D}=pu^{n^{p-1}}*u^n_x*{\bf A}_d,{\bf E}=u^{n^p}*{\bf B},\\ &{\bf F}^{n+1}= [g_a(t^{n+1})\ 0\ \cdots\ 0\ g_b(t^{n+1})]^T, \end{aligned}

where the symbol "*" means the componentwise multiplication.


Remark 1. Although the invertibility of matrix \bf{M} cannot be proved in general [18], it has been proved for \bf{A} by Micchelli for distinct collocation points [25]. However, Kansa argued that Micchelli's proof still holds in this case as well, since solving differential equations is a special case of generalized interpolation problem [19]. Furthermore, the singularities are rare in practice.


Remark 2. Eq. (2.10) represents a system of N linear equations in N unknown parameters $\lambda_j.$ This system can be solved by the Gaussian elimination method with partial pivoting.


Remark 3. The collocation matrix corresponding to TPS function becomes highly ill-conditioned because of a singularity at $r_{ij}$ where the sets of centers and collocation points coincide. In this case, we use the limiting value $\lim_{r\rightarrow 0}{r^4\log{(r)}}=0,$ to overcome the difficulty.


Now we are able to extract the following algorithm.


Algorithm of the method


  1. step(1): Choose the parameters $\delta t,~T$ and \theta such that $0\leq\theta\leq1.$

  2. step(2): Choose N collocation points in \overline\Omega.

  3. step(3): Set n:=0.

  4. step(4): Calculate the initial solution u^n using Eq. (2.2).

  5. step(5): Solve the linear system \bf{A}\lambda^n=u^n.

  6. step(6): Calculate the Matrix \bf{M} and vector \bf{R}.

  7. step(7): Set n:=n+1.

  8. step(8): Calculate the parameters \lambda_j^{n} using Eq. (2.10)

  9. step(9): Calculate the approximate solution u^{n} using Eq. (2.9).

  10. step(10): If n\,\delta t< T go to step 6 else stop.

Remark 4. The proposed method is valid for any value of $\theta \in [0,1],$ but we use $\theta=\frac{1}{2}.$ In this case, by taking into account of Eq. (2.5) and the second order accuracy of the Crank-Nicholson scheme in time, it can be deduced that the proposed method is second order accurate in time.


3 Stability analysis

In this section, following [34], the stability of the RBF approximation (2.10) is investigated by using the matrix method. Eq. (2.1) can be linearized by assuming the quantity u^p in the nonlinear term $u^pu_x$ as locally constant \hat{u}. The error e^n at the nth time level is given by


\begin{aligned} e^n=u^n_{\rm{exact}}-u^n_{\rm{app}}, \end{aligned}

where $u^n_{\rm{exact}},~u^n_{\rm{app}}$ are the exact and the numerical solutions at the n-th time level, respectively. The error equation for the linearized GEW equation can be written as


[{\bf H}+\theta\delta t{\bf K}]e^{n+1}=[{\bf H}-\delta t(1-\theta){\bf K}]e^{n},

(3.11)


where $\bf{H}=[\bf{A}-\mu\bf{C}]\bf{A}^{-1}$ and \bf{K}=\epsilon\bf{E}\bf{A}^{-1}. Eq. (3.11) can be written as


\begin{aligned} e^{n+1}={\bf P}e^n, \end{aligned}

where $\bf{P}=[\bf{H}+\theta\delta t\bf{K}]^{-1}[\bf{H}-\delta t(1-\theta)\bf{K}].$ The numerical scheme is stable if $\|\bf{P}\|_2\leq 1,$ which is equivalent to $\rho(\bf{P})\leq 1,$ where $\rho(\bf{P})$ denotes the spectral radius of the matrix \bf{P}. From Eq. (3.11), it can be seen that the stability is assured if all the eigenvalues of the matrix $[\bf{H}+\theta\delta t\bf{K}]^{-1}[\bf{H}-\delta t(1-\theta)\bf{K}$ satisfy the following condition


\begin{aligned} \bigg{|}\frac{\lambda_H-\delta t(1-\theta)\lambda_K}{\lambda_H+\delta t\theta\lambda_K}\bigg{|}\leq 1, \end{aligned}

(3.12)


where $\lambda_H$ and $\lambda_K$ are eigenvalues of the matrices \bf{H} and $\bf{K},$ respectively. When $\theta=0.5,$ the inequality (3.12) becomes


\begin{aligned} \bigg{|}\frac{\lambda_H-0.5\delta t\lambda_K}{\lambda_H+0.5\delta t\lambda_K}\bigg{|}\leq 1. \end{aligned}

(3.13)


In the case of complex eigenvalues $\lambda_{H}=a_h+ib_h$ and $\lambda_K=a_k+ib_k,$ where $a_h, a_k, b_h$ and $b_k$ are any real numbers, the inequality (3.13) takes the following form,


\begin{aligned} \bigg{|}\frac{(a_h-0.5\delta ta_k)+\rm{i}(b_h-0.5\delta tb_k)}{(a_h+0.5\delta ta_k)+\rm{i}(b_h+0.5\delta tb_k)}\bigg{|}\leq 1. \end{aligned}

(3.14)


The inequality (3.14) is satisfied if $a_ha_k+b_hb_k\geq0.$ For real eigenvalues, the inequality (3.13) holds true if either $(\lambda_h\geq 0$ and $\lambda_k\geq 0)$ or $(\lambda_h\leq 0$ and $\lambda_k\leq 0).$ This shows that the scheme (2.10) is unconditionally stable if $a_ha_k+b_hb_k\geq0,$ for complex eigenvalues and if either $(\lambda_h\geq 0$ and $\lambda_k\geq 0)$ or $(\lambda_h\leq 0$ and $\lambda_k\leq 0),$ for real eigenvalues. When $\theta=0,$ the inequality (3.12) becomes


\begin{aligned} \bigg{|}1-\frac{\delta t\lambda _K}{\lambda_H}\bigg{|}\leq 1, \end{aligned}

i.e.


\begin{aligned} \delta t\leq\frac{2\lambda_H}{\lambda_K}\qquad \rm{and}\qquad\frac{\lambda_H}{\lambda_K}\geq 0. \end{aligned}

Thus for $\theta=0,$ the scheme is conditionally stable.


The stability of the scheme (2.10) and conditioning of the component matrices $\bf{H},~\bf{K}$ of the matrix \bf{P} depend on the weight parameter \theta, the minimum distance between any two collocation points h in the domain set $[a,b],$ and the local shape parameter \varepsilon. Fornberg et al. [14] suggested that it can be advantageous to let the shape parameter vary spatially, rather than assigning a single value to it by considering the Runge phenomenon as a key error mechanism. Benefits typically include improvements in both accuracy and numerical conditioning. Still another benefit arises if one wishes to improve local accuracy by clustering nodes in selected areas. Cheng et al. [7] showed that when \varepsilon is very small then the RBFs system error is of exponential order. But there is a certain limit for the value \varepsilon after which the solution breaks down. For the limiting value of \varepsilon the condition number of the RBFs system becomes so large that the system leads to ill-conditioning. In the case of an ill-conditioned system, the numerical solution produced is not stable. Relation between the condition number of the matrix \bf{P} and the different values of the shape parameter \varepsilon is shown in Table 2 for EW equation and MQ case only. The critical value of the shape parameter \varepsilon in this case is $0.4$ and the condition number of the matrix \bf{P} is 5.046\times 10^{15}. It is clear that if the values of the shape parameter \varepsilon are smaller than the critical value, then the solution breaks down and hence the method becomes unstable. It can be seen from Table 2 that the RBFs approximation is not very sensitive to the values of the shape parameter \varepsilon. In particular, the method can tolerate a rather wide range of values of \varepsilon. The interval of stability in this case is $(0.4,10).$ In case of the parameter free RBFs such as TPS, the stability and conditioning depend on the weight parameter \theta, the eigenvalues $\lambda_H,~\lambda_K$ and h.


Table 2: Condition number versus shape parameter at $t=20,~N=120,~\delta t=0.25,~\nu=0.1,~0\leq x\leq 30.$


4 Numerical results

GEW equation (2.1) has the analytical solution given by


\begin{aligned} u(x,t)=\sqrt[p] {\frac{\nu(p+1)(p+2)}{2\epsilon}\mbox{sech}^2\bigg{[}\frac{p}{2\sqrt\mu}(x-\nu t-x_0)\bigg{]}}, \end{aligned}

which corresponds to a solitary wave of amplitude \sqrt[p]{\frac{\nu(p+1)(p+2)}{2\epsilon}}, speed $\nu,$ width $\frac{p}{2\sqrt\mu}$ and initially centered at $x_0.$ With the zero boundary conditions, solutions of GEW equation possess three invariants of the motion given by


\begin{aligned} &I_1=\int_a^budx,\qquad I_2=\int_a^b{(u^2+\mu u^2_x)}dx,\qquad I_3=\int_a^bu^{p+2}dx. \end{aligned}

Corresponding to mass, momentum and energy, respectively. In this section the proposed algorithm using different test problems related to the propagation of single solitary waves, the interaction of solitons and the Maxwellian initial condition into solitary waves for the EW $(p=1)$ and MEW $(p=2)$ equations by using the homogenous boundary conditions, are examined. Furthermore, to model the effect of non-homogenous boundary conditions, the development of the undular bore and the effect of forcing boundary conditions on the generation of solitary waves for EW equation, is examined. We use the following error norms to compare our numerical solution with the exact solution,


\begin{aligned} &L_2=\parallel u-\tilde{u} \parallel_2={\sqrt {h \sum_{j=1}^{N}{\mid u_j-\tilde{u}_j \mid}^2}},\\ &L_{\infty}=\parallel u-\tilde{u} \parallel_{\infty}=\max_{1\leq j\leq N}{\mid u_j-\tilde{u}_j\mid}, \end{aligned}

where u and \tilde{u} represent the exact and approximate solutions, respectively, and h is the minimum distance between any two collocation points of the domain \overline\Omega. The quantities $I_1,~I_2$ and $I_3$ are also applied to measure the conservation properties of the collocation scheme. The pointwise rate of convergence in time is also calculated by using the following formulae:


\begin{aligned} \rm{order}=\frac{\log_{10}(||u_{exact}-u_{{\delta t}_j}||/||u_{exact}-u_{{\delta t}_{j+1}}||)}{\log_{10}({\delta t}_{j}/{\delta t}_{j+1})} \end{aligned}

where u_{exact} is exact solution and $u_{{\delta t}_j}$ is the numerical solution with time step size ${\delta t}_j.$


4.1 Propagation of single solitary waves


4.1.1 EW equation


In our computational work in this case, we model the motion of a single solitary wave of EW equation with the three different amplitudes $0.3,~0.09$ and $0.03$ to make possible the comparisons between the collocation method using quartic spline finite elements [30], the least square method using linear spline elements [35] and the Petrov-Galerkin method using cubic spline finite elements [16] with the parameters $\mu=1,~x_0=10,~\epsilon=1$ and \delta t=0.25 through the interval $[0,30]$ with N=120 cells. The analytical values of conservation quantities can be found as


\begin{aligned} I_1=12\frac{\nu}{\epsilon}\sqrt{\mu},\qquad I_2=\frac{144}{5}\frac{\nu^2}{\epsilon^2}\sqrt{\mu},\qquad I_3=\frac{288}{5}\frac{\nu^3}{\epsilon^3}\sqrt{\mu}. \end{aligned}

The values of $L_\infty,~L_{2}$ and the invariants $I_1,~I_2$ and $I_3$ using five RBFs, i.e., MQ, IMQ, IQ, GA and TPS are given in Tables 3-4 up to time T=80. We find that our scheme is more accurate than the methods mentioned below and is in agreement with the method in [3].


  1. (i) The case with amplitude $0.3,$ the analytic values of the invariants are $I_1=1.2,~I_2=0.288$ and $I_3=0.0576.$ Referring to Table 3 at the time t=40 the value of L_2 is less than 4.914493\times 10^{-5} for the MQ, IMQ, IQ and GA types of RBFs, which is smaller than the results in [30, 35]. At the time t=80 the value of L_2 is less than 1.947625\times 10^{-4} for the all types of RBFs, which is smaller than the results in [35, 16]. Whereas the changes of the invariant $I_1$ is less than 3.8\times 10^{-4} and the invariants $I_2$ and $I_3$ are conserved during the experiment and they are all very close to the analytic values, our scheme is satisfactorily conservative.

  2. (ii) The case with amplitude $0.09,$ the analytic values of the invariants are $I_1=0.36,~I_2=0.02592,$ and $I_3=0.001555.$ Referring to Table 4 at the time t=80 the values of L_\infty and L_2 are less than 1.790504\times 10^{-5} and $2.235945\times 10^{-5},$ respectively, for the all types of RBFs, which are smaller than the results in [35]. Whereas $I_1,~I_2$ and $I_3$ are constants during the experiment and very close to the analytic values, our scheme is satisfactorily conservative.

  3. (iii) The case with amplitude $0.03,$ the analytic values of the invariants are $I_1=0.12,~I_2=0.00288$ and $I_3=0.000058.$ Referring to Table 4 at the time t=80 the values of L_\infty and L_2 are less than 2.447841\times 10^{-6} and $3.536053\times 10^{-6},$ respectively, for the all types of RBFs, which are smaller than the results in [35]. Whereas $I_1,~I_2$ and $I_3$ are constants during the experiment and very close to the analytic values, our scheme is satisfactorily conservative.

Table 3: Invariants and errors for single soliton of EW equation; $h=0.25,$ amplitude $0.3,~\delta t=0.25,~0\leq x\leq 30.$


Table 4: Invariants and errors for single soliton of EW equation; $h=0.25,~\delta t=0.25,~0\leq x\leq 30.$


The numerical solutions with five RBFs and also the exact solution at t=20 are graphed in Figure 1 for amplitude $3,~h=\delta t=0.2$ and $x_0=15$ over the region $[0,80].$ The solution curves are indistinguishable. In Table 5, the number of the collocation points is kept fixed at N=200 and the time step size ${\delta t}_i=0.8, 0.4, 0.2, 0.1, 0.05$ is varied to compute the time rate of convergence for each of the RBFs approximation. It can be noted from Table 5, that the rate of convergence increases with the smaller time step size and the method is convergent of order 2 in time which confirms the assertion in Remark 4.


Figure 1: Solitonic solutions of EW (left) and MEW (right) equations; $0\leq x\leq 80,$ time=$20.$


Table 5: Time rate of convergence for single soliton of EW equation at $t=20,~N=200,~\nu=0.1,~\varepsilon=1,~0\leq x\leq 30.$


4.1.2 MEW equation


For the numerical work in this case, we put $\nu=\frac{1}{32},~h=0.1,~\mu=1,~x_0=30,~\epsilon=3$ and $\delta t=0.05,$ through the interval $[0,70]$ to make possible the comparisons between the collocation method based on quadratic B-splines [13], Petrov-Galerkin method using quintic B-spline finite elements [37] and the quintic B-spline collocation algorithms [33], which are found in [33]. The analytical values of conservation quantities can be found as


\begin{aligned} &I_1=\pi\sqrt{\frac{6\nu}{\epsilon}}=0.785398,\qquad I_2=\frac{16\nu}{\epsilon}=0.166667,\\ &I_3=\frac{48\nu^2}{\epsilon^2}=0.00520833. \end{aligned}

The values of $L_\infty,~L_{2}$ and the invariants $I_1,~I_2$ and $I_3$ for this case using five RBFs, i.e., MQ, IMQ, IQ, GA and TPS are given in Table 6. Referring to Table 6, it can be noted that the we obtained smaller errors than the methods mentioned and the conservation quantities are all exactly equal to the analytical values, throughout the simulation. The numerical solutions with five RBFs and also the exact solution at t=20 are graphed in Figure 1 for h=0.1 and \delta t=0.2 over the region $[0,80].$ The solution curves are indistinguishable. In Table 7, the number of the collocation points is kept fixed at N=700 and the time step size ${\delta t}_i=0.8, 0.4, 0.2, 0.1, 0.05$ is varied to compute the time rate of convergence for each of the RBFs approximation. It can be noted from Table 7, that the rate of convergence increases with the smaller time step size and the method is convergent of order 2 in time which confirms the assertion in Remark 4.


Table 6: Invariants and errors for single soliton solution of MEW equation; $\nu=\frac{1}{32},~h=0.1,~\delta t=0.05,~0\leq x\leq 70$


Table 7: Time rate of convergence for single soliton of MEW equation at $t=20,~N=700,~\nu=\frac{1}{32},~\varepsilon=2.65,~0\leq x\leq 70.$


4.2 Interaction of two solitary waves


Our second experiment pertains to the interaction of two solitary solutions of GEW equation having different amplitudes and travelling in the same direction. We consider GEW equation with initial conditions given by the linear sum of two well separated solitary waves of various amplitudes as follows


\begin{aligned} u(x,0)=\displaystyle\sum_{i=1}^2 \sqrt[p] {\frac{\nu_i(p+1)(p+2)}{2\epsilon}\mbox{sech}^2\bigg{[}\frac{p}{2\sqrt\mu}(x-x_i)\bigg{]}}. \end{aligned}

4.2.1 EW equation


We study the interaction of two solitary solutions of EW equation with amplitude ratio of $4\!:\!1,$ so we take $\nu_1=2,~\nu_2=0.5,~x_1=10,~x_2=25,~h=0.1,~\delta t=0.1,~\mu=1$ and \epsilon=1 over the region $[0,80].$ The analytical values of the conservation quantities of this case can be found as


\begin{aligned} &I_1=12(\nu_1+\nu_2)=30,\qquad I_2=28.8(\nu_1^2+\nu_2^2)=122.4,\\ &I_3=57.6(\nu_1^3+\nu_2^3)=468. \end{aligned}

The values of $I_1,~I_2,$ and $I_3$ throughout the simulation are shown in Table 8. It can be seen that these quantities are very close to the corresponding analytical values. Figure 2 shows the interaction of these two solitary waves at different times for all types of RBFs. The curves are indistinguishable.


Table 8: Invariants for interaction of two and three solitary solutions of EW equation; $h=0.1,~\delta t=0.1,~0\leq x\leq 80$


Figure 2: Interaction of two solitary solutions of EW equation at selected time levels.


4.2.2 MEW equation


We study the interaction of two solitary solutions of MEW equation with amplitude ratio of $2\!:\!1,$ so we take $\nu_1=\frac{1}{2},~\nu_2=\frac{1}{8},~h=0.1,~\delta t=0.2,~x_1=15,~x_2=30,~\mu=1$ and \epsilon=3 over the region $[0,80].$ The analytical values of the conservation quantities can be found as


\begin{aligned} &I_1=\pi\Bigg({}\sqrt{\frac{6\nu_1}{\epsilon}}+\sqrt{\frac{6\nu_2}{\epsilon}}\Bigg{)}= 4.71239,\qquad I_2=\frac{16(\nu_1+\nu_2)}{\epsilon}=3.33333,\\ &I_3=\frac{48(\nu_1^2+\nu_2^2)}{\epsilon^2}=1.41667. \end{aligned}

The three invariants in this case are recorded in Table 9. It is clear that the quantities are reasonably constant, since the changes in $I_2$ and $I_3$ are less than 2.07\times 10^{-3} and $2.24\times 10^{-3},$ respectively, for all types of RBFs and the changes in $I_1$ approach to zero. In Figure 4, we show the interaction of these two solitary waves at different times for all types of RBFs. The curves are indistinguishable.


Table 9: Invariants for interaction of two and three solitary solutions of MEW equation; $h=0.1,~\delta t=0.2,~0\leq x\leq 80$


4.3 Interaction of three solitary waves


Now the interaction of three solitary solutions of GEW equation with different amplitudes and travelling in the same direction is presented. We consider GEW equation with initial conditions given by the linear sum of three well separated solitary waves of various amplitudes as follows


\begin{aligned} u(x,0)=\displaystyle\sum_{i=1}^3 \sqrt[p] {\frac{\nu_i(p+1)(p+2)}{2\epsilon}\mbox{sech}^2\bigg{[}\frac{p}{2\sqrt\mu}(x-x_i)\bigg{]}}. \end{aligned}

4.3.1 EW equation


We study the interaction of three solitary solutions of EW equation with amplitude ratio of $9\!:\!3\!:\!1,$ so we take $\nu_1=4.5,$ $\nu_2=1.5,$ $\nu_3=0.5,$ $h=0.1,$ $\delta t=0.1,$ $x_1=10,$ $x_2=25,$ $x_3=35,$ $\mu=1$ and \epsilon=1 over the region $[0,80].$ The analytical values of the conservation quantities can be found as follows


\begin{aligned} &I_1=12(\nu_1+\nu_2+\nu_3)=78,\qquad I_2=28.8(\nu_1^2+\nu_2^2+\nu_3^2)=655.2,\\ &I_3=57.6(\nu_1^3+\nu_2^3)=5450.4. \end{aligned}

The three invariants in this case are shown in Table 8. We find from our numerical scheme, that the invariants $I_1,~I_2$ and $I_3$ for the interaction of these solitary waves are sensible constants. In Figure 3, we show the interaction of these three solitary waves at different times for all types of RBFs. The curves are indistinguishable.


Figure 3: Interaction of three solitary solutions of EW equation at selected time levels.


Figure 4: Interaction of two solitary solutions of MEW equation at selected time levels.


4.3.2 MEW equation


We study the interaction of three solitary solutions of MEW equation with amplitude ratio of $4\!:\!2\!:\!1,$ so we take $\nu_1=\frac{1}{2},$ $\nu_2=\frac{1}{8},$ $\nu_3=\frac{1}{32},$ $h=0.1,$ $\delta t=0.2,$ $x_1=15,$ $x_2=30,$ $x_3=45,$ $\mu=1$ and \epsilon=3 over the region $[0,120].$ The analytical values of the conservation quantities can be found as follows


\begin{aligned} &I_1=\pi\Bigg({}\sqrt{\frac{6\nu_1}{\epsilon}}+\sqrt{\frac{6\nu_2}{\epsilon}}+\sqrt{\frac{6\nu_3}{\epsilon}}\Bigg{)}= 5.49779,\qquad I_2=\frac{16(\nu_1+\nu_2+\nu_3)}{\epsilon}=3.5,\\ &I_3=\frac{48(\nu_1^2+\nu_2^2+\nu_3^2)}{\epsilon^2}=1.42188. \end{aligned}

The three invariants in this case are shown in Table 9. We find from our numerical scheme, that the invariants $I_1,~I_2$ and $I_3$ for the interaction of these solitary waves are sensible constants, since the changes in $I_1,~I_2$ and $I_3$ are less than $1.4 \times 10^{-5},~6.81 \times 10^{-3}$ and $6.55 \times 10^{-3},$ respectively, for all types of RBFs.


4.4 Maxwellian initial condition


4.4.1 EW equation


In final series of numerical experiments of EW equation the development of the Maxwellian initial condition


\begin{aligned} u(x,0)=e^{-(x-7)^2}, \end{aligned}

with boundary conditions


\begin{aligned} u(x,t)=0,~ \mbox{as}~~ x\longrightarrow\pm \infty, ~t>0, \end{aligned}

into a sequence of solitary waves is examined. Since the behavior of the solution depends on the values of $\mu,$ we consider different values for \mu and run the simulation up to time t=4 with $h=0.1,~\delta t=0.1,~\epsilon=1$ and interval $[0,80].$ The recorded values of the quantities $I_1,~I_2$ and $I_3$ for $\mu=0.04$ and $\mu=0.2$ are given in Table 10. We find that the changes in these invariants for $\mu=0.04$ are less than $0.5\times 10^{-5},~2.43\times 10^{-3}$ and $3.87\times 10^{-3},$ respectively, and the changes for $\mu=0.2,$ are less than $0.2\times 10^{-5},~8.3\times 10^{-5}$ and $7.6\times 10^{-5},$ respectively, for all types of RBFs.


Table 10: Invariants of EW equation using Maxwellian initial condition; $h=0.1,~\delta t=0.1,~0\leq x\leq 80$


4.4.2 MEW equation


We have examined the evolution of an initial Maxwellian pulse into solitary waves of MEW equation, using initial condition of the form


\begin{aligned} u(x,0)=e^{-x^2}, \end{aligned}

with boundary conditions


\begin{aligned} u(x,t)=0,~ \mbox{as}~~ x\longrightarrow\pm \infty, ~t>0, \end{aligned}

into a sequence of solitary waves is examined. We consider different values for \mu and run the simulation up to time t=12.51 with $h=0.05,~\delta t=0.01,~\epsilon=3$ and interval $[-20,20].$ When \mu is reduced, more and more solitary waves are formed. For $\mu\!=\!1$ the Maxwellian develops into a pair of waves such that, one wave has negative amplitude while the other one has positive as indicated in Figure 5. Clearly, the additional small waves may occur between the two main solitary waves. For $\mu=0.5,$ in a different behavior, we have an oscillatory form and no clean solitons are formed. For $\mu=0.1,$ a clean solitary wave is observed. For $\mu=0.05$ two solitons are generated. We have three solitary waves for the case $\mu=0.2$ and finally, for the case $\mu=0.005,$ the Maxwellian initial condition has decayed into seven solitary waves. The recorded values of the quantities $I_1,~I_2$ and $I_3$ for $\mu=1$ and $\mu=0.5$ are given in Table 11, we find that the changes in these invariants for $\mu=1$ are less than $0.1\times 10^{-5},~0.5\times 10^{-5}$ and $0.13\times 10^{-5},$ respectively, and the changes in $I_2$ and $I_3$ for $\mu=0.5,$ are less than 0.2\times 10^{-5} and $0.12\times 10^{-5},$ respectively, whereas, the changes in $I_1$ approach to zero. Invariants in the cases $\mu=0.1$ and $\mu=0.05$ are reported in Table 11, in this table the changes in $I_2$ and $I_3$ for $\mu=0.1$ are less than 1.6\times 10^{-5} and 1.7\times 10^{-5}, respectively, whereas, the changes in $I_1$ approach to zero and the changes in these values for $\mu=0.05$ are less than $0.2\times 10^{-5},~7.8\times 10^{-5}$ and $1.19\times 10^{-5}$ , respectively.


Table 11: Invariants of MEW equation using Maxwellian initial condition; $h=0.05,~\delta t=0.01,~-20\leq x\leq 20$


Figure 5: Maxwellian initial condition of MEW equation with (a) $\mu=1,~t=12.51,$ (b) $\mu=0.5,~t=12.51,$ (c) $\mu=0.1,~t=12.51,$ (d) $\mu=0.05,~t=12.51,$ (e) $\mu=0.02,~t=12.51$ and (f) $\mu=0.005,~t=6.51.$


4.5 Wave undulation


The development of the undular is studied by using the initial condition


\begin{aligned} u(x,0)=\frac{U_0}{2}\left(1-\tanh\left(\frac{x-x_c}{d}\right)\right) \end{aligned}

and the boundary conditions $u(a,t)=U_0,~u(b,t)=0,$ where $u(x,0)$ shows the elevation of the water above the equilibrium surface at time $t=0,~d$ represents the slope between the still water and deeper water, and the change in water level of magnitude $U_0$ is centered on $x=x_c.$ Using the present scheme, the simulation is run until time t=800 in the range -20\leq x\leq 50 with the parameters $\mu=0.1666667,~U_0=0.1,~x_c=0,$ the gentle slope d=5 and the steep slope d=2. Variation of the invariants $I_1,~I_2~I_3,$ position and magnitude of the leading undulation agree with that referenced in the study [32]. The undulation profiles at times $t=200,~400,~600$ and 800 are graphed in Fig. 6 for d=5 and in Fig. 7 for d=2 with the TPS type of RBFs.


Figure 6: Undulation profiles for d=5 at different times.


Figure 7: Undulation profiles for d=2 at different times.


4.6 Wave generation


To model the effect of a wave maker, we study the generation of solitary waves within a semi-infinite region and using boundary conditions


\begin{aligned} u(a,t)=\left\{\begin{array}{ll} U_0\frac{t}{\tau},&\qquad 0\leq t\leq\tau,\\ U_0, &\qquad \tau< t< t_0-\tau,\\ U_0\frac{t_0-t}{\tau}, &\qquad t_0-\tau\leq t\leq t_0,\\ 0& \qquad\rm{otherwise},\end{array}\right. \end{aligned}

and $u(b,t)=0.$ Solitary waves will be continually generated until the forcing is off. The parameters are chosen as $h=0.4,~\delta t=0.1,~t_0=20,~\tau= 0.1$ and $U_0=2$ over the range $0\leq x\leq 260.$ During run time of the algorithm, magnitude of the waves agree with that referenced in the study [32]. With no forcing, wave generation stopped and these waves moved to the right. The numerical solutions at times t=100 and 200 are depicted in Fig. 8 with the TPS type of RBFs.


Figure 8: Solitary waves produced by boundary forcing of duration $t_0=20$ and $h=0.4,~\delta t=0.1,~\tau=0.1.$


5 Conclusion

It has been shown that the collocation method using five different types of RBFs, i.e., MQ, IMQ, IQ, GA and TPS, leads to a more accurate method than the rather reported methods for solving GEW equation. We tested our scheme through single solitary wave, the interaction of solitons, and the Maxwellian initial condition. Wave undulation and wave generation induced by the boundary forcing were also used. It is worthy of mention that the mesh free feature of technique is reserved due to capability of using non-uniform meshes. It should be noted that the technique is rather computationally costly, since it is $O(N^3)$ per time step due to recomputing the (non-sparse) matrix at each time step. Nevertheless, we can conclude with confidence, that the collocation method using all types of RBFs can be considered as a highly accurate numerical method for solving these kinds of nonlinear partial differential equations.


Acknowledgments

The author would like to thank the unknown referees for their careful reading and helpful comments.



References

  1. Hamid Panahipour, Application of Extended Tanh Method to Generalized Burgers-type Equations, Communications in Numerical Analysis, Volume 2012 (2012), 1-14.

  2. P. Donald Ariel, On a second parameter in the solution of the flow near a rotating disk by homotopy analysis method, Communication in Numerical Analysis, Volume 2012 (2012) , 1-13.

  3. A. H. A. Ali, Spectral method for solving the equal width equation based on Chebyshev polynomials, Nonlinear Dyn. 51 (2008), 59-70.

  4. T.B. Benjamin, J.L. Bona, J.J. Mahony, Model equations for waves in nonlinear dispersive systems, Phil. Trans. Roy. Soc. (London). A 227 (1972), 47-78.

  5. J.L. Bona, W.G. Pritchard, L.R. Scott, A comparison of solutions of two model equations for long waves, In: Norman R. Lebovitz (Ed.), Fluid Dynamics in Astrophysics and Geophysics, Lectures in Applied Mathematics, (1983), 235-267.

  6. J.L. Bona, W.G. Pritchard, L.R. Scott, An evaluation for water waves, Phil. Trans. Roy. Soc. (London). A 302 (1981), 457-510.

  7. A.H.D. Cheng, M.A. Golberg, E.J. Kansa, G. Zammito, Exponential convergence and H-c multiquadric collocation method for partial differential equations, Numer. Methods Partial Differential Equations. 19 (2003), 571-5-94.

  8. İ. Daǧ, Y. Dereli, Numericl Solution of the KDV equation using Radial basis Functions, Appl. Math. Modelling. 32 (2008), 535-546.

  9. M. Dehghan, A. Shokri, A numerical method for KdV equation using collocation and radial basis functions, Nonlinear Dyn. 50 (2007), 111-120.

  10. M. Dehghan, A. Shokri, A numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions, Mathematics and Computers in Simulation. 79 (2008), 700-715.

  11. K. Djidjeli, W.G. Price, E.H. Twizell, Q. Cao, A linearized implicit pseudo-spectral method for some model equations: the regularized long wave equations, Commun. Numer. Meth. Engng. 19 (2003), 847-863.

  12. R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, H.C. Morris, Solitons and Nonlinear Wave Equations, Academic Press, New York, 1982.

  13. D.J. Evans, K.R. Raslan, Solitary waves for the generalized equal width (GEW) equation, Int. J. Comput. Math. 82 (4) (2005), 445-455.

  14. B. Fornberg, J. Zuev, The Runge phenomenon and spatially variable shape parameters in RBF interpolation, Comput. Math. Appl. 54 (2007), 379-398.

  15. L.R.T. Gardner, G.A. Gardner, Solitary waves of the equal width wave equation, Comput. Phys. 101 (1992), 218-223.

  16. L.R.T. Gardner, G.A. Gardner, F.A. Ayoub, N.K. Amein, Simulations of the EWE undular bore, Commun. Num. Methods Eng. 13 (1997), 583-592.

  17. L.R.T. Gardner, G. Gardner, T. Geyikli, The boundary forced MKdV equation, J. Comput. Phys. 11 (1994), 5-12.

  18. Y.C. Hon, R. Schaback, On unsymmetric collocation by radial basis functions, Appl. Math. Comput. 119 (2001), 177-186.

  19. E.J. Kansa, Multiquadrics, A scattered data approximation scheme with applications to computational fluid-dynamics, I- Surface approximations and partial derivative estimates, Comput. Math. Appl. 19 (1990), 127-145.

  20. E.J. Kansa, Multiquadrics, A scattered data approximation scheme with applications to computational fluid-dynamics, II-Solutions to parabolic, hyperbolic and elliptic partial differential equations, Comput. Math. Appl. 19 (1990), 147-161.

  21. A. Karageorghis, C.S. Chen and Y.S. Smyrlis, A matrix decomposition RBF algorithm: Approximation of functions and their derivatives, Appl. Numer. Math. 57 (2007), 304-319.

  22. D. Kaya, A numerical simulation of solitary-wave solutions of the generalized long-wave equation, Appl. Math. Comput. 88 (1997), 153-175.

  23. D. Kaya, S.M. El-Sayed, An application of the decomposition method for the generalized KdV and RLW equations, Chaos Solitons Fractals. 17 (2003), 869-877.

  24. J.C. Lewis, J.A. Tjon, Resonant production of solitons in the RLW equation, Phys. Lett. A. 73 (1979), 275-279.

  25. C.A. Micchelli, Interpolation of scattered data: distance matrix and conditionally positive functions, Construct. Approx. 2 (1986), 11-22.

  26. D.H. Peregrine, Calculations of the development of an undular bore, J. Fluid Mech. 25 (1966), 321-330.

  27. D.H. Peregrine, Long waves on a beach, J. Fluid Mech. 27 (1967), 815-827.

  28. J.I. Ramos, Explicit finite difference methods for the EW and RLW equations, Appl. Math. Comput. 179 (2006), 622-638.

  29. K.R. Raslan, A computational method for the equal width equation, Int. J. Comput. Math. 81(1) (2004), 63-72.

  30. K.R. Raslan, Collocation method using quartic B-spline for the equal width (EW) equation, Appl. Math. Comput. 168 (2005), 795-805.

  31. K.R. Raslan, Numerical methods for partial differential equations, Dissertation, Univ. of Al-Azhar, Cairo, Egypt, (1999).

  32. B. Saka, A finite element method for equal width equation, Appl. Math. Comput. 175 (2006), 730-747.

  33. B. Saka, Algorithms for numerical solution of the modified equal width wave equation using collocation method, Math. Comput. Modeling. 45 (2007), 1096-1117.

  34. Siraj-ul-Islam, Sirajul Haq, Arshed Ali, A meshfree method for the numerical solution of the RLW equation, J. Comput. Appl. Math. 223(2) (2009), 997-1012.

  35. S.I. Zaki, A least-squares finite element scheme for the EW equation, Comput. Methods Appl. Mech. Engrg. 189 (2000), 587-594.

  36. S.I. Zaki, Solitary waves induced by the boundary forced EW equation, Comput. Methods Appl. Mech. Engrg. 190 (2001), 4881-4887.

  37. S.I. Zaki, Solitary wave interactions for the modified equal width equation, Comput. Phys. 126 (2000), 219-231.