# A Fourth Order Finite Difference Method for the Good Boussinesq Equation

Authors: M. S. Ismail; Farida Mosally
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The “good” Boussinesq equation is transformed into a first order differential system. A fourth order finite difference scheme is derived for this system. The resulting scheme is analyzed for accuracy and stability. Newton’s method and linearization techniques are used to solve the resulting nonlinear system. The exact solution and the conserved quantity are used to assess the accuracy and the efficiency of the derived method. Head-on and overtaking interactions of two solitons are also considered. The numerical results reveal the good performance of the derived method.


---

## Body

## 1. Introduction

In recent years, remarkable developments have taken place in the study of nonlinear evolutionary partial differential equations. It is realized that many such equations possess special solutions in the form of pulses which retain their shapes and velocities after interacting with each other. Such solutions are called solitons. Many equations admitting soliton solutions are as follows: sine Gordon and double sine Gordon equations, Schrodinger equation, and KdV, MKdV, and complex modified KdV equations; many research works have been done on these equations [1–8]. Most of the current research is directed to solve coupled nonlinear systems analytically and numerically [9–21]. Solitons are of great interest in many physical areas, as, for example, in dislocation theory of crystals, plasma and fluid dynamics, magnetohydrodynamics, laser and fiber optics, and the study of the water waves.In this work we will study numerically the “good” Boussinesq (GB) nonlinear equation(1)∂2u∂t2-∂2u∂x2+∂4u∂x4-∂2(u2)∂x2=0;2222222222L0<x<L1,t>t0,
which provides a balance between dispersion and nonlinearity, that leads to the existence of soliton solutions, similar to the Korteweg-de Vries (KdV) and cubic nonlinear Schrödinger equation [1, 8, 22].The initial displacement associated with the partial differential equation given in (1) is assumed to take the form [23–25]
(2)u(x,t0)=f(x);L0≤x≤L1,∂u(x,t0)∂t=g(x);L0≤x≤L1,
and the boundary conditions
(3)u(L0,t)=0,u(L1,t)=0;t>t0,∂u(L0,t)∂x=0,∂u(L1,t)∂x=0;t>t0.The solitary wave solution of (1) is given by [3, 14]
(4)u(x,t)=-Asech2[A6(x-ct+x0)+b+12],c=±-2(b+A3),
where A is amplitude of the pulse and b is an arbitrary parameter. In addition, the double soliton solution [3, 14, 24, 26] can be written as
(5)u(x,t)=∑i=12{-Aisech2[Ai6(x-cit+xi0)+bi+12]},(6)ci=±-2(bi+Ai3),i=1,2,
where ci is the velocity of the ith soliton, Ai is the corresponding amplitude, bi is an arbitrary parameter, and xi0 is the initial position. The exact solution (5) represents two solitary waves located initially at the positions x=x10 and x=x20 moving to the right or left according to their velocities.Many research works on Boussinesq equation have been developed. Analytical solution of this equation was studied by many authors, such as the construction ofN-soliton solutions using the bilinear form [4], multiple soliton solutions for the GB equation using a simplified version of Hirota method [27] and Adomian decomposition method [28]. Construction of soliton solutions and periodic solution of Boussinesq equation by modified decomposition method are given in [29, 30]. A variational iteration method is developed for GB equation [31]. A solitary wave solution of the Boussinesq equation with power law nonlinearity is derived in [32].Many numerical methods have been developed for solving the Boussinesq equation, such as Petrov-Galerkin method [19]. Bratsos et al. [23, 24, 33] have developed finite difference schemes and have considered the t Bad Boussinesq (BB) and GB equations. El-Zhoheiry has developed an implicit finite difference scheme [34]. Aydin and Karasözen [35] constructed second order symplectic and multisymplectic integrators for the GB equation using the two-stage Lobatto IIIA-IIIB partitioned Runge-Kutta method. Daripa and Hua [36] developed finite difference method for BB equation of second order in time and space, where they have used the filtration and regularization techniques to control the growth of the errors arising from the instabilities and provide better approximate solutions of this equation. Ismail and Bratsos [25] have derived a predictor corrector scheme for the GB and BB equations; the scheme is fourth order in time and second order in space, and it is conditionally stable. Matsuo [37] also has derived conservative finite difference schemes for certain classes of nonlinear wave equations with some applications for the nonlinear Klein-Gordon and GB equations. Mohebbi and Asgari [26] also have solved the GB equation using a fourth order time stepping schemes with combination of discrete Fourier transform. Split step Fourier method is also used to solve Boussinesq-type equations. Dehghan and Salehi [38] have derived a mesh free method for the classical Boussinesq equation.The paper is laid out as follows. In Section2, the GB equation is transformed into a first order differential system in time; a finite difference scheme is derived for this system. The method is analyzed for accuracy and stability. In Section 3, a linearization technique is used to convert the nonlinear system obtained into a linear block tridiagonal system. Numerical tests are presented in Section 4. Concluding remarks are given in Section 5.


