- Email: [email protected]
Galerkin methods for the Davey–Stewartson equations
Applied Mathematics and Computation 328 (2018) 144–161
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Galerkin methods for the Davey–Stewartson equationsR Yali Gao a, Liquan Mei a,∗, Rui Li b a b
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China School of Mathematics and Information Science, Shaanxi Normal University, Xi’an, 710062, China
a r t i c l e
i n f o
Keywords: Davey–Stewartson equations Galerkin finite element method Extrapolated Crank–Nicolson method Explicit multistep method
a b s t r a c t In this paper, we propose two Galerkin methods to investigate the evolution of the Davey– Stewartson equations. The extrapolated Crank–Nicolson scheme and decoupled semiimplicit multistep scheme are employed to increase the order of the time discrete accuracy, which only requires the solutions of a linear system at each time step. Four numerical experiments are presented to illustrate the features of the proposed numerical methods, such as the optimal convergence order, the conservation variable and the application in rogue waves. © 2018 Elsevier Inc. All rights reserved.
1. Introduction A variety of phenomena across a range of disciplines, which in applied mathematics and physics from fluid dynamics, solid state physics, quantum machines, plasmas physics to nonlinear optics and so forth, and in chemistry and biology as well [1,2], are described by nonlinear dispersive partial differential equations. For instance, the Schrödinger equation is the fundamental governing equation in quantum mechanics and quantum field theory [3], which is used to describe many-body theory and condensed matter physics like the Bose–Einstein condensate. It also has extensive applications in nonlinear optics [4] and water waves [5]. In view of the importance of these equations, many researchers tried their best to find analytical solutions or numerical solutions to the nonlinear dispersive equations by using various methods [6–14]. The Davey–Stewartson (DS) equations were originally introduced by Davey and Stewartson [15] to describe the evolution of a three-dimensional weakly nonlinear wave-packet on water in finite depth, in which the water waves have one direction of travel but the amplitude of waves is modulated in two spatial directions. The equations appear in studying the longwave-short-wave resonances [16] and describing gravity-capillarity surface wave packets in the limit of shallow water [17]. The evolution of a plasma under the action of a magnetic field in plasma physics [18] and certain (2+1)-dimensional nonlinear localized excitations can be described by DS equations. Furthermore, the (2+1)-dimensional Schrödinger equation with cubic nonlinearity fails to be integrable. However, Davey–Stewartson equations, a pair of coupled nonlinear equations in two dependent variables, were shown to be integrable [19–21] as the extension of the (1+1)-dimensional nonlinear Schrödinger equation. In fact, when the DS system is considered with various parameters in it, there are only two integrable cases which are known as DS I and DS II [21]. These two integrable cases are of elliptic-hyperbolic and hyperbolic-elliptic type, depending on the strength of surface tension. The DS II equations are thus the hyperbolic-elliptic ones, which are in particular relevant for surface gravity waves without surface tension. In addition, for DS I and DS II equations, special localized solutions, called respectively dromions and lumps, can be explicitly found. R ∗
The project is supported by China Scholarship Council (201506280175) and NSF of China (11371289, 11501441). Corresponding author. E-mail addresses: [email protected] (Y. Gao), [email protected] (L. Mei), [email protected] (R. Li).
https://doi.org/10.1016/j.amc.2018.01.044 0 096-30 03/© 2018 Elsevier Inc. All rights reserved.
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
145
Many works for Davey–Stewartson equations have been done in recent years. A vast method finding analytical solutions has been proposed [6,8,21–28]. Biswas [22] considered G’/G method and the solitary wave ansatz method to study traveling wave solution and 1-soliton solutions, respectively. El-Kalaawy and Ibrahim [23] obtained the solitary wave solutions by using extended mapping method technique. Ohta and Yang [27,28] derived the general rogue waves by the bilinear method. In parallel with the analytical studies, a surge of efforts have been devoted to the numerics of this system [11,13,29–34]. Stimming and Zhang [29] presented an efficient and accurate solver for the nonlocal potential in the Davey–Stewartson equations using nonuniform Fast Fourier Transform. Besse et al. [30] proposed a time splitting spectral method for numerical study of the DS equations, including exact soliton type solutions of the hyperbolic-elliptic equations and the blow up of focusing elliptic-elliptic equations. Klein and Roidot [13] investigated the semi-classical Davey–Stewartson equations by both fourth order time-stepping and spectral method. Gao et al. [31] proposed a Galerkin method with first-order time discretization. To the best of our knowledge, there are few results on numerics of the DS equations by finite element method, though there are numerous investigations for nonlinear dispersive equations [31,35]. In this paper, we study the numerical solutions of Davey–Stewartson equations by Galerkin finite element method with extrapolated Crank–Nicolson Scheme and decoupled multistep method. Due to the nonlinearity and coupling and complex structure system of DS equations, linearized semi-implicit schemes should be considered approximations in time direction since their implementation only requires the solution of a linear system. Among linearized schemes, the extrapolated Crank– Nicolson method is a popular one since it provides a second order accuracy in the time direction. However, linearized Crank–Nicolson method only solved the nonlinearity of the DS equations, we also need solve a coupled algebraic system. To overcome this shortcoming of coupled Crank–Nicolson method, we develop the decoupled scheme by using semi-implicit multistep method, which improves the efficiency of algorithm and remains the second order convergence rate in time as Crank–Nicolson time discretization. The outline of this paper is as follows. In the next section, we introduce the governing equations for DS equations. In Section 3, numerical algorithms of the DS equations are given. In Section 3.1, we study the weak form and the proof of L2 norm conservations. The full discretization based on linearized Crank–Nicolson method and decoupled multistep method are presented in Section 3.2 and Section 3.3, respectively. In Section 4, the L2 error estimates for extrapolated Crank–Nicolson and decoupled implicit-explicit two-step schemes are presented and the proof are given. In Section 5, some details during the progress of implementation are presented. In Section 6, we propose four numerical experiments for DS equations, which illustrates the efficiency of the presented numerical algorithms by the L2 , L∞ , H1 error accuracy and the relative error of L2 norm conservation laws. Furthermore, we apply the proposed method to the rogue waves in DS equations. Finally, a brief conclusion is given in the last section.
2. Governing equations Let = [L1 , L2 ] × [L3 , L4 ] be a bounded subset of R2 and T > 0 be a final time, J = (0, T ], T = × (0, T ). We consider the following (2+1)-dimensional Davey–Stewartson equations with power law nonlinearity
⎧ ∂ ∂ 2 ∂ 2 ⎪ ⎪ i +α +β + γ | |2 + δ x = 0, ⎪ 2 ⎪ ∂t ∂x ∂ y2 ⎪ ⎪ ⎨ ∂ 2 ∂ 2 ∂| |2 + +θ = 0, 2 2 ∂x ∂x ∂y ⎪ ⎪ ⎪ ⎪ = 0, = 0, ⎪ ⎪ ⎩ = 0 , = 0 ,
in T , in T , on
(2.1)
∂ × [0, T ],
on × {t = 0},
where (x, y, t) is a complex valued function which denotes the complex amplitude,while (x, y, t) is a real valued function which denotes the velocity of a underlying mean flow. Using the terminology in [31], sgnα = 1 and sgnα = −1 denote elliptic-elliptic DS and hyperbolic-elliptic DS, respectively. For both there is a focusing (sgnγ = 1) and a defocusing (sgn √γ = −1) version [29]. β , δ , θ are real parameters. 0 : → C and 0 : → R are the given initial conditions, and i = −1. In absence of y dependence, the DS equations can be viewed as a generalization of the (1+1)-dimensional cubic nonlinear Schrödinger equation,
i
∂ ∂ 2 +α + γ | |2 = 0. ∂t ∂ x2
such that = x instead of Moreover, the Davey–Stewartson system has in fact another form (when using the variable in (2.1), we still denote it by for simplicity)
⎧ ∂ ∂ 2 ∂ 2 ⎪ +α +β + γ | |2 + δ = 0, ⎨i 2 ∂t ∂x ∂ y2 2 2 2 2 ⎪ ⎩ ∂ + ∂ + θ ∂ | | = 0. ∂ x2 ∂ y2 ∂ x2
(2.2)
146
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Furthermore, the DS Eq. (2.1) possess two formal conservation laws [30] corresponding to the L2 norm
M (t ) =
| |2 dx dy,
(2.3)
and the energy (Hamiltonian)
E (t ) =
1 α|x |2 + β|y |2 + γ | |4 − δθ (∂x 2 + ∂y 2 ) dx dy. 2
3. Numerical methods We use the standard Sobolev space Wm,k (), where m is a nonnegative integer with 1 ≤ k ≤ ∞, and the space involving time Lk (0, T; Wm,k ) with norm · Lk (0,T ;W m,k ) . Let W m,2 () = H m () with the norm · m and ( · , · ) be the inner product on L2 () = H 0 () defined by
(u, v ) =
u(x )(v(x ))∗ dx dy,
∀ u, v ∈ L2 (),
where v∗ denotes the conjugate of v. The norm · ∞ denotes the essential supremum. Set V = H01 () = {v ∈ H 1 () : v|∂ = 0}. Applying Greens formula to problem (2.1) and using the boundary condition in the definition of V, the weak formulation of problem (2.1) can be written as follows: Find ( , ): J → V × V, such that, for ∀ (v, s ) ∈ V × V,
⎧ ∂ ∂v ∂ ∂v ∂ ⎪ ⎪ i , v − α , − β , + γ (| |2 , v ) + δ ( x , v ) = 0, ⎪ ⎪ ∂x ∂x ∂y ∂y ⎨ ∂ t
∂s ∂ ∂s ∂ ∂s , + , + θ | |2 , = 0, ⎪ ⎪ ∂ x ∂ x ∂ y ∂ y ∂ x ⎪ ⎪ ⎩ ( , v ) = (0 , v ), ( , s ) = ( 0 , s ).
(3.1)
3.1. Semi-discretization in space Let us consider a uniform regular partition h of triangular element e with = maxe ∈ h diam(e ). Set
e
e . We denote mesh size by h =
¯ , v|e ∈ Pr (e ), ∀ e ∈ h } ⊂ H 1 (), Vh = {v : v ∈ V ∩ C 0 () 0 where Pr (e ) denotes the polynomial space with degree less than or equal to r defined on e . Then from Eq. (3.1), we obtain the following semi-discrete form for problem (2.1): Find ( h , h ): J → Vh × Vh , such that, for ∀ (v, s ) ∈ Vh × Vh ,
⎧ ∂ h ∂ h ∂v ∂ h ∂v ⎪ ⎪ i , v − α , − β , + γ (|h |2 h , v ) + δ ( hx h , v ) = 0, ⎪ ⎪ ∂x ∂x ∂y ∂y ⎨ ∂ t
∂s ∂ h ∂ s ∂ h ∂ s , + , + θ | h | 2 , = 0, ⎪ ⎪ ∂ x ∂ x ∂ y ∂ y ∂ x ⎪ ⎪ ⎩ 0 (h , v ) = (0 , v ), ( 0h , s ) = ( 0 , s ).
(3.2)
Proposition 3.1. (L2 norm conservation) The solutions to the scheme (3.2) satisfy the L2 norm conservation
d d M= | |2 d x d y = 0. dt dt h
(3.3)
Proof. First, the complex conjugate of first equation of Eq. (3.2) is
∗ ∗
∂ h ∂v ∂ h ∂v ∂ h∗ −i ,v − α , −β , + γ (|h |2 h∗ , v ) + δ ( hx h∗ , v ) = 0. ∂t ∂x ∂x ∂y ∂y
(3.4)
Choosing v = h in the first equation of Eq. (3.2) and v = h∗ in Eq. (3.4), we obtain by subtracting
∗
∂ h ∗ ∂ h , h + , h = 0 , ∂t ∂t
(3.5)
namely,
d | |2 d x d y = 0. dt h
Then, we obtain the L2 norm conservation (3.3) of Eq. (3.2).
(3.6)
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
147
We now consider full discrete scheme of Eq. (3.2). For time discretization, we discuss the second order schemes: coupled Crank–Nicolson (CCN) method and decoupled BDF2 (DBDF2) method. Let 0 = t 0 < t 1 < · · · < t M = T be a uniform partition of J into subintervals J n = (t n , t n+1 ), n = 0, 1, . . . , M − 1, with length t = t n+1 − t n = MT . Then tn = n t, n = 0, 1, . . . , M denote the discrete time levels. (Un , Zn ) ∈ Vh × Vh define the fully-discrete approximations of the exact solutions ( (tn ), (tn )) ∈ V × V of DS problem (2.1). 3.2. Full discretization with extrapolated Crank–Nicolson method The fully discrete Galerkin approximation of Eq. (3.2) is given as follows: Find (hn+1 , nh+1 ) ∈ Vh × Vh , such that, for all 0 ≤ n ≤ M − 1, ∀ (v, s ) ∈ Vh × Vh ,
n+1 n+1
⎧ n+1 h − hn ∂ ¯ h ∂ ¯ h ∂v ∂v ⎪ ¯ n+1 |2 ¯ n+1 , v ) + δ ( ¯ n+1 ¯ n+1 , v ) = 0, ⎪ i , v − α , − β , + γ ( | ⎨ h h hx h t ∂x ∂x ∂y ∂y n+1
n+1
¯ ¯ ∂ ∂ ⎪ ∂s ∂s ⎪ h h ¯ n+1 |2 , ∂ s = 0, ⎩ , + , + θ | h ∂x ∂x ∂y ∂y ∂x
(3.7)
where
¯ hn+1 =
hn+1 + hn 2
,
¯ n+1 = h
nh+1 + nh 2
.
However, we should solve a nonlinear system at each time level and this progress decreases the accuracy of time. We can choose an approximation of second order accuracy to remain the same order of accuracy in the truncation error. This can be done by exploiting extrapolation technique
3 2
1 2
ˆ n+1 = n − n−1 , h h h
3 2
1 2
ˆ n+1 = n − n−1 , h h h
n ≥ 1.
Let denote hn+1 = 1nh+1 + i2nh+1 . We obtain the following full-discrete scheme with the extrapolated Crank–Nicolson method: Algorithm 3.1. (Coupled Crank–Nicolson (CCN) scheme) Find (hn+1 , nh+1 ) ∈ Vh × Vh , such that, for all 1 ≤ n ≤ M − 1, ∀ (v, s ) ∈ Vh × Vh ,
n+1 n+1 ⎧ n+1 h − hn ∂ ¯ h ∂ ¯ h ∂v ∂v ⎪ ˆ n+1 |2 ¯ n+1 , v ) + δ ( ¯ n+1 ˆ n+1 , v ) = 0, ⎪ ,v − α , −β , + γ ( | ⎨i h h hx h t ∂x ∂x ∂y ∂y n+1
n+1
¯ ¯ ∂ ∂ ⎪ ∂s ∂s ⎪ h h ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 , ∂ s = 0. ⎩ , + , +θ 1h 1h 2h 2h ∂x ∂x ∂y ∂y ∂x
(3.8)
Remark 3.1. The CCN scheme for DS problem is a linear scheme and satisfies the L2 norm conservation (2.3). Therefore, one only need to solve the coupled linear equations which can be solved very efficiently at each time level. However, we should solve a fully coupled system (3.8) at each time level which is time consuming in practical computation, especially in three dimensions. This motivates us to consider the following decoupled scheme for Eq. (3.2) by semi-implicit multistep scheme in time discretization. 3.3. Full discretization with semi-implicit multistep scheme In this paper, we consider 2-steps full discretization method for DS Eq. (2.1), which shall be referred to as two-step backward differentiation formula (BDF2). Thus, the decoupled numerical method can be written as follows: Algorithm 3.2. (Decoupled BDF2 (DBDF2) method) Find (hn+1 , nh+1 ) ∈ Vh × Vh , such that, for all 1 ≤ n ≤ M − 1, ∀ (v, s ) ∈ Vh × Vh ,
⎧ n+1 n+1 n+1
3h − 4hn + hn−1 ∂ h ∂ h ∂v ∂v ⎪ ⎪ ,v − α , −β , ⎪i ⎪ 2 t ∂x ∂x ∂y ∂y ⎪ ⎨ +γ (2|hn |2 hn − |hn−1 |2 hn−1 , v ) + δ ( nhx+1 (2hn − hn−1 ), v ) = 0, ⎪ n+1
n+1
⎪ ⎪ ∂ h ∂ h ∂s ∂s ∂s ⎪ ⎪ , + , + θ 2|hn |2 − |hn−1 |2 , = 0. ⎩ ∂x ∂x ∂y ∂y ∂x
(3.9)
Remark 3.2. The DBDF2 scheme for DS problem is totally decoupled, linear scheme and satisfies the L2 norm conservation (2.3). Therefore, the sequence linear algebraic system can be efficiently solved at each time level.
148
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
As both Eqs. (3.8) and (3.9) are used only for n ≥ 1, we must need start approximation (h1 , 1h ). A predictor-corrector
method is used for this purpose. We begin with an approximation (h1,0 , 1h,0 ), computed by the implicit-explicit Euler scheme,
1,0 1,0 ⎧ 1,0 h − h0 ∂ ¯ h ∂v ∂ ¯ h ∂v ⎪ ¯ 1,0 , v ) + δ ( ¯ 1,0 0 , v ) = 0 , ⎪ ,v − α , −β , + γ (|h0 |2 ⎨i h h hx t ∂x ∂x ∂y ∂y 1,0
1,0
¯ ¯ ∂ ∂ ⎪ ∂s ∂s ⎪ h h ¯ 1,0 + 0 ¯ 1,0 , ∂ s = 0 , ⎩ , + , − θ 10h 2h 2h ∂ x 1h ∂x ∂x ∂y ∂y
(3.10)
where
¯ h1,0 =
h1,0 + h0 2
¯ 1,0 = h
,
1h,0 + 0h 2
.
and correct it by
1 1
⎧ 1 h − h0 ∂ ¯ h ∂v ∂ ¯ h ∂v ⎪ ¯ 1,0 |2 ¯ 1 , v ) + δ ( ¯1 ¯ 1,0 , v ) = 0 , ⎪ ,v − α , −β , + γ ( | ⎨i h hx h t ∂x ∂x ∂y ∂y 1
1
¯ ¯ ∂ ∂ ⎪ ⎪ h ∂s h ∂s ¯ 1,0 ¯ 1 + ¯ 1,0 ¯ 1 , ∂ s = 0. ⎩ , + , −θ 1h 2h ∂ x 1h 2h ∂x ∂x ∂y ∂y
(3.11)
4. Error estimates In this section, we shall analyze the optimal error estimates for the full discrete scheme presented in Algorithm 3.1 with Crank–Nicolson method and Algorithm 3.2 with two-step backward differentiation formula. We introduce the projection : V → Vh , : V → Vh defined by, h h
∂χ ∂ α ( v − v ), ∂x h ∂x
∂χ ∂ +β ( v − v ), ∂y h ∂y
(∇ ( h v − v ), ∇χ ) = 0,
∀χ ∈ Vh , v ∈ V.
= 0,
(4.1)
∀χ ∈ Vh , v ∈ V.
(4.2)
r+1 v − vr+1 , ∀ v ∈ H r+1 () ∩ H01 (), h v + h (v − h v ) ≤ Ch
(4.3)
r+1 v − vr+1 , ∀ v ∈ H r+1 () ∩ H01 (). h v + h (v − h v ) ≤ Ch
(4.4)
Then, it may be shown that
Let denote n n hn − (tn ) = hn − h (tn ) + h (tn ) − (tn ) = e + E , n n nh − (tn ) = nh − h (tn ) + h (tn ) − (tn ) = e + E .
Theorem 4.1. Let {hn+1 , nh+1 }M satisfy the Eq. (3.8), the solutions ((t n+1 ), (t n+1 )) of Eq. (2.1) be sufficiently smooth and n=0
2 h0 − 0 2 + h1 − (t 1 )2 + t ( 0h − 0 2 + 1h − (t 1 )2 ) ≤ C hr+1 + t 2 ,
Then, for h and k sufficiently small, there exists a constant C, such that
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 .
0≤n≤M
(4.5)
(4.6)
where the constant C is independent of h and t. n , E n ) is bounded as claimed in (4.3), we only need to consider (en , en ). We consider the error for n ≥ 1. Proof. Since (E Using the projection property as presented in (4.1) and (4.2), and subtracting (3.8) from (3.1) results in the following error estimates in (en , en ), for n ≥ 1,
∂ n+1 ∂v ∂ n+1 ∂v −α −β + γ (J1n+1 , v ) + δ (J2n+1 , v ) e¯ , e¯ , ∂x ∂x ∂y ∂y
n+1 n E − E (t n+1 ) − (t n ) n+ 12 = −i − t (t )+ , v + (J3n+1 , v ), t t
en+1 − en i ,v t
(e¯n +1 , s ) + θ (J4n+1 , sx ) = −θ (J5n+1 , sx ) − (J6n+1 , s ),
(4.7) (4.8)
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
149
where
ˆ n+1 |2 ¯ n+1 − | ˆ n+1 |2 ¯ n+1 , J1n+1 = | h h 1
n+1 ¯ n+1 − (t n+ 2 ), J31 =
n+1 n+1 n+1 J3n+1 = α J31 − γ J32 − δ J33 ,
¯ nx +1 ¯ n+1 ˆ n+1 − ˆ n+1 , J2n+1 = hx h
n+1 ˆ n+1 |2 ¯ n+1 − |(t n+ 2 )|2 (t n+ 2 ), J32 = | 1
1
1
1
n+1 ¯ nx +1 ˆ n+1 − x (t n+ 2 )(t n+ 2 ), J33 =
ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 ) − ( ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 ), J4n+1 = ( 1 1 2 2 1h 1h 2h 2h 1
1
ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 ) − ( 2 (t n+ 2 ) + 2 (t n+ 2 )), J5n+1 = ( 1 2 1 1 2 2
1
¯ n+1 − (t n+ 2 ). J6n+1 =
Choosing s = e¯n +1 in (4.8), we obtain
e¯n +1 ≤ C (J4n+1 + J5n+1 + J6n+1 ). n+1
Choosing v = e¯
(4.9)
in (4.7), we get
i (en+1 2 − en 2 ) − αe¯n+1 2 + γ (J1n+1 , e¯n+1 ) + δ (J2n+1 , e¯n+1 ) 2 t
n E n+1 − E (t n+1 ) − (t n ) 1 = −i − t (t n+ 2 ) + , e¯n+1 + (J3n+1 , e¯n+1 ). t t
(4.10)
Using Taylor expansions, the first term on the right-hand of (4.10) satisfies
2 t n+1 (t n+1 ) − (t n ) n+ 12 3 − ( t ) ≤ C ( t ) t t t (s )2 ds. t t tn
(4.11)
Taking the imaginary part of Eq. (4.10), and using Eq. (4.3), we get
2 1 (en+1 2 − en 2 ) ≤C (J1n+1 2 + J2n+1 2 + J3n+1 2 + J4n+1 2 + J5n+1 2 + J6n+1 2 ) + C e¯n+1 2 +C hr+1 + t 2 . 2 t (4.12) Now, we estimate terms on the right-hand of Eq. (4.12). Exploiting inequality (4.3), we have
J1n+1 = (¯ hn+1 − ¯ n+1 )|ˆ hn+1 |2 + (|ˆ n+1 |2 − |ˆ hn+1 |2 )¯ n+1 ≤ C (e¯n+1 + eˆn+1 ) + Chr+1 ,
(4.13)
ˆ n+1 ) + ( ¯ nx +1 ) ¯ n+1 ( ˆ n+1 − ¯ n+1 − ˆ n+1 ≤ C (eˆn+1 + e¯n+1 ) + Chr+1 , J2n+1 = hx h hx
(4.14)
J4n+1 = ˆ 1nh+1 (¯ 1nh+1 − ¯ 1n+1 ) + ˆ 2nh+1 (¯ 2nh+1 − ¯ 2n+1 ) + (ˆ 1nh+1 − ˆ 1n+1 )¯ 1n+1 + (ˆ 2nh+1 − ˆ 2n+1 )¯ 2n+1 ≤ C (e¯n+1 + eˆn+1 ) + Chr+1 .
(4.15)
By Taylor formula, it is easy to see that, for j = 3, 5, 6
Jnj +1 2 ≤ C ( t )3
t n+1
tn
tt (s )2 ds.
(4.16)
Using Eqs. (4.13)–(4.16), we obtain
2 1 (en+1 2 − en 2 ) ≤ C (en+1 2 + en 2 + en−1 2 ) + C hr+1 + t 2 , t
(4.17)
There exists a constant k0 , for t ∈ [0, k0 ], such that 1 − C t > 0, thus, Eq. (4.16) can be written as
en+1 2 ≤
2 1 + C t n 2 e + C t en−1 2 + C t hr+1 + t 2 , 1 − C t
namely,
e + C t e ≤ n+1 2
n
2
(4.18)
2 1 + C t + C t (en 2 + en−1 2 ) + C t hr+1 + t 2 . 1 − C t
(4.19)
For n t ≤ T, n ≥ 1, we obtain
2 en+1 2 ≤ C e1 2 + t e0 2 + C hr+1 + t 2 .
(4.20)
Combining Eqs. (4.9), (4.20) and (4.5), we get
2 en+1 2 + t en +1 2 ≤ C hr+1 + t 2 .
(4.21)
150
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
From Eqs. (4.3) and (4.21), we obtain by triangle inequality
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 ,
0≤n≤M
We complete the proof of Theorem 4.1. Theorem 4.2. Let
{hn+1 , nh+1 }M n=0
satisfy the Eq. (3.9), the solutions ((t n+1 ), (t n+1 )) of Eq. (2.1) be sufficiently smooth and
2 h0 − 0 2 + h1 − (t 1 )2 + t ( 0h − 0 2 + 1h − (t 1 )2 ) ≤ C hr+1 + t 2 . Then, for h and k sufficiently small, there exists a constant C, such that
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 ,
(4.22)
(4.23)
0≤n≤M
where the constant C is independent of h and t. Proof. Subtracting Eq. (3.8) from Eq. (3.1) by the projection definition (4.1) and (4.2), we obtain the following error estimates in (en , en ), for n ≥ 1,
∂ n+1 ∂v ∂ n+1 ∂v e , −β e , + γ (J1n+1 , v ) + δ (J2n+1 , v ) ∂x ∂x ∂y ∂y
n+1
n−1 n 3E − 4E + E 3(t n+1 ) − 4(t n ) + (t n−1 ) 1 = −i − t (t n+ 2 ), v − i , v + (J3n+1 , v ), 2 t 2 t
i
3en+1 − 4en + en−1 ,v 2 t
−α
(en +1 , s ) + θ (J4n+1 , sx ) = (J5n+1 , sx ),
(4.24) (4.25)
where
J1n+1 = (2|hn |2 hn − |hn−1 |2 hn−1 ) − (2|(t n )|2 (t n ) − |(t n−1 )|2 (t n−1 )), J2n+1 = nhx+1 (2hn − hn−1 ) − x (t n+1 )(2(t n ) − (t n+1 )), J3n+1 = −γ [(2|(t n )|2 (t n ) − |(t n−1 )|2 (t n−1 )) − |(t n+1 )|2 (t n+1 )] −δ [ x (t n+1 )(2(t n ) − (t n−1 )) − x (t n+1 )(t n+1 )], J4n+1
= (2|hn |2 − |hn−1 |2 ) − (2|(t n )|2 − |(t n−1 )|2 ),
J5n+1 = θ (2|(t n )|2 − |(t n−1 )|2 − |(t n+1 )|2 ).
Choosing s = en +1 in Eq. (4.25), we obtain
en +1 ≤ C (J4n+1 + J5n+1 ).
(4.26)
Choosing c = en+1 in Eq. (4.24), we obtain
i t
n
j=1
= −i t
j+1 3e − 4en + en−1 j+1 , e 2 t
n j=1
−γ t
n j=1
− t
n
α (en+1 , ej+1 )
j=1
n−1 n 3E n+1 − 4E + E 3(t j+1 ) − 4(t j ) + (t j−1 ) j+1 − t (t n+1 ) + , e 2 t 2 t
(J1j+1 , ej+1 ) − δ t
n
(J2j+1 , ej+1 ) + t
j=1
n
(J3n+1 , ej+1 ).
(4.27)
j=1
First, we estimate the first term on the left-hand of Eq. (4.27)
t
n j=1
j+1 j j−1 3e − 4e + e j+1 , e 2 t
=
n
ej+1 2 − ej 2 −
j=1
3 ≥ en+1 2 − 4 3 ≥ en+1 2 − 4
1 n 2 e − 4 1 n 2 e − 4
1 j+1 2 1 e − ej−1 2 + ej+1 −ej 2 − ej+1 −ej−1 2 4 4 3 1 2 e + 4 3 1 2 e − 4
1 0 2 1 n+1 e + e − en 2 − e1 − e0 2 4 2 1 0 2 e . (4.28) 4
Next, using Taylor expansions and property (4.3), the right-hand side of Eqs. (4.25) and (4.27) have the following estimates, for j ≥ 1
tj 3(t j+1 ) − 4(t j ) + (t j−1 ) n+1 − ( t ) ≤ C t t t t ds, t 2 t t j−2
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
151
j+1 3E − 4Ej + Ej−1 1 = ( h − I )(3(t j+1 ) − 4(t j ) + (t j−1 )) t 2 2 t t j+1 t j+1 r+1 ≤ C ( − I ) ds ≤ Ch t r+1 ds, t h t j−1 t j−1 and n
Jlj+1 ≤ C t, l = 3, 5, Jlj+1 ≤ C (ej + ej−1 + hr+1 ), l = 1, 4,
j=1
J2j+1 ≤ C (ej + ej−1 + e j+1 ) + Chr+1 . Combining the above estimates, taking the imaginary part of Eq. (4.27), we obtain
en+1 2 ≤
4 1 n 2 e + e1 2 + e0 2 + t (J1j+1 + J2j+1 + J3j+1 + J4j+1 + J5j+1 )ej+1 . 3 3 n
(4.29)
j=1
Suppose m be chosen so that em = max0≤ j≤n e . Then it yields j+1
1 0 r+1 en+1 ≤ em + t 2 ). ≤ C e + e + C ( h
(4.30)
Hence, from Eqs. (4.3), (4.26), (4.22) and (4.30), we get by the triangle inequality
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 .
0≤n≤M
We complete the proof of Theorem 4.2.
5. Basis functions and implementation issues Let {ϕ j }2j=0 be the linear basis functions on a triangle element e , namely, ϕ j ∈ P1 (e ), j = 0, 1, 2, Vh |e = span{ϕ0 , ϕ1 , ϕ2 }. Assume (Un , Zn ) be the approximation of the exact solutions ( n , n ) = ((t n ), (t n )). To overcome the deficiency of solving complex algebraical system (3.8), we decompose the complex value function into real and imaginary parts, namely, (tn ) = u1 (tn ) + iu2 (tn ), and approximate u j (tn ) = unj by U nj ∈ Vh , namely, hn = U1n + iU2n . We denote Uj and Z on an element by 2
U je =
e U jm (t )ϕm (x, y ),
Ze =
m=0
2
e Zm (t )ϕm (x, y ).
m=0
Then, the CCN scheme (3.8), namely, Algorithm 3.1 has the following matrix form on an element:
Me
X e,n+1 − X e,n + Ae (Xˆ e,n+1 )X e,n+1 + Be X e,n+1 = 0, t
(5.1)
where X e,n = (U e,n , Z e,n ) = (U1e,n , U2e,n , Z e,n ) , the element matrices are, for i, s = 0, 1, 2,
˜e = M is Ae (Uˆ j
D (Uˆ j e
e,n
e,n
)is = )is =
A˜ e (Uˆ e,n )is =
ϕi ϕs dx dy, Be1,is = e
2
e
e
e
˜e M 0
M˜ e e
M =
2
e
ϕix ϕsx dx dy, Be2,is =
ϕi ϕs dx dy,
e
ϕiy ϕsy dx dy,
j = 1, 2
e,n Uˆ jm
ϕm ϕix ϕsx dx dy,
j = 1, 2
m=0
⎛ 2 2 ⎞ 2 2 ⎝ ⎠ϕi ϕs dx dy, Uˆ1e,n ϕ + Uˆ2e,n ϕ m m m m
Ee =
e,n Uˆ jm ϕm
m=0
0 0
m=0
0 , ˜e −M 0 ˜e −M 0
m=0
F e (Uˆ e,n ) = 0 0 , 0
0 γ A˜ e (Uˆ e,n )
A (U e
ˆ e,n
)=
γ A˜ e (Uˆ e,n ) 0
0 e ˆ e,n ˜ γ A (U ) e,n −2θ De (Uˆ1 )
,
Ge =
γ A˜ e (Uˆ e,n ) 0 e,n −2θ De (Uˆ2 )
0 −α Be1 − β Be2 e,n δ Ae (Uˆ2 ) e,n δ Ae (Uˆ1 ) ,
0
−α Be1 − β Be2 , 0
152
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
e
B =
−α
Be1
0 − β Be2 0
−α Be1 − β Be2 0 0
0 0 . −Be1 − Be2
˜ e , Be , Be , Ae , De and A˜ e is 3 × 3. Assembling contributions from all elements leads to the The size of element matrices M 1 2 following matrix form of the linear algebraic equations, from Eq. (5.1),
M
X n+1 − X n + A(Xˆ n+1 )X n+1 + BX n+1 = 0, t
(5.2)
for n ≥ 0, where Xˆ n+1 = 32 Xˆ n − 12 Xˆ n−1 and X = (U, Z ) contains all the nodal parameters. Similarly, the matrix form on an element e of the DBDF2 method (3.9), namely, Algorithm 3.2 can be written as follows:
⎧ ⎨ 3U n+1 − 4U n + U n−1 E + E (U n )U n+1 + GU n+1 = bn , 2 t ⎩(−B − B )Z n+1 = 2θ DU n+1 U n+1 + DU n+1 U n+1 , 1 2 1 1 2 2
(5.3)
where be,n = (−δ Ae (U2 e,n )Z n , −δ Ae (U2 e,n )Z n ) . ˜ , A, B, E, F and F with the total domain are sparse block matrices, the linear As mentioned above, the matrices M, M algebraic system (5.2) and (5.3) can be respectively implemented efficiently. 6. Numerical experiments In this section, four experiences: the propagation of elliptic-elliptic and hyperbolic-elliptic DS equations are presented to illustrate the accuracy and L2 norm conservation laws of the numerical schemes based on CCN method and DBDF2 method. Moreover, we extend the proposed numerical methods to study the dynamic of rogue waves including line rogue wave for elliptic-elliptic and hyperbolic–elliptic DS equations, two-rogue waves for elliptic–elliptic DS equations. To illustrate the accuracy of the numerical methods, we denote the error norms between numerical solution Uhn and accuracy solution U(tn ) as follows, for j = 1, 2, U = ,
1 Uhn − U (t n )2 dx dy 2 , L∞ = Uhn − U (t n )∞ = max Uhn − U (t n ), L2 = Uhn − U (t n )2 =
H1 =
(x,y )∈
1 n n − Ux (t n )22 + Uhy − Uy (t n )22 2 , Uhx
(6.1)
where the nine points Gauss numerical quadrature formula for the integration is used. 6.1. Convergence and accuracy Example 1. Elliptic–elliptic DS equations The exact solutions of Eq. (2.1) are taken as [36]
⎧ ⎪ ⎨(x, y, t ) = 2 2
1 exp(iη ), 3 coth(ξ )
⎪ ⎩ (x, y, t ) = 4 (ξ − tanh(ξ )),
(6.2)
3
where
ξ = x + y − 4t, η = x + y − 6t. We choose the parameter α = 1, β = −1, γ = −1, δ = −1 and θ = −1 with the domain [−4, 4]2 . The evolution of traveling wave is done up to t = 1. To illustrate the accuracy of our algorithms, we compute the pointwise convergence rate in space. Setting the temporal h 1 1 1 step size t = 16 , we list the L2 , L∞ and H1 error accuracy of space in Tables 1 and 2 with varying spacing h = 10 , 20 , 40 at time t = 1. From the second and fourth columns in Table 1 based on CCN method and Table 2 based on DBDF2 method, the accuracy of the L2 and L∞ norms for and are closed to second order in space and the H1 convergence rate is O(h) as the sixth column shown. Therefore, the numerical experiment strongly suggests that the orders of convergence in space in error estimates for the L2 , L∞ and H1 norms of and are optimal and confirm with the theoretical results. To investigate the order of convergence in time, we take a more accurate, but slightly more sophisticated, approach to examine the orders of convergence with respect to the time step t because the approximation errors O( t γ ) + O(hμ ) are dominated by the h-terms O(hμ ) [37]. For instance, assuming
vh, t ≈ v(x, t m ) + C1 (x, t m ) t γ + C2 (x, t m )hμ ,
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161 Table 1 The order of convergence in space for error norms with t =
| |
by CCN method.
h
L∞
Order
L2
Order
H1
Order
1/10 1/20 1/40
1.257606E−1 3.089576E−2 7.673786E−3
2.0252 2.0093
3.310235E−1 8.054239E−2 1.987874E−2
2.0391 2.0185
8.686715E−1 3.506051E−1 1.639445E−1
1.3089 1.0966
1/10 1/20 1/40
2.025977E−2 5.158323E−3 1.300563E−3
1.9736 1.9877
4.412639E−2 1.157752E−2 2.952151E−3
1.9303 1.9715
1.439311E−1 6.069855E−2 2.869494E−2
1.2456 1.0808
Order
H1
Order 1.3053 1.1072 1.2504 1.0768
Table 2 The order of convergence in space for error norms with t =
| |
h 16
153
∞
h 16
by DBDF2 method.
Order
L2
1.337010E−1 3.243466E−2 8.109017E−3
2.0434 1.9999
3.475989E−1 8.491867E−2 2.099681E−2
2.0332 2.0159
8.852042E−1 3.581702E−1 1.662515E−1
1.935581E−2 5.015085E−3 1.248773E−3
1.9484 2.0057
4.571451E−2 1.202377E−2 3.069377E−3
1.9267 1.9698
1.437976E−1 6.044263E−2 2.865467E−2
h
L
1/10 1/20 1/40 1/10 1/20 1/40
it follows that
1 t vh, t − vh, 2 ≈ C1 (x, t m ) 1 − γ t γ . 2
Thus, t
ρv, t =
vh, t − vh, 2 4γ − 2γ ≈ γ = 2γ . t t 2 −1 vh, 2 − vh, 4
Here, v can be and . In particular, ρv, t approaches to 4 for γ = 2 when the corresponding order of convergence in 1 time is of O( t2 ). Fig. 1 shows the temporal accuracy in L2 , L∞ and H1 norms with fixed the spatial step size h = and 10 0.025 different temporal partition t = k , k = 0, 1, 2, 3. These results match the regular expectations of accuracy order O( t2 ) 2 arising from CCN method and DBDF2 method. Moreover, we compare CPU time needed for solving the DS problem (2.1) using the proposed CCN and DBDF2 numerical schemes and run simulations done up to t = 1. The CPU time as listed in Table 3 is given in seconds and includes the time cost by generating mesh structure, assembling the stiffness matrix and load vector, solving the corresponding algebraic system and computing errors on all the time steps. The coupled linear system is needed to solve for CCN method at each time, which causes that the DBDF2 numerical method needs much less CPU time than CCN method as presented in Table 3. Furthermore, we discuss numerical L2 norm conservation as presented in Eq. (2.3). The changes of the invariant from the initial variant approach zero by proposed Algorithms 3.1 and 3.2 as shown in Fig. 2. These results display that both the precision and the conservation are well in the proposed two numerical schemes. Example 2. Hyperbolic–elliptic DS equations Consider the following exact solutions in [24]
⎧ 2 2 2 2 ⎪ ⎪(x, y, t ) = 2(B1 + B2 )(ω + α k1 + β k2 ) sech(ξ ) exp(iη ), ⎨ 2 2 B1 (γ + δθ ) + γ B2 2 2 ⎪ −2 θ B 1 ( ω + α k1 + β k2 ) ⎪ ⎩ (x, y, t ) = − tanh(ξ ), B21 (γ + δθ ) + γ B22
(6.3)
where
ξ = B1 x + B2 y − ν t, η = −k1 x − k2 y − ωt, ν = −2(α k1 B1 + β k2 B2 ). Choosing parameters α = β = 12 , γ = −1, δ = 2 and θ = −2 in DS system (2.1), and B1 = B2 = 1, k1 = −0.1, k2 = −0.05 and ω = 0.05 in system (6.3) with the domain [−10, 10]2 and terminated time t = 4 of the evolution. The L2 , L∞ and H1 error accuracy of space are reported in Tables 4 and 5 with varying spacing h = 12 , 14 , 18 at time t = 4, which are almost completely closed to spatial order O(h2 ), O(h2 ) and O(h), respectively, and confirm with the theoretical results. Fig. 3 shows the temporal order of error norms can reach the theoretical result O( t2 ), which illustrates the efficiency and accuracy of the presented algorithms in this paper. In order to compare the proposed CCN and DBDF2 methods, we choose step sizes h = 16 and t = 0.1. Tables 6 illustrates that both CCN method and DBDF2 method have higher accuracy than the Lie-Trotter time-splitting (LTTS) scheme in [31],
154
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 1. Error norms at time t = 1.0 with different temporal step size t =
0.025 ,k 2k
= 0, 1, 2, 3 and the fixed mesh h =
1 . 10
due to the fact that the CCN and DBDF2 methods can derive the order convergence in time O( t2 ) while LTTS method only has first-order temporal convergence. Moreover, our numerical methods for both CCN method and DBDF2 method conserve the L2 norm (2.3) as presented in Fig. 4. The accuracy, efficiency and conservation of our proposed two numerical methods are verified again.
6.2. Application Example 3. Dynamics of fundamental rogue waves In this example, we consider the dynamics of the one line rogue waves for hyperbolic–elliptic and elliptic–elliptic DS Eq. (2.2) as presented in [26,27]. In the first case, we study the elliptic-elliptic type equation. Taking parameters α = 1,
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161 Table 3 The CPU time of different numerical methods with t =
155
h . 16
Method
h = 1/10
h = 1/20
h = 1/40
CCN DBDF2
203 186
2311 2144
30,576 28,107
Fig. 2. The L2 norm conservation M of numerical solutions by different numerical methods with h = Table 4 The order of convergence in space for error norms with t =
| |
∞
by CCN method.
Order
L2
8.768628E−4 2.223408E−4 5.579751E−5 4.424144E−4 1.111999E−4 2.806503E−5
h
L
1/2 1/4 1/8 1/2 1/4 1/8
Order
H1
Order
1.9795 1.9945
4.551711E−3 1.140348E−3 2.853134E−4
1.9969 1.9988
3.863076E−2 1.932430E−2 9.663248E−3
0.9933 0.9998
1.9819 1.9962
2.380935E−3 5.968542E−3 1.493584E−4
1.9961 1.9986
1.968664E−2 9.853213E−3 4.927838E−3
0.9986 0.9996
Table 5 The order of convergence in space for error norms with t =
| |
h 16
h
L∞
1/2 1/4 1/8 1/2 1/4 1/8
1 . 40
h 16
by DBDF2 method.
Order
L2
Order
H1
Order
8.749512E−4 2.220580E−4 5.575367E−5
1.9782 1.9938
4.547120E−3 1.139739E−3 2.852320E−4
1.9969 1.9985
3.863045E−2 1.932428E−2 9.663247E−3
0.9933 0.9998
4.389001E−4 1.112386E−4 2.789141E−5
1.9802 1.9957
2.332774E−3 5.850946E−3 1.464577E−4
1.9953 1.9981
1.968638E−2 9.853186E−3 4.927835E−03
0.9985 0.9996
β = −1, γ = 2, δ = 2 and θ = 1 in DS system (2.3), the initial conditions are given as √
(x, y, 0 ) = 2 1 −
√
4
1 + ( 3x − y )2
,
√
(x, y, 0 ) = −1 − 4 3
√ 1 − ( 3x − y ) 2 . √ [1 + ( 3x − y )2 ]2
(6.4)
with domain [0, 0.5] × [−10, 10] and time interval t ∈ [−1, 0.25]. The solutions describe a line wave with the line oriented √ in the ( 12 , 23 ) direction of the (x, y) plane, and the orientation angle is π3 . This fundamental rogue wave is illustrated in Fig. 5. At time t = 0 in Fig. 5(e), reaches the maximum amplitude at the center of the line.
156
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 3. Error norms at time t = 4.0 with different temporal step size t =
0.4 ,k 2k
= 0, 1, 2, 3 and the fixed mesh h =
1 . 8
Table 6 Comparison error norms of different methods at time t = 4. Method
CCN DBDF2 LTTS [31]
| |
L∞
L2
H1
L∞
L2
H1
2.095568E−4 2.056871E−4 6.535216E−4
1.016200E−3 1.006211E−3 4.653657E−3
1.306883E−2 1.306835E−2 1.379714E−2
9.637222E−5 9.824542E−5 3.218093E−4
4.727967E−4 4.955172E−4 1.820333E−3
6.554249E−3 6.5550 0 0E−3 6.785512E−3
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 4. The L2 norm conservation M of numerical solutions by different numerical methods with h =
157
1 . 6
Fig. 5. The profiles of numerical solutions | | at different time.
In the second case, the evolution of line rogue waves of hyperbolic-elliptic type DS equations is studied. Most of the parameters used in this simulation are the same as those in the first case, except β = γ = θ = −1. We take the following initial conditions
√ (x, y, 0 ) = 2 1 −
4 , 1 + ( k1 x + k2 y )2
(x, y, 0 ) = 1 − 4k21
1 − ( k1 x + k2 y )2 . [1 + (k1 x + k2 y )2 ]2
(6.5)
158
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 6. The profiles of numerical solutions | | at different time.
Fig. 7. The profiles of numerical solutions | | at different time.
with k1 = 56 and k2 = 13 6 . The computational domain is [−4, 4] × [0, 1] and time interval is t ∈ [−1, 1]. The solutions describe 5 a line rogue wave with the line oriented in the ( 13 6 , − 6 ) direction of the (x, y) plane, thus the fundamental rogue waves in the DS equations are line rogue waves. Fig. 6 shows the evolution of fundamental rogue wave. The maximum amplitude can be reached at time t = 0 in Fig. 6(d). Example 4. Dynamics of two-rogue waves In this example, we consider the dynamics of two-rogue waves in DS equations, which describes the interaction between two fundamental line rogue waves. Choosing parameters α = −1, β = 1, γ = 1, δ = 2, θ = 1 and m = 1 in DS system (2.2).
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
159
Fig. 8. The profiles of numerical solutions at different time.
Owing to the works in [27], the initial conditions can be given
√ τ 2 = 2 0 , = −1 − 2
τ1
τ1
2 ∂ 2 τ1 ∂τ1 τ1 2 − , ∂x ∂x
(6.6)
where τs , s = 0, 1 is the determinant of 2 × 2 block matrix with each block 2 × 2,
m (s ) ij τs = (−s) −m i j
and (s )
mi j
1 = pi + q j
ξi =
, mi(j−s ) i(js ) m
ξi + s −
(6.7)
pi + ci1 , pi + q j
pi + 1/pi pi − 1/pi x+ −1y. 2 2
i(js ) = m
1 pi q j − 1
ξi + s −
pi q j pi q j − 1
+ ci1 ,
√ √ with the parameters p1 = 1, p2 = −1, q1 = 0, q2 = −3, c11 = 0 and c21 = 2−1 on the domain [−10, 10]2 and time interval 1 t ∈ [−0.3, 0.6]. The corresponding numerical solutions with step size h = 30 and t = 0.001 is displayed in Figs. 7–10. It is seen that a cross-shape rogue wave appears, which describes the interaction between two fundamental line rogue waves: one oriented along the y-direction (corresponding to the parameters p1 ) and the other one oriented along the x-direction √ (corresponding to the parameters p2 ). These two individual line waves reach their peak amplitude 3 2 at different time with the x-direction one peaking at t ≈ − 41 as shown in Fig. 7(a) and the y-direction one peaking at t ≈ 0 in Fig. 7(c). This phenomenon is similar to the results as presented in [27], which validates the efficiency of our numerical scheme.
7. Conclusion In this paper, numerical algorithms for the numerical solutions of the Davey–Stewartson equations are proposed using Galerkin finite element with second order time discretization, which employs linearized Crank–Nicolson scheme and decoupled multistep scheme. The high efficiency and accuracy of our numerical methods are demonstrated by the numerical experiments: the propagation of elliptic-elliptic and hyperbolic-elliptic DS equations. The numerical results recorded in
160
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 9. Contours of numerical solutions | | at different time.
Fig. 10. Contours of numerical solutions at different time.
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
161
Tables 1–4 and Figs. 1 and 3 show that the proposed numerical scheme have good accuracy. Moreover, our numerical methods satisfy L2 norm conservation laws as shown in Figs. 2 and 4. We extend the proposed numerical methods to study the dynamics of rogue waves, including fundamental rogue wave and interaction of two line rogue waves. What is more, the proposed methods are very efficient since their implementation requires at each time step the solutions of a linear system with the same matrix for all time levels. References [1] F. Linares, G. Ponce, Introduction to Nonlinear Dispersive Equations, Springer, New York, 2009. [2] J.A. Pava, Nonlinear Dispersive Equations: Existence and Stability of Solitary and Periodic Travelling Wave Solutions, American Mathematical Society, Providence, R.I., 2009. [3] Y.V. Kartashov, B.A. Malomed, L. Torner, Solitons in nonlinear lattices, Rev. Mod. Phys. 83 (1) (2011) 247–305. [4] R. Carles, E. Dumas, C. Sparber, Multiphase weakly nonlinear geometric optics for Schrödinger equations, SIAM J. Math. Anal. 42 (2010) 489–518. [5] R. Thomas, C. Kharif, M. Manna, A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity, Phys. Fluids 24 (12) (2012) 395–416. [6] L. Chang, Y. Pan, X. Ma, New exact travelling wave solutions of Davey–Stewartson equation, J. Comput. Inf. Syst. 9 (2013) 1687–1693. [7] M. Dehghan, M. Abbaszadeh, A. Mohebbi, Analysis of two methods based on Galerkin weak form for fractional diffusion-wave: meshless interpolating element free Galerkin (IEFG) and finite element methods, Eng. Anal. Bound Elem. 64 (4) (2016) 205–221. [8] M. Song, A. Biswas, Topological defects and bifurcation analysis of the DS equation with power law nonlinearity, Appl. Math. Inf. Sci. 9 (4) (2015) 1719–1724. [9] M. Dehghan, M. Abbaszadeh, A finite element method for the numerical solution of Rayleigh–Stokes problem for a heated generalized second grade fluid with fractional derivatives, Eng. Comput. Ger. 33 (3) (2017) 587–605. [10] M. Dehghan, A. Shokri, Numerical solution of the nonlinear Klein–Gordon equation using radial basis functions, J. Comput. Appl. Math. 230 (2) (2009) 400–410. [11] C. Klein, B. Muite, K. Roidot, Numerical study of blowup in the Davey–Stewartson system, Discr. Cont. Dyn. Syst. B 18 (2013) 1361–1387. [12] M. Dehghan, Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices, Math. Comput. Simul. 71 (1) (2006) 16–30. [13] C. Klein, K. Roidot, Fourth order time-stepping for Kadomtsev–Petviashvili and Davey–Stewartson equations, SIAM J. Sci. Comput. 33 (2011) 3333–3356. [14] M. Dehghan, F. Fakhar-Izadi, The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves, Math. Comput. Model. 53 (9) (2011) 1865–1877. [15] A. Davey, K. Stewartson, On three-dimensional packets of surface waves, Proc. R. Soc. Lond. A 338 (1974) 101–110. [16] C. Babaoglu, Long-wave short-wave resonance case for a generalized Davey–Stewartson system, Chaos Soliton Fract. 38 (2008) 48–54. [17] V.D. Djordjevic´ , L.G. Redekopp, On two-dimensional packets of capillary-gravity waves, J. Fluid Mech. 79 (1977) 703–714. [18] K. Nishinari, K. Abe, J. Satsuma, A new-type of soliton behavior of the Davey–Stewartson equations in a plasma system, Teoret. Mat. Fiz. 99 (1994) 487–498. [19] J.M. Ablowitz, R. Haberman, Nonlinear evolution equations in two and three dimensions, Phys. Rev. Lett. 35 (1975) 1185–1188. [20] A. Eden, T.B. Gürel, On the integrability of a generalized Davey–Stewartson system, Phys. D 259 (2013) 1–7. [21] H. Jafari, A. Sooraki, Y. Talebi, A. Biswas, The first integral method and traveling wave solutions to Davey–Stewartson equation, Nonlinear Anal. Model. Control 17 (2) (2012) 182–193. [22] G. Ebadi, A. Biswas, The G’/G method and 1-soliton solution of the Davey–Stewartson equation, Math. Comput. Model. 53 (2011) 694–698. [23] O.H. El-Kalaawy, R.S. Ibrahim, Solitary wave solution of the two-dimensional regularized long-wave and Davey–Stewartson equations in fluids and plasmas, Appl. Math. Ser. B 3 (8) (2012) 833–843. [24] G. Ebadi, E.V. Krishnan, M. Labidi, E. Zerrad, A. Biswas, Analytical and numerical solutions to the Davey–Stewartson equation with power-law nonlinearity, Wave Random Complex 21 (2011) 559–590. [25] J. Shi, J. Li, S. Li, Analytical travelling wave solutions and parameter analysis for the (2+1)-dimensional Davey–Stewartson-type equations, Pramana J. Phys. 81 (2013) 747–762. [26] V.A. Arkadiev, A.K. Pogrebkov, M.C. Polivanov, Inverse scattering transform method and soliton solutions for the Davey–Stewartson II equation, Phys. D 36 (1989) 189–196. [27] Y. Ohta, J. Yang, Dynamics of rogue waves in the Davey–Stewartson II equation, J. Phys. A Math. Theor. 46 (10) (2013) 105202. [28] Y. Ohta, J. Yang, Rogue waves in the Davey–Stewartson I equation, Phys. Rev. A Math. Theor. 86 (2012) 036604. [29] Y. Zhang, H.P. Stimming, N.J. Mauser, A novel nonlocal potential solver based on nonuniform FFT for efficient simulation of the Davey–Stewartson equations, ESAIM Math. Model. Num. 51 (4) (2017) 1527–1538. [30] C. Besse, N.J. Mauser, H.P. Stimming, Numerical study of the Davey–Stewartson system, Math. Model. Numer. Anal. 38 (2004) 1035–1054. [31] Y. Gao, L. Mei, R. Li, A time-splitting galerkin finite element method for the Davey–Stewartson equations, Comput. Phys. Commun. 197 (2015) 35–42. [32] J.M. Ghidaglia, J.C. Saut, On the initial value problem for the Davey–Stewartson systems, Nonlinearity 3 (1990) 475–506. [33] C. Klein, K. Roidot, Numerical study of the semiclassical limit of the Davey–Stewartson II equations, Nonlinearity 27 (2014) 2177–2214. [34] G.M. Muslu, H.A. Erbay, Numerical simulation of blow-up solutions for the generalized Davey–Stewartson system, Int. J. Comput. Math. 88 (2011) 805–815. [35] Y. Gao, L. Mei, Implicit-explicit multistep methods for general two-dimenisonal nonlinear Schrödinger equations, Appl. Numer. Math. 109 (2016) 41–60. [36] M.A. Bashir, A.A. Moussa, The cotha (ξ ) expansion method and its application to the Davey–Stewartson equation, Appl. Math. Sci. 8 (2014) 3851–3868. [37] L. Shan, H. Zheng, Partitioned time stepping method for fully evolutionary Stokes–Darcy flow with the Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 51 (2013) 813–839.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Galerkin methods for the Davey–Stewartson equationsR Yali Gao a, Liquan Mei a,∗, Rui Li b a b
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China School of Mathematics and Information Science, Shaanxi Normal University, Xi’an, 710062, China
a r t i c l e
i n f o
Keywords: Davey–Stewartson equations Galerkin finite element method Extrapolated Crank–Nicolson method Explicit multistep method
a b s t r a c t In this paper, we propose two Galerkin methods to investigate the evolution of the Davey– Stewartson equations. The extrapolated Crank–Nicolson scheme and decoupled semiimplicit multistep scheme are employed to increase the order of the time discrete accuracy, which only requires the solutions of a linear system at each time step. Four numerical experiments are presented to illustrate the features of the proposed numerical methods, such as the optimal convergence order, the conservation variable and the application in rogue waves. © 2018 Elsevier Inc. All rights reserved.
1. Introduction A variety of phenomena across a range of disciplines, which in applied mathematics and physics from fluid dynamics, solid state physics, quantum machines, plasmas physics to nonlinear optics and so forth, and in chemistry and biology as well [1,2], are described by nonlinear dispersive partial differential equations. For instance, the Schrödinger equation is the fundamental governing equation in quantum mechanics and quantum field theory [3], which is used to describe many-body theory and condensed matter physics like the Bose–Einstein condensate. It also has extensive applications in nonlinear optics [4] and water waves [5]. In view of the importance of these equations, many researchers tried their best to find analytical solutions or numerical solutions to the nonlinear dispersive equations by using various methods [6–14]. The Davey–Stewartson (DS) equations were originally introduced by Davey and Stewartson [15] to describe the evolution of a three-dimensional weakly nonlinear wave-packet on water in finite depth, in which the water waves have one direction of travel but the amplitude of waves is modulated in two spatial directions. The equations appear in studying the longwave-short-wave resonances [16] and describing gravity-capillarity surface wave packets in the limit of shallow water [17]. The evolution of a plasma under the action of a magnetic field in plasma physics [18] and certain (2+1)-dimensional nonlinear localized excitations can be described by DS equations. Furthermore, the (2+1)-dimensional Schrödinger equation with cubic nonlinearity fails to be integrable. However, Davey–Stewartson equations, a pair of coupled nonlinear equations in two dependent variables, were shown to be integrable [19–21] as the extension of the (1+1)-dimensional nonlinear Schrödinger equation. In fact, when the DS system is considered with various parameters in it, there are only two integrable cases which are known as DS I and DS II [21]. These two integrable cases are of elliptic-hyperbolic and hyperbolic-elliptic type, depending on the strength of surface tension. The DS II equations are thus the hyperbolic-elliptic ones, which are in particular relevant for surface gravity waves without surface tension. In addition, for DS I and DS II equations, special localized solutions, called respectively dromions and lumps, can be explicitly found. R ∗
The project is supported by China Scholarship Council (201506280175) and NSF of China (11371289, 11501441). Corresponding author. E-mail addresses: [email protected] (Y. Gao), [email protected] (L. Mei), [email protected] (R. Li).
https://doi.org/10.1016/j.amc.2018.01.044 0 096-30 03/© 2018 Elsevier Inc. All rights reserved.
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
145
Many works for Davey–Stewartson equations have been done in recent years. A vast method finding analytical solutions has been proposed [6,8,21–28]. Biswas [22] considered G’/G method and the solitary wave ansatz method to study traveling wave solution and 1-soliton solutions, respectively. El-Kalaawy and Ibrahim [23] obtained the solitary wave solutions by using extended mapping method technique. Ohta and Yang [27,28] derived the general rogue waves by the bilinear method. In parallel with the analytical studies, a surge of efforts have been devoted to the numerics of this system [11,13,29–34]. Stimming and Zhang [29] presented an efficient and accurate solver for the nonlocal potential in the Davey–Stewartson equations using nonuniform Fast Fourier Transform. Besse et al. [30] proposed a time splitting spectral method for numerical study of the DS equations, including exact soliton type solutions of the hyperbolic-elliptic equations and the blow up of focusing elliptic-elliptic equations. Klein and Roidot [13] investigated the semi-classical Davey–Stewartson equations by both fourth order time-stepping and spectral method. Gao et al. [31] proposed a Galerkin method with first-order time discretization. To the best of our knowledge, there are few results on numerics of the DS equations by finite element method, though there are numerous investigations for nonlinear dispersive equations [31,35]. In this paper, we study the numerical solutions of Davey–Stewartson equations by Galerkin finite element method with extrapolated Crank–Nicolson Scheme and decoupled multistep method. Due to the nonlinearity and coupling and complex structure system of DS equations, linearized semi-implicit schemes should be considered approximations in time direction since their implementation only requires the solution of a linear system. Among linearized schemes, the extrapolated Crank– Nicolson method is a popular one since it provides a second order accuracy in the time direction. However, linearized Crank–Nicolson method only solved the nonlinearity of the DS equations, we also need solve a coupled algebraic system. To overcome this shortcoming of coupled Crank–Nicolson method, we develop the decoupled scheme by using semi-implicit multistep method, which improves the efficiency of algorithm and remains the second order convergence rate in time as Crank–Nicolson time discretization. The outline of this paper is as follows. In the next section, we introduce the governing equations for DS equations. In Section 3, numerical algorithms of the DS equations are given. In Section 3.1, we study the weak form and the proof of L2 norm conservations. The full discretization based on linearized Crank–Nicolson method and decoupled multistep method are presented in Section 3.2 and Section 3.3, respectively. In Section 4, the L2 error estimates for extrapolated Crank–Nicolson and decoupled implicit-explicit two-step schemes are presented and the proof are given. In Section 5, some details during the progress of implementation are presented. In Section 6, we propose four numerical experiments for DS equations, which illustrates the efficiency of the presented numerical algorithms by the L2 , L∞ , H1 error accuracy and the relative error of L2 norm conservation laws. Furthermore, we apply the proposed method to the rogue waves in DS equations. Finally, a brief conclusion is given in the last section.
2. Governing equations Let = [L1 , L2 ] × [L3 , L4 ] be a bounded subset of R2 and T > 0 be a final time, J = (0, T ], T = × (0, T ). We consider the following (2+1)-dimensional Davey–Stewartson equations with power law nonlinearity
⎧ ∂ ∂ 2 ∂ 2 ⎪ ⎪ i +α +β + γ | |2 + δ x = 0, ⎪ 2 ⎪ ∂t ∂x ∂ y2 ⎪ ⎪ ⎨ ∂ 2 ∂ 2 ∂| |2 + +θ = 0, 2 2 ∂x ∂x ∂y ⎪ ⎪ ⎪ ⎪ = 0, = 0, ⎪ ⎪ ⎩ = 0 , = 0 ,
in T , in T , on
(2.1)
∂ × [0, T ],
on × {t = 0},
where (x, y, t) is a complex valued function which denotes the complex amplitude,while (x, y, t) is a real valued function which denotes the velocity of a underlying mean flow. Using the terminology in [31], sgnα = 1 and sgnα = −1 denote elliptic-elliptic DS and hyperbolic-elliptic DS, respectively. For both there is a focusing (sgnγ = 1) and a defocusing (sgn √γ = −1) version [29]. β , δ , θ are real parameters. 0 : → C and 0 : → R are the given initial conditions, and i = −1. In absence of y dependence, the DS equations can be viewed as a generalization of the (1+1)-dimensional cubic nonlinear Schrödinger equation,
i
∂ ∂ 2 +α + γ | |2 = 0. ∂t ∂ x2
such that = x instead of Moreover, the Davey–Stewartson system has in fact another form (when using the variable in (2.1), we still denote it by for simplicity)
⎧ ∂ ∂ 2 ∂ 2 ⎪ +α +β + γ | |2 + δ = 0, ⎨i 2 ∂t ∂x ∂ y2 2 2 2 2 ⎪ ⎩ ∂ + ∂ + θ ∂ | | = 0. ∂ x2 ∂ y2 ∂ x2
(2.2)
146
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Furthermore, the DS Eq. (2.1) possess two formal conservation laws [30] corresponding to the L2 norm
M (t ) =
| |2 dx dy,
(2.3)
and the energy (Hamiltonian)
E (t ) =
1 α|x |2 + β|y |2 + γ | |4 − δθ (∂x 2 + ∂y 2 ) dx dy. 2
3. Numerical methods We use the standard Sobolev space Wm,k (), where m is a nonnegative integer with 1 ≤ k ≤ ∞, and the space involving time Lk (0, T; Wm,k ) with norm · Lk (0,T ;W m,k ) . Let W m,2 () = H m () with the norm · m and ( · , · ) be the inner product on L2 () = H 0 () defined by
(u, v ) =
u(x )(v(x ))∗ dx dy,
∀ u, v ∈ L2 (),
where v∗ denotes the conjugate of v. The norm · ∞ denotes the essential supremum. Set V = H01 () = {v ∈ H 1 () : v|∂ = 0}. Applying Greens formula to problem (2.1) and using the boundary condition in the definition of V, the weak formulation of problem (2.1) can be written as follows: Find ( , ): J → V × V, such that, for ∀ (v, s ) ∈ V × V,
⎧ ∂ ∂v ∂ ∂v ∂ ⎪ ⎪ i , v − α , − β , + γ (| |2 , v ) + δ ( x , v ) = 0, ⎪ ⎪ ∂x ∂x ∂y ∂y ⎨ ∂ t
∂s ∂ ∂s ∂ ∂s , + , + θ | |2 , = 0, ⎪ ⎪ ∂ x ∂ x ∂ y ∂ y ∂ x ⎪ ⎪ ⎩ ( , v ) = (0 , v ), ( , s ) = ( 0 , s ).
(3.1)
3.1. Semi-discretization in space Let us consider a uniform regular partition h of triangular element e with = maxe ∈ h diam(e ). Set
e
e . We denote mesh size by h =
¯ , v|e ∈ Pr (e ), ∀ e ∈ h } ⊂ H 1 (), Vh = {v : v ∈ V ∩ C 0 () 0 where Pr (e ) denotes the polynomial space with degree less than or equal to r defined on e . Then from Eq. (3.1), we obtain the following semi-discrete form for problem (2.1): Find ( h , h ): J → Vh × Vh , such that, for ∀ (v, s ) ∈ Vh × Vh ,
⎧ ∂ h ∂ h ∂v ∂ h ∂v ⎪ ⎪ i , v − α , − β , + γ (|h |2 h , v ) + δ ( hx h , v ) = 0, ⎪ ⎪ ∂x ∂x ∂y ∂y ⎨ ∂ t
∂s ∂ h ∂ s ∂ h ∂ s , + , + θ | h | 2 , = 0, ⎪ ⎪ ∂ x ∂ x ∂ y ∂ y ∂ x ⎪ ⎪ ⎩ 0 (h , v ) = (0 , v ), ( 0h , s ) = ( 0 , s ).
(3.2)
Proposition 3.1. (L2 norm conservation) The solutions to the scheme (3.2) satisfy the L2 norm conservation
d d M= | |2 d x d y = 0. dt dt h
(3.3)
Proof. First, the complex conjugate of first equation of Eq. (3.2) is
∗ ∗
∂ h ∂v ∂ h ∂v ∂ h∗ −i ,v − α , −β , + γ (|h |2 h∗ , v ) + δ ( hx h∗ , v ) = 0. ∂t ∂x ∂x ∂y ∂y
(3.4)
Choosing v = h in the first equation of Eq. (3.2) and v = h∗ in Eq. (3.4), we obtain by subtracting
∗
∂ h ∗ ∂ h , h + , h = 0 , ∂t ∂t
(3.5)
namely,
d | |2 d x d y = 0. dt h
Then, we obtain the L2 norm conservation (3.3) of Eq. (3.2).
(3.6)
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
147
We now consider full discrete scheme of Eq. (3.2). For time discretization, we discuss the second order schemes: coupled Crank–Nicolson (CCN) method and decoupled BDF2 (DBDF2) method. Let 0 = t 0 < t 1 < · · · < t M = T be a uniform partition of J into subintervals J n = (t n , t n+1 ), n = 0, 1, . . . , M − 1, with length t = t n+1 − t n = MT . Then tn = n t, n = 0, 1, . . . , M denote the discrete time levels. (Un , Zn ) ∈ Vh × Vh define the fully-discrete approximations of the exact solutions ( (tn ), (tn )) ∈ V × V of DS problem (2.1). 3.2. Full discretization with extrapolated Crank–Nicolson method The fully discrete Galerkin approximation of Eq. (3.2) is given as follows: Find (hn+1 , nh+1 ) ∈ Vh × Vh , such that, for all 0 ≤ n ≤ M − 1, ∀ (v, s ) ∈ Vh × Vh ,
n+1 n+1
⎧ n+1 h − hn ∂ ¯ h ∂ ¯ h ∂v ∂v ⎪ ¯ n+1 |2 ¯ n+1 , v ) + δ ( ¯ n+1 ¯ n+1 , v ) = 0, ⎪ i , v − α , − β , + γ ( | ⎨ h h hx h t ∂x ∂x ∂y ∂y n+1
n+1
¯ ¯ ∂ ∂ ⎪ ∂s ∂s ⎪ h h ¯ n+1 |2 , ∂ s = 0, ⎩ , + , + θ | h ∂x ∂x ∂y ∂y ∂x
(3.7)
where
¯ hn+1 =
hn+1 + hn 2
,
¯ n+1 = h
nh+1 + nh 2
.
However, we should solve a nonlinear system at each time level and this progress decreases the accuracy of time. We can choose an approximation of second order accuracy to remain the same order of accuracy in the truncation error. This can be done by exploiting extrapolation technique
3 2
1 2
ˆ n+1 = n − n−1 , h h h
3 2
1 2
ˆ n+1 = n − n−1 , h h h
n ≥ 1.
Let denote hn+1 = 1nh+1 + i2nh+1 . We obtain the following full-discrete scheme with the extrapolated Crank–Nicolson method: Algorithm 3.1. (Coupled Crank–Nicolson (CCN) scheme) Find (hn+1 , nh+1 ) ∈ Vh × Vh , such that, for all 1 ≤ n ≤ M − 1, ∀ (v, s ) ∈ Vh × Vh ,
n+1 n+1 ⎧ n+1 h − hn ∂ ¯ h ∂ ¯ h ∂v ∂v ⎪ ˆ n+1 |2 ¯ n+1 , v ) + δ ( ¯ n+1 ˆ n+1 , v ) = 0, ⎪ ,v − α , −β , + γ ( | ⎨i h h hx h t ∂x ∂x ∂y ∂y n+1
n+1
¯ ¯ ∂ ∂ ⎪ ∂s ∂s ⎪ h h ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 , ∂ s = 0. ⎩ , + , +θ 1h 1h 2h 2h ∂x ∂x ∂y ∂y ∂x
(3.8)
Remark 3.1. The CCN scheme for DS problem is a linear scheme and satisfies the L2 norm conservation (2.3). Therefore, one only need to solve the coupled linear equations which can be solved very efficiently at each time level. However, we should solve a fully coupled system (3.8) at each time level which is time consuming in practical computation, especially in three dimensions. This motivates us to consider the following decoupled scheme for Eq. (3.2) by semi-implicit multistep scheme in time discretization. 3.3. Full discretization with semi-implicit multistep scheme In this paper, we consider 2-steps full discretization method for DS Eq. (2.1), which shall be referred to as two-step backward differentiation formula (BDF2). Thus, the decoupled numerical method can be written as follows: Algorithm 3.2. (Decoupled BDF2 (DBDF2) method) Find (hn+1 , nh+1 ) ∈ Vh × Vh , such that, for all 1 ≤ n ≤ M − 1, ∀ (v, s ) ∈ Vh × Vh ,
⎧ n+1 n+1 n+1
3h − 4hn + hn−1 ∂ h ∂ h ∂v ∂v ⎪ ⎪ ,v − α , −β , ⎪i ⎪ 2 t ∂x ∂x ∂y ∂y ⎪ ⎨ +γ (2|hn |2 hn − |hn−1 |2 hn−1 , v ) + δ ( nhx+1 (2hn − hn−1 ), v ) = 0, ⎪ n+1
n+1
⎪ ⎪ ∂ h ∂ h ∂s ∂s ∂s ⎪ ⎪ , + , + θ 2|hn |2 − |hn−1 |2 , = 0. ⎩ ∂x ∂x ∂y ∂y ∂x
(3.9)
Remark 3.2. The DBDF2 scheme for DS problem is totally decoupled, linear scheme and satisfies the L2 norm conservation (2.3). Therefore, the sequence linear algebraic system can be efficiently solved at each time level.
148
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
As both Eqs. (3.8) and (3.9) are used only for n ≥ 1, we must need start approximation (h1 , 1h ). A predictor-corrector
method is used for this purpose. We begin with an approximation (h1,0 , 1h,0 ), computed by the implicit-explicit Euler scheme,
1,0 1,0 ⎧ 1,0 h − h0 ∂ ¯ h ∂v ∂ ¯ h ∂v ⎪ ¯ 1,0 , v ) + δ ( ¯ 1,0 0 , v ) = 0 , ⎪ ,v − α , −β , + γ (|h0 |2 ⎨i h h hx t ∂x ∂x ∂y ∂y 1,0
1,0
¯ ¯ ∂ ∂ ⎪ ∂s ∂s ⎪ h h ¯ 1,0 + 0 ¯ 1,0 , ∂ s = 0 , ⎩ , + , − θ 10h 2h 2h ∂ x 1h ∂x ∂x ∂y ∂y
(3.10)
where
¯ h1,0 =
h1,0 + h0 2
¯ 1,0 = h
,
1h,0 + 0h 2
.
and correct it by
1 1
⎧ 1 h − h0 ∂ ¯ h ∂v ∂ ¯ h ∂v ⎪ ¯ 1,0 |2 ¯ 1 , v ) + δ ( ¯1 ¯ 1,0 , v ) = 0 , ⎪ ,v − α , −β , + γ ( | ⎨i h hx h t ∂x ∂x ∂y ∂y 1
1
¯ ¯ ∂ ∂ ⎪ ⎪ h ∂s h ∂s ¯ 1,0 ¯ 1 + ¯ 1,0 ¯ 1 , ∂ s = 0. ⎩ , + , −θ 1h 2h ∂ x 1h 2h ∂x ∂x ∂y ∂y
(3.11)
4. Error estimates In this section, we shall analyze the optimal error estimates for the full discrete scheme presented in Algorithm 3.1 with Crank–Nicolson method and Algorithm 3.2 with two-step backward differentiation formula. We introduce the projection : V → Vh , : V → Vh defined by, h h
∂χ ∂ α ( v − v ), ∂x h ∂x
∂χ ∂ +β ( v − v ), ∂y h ∂y
(∇ ( h v − v ), ∇χ ) = 0,
∀χ ∈ Vh , v ∈ V.
= 0,
(4.1)
∀χ ∈ Vh , v ∈ V.
(4.2)
r+1 v − vr+1 , ∀ v ∈ H r+1 () ∩ H01 (), h v + h (v − h v ) ≤ Ch
(4.3)
r+1 v − vr+1 , ∀ v ∈ H r+1 () ∩ H01 (). h v + h (v − h v ) ≤ Ch
(4.4)
Then, it may be shown that
Let denote n n hn − (tn ) = hn − h (tn ) + h (tn ) − (tn ) = e + E , n n nh − (tn ) = nh − h (tn ) + h (tn ) − (tn ) = e + E .
Theorem 4.1. Let {hn+1 , nh+1 }M satisfy the Eq. (3.8), the solutions ((t n+1 ), (t n+1 )) of Eq. (2.1) be sufficiently smooth and n=0
2 h0 − 0 2 + h1 − (t 1 )2 + t ( 0h − 0 2 + 1h − (t 1 )2 ) ≤ C hr+1 + t 2 ,
Then, for h and k sufficiently small, there exists a constant C, such that
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 .
0≤n≤M
(4.5)
(4.6)
where the constant C is independent of h and t. n , E n ) is bounded as claimed in (4.3), we only need to consider (en , en ). We consider the error for n ≥ 1. Proof. Since (E Using the projection property as presented in (4.1) and (4.2), and subtracting (3.8) from (3.1) results in the following error estimates in (en , en ), for n ≥ 1,
∂ n+1 ∂v ∂ n+1 ∂v −α −β + γ (J1n+1 , v ) + δ (J2n+1 , v ) e¯ , e¯ , ∂x ∂x ∂y ∂y
n+1 n E − E (t n+1 ) − (t n ) n+ 12 = −i − t (t )+ , v + (J3n+1 , v ), t t
en+1 − en i ,v t
(e¯n +1 , s ) + θ (J4n+1 , sx ) = −θ (J5n+1 , sx ) − (J6n+1 , s ),
(4.7) (4.8)
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
149
where
ˆ n+1 |2 ¯ n+1 − | ˆ n+1 |2 ¯ n+1 , J1n+1 = | h h 1
n+1 ¯ n+1 − (t n+ 2 ), J31 =
n+1 n+1 n+1 J3n+1 = α J31 − γ J32 − δ J33 ,
¯ nx +1 ¯ n+1 ˆ n+1 − ˆ n+1 , J2n+1 = hx h
n+1 ˆ n+1 |2 ¯ n+1 − |(t n+ 2 )|2 (t n+ 2 ), J32 = | 1
1
1
1
n+1 ¯ nx +1 ˆ n+1 − x (t n+ 2 )(t n+ 2 ), J33 =
ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 ) − ( ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 ), J4n+1 = ( 1 1 2 2 1h 1h 2h 2h 1
1
ˆ n+1 ¯ n+1 + ˆ n+1 ¯ n+1 ) − ( 2 (t n+ 2 ) + 2 (t n+ 2 )), J5n+1 = ( 1 2 1 1 2 2
1
¯ n+1 − (t n+ 2 ). J6n+1 =
Choosing s = e¯n +1 in (4.8), we obtain
e¯n +1 ≤ C (J4n+1 + J5n+1 + J6n+1 ). n+1
Choosing v = e¯
(4.9)
in (4.7), we get
i (en+1 2 − en 2 ) − αe¯n+1 2 + γ (J1n+1 , e¯n+1 ) + δ (J2n+1 , e¯n+1 ) 2 t
n E n+1 − E (t n+1 ) − (t n ) 1 = −i − t (t n+ 2 ) + , e¯n+1 + (J3n+1 , e¯n+1 ). t t
(4.10)
Using Taylor expansions, the first term on the right-hand of (4.10) satisfies
2 t n+1 (t n+1 ) − (t n ) n+ 12 3 − ( t ) ≤ C ( t ) t t t (s )2 ds. t t tn
(4.11)
Taking the imaginary part of Eq. (4.10), and using Eq. (4.3), we get
2 1 (en+1 2 − en 2 ) ≤C (J1n+1 2 + J2n+1 2 + J3n+1 2 + J4n+1 2 + J5n+1 2 + J6n+1 2 ) + C e¯n+1 2 +C hr+1 + t 2 . 2 t (4.12) Now, we estimate terms on the right-hand of Eq. (4.12). Exploiting inequality (4.3), we have
J1n+1 = (¯ hn+1 − ¯ n+1 )|ˆ hn+1 |2 + (|ˆ n+1 |2 − |ˆ hn+1 |2 )¯ n+1 ≤ C (e¯n+1 + eˆn+1 ) + Chr+1 ,
(4.13)
ˆ n+1 ) + ( ¯ nx +1 ) ¯ n+1 ( ˆ n+1 − ¯ n+1 − ˆ n+1 ≤ C (eˆn+1 + e¯n+1 ) + Chr+1 , J2n+1 = hx h hx
(4.14)
J4n+1 = ˆ 1nh+1 (¯ 1nh+1 − ¯ 1n+1 ) + ˆ 2nh+1 (¯ 2nh+1 − ¯ 2n+1 ) + (ˆ 1nh+1 − ˆ 1n+1 )¯ 1n+1 + (ˆ 2nh+1 − ˆ 2n+1 )¯ 2n+1 ≤ C (e¯n+1 + eˆn+1 ) + Chr+1 .
(4.15)
By Taylor formula, it is easy to see that, for j = 3, 5, 6
Jnj +1 2 ≤ C ( t )3
t n+1
tn
tt (s )2 ds.
(4.16)
Using Eqs. (4.13)–(4.16), we obtain
2 1 (en+1 2 − en 2 ) ≤ C (en+1 2 + en 2 + en−1 2 ) + C hr+1 + t 2 , t
(4.17)
There exists a constant k0 , for t ∈ [0, k0 ], such that 1 − C t > 0, thus, Eq. (4.16) can be written as
en+1 2 ≤
2 1 + C t n 2 e + C t en−1 2 + C t hr+1 + t 2 , 1 − C t
namely,
e + C t e ≤ n+1 2
n
2
(4.18)
2 1 + C t + C t (en 2 + en−1 2 ) + C t hr+1 + t 2 . 1 − C t
(4.19)
For n t ≤ T, n ≥ 1, we obtain
2 en+1 2 ≤ C e1 2 + t e0 2 + C hr+1 + t 2 .
(4.20)
Combining Eqs. (4.9), (4.20) and (4.5), we get
2 en+1 2 + t en +1 2 ≤ C hr+1 + t 2 .
(4.21)
150
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
From Eqs. (4.3) and (4.21), we obtain by triangle inequality
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 ,
0≤n≤M
We complete the proof of Theorem 4.1. Theorem 4.2. Let
{hn+1 , nh+1 }M n=0
satisfy the Eq. (3.9), the solutions ((t n+1 ), (t n+1 )) of Eq. (2.1) be sufficiently smooth and
2 h0 − 0 2 + h1 − (t 1 )2 + t ( 0h − 0 2 + 1h − (t 1 )2 ) ≤ C hr+1 + t 2 . Then, for h and k sufficiently small, there exists a constant C, such that
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 ,
(4.22)
(4.23)
0≤n≤M
where the constant C is independent of h and t. Proof. Subtracting Eq. (3.8) from Eq. (3.1) by the projection definition (4.1) and (4.2), we obtain the following error estimates in (en , en ), for n ≥ 1,
∂ n+1 ∂v ∂ n+1 ∂v e , −β e , + γ (J1n+1 , v ) + δ (J2n+1 , v ) ∂x ∂x ∂y ∂y
n+1
n−1 n 3E − 4E + E 3(t n+1 ) − 4(t n ) + (t n−1 ) 1 = −i − t (t n+ 2 ), v − i , v + (J3n+1 , v ), 2 t 2 t
i
3en+1 − 4en + en−1 ,v 2 t
−α
(en +1 , s ) + θ (J4n+1 , sx ) = (J5n+1 , sx ),
(4.24) (4.25)
where
J1n+1 = (2|hn |2 hn − |hn−1 |2 hn−1 ) − (2|(t n )|2 (t n ) − |(t n−1 )|2 (t n−1 )), J2n+1 = nhx+1 (2hn − hn−1 ) − x (t n+1 )(2(t n ) − (t n+1 )), J3n+1 = −γ [(2|(t n )|2 (t n ) − |(t n−1 )|2 (t n−1 )) − |(t n+1 )|2 (t n+1 )] −δ [ x (t n+1 )(2(t n ) − (t n−1 )) − x (t n+1 )(t n+1 )], J4n+1
= (2|hn |2 − |hn−1 |2 ) − (2|(t n )|2 − |(t n−1 )|2 ),
J5n+1 = θ (2|(t n )|2 − |(t n−1 )|2 − |(t n+1 )|2 ).
Choosing s = en +1 in Eq. (4.25), we obtain
en +1 ≤ C (J4n+1 + J5n+1 ).
(4.26)
Choosing c = en+1 in Eq. (4.24), we obtain
i t
n
j=1
= −i t
j+1 3e − 4en + en−1 j+1 , e 2 t
n j=1
−γ t
n j=1
− t
n
α (en+1 , ej+1 )
j=1
n−1 n 3E n+1 − 4E + E 3(t j+1 ) − 4(t j ) + (t j−1 ) j+1 − t (t n+1 ) + , e 2 t 2 t
(J1j+1 , ej+1 ) − δ t
n
(J2j+1 , ej+1 ) + t
j=1
n
(J3n+1 , ej+1 ).
(4.27)
j=1
First, we estimate the first term on the left-hand of Eq. (4.27)
t
n j=1
j+1 j j−1 3e − 4e + e j+1 , e 2 t
=
n
ej+1 2 − ej 2 −
j=1
3 ≥ en+1 2 − 4 3 ≥ en+1 2 − 4
1 n 2 e − 4 1 n 2 e − 4
1 j+1 2 1 e − ej−1 2 + ej+1 −ej 2 − ej+1 −ej−1 2 4 4 3 1 2 e + 4 3 1 2 e − 4
1 0 2 1 n+1 e + e − en 2 − e1 − e0 2 4 2 1 0 2 e . (4.28) 4
Next, using Taylor expansions and property (4.3), the right-hand side of Eqs. (4.25) and (4.27) have the following estimates, for j ≥ 1
tj 3(t j+1 ) − 4(t j ) + (t j−1 ) n+1 − ( t ) ≤ C t t t t ds, t 2 t t j−2
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
151
j+1 3E − 4Ej + Ej−1 1 = ( h − I )(3(t j+1 ) − 4(t j ) + (t j−1 )) t 2 2 t t j+1 t j+1 r+1 ≤ C ( − I ) ds ≤ Ch t r+1 ds, t h t j−1 t j−1 and n
Jlj+1 ≤ C t, l = 3, 5, Jlj+1 ≤ C (ej + ej−1 + hr+1 ), l = 1, 4,
j=1
J2j+1 ≤ C (ej + ej−1 + e j+1 ) + Chr+1 . Combining the above estimates, taking the imaginary part of Eq. (4.27), we obtain
en+1 2 ≤
4 1 n 2 e + e1 2 + e0 2 + t (J1j+1 + J2j+1 + J3j+1 + J4j+1 + J5j+1 )ej+1 . 3 3 n
(4.29)
j=1
Suppose m be chosen so that em = max0≤ j≤n e . Then it yields j+1
1 0 r+1 en+1 ≤ em + t 2 ). ≤ C e + e + C ( h
(4.30)
Hence, from Eqs. (4.3), (4.26), (4.22) and (4.30), we get by the triangle inequality
2
max (hn+1 − (t n+1 )2 + t nh+1 − (t n+1 )2 ) ≤ C hr+1 + t 2 .
0≤n≤M
We complete the proof of Theorem 4.2.
5. Basis functions and implementation issues Let {ϕ j }2j=0 be the linear basis functions on a triangle element e , namely, ϕ j ∈ P1 (e ), j = 0, 1, 2, Vh |e = span{ϕ0 , ϕ1 , ϕ2 }. Assume (Un , Zn ) be the approximation of the exact solutions ( n , n ) = ((t n ), (t n )). To overcome the deficiency of solving complex algebraical system (3.8), we decompose the complex value function into real and imaginary parts, namely, (tn ) = u1 (tn ) + iu2 (tn ), and approximate u j (tn ) = unj by U nj ∈ Vh , namely, hn = U1n + iU2n . We denote Uj and Z on an element by 2
U je =
e U jm (t )ϕm (x, y ),
Ze =
m=0
2
e Zm (t )ϕm (x, y ).
m=0
Then, the CCN scheme (3.8), namely, Algorithm 3.1 has the following matrix form on an element:
Me
X e,n+1 − X e,n + Ae (Xˆ e,n+1 )X e,n+1 + Be X e,n+1 = 0, t
(5.1)
where X e,n = (U e,n , Z e,n ) = (U1e,n , U2e,n , Z e,n ) , the element matrices are, for i, s = 0, 1, 2,
˜e = M is Ae (Uˆ j
D (Uˆ j e
e,n
e,n
)is = )is =
A˜ e (Uˆ e,n )is =
ϕi ϕs dx dy, Be1,is = e
2
e
e
e
˜e M 0
M˜ e e
M =
2
e
ϕix ϕsx dx dy, Be2,is =
ϕi ϕs dx dy,
e
ϕiy ϕsy dx dy,
j = 1, 2
e,n Uˆ jm
ϕm ϕix ϕsx dx dy,
j = 1, 2
m=0
⎛ 2 2 ⎞ 2 2 ⎝ ⎠ϕi ϕs dx dy, Uˆ1e,n ϕ + Uˆ2e,n ϕ m m m m
Ee =
e,n Uˆ jm ϕm
m=0
0 0
m=0
0 , ˜e −M 0 ˜e −M 0
m=0
F e (Uˆ e,n ) = 0 0 , 0
0 γ A˜ e (Uˆ e,n )
A (U e
ˆ e,n
)=
γ A˜ e (Uˆ e,n ) 0
0 e ˆ e,n ˜ γ A (U ) e,n −2θ De (Uˆ1 )
,
Ge =
γ A˜ e (Uˆ e,n ) 0 e,n −2θ De (Uˆ2 )
0 −α Be1 − β Be2 e,n δ Ae (Uˆ2 ) e,n δ Ae (Uˆ1 ) ,
0
−α Be1 − β Be2 , 0
152
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
e
B =
−α
Be1
0 − β Be2 0
−α Be1 − β Be2 0 0
0 0 . −Be1 − Be2
˜ e , Be , Be , Ae , De and A˜ e is 3 × 3. Assembling contributions from all elements leads to the The size of element matrices M 1 2 following matrix form of the linear algebraic equations, from Eq. (5.1),
M
X n+1 − X n + A(Xˆ n+1 )X n+1 + BX n+1 = 0, t
(5.2)
for n ≥ 0, where Xˆ n+1 = 32 Xˆ n − 12 Xˆ n−1 and X = (U, Z ) contains all the nodal parameters. Similarly, the matrix form on an element e of the DBDF2 method (3.9), namely, Algorithm 3.2 can be written as follows:
⎧ ⎨ 3U n+1 − 4U n + U n−1 E + E (U n )U n+1 + GU n+1 = bn , 2 t ⎩(−B − B )Z n+1 = 2θ DU n+1 U n+1 + DU n+1 U n+1 , 1 2 1 1 2 2
(5.3)
where be,n = (−δ Ae (U2 e,n )Z n , −δ Ae (U2 e,n )Z n ) . ˜ , A, B, E, F and F with the total domain are sparse block matrices, the linear As mentioned above, the matrices M, M algebraic system (5.2) and (5.3) can be respectively implemented efficiently. 6. Numerical experiments In this section, four experiences: the propagation of elliptic-elliptic and hyperbolic-elliptic DS equations are presented to illustrate the accuracy and L2 norm conservation laws of the numerical schemes based on CCN method and DBDF2 method. Moreover, we extend the proposed numerical methods to study the dynamic of rogue waves including line rogue wave for elliptic-elliptic and hyperbolic–elliptic DS equations, two-rogue waves for elliptic–elliptic DS equations. To illustrate the accuracy of the numerical methods, we denote the error norms between numerical solution Uhn and accuracy solution U(tn ) as follows, for j = 1, 2, U = ,
1 Uhn − U (t n )2 dx dy 2 , L∞ = Uhn − U (t n )∞ = max Uhn − U (t n ), L2 = Uhn − U (t n )2 =
H1 =
(x,y )∈
1 n n − Ux (t n )22 + Uhy − Uy (t n )22 2 , Uhx
(6.1)
where the nine points Gauss numerical quadrature formula for the integration is used. 6.1. Convergence and accuracy Example 1. Elliptic–elliptic DS equations The exact solutions of Eq. (2.1) are taken as [36]
⎧ ⎪ ⎨(x, y, t ) = 2 2
1 exp(iη ), 3 coth(ξ )
⎪ ⎩ (x, y, t ) = 4 (ξ − tanh(ξ )),
(6.2)
3
where
ξ = x + y − 4t, η = x + y − 6t. We choose the parameter α = 1, β = −1, γ = −1, δ = −1 and θ = −1 with the domain [−4, 4]2 . The evolution of traveling wave is done up to t = 1. To illustrate the accuracy of our algorithms, we compute the pointwise convergence rate in space. Setting the temporal h 1 1 1 step size t = 16 , we list the L2 , L∞ and H1 error accuracy of space in Tables 1 and 2 with varying spacing h = 10 , 20 , 40 at time t = 1. From the second and fourth columns in Table 1 based on CCN method and Table 2 based on DBDF2 method, the accuracy of the L2 and L∞ norms for and are closed to second order in space and the H1 convergence rate is O(h) as the sixth column shown. Therefore, the numerical experiment strongly suggests that the orders of convergence in space in error estimates for the L2 , L∞ and H1 norms of and are optimal and confirm with the theoretical results. To investigate the order of convergence in time, we take a more accurate, but slightly more sophisticated, approach to examine the orders of convergence with respect to the time step t because the approximation errors O( t γ ) + O(hμ ) are dominated by the h-terms O(hμ ) [37]. For instance, assuming
vh, t ≈ v(x, t m ) + C1 (x, t m ) t γ + C2 (x, t m )hμ ,
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161 Table 1 The order of convergence in space for error norms with t =
| |
by CCN method.
h
L∞
Order
L2
Order
H1
Order
1/10 1/20 1/40
1.257606E−1 3.089576E−2 7.673786E−3
2.0252 2.0093
3.310235E−1 8.054239E−2 1.987874E−2
2.0391 2.0185
8.686715E−1 3.506051E−1 1.639445E−1
1.3089 1.0966
1/10 1/20 1/40
2.025977E−2 5.158323E−3 1.300563E−3
1.9736 1.9877
4.412639E−2 1.157752E−2 2.952151E−3
1.9303 1.9715
1.439311E−1 6.069855E−2 2.869494E−2
1.2456 1.0808
Order
H1
Order 1.3053 1.1072 1.2504 1.0768
Table 2 The order of convergence in space for error norms with t =
| |
h 16
153
∞
h 16
by DBDF2 method.
Order
L2
1.337010E−1 3.243466E−2 8.109017E−3
2.0434 1.9999
3.475989E−1 8.491867E−2 2.099681E−2
2.0332 2.0159
8.852042E−1 3.581702E−1 1.662515E−1
1.935581E−2 5.015085E−3 1.248773E−3
1.9484 2.0057
4.571451E−2 1.202377E−2 3.069377E−3
1.9267 1.9698
1.437976E−1 6.044263E−2 2.865467E−2
h
L
1/10 1/20 1/40 1/10 1/20 1/40
it follows that
1 t vh, t − vh, 2 ≈ C1 (x, t m ) 1 − γ t γ . 2
Thus, t
ρv, t =
vh, t − vh, 2 4γ − 2γ ≈ γ = 2γ . t t 2 −1 vh, 2 − vh, 4
Here, v can be and . In particular, ρv, t approaches to 4 for γ = 2 when the corresponding order of convergence in 1 time is of O( t2 ). Fig. 1 shows the temporal accuracy in L2 , L∞ and H1 norms with fixed the spatial step size h = and 10 0.025 different temporal partition t = k , k = 0, 1, 2, 3. These results match the regular expectations of accuracy order O( t2 ) 2 arising from CCN method and DBDF2 method. Moreover, we compare CPU time needed for solving the DS problem (2.1) using the proposed CCN and DBDF2 numerical schemes and run simulations done up to t = 1. The CPU time as listed in Table 3 is given in seconds and includes the time cost by generating mesh structure, assembling the stiffness matrix and load vector, solving the corresponding algebraic system and computing errors on all the time steps. The coupled linear system is needed to solve for CCN method at each time, which causes that the DBDF2 numerical method needs much less CPU time than CCN method as presented in Table 3. Furthermore, we discuss numerical L2 norm conservation as presented in Eq. (2.3). The changes of the invariant from the initial variant approach zero by proposed Algorithms 3.1 and 3.2 as shown in Fig. 2. These results display that both the precision and the conservation are well in the proposed two numerical schemes. Example 2. Hyperbolic–elliptic DS equations Consider the following exact solutions in [24]
⎧ 2 2 2 2 ⎪ ⎪(x, y, t ) = 2(B1 + B2 )(ω + α k1 + β k2 ) sech(ξ ) exp(iη ), ⎨ 2 2 B1 (γ + δθ ) + γ B2 2 2 ⎪ −2 θ B 1 ( ω + α k1 + β k2 ) ⎪ ⎩ (x, y, t ) = − tanh(ξ ), B21 (γ + δθ ) + γ B22
(6.3)
where
ξ = B1 x + B2 y − ν t, η = −k1 x − k2 y − ωt, ν = −2(α k1 B1 + β k2 B2 ). Choosing parameters α = β = 12 , γ = −1, δ = 2 and θ = −2 in DS system (2.1), and B1 = B2 = 1, k1 = −0.1, k2 = −0.05 and ω = 0.05 in system (6.3) with the domain [−10, 10]2 and terminated time t = 4 of the evolution. The L2 , L∞ and H1 error accuracy of space are reported in Tables 4 and 5 with varying spacing h = 12 , 14 , 18 at time t = 4, which are almost completely closed to spatial order O(h2 ), O(h2 ) and O(h), respectively, and confirm with the theoretical results. Fig. 3 shows the temporal order of error norms can reach the theoretical result O( t2 ), which illustrates the efficiency and accuracy of the presented algorithms in this paper. In order to compare the proposed CCN and DBDF2 methods, we choose step sizes h = 16 and t = 0.1. Tables 6 illustrates that both CCN method and DBDF2 method have higher accuracy than the Lie-Trotter time-splitting (LTTS) scheme in [31],
154
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 1. Error norms at time t = 1.0 with different temporal step size t =
0.025 ,k 2k
= 0, 1, 2, 3 and the fixed mesh h =
1 . 10
due to the fact that the CCN and DBDF2 methods can derive the order convergence in time O( t2 ) while LTTS method only has first-order temporal convergence. Moreover, our numerical methods for both CCN method and DBDF2 method conserve the L2 norm (2.3) as presented in Fig. 4. The accuracy, efficiency and conservation of our proposed two numerical methods are verified again.
6.2. Application Example 3. Dynamics of fundamental rogue waves In this example, we consider the dynamics of the one line rogue waves for hyperbolic–elliptic and elliptic–elliptic DS Eq. (2.2) as presented in [26,27]. In the first case, we study the elliptic-elliptic type equation. Taking parameters α = 1,
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161 Table 3 The CPU time of different numerical methods with t =
155
h . 16
Method
h = 1/10
h = 1/20
h = 1/40
CCN DBDF2
203 186
2311 2144
30,576 28,107
Fig. 2. The L2 norm conservation M of numerical solutions by different numerical methods with h = Table 4 The order of convergence in space for error norms with t =
| |
∞
by CCN method.
Order
L2
8.768628E−4 2.223408E−4 5.579751E−5 4.424144E−4 1.111999E−4 2.806503E−5
h
L
1/2 1/4 1/8 1/2 1/4 1/8
Order
H1
Order
1.9795 1.9945
4.551711E−3 1.140348E−3 2.853134E−4
1.9969 1.9988
3.863076E−2 1.932430E−2 9.663248E−3
0.9933 0.9998
1.9819 1.9962
2.380935E−3 5.968542E−3 1.493584E−4
1.9961 1.9986
1.968664E−2 9.853213E−3 4.927838E−3
0.9986 0.9996
Table 5 The order of convergence in space for error norms with t =
| |
h 16
h
L∞
1/2 1/4 1/8 1/2 1/4 1/8
1 . 40
h 16
by DBDF2 method.
Order
L2
Order
H1
Order
8.749512E−4 2.220580E−4 5.575367E−5
1.9782 1.9938
4.547120E−3 1.139739E−3 2.852320E−4
1.9969 1.9985
3.863045E−2 1.932428E−2 9.663247E−3
0.9933 0.9998
4.389001E−4 1.112386E−4 2.789141E−5
1.9802 1.9957
2.332774E−3 5.850946E−3 1.464577E−4
1.9953 1.9981
1.968638E−2 9.853186E−3 4.927835E−03
0.9985 0.9996
β = −1, γ = 2, δ = 2 and θ = 1 in DS system (2.3), the initial conditions are given as √
(x, y, 0 ) = 2 1 −
√
4
1 + ( 3x − y )2
,
√
(x, y, 0 ) = −1 − 4 3
√ 1 − ( 3x − y ) 2 . √ [1 + ( 3x − y )2 ]2
(6.4)
with domain [0, 0.5] × [−10, 10] and time interval t ∈ [−1, 0.25]. The solutions describe a line wave with the line oriented √ in the ( 12 , 23 ) direction of the (x, y) plane, and the orientation angle is π3 . This fundamental rogue wave is illustrated in Fig. 5. At time t = 0 in Fig. 5(e), reaches the maximum amplitude at the center of the line.
156
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 3. Error norms at time t = 4.0 with different temporal step size t =
0.4 ,k 2k
= 0, 1, 2, 3 and the fixed mesh h =
1 . 8
Table 6 Comparison error norms of different methods at time t = 4. Method
CCN DBDF2 LTTS [31]
| |
L∞
L2
H1
L∞
L2
H1
2.095568E−4 2.056871E−4 6.535216E−4
1.016200E−3 1.006211E−3 4.653657E−3
1.306883E−2 1.306835E−2 1.379714E−2
9.637222E−5 9.824542E−5 3.218093E−4
4.727967E−4 4.955172E−4 1.820333E−3
6.554249E−3 6.5550 0 0E−3 6.785512E−3
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 4. The L2 norm conservation M of numerical solutions by different numerical methods with h =
157
1 . 6
Fig. 5. The profiles of numerical solutions | | at different time.
In the second case, the evolution of line rogue waves of hyperbolic-elliptic type DS equations is studied. Most of the parameters used in this simulation are the same as those in the first case, except β = γ = θ = −1. We take the following initial conditions
√ (x, y, 0 ) = 2 1 −
4 , 1 + ( k1 x + k2 y )2
(x, y, 0 ) = 1 − 4k21
1 − ( k1 x + k2 y )2 . [1 + (k1 x + k2 y )2 ]2
(6.5)
158
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 6. The profiles of numerical solutions | | at different time.
Fig. 7. The profiles of numerical solutions | | at different time.
with k1 = 56 and k2 = 13 6 . The computational domain is [−4, 4] × [0, 1] and time interval is t ∈ [−1, 1]. The solutions describe 5 a line rogue wave with the line oriented in the ( 13 6 , − 6 ) direction of the (x, y) plane, thus the fundamental rogue waves in the DS equations are line rogue waves. Fig. 6 shows the evolution of fundamental rogue wave. The maximum amplitude can be reached at time t = 0 in Fig. 6(d). Example 4. Dynamics of two-rogue waves In this example, we consider the dynamics of two-rogue waves in DS equations, which describes the interaction between two fundamental line rogue waves. Choosing parameters α = −1, β = 1, γ = 1, δ = 2, θ = 1 and m = 1 in DS system (2.2).
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
159
Fig. 8. The profiles of numerical solutions at different time.
Owing to the works in [27], the initial conditions can be given
√ τ 2 = 2 0 , = −1 − 2
τ1
τ1
2 ∂ 2 τ1 ∂τ1 τ1 2 − , ∂x ∂x
(6.6)
where τs , s = 0, 1 is the determinant of 2 × 2 block matrix with each block 2 × 2,
m (s ) ij τs = (−s) −m i j
and (s )
mi j
1 = pi + q j
ξi =
, mi(j−s ) i(js ) m
ξi + s −
(6.7)
pi + ci1 , pi + q j
pi + 1/pi pi − 1/pi x+ −1y. 2 2
i(js ) = m
1 pi q j − 1
ξi + s −
pi q j pi q j − 1
+ ci1 ,
√ √ with the parameters p1 = 1, p2 = −1, q1 = 0, q2 = −3, c11 = 0 and c21 = 2−1 on the domain [−10, 10]2 and time interval 1 t ∈ [−0.3, 0.6]. The corresponding numerical solutions with step size h = 30 and t = 0.001 is displayed in Figs. 7–10. It is seen that a cross-shape rogue wave appears, which describes the interaction between two fundamental line rogue waves: one oriented along the y-direction (corresponding to the parameters p1 ) and the other one oriented along the x-direction √ (corresponding to the parameters p2 ). These two individual line waves reach their peak amplitude 3 2 at different time with the x-direction one peaking at t ≈ − 41 as shown in Fig. 7(a) and the y-direction one peaking at t ≈ 0 in Fig. 7(c). This phenomenon is similar to the results as presented in [27], which validates the efficiency of our numerical scheme.
7. Conclusion In this paper, numerical algorithms for the numerical solutions of the Davey–Stewartson equations are proposed using Galerkin finite element with second order time discretization, which employs linearized Crank–Nicolson scheme and decoupled multistep scheme. The high efficiency and accuracy of our numerical methods are demonstrated by the numerical experiments: the propagation of elliptic-elliptic and hyperbolic-elliptic DS equations. The numerical results recorded in
160
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
Fig. 9. Contours of numerical solutions | | at different time.
Fig. 10. Contours of numerical solutions at different time.
Y. Gao et al. / Applied Mathematics and Computation 328 (2018) 144–161
161
Tables 1–4 and Figs. 1 and 3 show that the proposed numerical scheme have good accuracy. Moreover, our numerical methods satisfy L2 norm conservation laws as shown in Figs. 2 and 4. We extend the proposed numerical methods to study the dynamics of rogue waves, including fundamental rogue wave and interaction of two line rogue waves. What is more, the proposed methods are very efficient since their implementation requires at each time step the solutions of a linear system with the same matrix for all time levels. References [1] F. Linares, G. Ponce, Introduction to Nonlinear Dispersive Equations, Springer, New York, 2009. [2] J.A. Pava, Nonlinear Dispersive Equations: Existence and Stability of Solitary and Periodic Travelling Wave Solutions, American Mathematical Society, Providence, R.I., 2009. [3] Y.V. Kartashov, B.A. Malomed, L. Torner, Solitons in nonlinear lattices, Rev. Mod. Phys. 83 (1) (2011) 247–305. [4] R. Carles, E. Dumas, C. Sparber, Multiphase weakly nonlinear geometric optics for Schrödinger equations, SIAM J. Math. Anal. 42 (2010) 489–518. [5] R. Thomas, C. Kharif, M. Manna, A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity, Phys. Fluids 24 (12) (2012) 395–416. [6] L. Chang, Y. Pan, X. Ma, New exact travelling wave solutions of Davey–Stewartson equation, J. Comput. Inf. Syst. 9 (2013) 1687–1693. [7] M. Dehghan, M. Abbaszadeh, A. Mohebbi, Analysis of two methods based on Galerkin weak form for fractional diffusion-wave: meshless interpolating element free Galerkin (IEFG) and finite element methods, Eng. Anal. Bound Elem. 64 (4) (2016) 205–221. [8] M. Song, A. Biswas, Topological defects and bifurcation analysis of the DS equation with power law nonlinearity, Appl. Math. Inf. Sci. 9 (4) (2015) 1719–1724. [9] M. Dehghan, M. Abbaszadeh, A finite element method for the numerical solution of Rayleigh–Stokes problem for a heated generalized second grade fluid with fractional derivatives, Eng. Comput. Ger. 33 (3) (2017) 587–605. [10] M. Dehghan, A. Shokri, Numerical solution of the nonlinear Klein–Gordon equation using radial basis functions, J. Comput. Appl. Math. 230 (2) (2009) 400–410. [11] C. Klein, B. Muite, K. Roidot, Numerical study of blowup in the Davey–Stewartson system, Discr. Cont. Dyn. Syst. B 18 (2013) 1361–1387. [12] M. Dehghan, Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices, Math. Comput. Simul. 71 (1) (2006) 16–30. [13] C. Klein, K. Roidot, Fourth order time-stepping for Kadomtsev–Petviashvili and Davey–Stewartson equations, SIAM J. Sci. Comput. 33 (2011) 3333–3356. [14] M. Dehghan, F. Fakhar-Izadi, The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves, Math. Comput. Model. 53 (9) (2011) 1865–1877. [15] A. Davey, K. Stewartson, On three-dimensional packets of surface waves, Proc. R. Soc. Lond. A 338 (1974) 101–110. [16] C. Babaoglu, Long-wave short-wave resonance case for a generalized Davey–Stewartson system, Chaos Soliton Fract. 38 (2008) 48–54. [17] V.D. Djordjevic´ , L.G. Redekopp, On two-dimensional packets of capillary-gravity waves, J. Fluid Mech. 79 (1977) 703–714. [18] K. Nishinari, K. Abe, J. Satsuma, A new-type of soliton behavior of the Davey–Stewartson equations in a plasma system, Teoret. Mat. Fiz. 99 (1994) 487–498. [19] J.M. Ablowitz, R. Haberman, Nonlinear evolution equations in two and three dimensions, Phys. Rev. Lett. 35 (1975) 1185–1188. [20] A. Eden, T.B. Gürel, On the integrability of a generalized Davey–Stewartson system, Phys. D 259 (2013) 1–7. [21] H. Jafari, A. Sooraki, Y. Talebi, A. Biswas, The first integral method and traveling wave solutions to Davey–Stewartson equation, Nonlinear Anal. Model. Control 17 (2) (2012) 182–193. [22] G. Ebadi, A. Biswas, The G’/G method and 1-soliton solution of the Davey–Stewartson equation, Math. Comput. Model. 53 (2011) 694–698. [23] O.H. El-Kalaawy, R.S. Ibrahim, Solitary wave solution of the two-dimensional regularized long-wave and Davey–Stewartson equations in fluids and plasmas, Appl. Math. Ser. B 3 (8) (2012) 833–843. [24] G. Ebadi, E.V. Krishnan, M. Labidi, E. Zerrad, A. Biswas, Analytical and numerical solutions to the Davey–Stewartson equation with power-law nonlinearity, Wave Random Complex 21 (2011) 559–590. [25] J. Shi, J. Li, S. Li, Analytical travelling wave solutions and parameter analysis for the (2+1)-dimensional Davey–Stewartson-type equations, Pramana J. Phys. 81 (2013) 747–762. [26] V.A. Arkadiev, A.K. Pogrebkov, M.C. Polivanov, Inverse scattering transform method and soliton solutions for the Davey–Stewartson II equation, Phys. D 36 (1989) 189–196. [27] Y. Ohta, J. Yang, Dynamics of rogue waves in the Davey–Stewartson II equation, J. Phys. A Math. Theor. 46 (10) (2013) 105202. [28] Y. Ohta, J. Yang, Rogue waves in the Davey–Stewartson I equation, Phys. Rev. A Math. Theor. 86 (2012) 036604. [29] Y. Zhang, H.P. Stimming, N.J. Mauser, A novel nonlocal potential solver based on nonuniform FFT for efficient simulation of the Davey–Stewartson equations, ESAIM Math. Model. Num. 51 (4) (2017) 1527–1538. [30] C. Besse, N.J. Mauser, H.P. Stimming, Numerical study of the Davey–Stewartson system, Math. Model. Numer. Anal. 38 (2004) 1035–1054. [31] Y. Gao, L. Mei, R. Li, A time-splitting galerkin finite element method for the Davey–Stewartson equations, Comput. Phys. Commun. 197 (2015) 35–42. [32] J.M. Ghidaglia, J.C. Saut, On the initial value problem for the Davey–Stewartson systems, Nonlinearity 3 (1990) 475–506. [33] C. Klein, K. Roidot, Numerical study of the semiclassical limit of the Davey–Stewartson II equations, Nonlinearity 27 (2014) 2177–2214. [34] G.M. Muslu, H.A. Erbay, Numerical simulation of blow-up solutions for the generalized Davey–Stewartson system, Int. J. Comput. Math. 88 (2011) 805–815. [35] Y. Gao, L. Mei, Implicit-explicit multistep methods for general two-dimenisonal nonlinear Schrödinger equations, Appl. Numer. Math. 109 (2016) 41–60. [36] M.A. Bashir, A.A. Moussa, The cotha (ξ ) expansion method and its application to the Davey–Stewartson equation, Appl. Math. Sci. 8 (2014) 3851–3868. [37] L. Shan, H. Zheng, Partitioned time stepping method for fully evolutionary Stokes–Darcy flow with the Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 51 (2013) 813–839.