## 2. The Numerical Method

In order to derive a highly accurate method, we transform the GB equation (1) into the first order differential system in time as [19]
(7)∂u∂t=∂2w∂x2,L0<x<L1,(8)∂w∂t=-∂2u∂x2+u+u2,
with initial conditions
(9)u(x,0)=f(x),w(x,0)=∬L0xut(y,0)dy=g(x)on⁡L0<x<L1
and boundary conditions
(10)∂u∂x=0,∂w∂x=0atx=L0,L1.By using the boundary conditions (10), the system (7) and (8) has the conserved quantity [24–26, 38]
(11)I=∫L0L1u(x,t)dx=constant.In order to derive a numerical method for solving the system (7) and (8), the region R=[t>0]×[L0<x<L1] with its boundary ∂R consisting of the lines x=L0, x=L1 and the axis t=0 is covered with a rectangular mesh of points with coordinates [25]
(12)x=xm=xL+mh,m=0,1,…,N,t=tn=nk,n=0,1,2,…,
where h and k are the space and time increments, respectively. We denote the theoretical solution of (7) and (8) at the typical mesh point (xm,tn) by umn and wmn, and the numerical solution by Umn and Vmn. The fourth finite difference formula [39]
(13)uxx=1h2(1+112δx2)-1δx2umn+𝒪(h4),(14)δx2umn=um+1n-2umn+um-1n
is used to approximate the second order space derivative which appears in our system. By making use of (13) and the implicit midpoint rule, the finite difference nonlinear implicit scheme
(15)Umn+1-Umnk=1h2(1+112δx2)-1δx2[Wmn+1+Wmn2],(16)Wmn+1-Wmnk=-1h2(1+112δx2)-1δx2[Umn+1+Umn2]+[Umn+1+Umn2]2+[Umn+1+Umn2]
is obtained. The schemes (15) and (16) can be written as
(17)(1+112δx2)(Umn+1-Umn)=pδx2[Wmn+1+Wmn2],(18)(1+112δx2)(Wmn+1-Wmn)=-pδx2[Umn+1+Umn2]+k(1+112δx2)[Umn+1+Umn2]2=+k(1+112δx2)[Umn+1+Umn2],
where p=k/h2. On expansion of the central differences using (14), this will lead us to the difference scheme
(19)112[(Um-1n+1-Um-1n)+10(Umn+1-Umn)+(Um+1n+1-Um+1n)]=12p[(Wm-1n+1+Wm-1n)-2(Wmn+1+Wmn)=12p+(Wm+1n+1+Wm+1n)],(20)112[(Wm-1n+1-Wm-1n)+10(Wmn+1-Wmn)+(Wm+1n+1-Wm+1n)]=-12p[(Um-1n+1+Um-1n)-2(Umn+1+Umn)+(Um+1n+1+Um+1n)]=+148k[(Um-1n+1+Um-1n)2+10(Umn+1+Umn)2+(Um+1n+1+Um+1n)2]=+124k[(Um-1n+1+Um-1n)+10(Umn+1+Umn)+(Um+1n+1+Um+1n)],m=1,2,…,N.
By imposing the boundary conditions, a coupled nonlinear system is obtained, which can be solved by any iterative scheme. This system can be solved by Newton's method or fixed point method. Linearization method can also be used.


### 2.1. Accuracy of the Scheme

To study the accuracy of the proposed schemes (19) and (20), the numerical solutions Umn and Wmn are replaced by the exact solutions umn and wmn, respectively. By making use of the following expansions about the point (xm,tn):
(21)(1+112δx2)[umn+1-umnk]=[∂u∂t+k2∂2u∂t2+h212∂3u∂x2∂t+kh224∂4u∂x2∂t2+k36∂3u∂t3+⋯]mn,12h2δx2[umn+1+umn]=[∂2u∂x2+k2∂3u∂x2∂t+h212∂4u∂x4+kh224∂5u∂x4∂t+k24∂4u∂x2∂t2+⋯]mn,(1+112δx2)G[umn+1+umn2]=[G(u)+k2∂G∂t+h22∂2G∂x2+kh224∂3G∂x2∂t+⋯]mn
into (19) to obtain
(22)(1+112δx2)[umn+1-umnk]-12h2δx2[wmn+1+wmn]=(1+k2∂∂t)(∂u∂t-∂w∂x2)+h212∂2∂x2×((1+k2∂∂t)(∂u∂t-∂w∂x2))+𝒪(k2+h4).By making use of the differential system (7) and (8), all terms inside the brackets in (22) are zero. Similar analysis can be done for (20). This will lead us to the conclusion that the derived scheme is of second order in time and fourth order in space, that is, 𝒪(k2+h4). The numerical results confirm this.


### 2.2. Stability of the Scheme

To study the stability of the derived scheme, von Neumann stability analysis will be used. We consider the linear version of the proposed schemes (17) and (18) which can be given as
(23)(1+112δx2)(Umn+1-Umn)=p1δx2[Wmn+1+Wmn2],(24)(1+112δx2)(Wmn+1-Wmn)=-p1δx2[Umn+1+Umn2]+12αk(1+112δx2)[Umn+1+Umn],
where α is a constant quantity defined by
(25)α=max⁡1≤m≤N{1+Umn}.To apply von Neumann method, we assume that the solutions of (22) and (23) are
(26)Umn=λnexp⁡⁡(iβmh),Wmn=γnexp⁡⁡(iβmh),
where λ, γ, and β are constants. By substituting these solutions into (23) and (24), and after some manipulation, we get the system
(27)(1-13s)(Umn+1-Umn)=-2p1s[Wmn+1+Wmn],(1-13s)(Wmn+1-Wmn)=2p1s[Umn+1+Umn]+12αk(1-13s)[Umn+1+Umn],
where s=sin2(βh/2). Equations (27) can be written in a matrix vector form as
(28)Ψn+1=MΨn,
where
(29)Ψ=[UW]t,M=[ab-ca]-1[a-bca],a=1-13s,b=2p1s,c=2p1s+12αk(1-13s).The von Neumann necessary condition for the stability of this system is the maximum eigenvalue of value of the amplification matrixM in (28) is to be less than or equal to one. By direct calculation, the eigenvalues of the matrix M are
(30)λ1,2=a±ibca±ibc,
with modulus equal to one, and hence the proposed scheme is unconditionally stable.


## 2.1. Accuracy of the Scheme

To study the accuracy of the proposed schemes (19) and (20), the numerical solutions Umn and Wmn are replaced by the exact solutions umn and wmn, respectively. By making use of the following expansions about the point (xm,tn):
(21)(1+112δx2)[umn+1-umnk]=[∂u∂t+k2∂2u∂t2+h212∂3u∂x2∂t+kh224∂4u∂x2∂t2+k36∂3u∂t3+⋯]mn,12h2δx2[umn+1+umn]=[∂2u∂x2+k2∂3u∂x2∂t+h212∂4u∂x4+kh224∂5u∂x4∂t+k24∂4u∂x2∂t2+⋯]mn,(1+112δx2)G[umn+1+umn2]=[G(u)+k2∂G∂t+h22∂2G∂x2+kh224∂3G∂x2∂t+⋯]mn
into (19) to obtain
(22)(1+112δx2)[umn+1-umnk]-12h2δx2[wmn+1+wmn]=(1+k2∂∂t)(∂u∂t-∂w∂x2)+h212∂2∂x2×((1+k2∂∂t)(∂u∂t-∂w∂x2))+𝒪(k2+h4).By making use of the differential system (7) and (8), all terms inside the brackets in (22) are zero. Similar analysis can be done for (20). This will lead us to the conclusion that the derived scheme is of second order in time and fourth order in space, that is, 𝒪(k2+h4). The numerical results confirm this.


## 2.2. Stability of the Scheme

To study the stability of the derived scheme, von Neumann stability analysis will be used. We consider the linear version of the proposed schemes (17) and (18) which can be given as
(23)(1+112δx2)(Umn+1-Umn)=p1δx2[Wmn+1+Wmn2],(24)(1+112δx2)(Wmn+1-Wmn)=-p1δx2[Umn+1+Umn2]+12αk(1+112δx2)[Umn+1+Umn],
where α is a constant quantity defined by
(25)α=max⁡1≤m≤N{1+Umn}.To apply von Neumann method, we assume that the solutions of (22) and (23) are
(26)Umn=λnexp⁡⁡(iβmh),Wmn=γnexp⁡⁡(iβmh),
where λ, γ, and β are constants. By substituting these solutions into (23) and (24), and after some manipulation, we get the system
(27)(1-13s)(Umn+1-Umn)=-2p1s[Wmn+1+Wmn],(1-13s)(Wmn+1-Wmn)=2p1s[Umn+1+Umn]+12αk(1-13s)[Umn+1+Umn],
where s=sin2(βh/2). Equations (27) can be written in a matrix vector form as
(28)Ψn+1=MΨn,
where
(29)Ψ=[UW]t,M=[ab-ca]-1[a-bca],a=1-13s,b=2p1s,c=2p1s+12αk(1-13s).The von Neumann necessary condition for the stability of this system is the maximum eigenvalue of value of the amplification matrixM in (28) is to be less than or equal to one. By direct calculation, the eigenvalues of the matrix M are
(30)λ1,2=a±ibca±ibc,
with modulus equal to one, and hence the proposed scheme is unconditionally stable.


## 3. Linearization Technique

In order to avoid the solution of the nonlinear systems (19) and (20), a linearization technique will be developed to overcome this difficulty. Using Taylor's series expansion of the nonlinear term about umn, we obtain
(31)(umn+1+umn2)2=(umn)2+2umn(umn+1+umn2-umn)+𝒪(k2),
and this will lead us to the following approximation:
(32)(umn+1+umn2)2=umn+1umn+𝒪(k2),
which preserves second order temporal accuracy. By replacing the exact solution by the approximating one, we get
(33)(Umn+1+Umn2)2=Umn+1Umn.By substituting (33) into (18), we obtain the linearly implicit finite difference scheme(34)(1+112δx2)(Umn+1-Umn)=pδx2[Wmn+1+Wmn2],(1+112δx2)(Wmn+1-Wmn)=-pδx2[Umn+1+Umn2]+k(1+112δx2)Umn+1Umn=+k(1+112δx2)[Umn+1+Umn2],m=1,2,…,N.By expanding the central difference operator in (34), we obtain the linearized scheme
(35)112(Um-1n+1+10Umn+1+Um+1n+1)-12p(Wm-1n+1-2Wmn+1+Wm+1n+1)=112(Um-1n+10Umn+Um+1n)=+12p(Wm-1n-2Wmn+Wm+1n),112(Wm-1n+1+10Wmn+1+Wm+1n+1)+12p(Um-1n+1-2Umn+1+Um+1n+1)-112k(Um-1nUm-1n+1+10UmnUmn+1+Um+1nUm+1n+1)-124k(Um-1n+1+10Umn+1+Um+1n+1)=112(Wm-1n+10Wmn+Wm+1n)=-12p(Um-1n-2Umn+Um+1n)=+124k(Um-1n+10Umn+Um+1n),m=1,2,…,N.By utilizing the boundary conditions, the finite difference scheme (35) produces a linear block tridiagonal system in the unknowns {Umn+1,Wmn+1}, which can be solved directly by Crout's method. The obtained difference scheme is still of second order in time and fourth order method in space, and it is unconditionally stable.


## 4. Numerical Results

In all experiments we use the following values:L0=-100, L1=100, h=0.1, k=0.01, and b=-1/2. We study the accuracy of the proposed method by calculating the infinity error norm
(36)L∞=max⁡1≤m≤N|Umn-umn|.
The numerical solution and trapezoidal rule are used to calculate the error and the conserved quantity [19, 25, 26, 37]. All numerical results in this section are obtained from the solution of the nonlinear schemes (19) and (20) using Newton's method. Linearization method is used for comparison purpose only.


### 4.1. Single Soliton

To test the derived method, we consider the initial condition(37)u(x,0)=-Asech2[A6(x+x0)],c=±(1-2A3)1/2,w(x,0)=6Actanh[A6(x+x0)],
where (37) represent soliton and kink solution, respectively. The parameters A=0.369, x0=0.0, and c=0.86833 are used in this test. The numerical results for the nonlinear and linearized schemes are given in Tables 1 and 2, respectively. The numerical results (amplitude and the conserved quantity) obtained display the high accuracy of the proposed method. The execution time required for the nonlinear scheme is 24.89 sec. compared to 11.34 sec. for the linearized scheme. The numerical results produced in Table 1 are more accurate than the one in Table 2, and this is due to the approximation of the nonlinear term in the linearization process. A comparison of some existing methods is given in Table 3, which indicates that our method is the most accurate one. In Figures 1 and 2, we display the numerical solutions Umn and Wmn moving to the right for t=0,1,…,50.Table 1
Nonlinear scheme of single solitonA=0.369,h=0.1,k=0.01.
Time
A
I
L
∞
L
2
0.0
0.369000
−2.975902
0.0
0.0
10.0
0.368994
−2.975901
1
E
-
06
8
E
-
06
20.0
0.368975
−2.975901
2
E
-
06
1.3
E
-
05
30.0
0.368944
−2.975901
2
E
-
06
1.7
E
-
05
40.0
0.368976
−2.975903
3
E
-
06
2.0
E
-
05
50.0
0.368995
−2.975902
3
E
-
06
2.3
E
-
05Table 2
Linearization method of single soliton withA=0.369,h=0.1,k=0.01.
Time
A
I
L
∞
L
2
0.0
0.369000
−2.975902
0.0
0.0
10.0
0.369000
−2.975902
8.0
E
-
06
5.2
E
-
05
20.0
0.368987
−2.975989
1.9
E
-
05
1.11
E
-
04
30.0
0.368967
−2.976034
2.9
E
-
05
1.83
E
-
04
40.0
0.369008
−2.976080
4.4
E
-
05
2.71
E
-
04
50.0
0.369034
−2.976125
6.0
E
-
05
3.76
E
-
04Table 3
Comparison of numerical solutions withh=0.1,k=0.001 at T=72.
Method
L
∞
Present method (nonlinear)
0.970
E
-
07
Present method (linearization)
0.975
E
-
03
Bratsos [23]
0.378
E
-
03
Dehghan and Salehi [38]
0.488
E
-
03Figure 1
Single soliton: numerical solution ofUmn.Figure 2
Single kink: numerical solution ofWmn.To test whether the proposed numerical scheme exhibits the expected convergence rate in space, we perform some numerical experiments for various values ofh and fixed value of k. In these experiments we choose k=0.0001 to ensure that the temporal error is negligible. The rate of convergence for the scheme is calculated using the formula
(38)Rate  of  convergence=ln⁡⁡[L∞(h)/L∞(h/2)]ln⁡⁡2,
where L∞(h) is the infinity error norm. The errors and the rate of convergence for T=5 are given in Table 4.Table 4
Rate of convergence withk=0.0001,T=5.
h
L
∞
Order
0.8
0.000051587
—
0.4
0.000003161
4.03
0.2
0.000000197
4.00
0.1
0.000000013
3.92


### 4.2. Head-On Soliton Interaction

To study the head-on collision of two solitons, we choose the initial conditions(39)u(x,0)=u1(x,0)+u2(x,0),w(x,0)=w1(x,0)+w2(x,0),
where
(40)ui(x,0)=-Aisech2[Ai6(x+xi0)],ci=±(1-2Ai3)1/2,wi(x,0)=6Aicitanh[Ai6(x+xi0)],i=1,2.In [19] it was reported that solution blew up numerically when A>0.3691. By choosing the parameters x10=-x20=50.0, A1=A2=0.369, and c1=-c2=0.86833. The initial conditions represent two solitons and two kinks located at x10=-50 and x20=50. In Figures 3 and 4, we display the interaction scenario for t=0 to t=120; we can easily observe that the two solitons and the two kinks have been separated completely after the interaction and restored their original shapes and velocities. The velocity and the conserved quantity are given in Table 5. It is noted that when the two waves overlap (Figures 3 and 4), the joint amplitude is greater than the sum of individual amplitudes; this is in full agreement with [37]. By choosing A1=A2=0.4, the situation dramatically changes; both Figures 5 and 6 end around t~60. Newton's method fails to find the solution there. This agrees with the blowup results in [19, 37].Table 5
Interaction of two solitonsA1=A2=0.369,c1=-c2=0.86833.
T
A
I
0.0
0.36900
−5.951801
20
0.368975
−5.951801
40
0.368976
−5.951802
60
1.223350
−5.591803
80
0.369014
−5.951803
100
0.368973
−5.951809
120
0.368978
−5.951802Figure 3
Head-on interaction: numerical solution ofUmn with A1=A2=0.369.Figure 4
Head-on interaction: numerical solution ofWmn with A1=A2=0.369.Figure 5
Blowup: numerical solution ofUmn with A1=A2=0.4.Figure 6
Blowup: numerical solution ofWmn with A1=A2=0.4.


### 4.3. Overtaking Soliton Interaction

In this test we will study the interaction of two solitons moving in the same direction. In this case we choose the initial conditions(41)u(x,0)=u1(x,0)+u2(x,0),w(x,0)=w1(x,0)+w2(x,0),
where
(42)ui(x,0)=-Aisech2[A6(x-xi0)],ci=±(1-2bA3)1/2,(43)wi(x,0)=6Aicitanh[A6(x-xi0)],i=1,2.The following parameters are selected:(44)L0=-L1=-100,h=0.1,k=0.01,A1=0.3,A2=1.0,x10=-80,x20=-50.The initial conditions represent two solitons and two kinks are initially located atx10=-80 and x20=-50, moving in the same direction with velocities c1=0.8944 and c2=0.5774. The usual interaction has taken place and the faster wave interacted and separated from the slower one and left it behind; the interaction scenario is given in Figures 7 and 8. In Table 6, we track the amplitude and the conserved quantity during the interaction scenario. In Figure 9, we display the contours of the numerical solution Umn. Many numerical tests have been conducted for different amplitudes. We have noticed that the interaction occurred for Ai≤1.15, i=1,2, which is completely different from the head-on interaction case, where the solution blew up numerically when Ai>0.3691, i=1,2. This type of interaction for the GB equation seems to be reported in the literature for the first time as far as the authors know. We have noticed that when the two waves overlap (Figures 7 and 8), the joint amplitude is less than the sum of individual amplitudes.Table 6
Overtaking interaction withA1=0.3,A2=0.3.
T
A
v
I
0.0
1.0
0.5774
−7.581908
40.0
0.997507
0.578788
−7.581930
80.0
0.652932
0.751473
−7.581911
120.0
0.997311
0.578901
−7.581908
160.0
0.999744
0.577498
−7.581910
180
0.999982
0.577361
−7.581907Figure 7
Overtaking soliton interaction: numerical solution ofUmn with A1=0.3, A2=1.0.Figure 8
Overtaking soliton interaction: numerical solution ofWmn with A1=0.3, A2=1.0.Figure 9
Contours of overtaking soliton interaction, whereA1=0.3, A2=1.0.


### 4.4. Birth of Solitons

To study the evolution of arbitrary pulse [19], we choose the initial conditions
(45)u(x,0)=-Asech2[A6(x+x0)],w(x,0)=0.This test is carried out for values ofA, mainly in the vicinity of A=1.5, the theoretical value of the amplitude for c=0. We considered two cases A=1.2 and A=1.499999. In the first case we choose A=1.2; we have noticed that as time evolves, the initial pulse split into two solitons moving in the opposite directions with equal amplitudes A1=A2=0.354663 and velocities c1=-c2=0.873818; see Table 7. In Figures 10 and 11, we display the splitting scenario and some dispersive oscillations which occurred after splitting. In the second case we choose A=1.499999; as time evolves, the initial pulse is split into two solitons moving in the opposite direction with amplitudes A1=A2=0.37500; see Table 8. Figures 12 and 13 display the splitting scenario; no dispersive oscillation is observed in this case (linearized scheme failed to produce the numerical solution in this case). Many numerical tests have been conducted for different amplitudes; we have noticed that the numerical solution blew up if A>1.499999.Table 7
Birth of solitons withA=1.2,h=0.1,k=0.01.
T
A
v
I
10
0.367859
0.868770
−5.366562
20
0.360967
0.871410
−5.366563
30
0.357963
0.872558
−5.366565
40
0.356061
0.873284
−5.366560
50
0.354663
0.873818
−5.366560Table 8
Birth of solitons withA=1.499999,  h=0.1,k=0.01.
T
A
v
I
20
1.498385
0.032808
−5.999997
40
0.384944
0.862189
−5.999997
60
0.374967
0.866038
−5.999997
80
0.374994
0.866028
−5.999996
100
0.375001
0.866025
−5.999987Figure 10
Birth of solitons, whereA=1.2.Figure 11
Birth of solitons, whereA=1.2.Figure 12
Birth of solitons, whereA=1.499999.Figure 13
Birth of solitons, whereA=1.499999.


## 4.1. Single Soliton

To test the derived method, we consider the initial condition(37)u(x,0)=-Asech2[A6(x+x0)],c=±(1-2A3)1/2,w(x,0)=6Actanh[A6(x+x0)],
where (37) represent soliton and kink solution, respectively. The parameters A=0.369, x0=0.0, and c=0.86833 are used in this test. The numerical results for the nonlinear and linearized schemes are given in Tables 1 and 2, respectively. The numerical results (amplitude and the conserved quantity) obtained display the high accuracy of the proposed method. The execution time required for the nonlinear scheme is 24.89 sec. compared to 11.34 sec. for the linearized scheme. The numerical results produced in Table 1 are more accurate than the one in Table 2, and this is due to the approximation of the nonlinear term in the linearization process. A comparison of some existing methods is given in Table 3, which indicates that our method is the most accurate one. In Figures 1 and 2, we display the numerical solutions Umn and Wmn moving to the right for t=0,1,…,50.Table 1
Nonlinear scheme of single solitonA=0.369,h=0.1,k=0.01.
Time
A
I
L
∞
L
2
0.0
0.369000
−2.975902
0.0
0.0
10.0
0.368994
−2.975901
1
E
-
06
8
E
-
06
20.0
0.368975
−2.975901
2
E
-
06
1.3
E
-
05
30.0
0.368944
−2.975901
2
E
-
06
1.7
E
-
05
40.0
0.368976
−2.975903
3
E
-
06
2.0
E
-
05
50.0
0.368995
−2.975902
3
E
-
06
2.3
E
-
05Table 2
Linearization method of single soliton withA=0.369,h=0.1,k=0.01.
Time
A
I
L
∞
L
2
0.0
0.369000
−2.975902
0.0
0.0
10.0
0.369000
−2.975902
8.0
E
-
06
5.2
E
-
05
20.0
0.368987
−2.975989
1.9
E
-
05
1.11
E
-
04
30.0
0.368967
−2.976034
2.9
E
-
05
1.83
E
-
04
40.0
0.369008
−2.976080
4.4
E
-
05
2.71
E
-
04
50.0
0.369034
−2.976125
6.0
E
-
05
3.76
E
-
04Table 3
Comparison of numerical solutions withh=0.1,k=0.001 at T=72.
Method
L
∞
Present method (nonlinear)
0.970
E
-
07
Present method (linearization)
0.975
E
-
03
Bratsos [23]
0.378
E
-
03
Dehghan and Salehi [38]
0.488
E
-
03Figure 1
Single soliton: numerical solution ofUmn.Figure 2
Single kink: numerical solution ofWmn.To test whether the proposed numerical scheme exhibits the expected convergence rate in space, we perform some numerical experiments for various values ofh and fixed value of k. In these experiments we choose k=0.0001 to ensure that the temporal error is negligible. The rate of convergence for the scheme is calculated using the formula
(38)Rate  of  convergence=ln⁡⁡[L∞(h)/L∞(h/2)]ln⁡⁡2,
where L∞(h) is the infinity error norm. The errors and the rate of convergence for T=5 are given in Table 4.Table 4
Rate of convergence withk=0.0001,T=5.
h
L
∞
Order
0.8
0.000051587
—
0.4
0.000003161
4.03
0.2
0.000000197
4.00
0.1
0.000000013
3.92


## 4.2. Head-On Soliton Interaction

To study the head-on collision of two solitons, we choose the initial conditions(39)u(x,0)=u1(x,0)+u2(x,0),w(x,0)=w1(x,0)+w2(x,0),
where
(40)ui(x,0)=-Aisech2[Ai6(x+xi0)],ci=±(1-2Ai3)1/2,wi(x,0)=6Aicitanh[Ai6(x+xi0)],i=1,2.In [19] it was reported that solution blew up numerically when A>0.3691. By choosing the parameters x10=-x20=50.0, A1=A2=0.369, and c1=-c2=0.86833. The initial conditions represent two solitons and two kinks located at x10=-50 and x20=50. In Figures 3 and 4, we display the interaction scenario for t=0 to t=120; we can easily observe that the two solitons and the two kinks have been separated completely after the interaction and restored their original shapes and velocities. The velocity and the conserved quantity are given in Table 5. It is noted that when the two waves overlap (Figures 3 and 4), the joint amplitude is greater than the sum of individual amplitudes; this is in full agreement with [37]. By choosing A1=A2=0.4, the situation dramatically changes; both Figures 5 and 6 end around t~60. Newton's method fails to find the solution there. This agrees with the blowup results in [19, 37].Table 5
Interaction of two solitonsA1=A2=0.369,c1=-c2=0.86833.
T
A
I
0.0
0.36900
−5.951801
20
0.368975
−5.951801
40
0.368976
−5.951802
60
1.223350
−5.591803
80
0.369014
−5.951803
100
0.368973
−5.951809
120
0.368978
−5.951802Figure 3
Head-on interaction: numerical solution ofUmn with A1=A2=0.369.Figure 4
Head-on interaction: numerical solution ofWmn with A1=A2=0.369.Figure 5
Blowup: numerical solution ofUmn with A1=A2=0.4.Figure 6
Blowup: numerical solution ofWmn with A1=A2=0.4.


## 4.3. Overtaking Soliton Interaction

In this test we will study the interaction of two solitons moving in the same direction. In this case we choose the initial conditions(41)u(x,0)=u1(x,0)+u2(x,0),w(x,0)=w1(x,0)+w2(x,0),
where
(42)ui(x,0)=-Aisech2[A6(x-xi0)],ci=±(1-2bA3)1/2,(43)wi(x,0)=6Aicitanh[A6(x-xi0)],i=1,2.The following parameters are selected:(44)L0=-L1=-100,h=0.1,k=0.01,A1=0.3,A2=1.0,x10=-80,x20=-50.The initial conditions represent two solitons and two kinks are initially located atx10=-80 and x20=-50, moving in the same direction with velocities c1=0.8944 and c2=0.5774. The usual interaction has taken place and the faster wave interacted and separated from the slower one and left it behind; the interaction scenario is given in Figures 7 and 8. In Table 6, we track the amplitude and the conserved quantity during the interaction scenario. In Figure 9, we display the contours of the numerical solution Umn. Many numerical tests have been conducted for different amplitudes. We have noticed that the interaction occurred for Ai≤1.15, i=1,2, which is completely different from the head-on interaction case, where the solution blew up numerically when Ai>0.3691, i=1,2. This type of interaction for the GB equation seems to be reported in the literature for the first time as far as the authors know. We have noticed that when the two waves overlap (Figures 7 and 8), the joint amplitude is less than the sum of individual amplitudes.Table 6
Overtaking interaction withA1=0.3,A2=0.3.
T
A
v
I
0.0
1.0
0.5774
−7.581908
40.0
0.997507
0.578788
−7.581930
80.0
0.652932
0.751473
−7.581911
120.0
0.997311
0.578901
−7.581908
160.0
0.999744
0.577498
−7.581910
180
0.999982
0.577361
−7.581907Figure 7
Overtaking soliton interaction: numerical solution ofUmn with A1=0.3, A2=1.0.Figure 8
Overtaking soliton interaction: numerical solution ofWmn with A1=0.3, A2=1.0.Figure 9
Contours of overtaking soliton interaction, whereA1=0.3, A2=1.0.


## 4.4. Birth of Solitons

To study the evolution of arbitrary pulse [19], we choose the initial conditions
(45)u(x,0)=-Asech2[A6(x+x0)],w(x,0)=0.This test is carried out for values ofA, mainly in the vicinity of A=1.5, the theoretical value of the amplitude for c=0. We considered two cases A=1.2 and A=1.499999. In the first case we choose A=1.2; we have noticed that as time evolves, the initial pulse split into two solitons moving in the opposite directions with equal amplitudes A1=A2=0.354663 and velocities c1=-c2=0.873818; see Table 7. In Figures 10 and 11, we display the splitting scenario and some dispersive oscillations which occurred after splitting. In the second case we choose A=1.499999; as time evolves, the initial pulse is split into two solitons moving in the opposite direction with amplitudes A1=A2=0.37500; see Table 8. Figures 12 and 13 display the splitting scenario; no dispersive oscillation is observed in this case (linearized scheme failed to produce the numerical solution in this case). Many numerical tests have been conducted for different amplitudes; we have noticed that the numerical solution blew up if A>1.499999.Table 7
Birth of solitons withA=1.2,h=0.1,k=0.01.
T
A
v
I
10
0.367859
0.868770
−5.366562
20
0.360967
0.871410
−5.366563
30
0.357963
0.872558
−5.366565
40
0.356061
0.873284
−5.366560
50
0.354663
0.873818
−5.366560Table 8
Birth of solitons withA=1.499999,  h=0.1,k=0.01.
T
A
v
I
20
1.498385
0.032808
−5.999997
40
0.384944
0.862189
−5.999997
60
0.374967
0.866038
−5.999997
80
0.374994
0.866028
−5.999996
100
0.375001
0.866025
−5.999987Figure 10
Birth of solitons, whereA=1.2.Figure 11
Birth of solitons, whereA=1.2.Figure 12
Birth of solitons, whereA=1.499999.Figure 13
Birth of solitons, whereA=1.499999.


## 5. Concluding Remarks

In this paper, a highly accurate finite difference method to solve GB equation is proposed. The main idea is to transform the GB equation into a first order differential system in time. The method is fourth order in space and second order in time, and it is unconditionally stable. The method produced a coupled nonlinear system; Newton's and linearization methods are used to solve it. The nonlinear scheme produced more accurate results than the linearized scheme. Numerical results for single soliton, head-on, and overtaking solitons interactions are given; the most interesting phenomena we observed were when two solitons collide; if both solitons are small enough, they pass through each other like the other usual solitons do, but when they exceed some limit, the solution blows up at the collision, even if both amplitudes are smaller than3/2 for being stable solitons [14, 15]. The overtaking solitons interaction is also discussed as a new numerical test. Elastic interaction occurred for large amplitudes (Ai≤1.15, i=1,2), and this is contrary to the head-on interaction, where blowup is observed for moderate values of the amplitudes (Ai>0.3691, i=1,2). Birth of solitons is also observed for initial pulse with amplitude A≤1.499999. A blowup was observed when A=1.5 and this is in full agreement with the existing methods [14, 15].



---

*Source: 323260-2014-02-26.xml*