Energy-preserving finite volume element method ... - Semantic Scholar
Journal of Computational Physics 270 (2014) 58–69
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Energy-preserving finite volume element method for the improved Boussinesq equation Quanxiang Wang a,b , Zhiyue Zhang a,∗ , Xinhua Zhang b , Quanyong Zhu a,c a b c
Jiangsu Key Laboratory for NSLSCS, School of Mathematical Sciences, Nanjing Normal University, Nanjing 210023, China College of Engineering, Nanjing Agricultural University, Nanjing 210031, China School of Science, Lishui University, Lishui 323000, China
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 2 January 2014 Received in revised form 23 March 2014 Accepted 27 March 2014 Available online 2 April 2014
In this paper, we design an energy-preserving finite volume element scheme for solving the initial boundary problems of the improved Boussinesq equation. Theoretical analysis shows that the proposed numerical schemes can conserve the energy and mass. Numerical experiments are performed to illustrate the efficiency of the scheme and theoretical analysis. While the results demonstrate that the proposed finite volume element scheme is second-order accuracy in space and time. Moreover, the new scheme can conserve mass and energy. © 2014 Elsevier Inc. All rights reserved.
Keywords: Conservative Improved Boussinesq equation Energy Finite volume element
1. Introduction The Korteweg–de Vries (KdV) equation is often used to describe nonlinear dispersive wave phenomena, such as the propagation of long waves on the surface of shallow water, the propagation of long waves in one-dimensional nonlinear lattices, and the propagation of ion-acoustic waves. In the past decades, much attention has been paid to investigate the KdV equation theoretically and numerically by many authors [1–4]. Boussinesq equation was first presented by Joseph Boussinesq in 1873, it is a fourth-order nonlinear partial differential equation and belongs to Korteweg–de Vries equation (KdV) family. The Boussinesq equation which describes motions of long waves in shallow water under gravity propagating in both directions, is given by
utt = u xx + qu xxxx + u 2
xx
(1)
,
where u = u (x, t ) is a sufficiently differentiable function, which for q = −1 gives the good Boussinesq or well-posed equation, while for q = 1 the bad or ill-posed Boussinesq equation [5,6]. In this study, we will take q = 1. If the term u xxxx in Eq. (1) is replaced with u xxtt , it gives the so-called improved Boussinesq equation
utt = u xx + u xxtt + u 2
xx
,
x ∈ [a, b].
(2)
Researchers do not observe any instabilities for the above improved Boussinesq equation [7]. In addition, the equation has in common the mass and energy conservation property [8]:
*
Corresponding author. E-mail address: [email protected] (Z. Zhang).
http://dx.doi.org/10.1016/j.jcp.2014.03.053 0021-9991/© 2014 Elsevier Inc. All rights reserved.
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
dM dt
59
b = 0,
M=
u dx,
(3)
a
dE dt
b = 0,
E=
v2 2
+
(ut )2 2
+ G (u ) dx,
(4)
a 2
3
under some suitable boundary conditions, where ut = v x , G (u ) = u2 + u3 . The improved Boussinesq equation has been the subject of extensive investigations and some numerical schemes have been designed to study it. For example, El-Zoheiry [9] designed a three-level iterative scheme based on the compact implicit method for solving the improved Boussinesq equation. Bratsos [5] proposed a second order implicit finite difference scheme by rational approximations. Recently, Bratsos [6] constructed a predictor–corrector scheme to avoid the solution of a nonlinear system. Lin et al. [10] developed a finite element scheme for solving initial boundary value problems for the improved Boussinesq equation. Irk and Da˘g [11] proposed two numerical methods based on classical five point finite difference approximation. Zhang and Lu [12] designed quadratic finite volume element schemes for the equation. Then, Wang et al. [13] numerically studied stochastic damped improved Boussinesq equation using the finite volume element schemes. However, these numerical methods are not concerned about the conservation properties of numerical schemes. The finite volume element method has been one of the most commonly used numerical methods for solving partial differential equations. One main attractive property of the method is that, the mass conservation law is maintained, which is fairly desirable for fluid and underground fluid computations. So it has been extensively used in computational fluid dynamics [14–17]. In this study, we propose an energy-preserving finite volume element method for the improved Boussinesq equation. The method conserves mass M and energy E of the improved Boussinesq equation. The numerical schemes are constructed using the discrete variational derivative method (DVDM). Furihata has used the DVDM to design a stable finite difference scheme for the Cahn–Hilliard equation [18]. Further, Matsuo et al. [19,20] have extended the method to solve nonlinear Klein–Gordon equation and Degasperis–Procesi equation. The organization of the paper is the following. In Section 2, preliminaries and notations are introduced. In Section 3, the proposed energy-preserving finite volume element schemes are presented. We provide some numerical experiments to illustrate the effectiveness and accuracy of the proposed schemes in Section 4. Conclusions are drawn in Section 5. 2. Preliminaries and notations In this section, we introduce notations and useful propositions employed in this paper. 2.1. Trial and test function spaces The region of interest [a, b] can be decomposed into a grid T h with nodes
a = x0 < x1 < x2 < · · · < x N −1 < x N = b, where xi − xi −1 = h and h = nodes
b−a . N
(5)
Denote T h = { I i : I i = [xi −1 , xi ], i = 1, 2, . . . , N }. Accordingly we place a dual grid T h∗ with
a = x0 < x1/2 < x3/2 < · · · < x N −3/2 < x N −1/2 < x N = b, where xi −1/2 = xi − h/2 (i = 1, 2, . . . , N ). Denote I 0∗
= [x0 , x1/2 ], I i∗
(6)
= [xi −1/2 , xi +1/2 ] and I ∗N
= [xN −1/2 , xN ]. Select the trial function space U h as the linear element space with respect to T h . The basis function of the node xi is
φi (x) =
1 − |x − xi |/h, 0,
x i −1 x x i +1 , elsewhere.
Then numerical solution U for Eq. (2) can be uniquely written as U = element I i , we have
U = U i −1 (1 − μ) + U i μ,
U = (U i − U i −1 )/h,
(7)
N
i =1
U i φi (x), where U i = U (xi , t ). So in the
(8) (9)
where μ = (x − xi −1 )/h. The test function space V h corresponding to T h∗ is taken as the piecewise constant function space. The test function of the nodes x j is
ψ j (x) =
1, x j −1/2 x x j +1/2 , 0, elsewhere.
(10)
60
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
The numerical solution is denoted by U kn u (kh, nt ), where t is the time mesh size. For the improved Boussinesq equation, we adopt the following periodic boundary conditions
∂ j u ∂ j u = ( j = 0, 1, 2). ∂ x j x=a ∂ x j x=b
(11)
So the discrete periodic boundary conditions are
U kn = U kn mod N
for ∀k ∈ Z.
(12)
2.2. Conserved quantities In order to see the conservation property, it is more convenient to work with the following system representation
ut = v x , v t − v xxt = ∂∂x δδGu ,
(13)
where δδGu is the Euler–Lagrange variational derivative defined by
∂G ∂ δG δG . = − δu ∂u ∂ x δux
So the numerical solution V for Eq. (13) can be uniquely written as V = Proposition 1. Let u be the solution of Eq. (13). Then the total mass
d
b a
N
i =1
V i φi (x), where V i = V (xi , t ).
u dx is independent of t. Namely,
b u dx = 0.
dt
(14)
a
The conservation of mass can be easily proved [21] and the standard finite volume element schemes can conserve the mass of IBq [13]. Proposition 2. Let u be the solution of Eq. (13). Then the total energy
d
b
dt
b a
G˜ (u , ut , v ) dx is independent of t. Namely,
G˜ (u , ut , v ) dx = 0,
(15)
a
where
v2
G˜ (u , ut , v ) =
2
+
(ut )2 2
+ G (u ).
Proof. According to the chain rule, we have
d
b
dt
G˜ dx =
a
b v v t + ut utt + a
δG ut dx δu
By the first equality of Eq. (13), we can get
b = a
δG v v t + v x v xt + v x dx δu
Using the integration by part formula, we have
b =
v v t − v xxt − a
∂ δG dx ∂ x δu
At last, we can complete the proof using the second equality of Eq. (13).
2
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
61
3. Design of schemes 3.1. Discrete variational derivative method In the following, we recall the discrete variational derivative method briefly. Firstly, we define the energy as follows:
b J (u ) =
G˜ (u , ut , v ) dx,
(16)
a
where we call G˜ the energy function throughout this paper. As the discretization of energy function, we employ the following finite volume element approximation
G˜ d (U , U t , V ) G˜ .
(17)
Then we define the total discrete energy as follows:
b J (u ) =
G˜d (U , U t , V ) dx.
(18)
a
˜ ˜ To obtain the approximation δ(Uδ G, V ) of Euler–Lagrange derivative δδGu , we need to compute the following difference:
b J (U ) − J ( V ) = a
δ G˜ (U − V ) dx. δ(U , V )
(19)
3.2. Energy-preserving finite volume element scheme Recall the definition G˜ (u , ut , v ) =
δ G δ G˜ δ u , δ ut
˜ = ut and δδGv
v2 2
+
(ut )2 2
+ G (u ). According the definition to Euler–Lagrange derivative, we have
δ G˜ δu
=
= v. So Eq. (13) can be rewritten as
δ G˜
ut = ∂∂x δ v , ˜ (1 − ∂x2 ) v t = ∂∂x δδGu ,
(20)
2 where ∂x2 = ∂∂x2 .
Using any ψ ∈ V h and ψ˜ ∈ V h to multiply both sides of Eq. (20) and integrating on the interval [a, b], we can obtain the following semi-discrete scheme
,ψ , ˜ = ∂ δ G˜ , ψ˜ , ((1 − ∂x2 ) v t , ψ) ∂ x δu (ut , ψ) =
∂ δ G˜ ∂x δv
(21)
where (·, ·) denotes the inner product of L 2 ([a, b]). Now we discretize the above semi-discrete scheme in time. To simplify the notation, we will use the difference operator 1
δn un = δ G˜ δv
un+1 −un−1 . 2t
˜ Additionally, we define the discrete version of Euler–Lagrange derivative δδGu =
u2 2
+
˜ u3 , δδuG 3 t
= ut and
= v by
δ G˜ d u n ( u n +1 + u n + u n −1 ) n u + , = 3 δ(un+1 , un , un−1 ) δ G˜ d 1 = δn un , 1 n δ(δn u ) δ G˜ d = vn. δ( v n )
(22) (23)
(24)
At last, we can get the following energy-preserving finite volume element scheme
⎧ ⎨ (δ 1 un , ψ) =
δ G˜ d ∂ ∂ x δ( v n ) , ψ , δ G˜ d ⎩ ((1 − ∂ 2 )δn1 v n , ψ) ˜ ˜ = ∂ x ∂ x δ(un+1 ,un ,un−1 ) , ψ . n
Eq. (25) has the following discrete mass and energy conservation law.
(25)
62
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Theorem 1 (Discrete mass conservation law). Under the discrete periodic boundary condition, then numerical solution by Eq. (25) conserves mass
b un dx = const.
(26)
a
Theorem 2 (Discrete energy conservation law). Under the discrete periodic boundary condition, then numerical solution by Eq. (25) conserves energy
b
G˜ nd dx = const.
(27)
a
Proof.
b
1
t
˜ 1 1 δn un+1 − δn un−1 δ G˜ d δ G d 1 n δ G˜ d 1 n + + , u v , δ , δ n n 1 2t δ( v n ) δ(un+1 , un , un−1 ) δ(δn un ) ˜ ˜ n+1 ˜ ∂ δ G˜ d δGd δGd δ G d 1 n ∂ v x − v nx −1 + , = + v , , δ n ∂ x δ( v n ) 2t δ( v n ) δ(un+1 , un , un−1 ) ∂ x δ( v n ) ˜ δ G˜ d δ G˜ d ∂ 2 1 n δ G d + δ 1 − ∂ v , , =− n x ∂ x δ(un+1 , un , un−1 ) δ( v n ) δ( v n ) = 0. 2
G˜ nd+1 − G˜ nd dx =
a
3.3. Finite volume element scheme For comparison, we employ the standard finite volume element scheme constructed based on the
u = v , t x (1 − ∂x2 ) v t = (u + u 2 )x .
(28)
Then we can obtain the following finite volume element scheme
1
(δn un , ψ) = ( v nx , ψ), 1 ˜ = ((un + (un )2 )x , ψ). ˜ ((1 − ∂x2 )δn v n , ψ)
(29)
Theorem 3 (Discrete mass conservation law). Under the discrete periodic boundary condition, then numerical solution by Eq. (29) conserves mass
b un dx = const.
(30)
a
4. Numerical experiments In this section, we present some numerical examples using the proposed energy-preserving finite volume element method. To show the efficiency of the present method, we calculate L 1 , L 2 and L ∞ error norms defined as follows:
L1 =
N uh (xi ) − u (xi ) h,
(31)
i =1
N uh (xi ) − u (xi ) 2 h, L2 = i =1
L ∞ = max u h (xi ) − u (xi ) . i
Define the rate of convergence r:
r = log2
u 2h − u
,
u h − u
(32) (33)
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
63
Table 1 The rate of convergence in space for the EFVEM. h
L ∞ error
r
L 2 error
r
L 1 error
r
0.40 0.20 0.10 0.05
9.78e−04 2.46e−04 6.23e−05 1.58e−05
– 1.99 1.98 1.98
1.95e−03 4.92e−04 1.24e−04 3.18e−05
– 1.98 1.98 1.96
6.35e−03 1.59e−03 4.04e−04 1.03e−04
– 2.00 1.98 1.97
Table 2 The rate of convergence in time for the EFVEM.
t
L ∞ error
r
L 2 error
r
L 1 error
r
0.040 0.020 0.010 0.005
1.42e−03 2.99e−04 7.47e−05 1.87e−05
– 2.24 2.00 2.00
2.53e−03 6.15e−04 1.54e−04 3.84e−05
– 2.04 2.00 2.00
8.24e−03 2.02e−03 5.06e−04 1.26e−04
– 2.03 2.00 2.01
where u h is the numerical solution with space step size h, u the analytical solution, N the number of node. The rate approaches the number 2 would indicate second-order accuracy in space and time. 4.1. Single solitary wave Eq. (2) has a special solution [22] which describes solitary wave:
u (x, t ) = α sech
2
α x − x0 − β t 6
β
,
β=
1+
2α 3
(34)
,
where α is amplitude of the solitary wave which is initially centered at x = x0 traveling with velocity β . Differentiating the above equation with respect to time and then setting t = 0, the initial conditions are obtained
u (x, 0) = α sech2
ut (x, 0) = 2α
α 6
α x − x0 6
sech2
β
α x − x0 6
(35)
,
β
tanh
α x − x0 6
β
.
(36)
To compare our results with those given in other references, we will take x0 = 0 and α = 0.5 in the following. In order to demonstrate the accuracy in space and time, we discretize the problem over the region [−60, 100] up to time t = 1. Firstly, we choose a relatively small time step t = 0.001 such that the error from the time discretization is negligible compared with those from the space. The L 1 , L 2 and L ∞ error norms for u have been computed by the energy-preserving finite volume element method (EFVEM). They are displayed in Table 1. Examination of the table shows that the error measures of the EFVEM diminish approximately quadratically as the space step is halved. In other words, the proposed energy-preserving scheme is second-order in space. To investigate the accuracy in time, the time and space step are chosen such that the ratio t /h = 0.1. We list the errors of the numerical solution computed by our new scheme in Table 2 in the L 1 , L 2 and L ∞ error norms. Second-order convergence can be observed. We perform a simulation with space step h = 0.2 and time step t = 0.1 for single solitary wave on the time interval [0, 40]. Fig. 1 presents that the solitary wave propagates towards the right without change in form. It is easy to find that the solitary wave continues to evolve with the initial velocity. We also perform a simulation with the resolution to check the conservation of schemes. Fig. 2 shows the error in the discrete mass for standard finite volume element scheme (FVEM) and energy-preserving finite volume element scheme, which well agree with the discrete mass conservation law. Fig. 3 presents the error in the discrete energy for FVEM and EFVEM. From the figure, we can find that standard finite volume scheme can conserve energy within 10−3 , which though very good is not strict in numerical sense. On the contrary, the proposed energy-preserving finite volume element scheme conserves energy within 10−14 , which is strict conservation in numerical sense and conforms the theoretical result. To make comparison with the earlier works [5,11,22], the numerical solution is given at t = 72 with α = 0.5 and x ∈ [−80, 100]. The results are presented in Table 3. We can find that our numerical method is much more efficient than those presented in the references. For example, L ∞ error = 4.1403e−04 at t = 72 is achieved with t = 0.1. Fig. 4 shows the single solitary wave at t = 72. Contrarily to [22], the amplitude of solitary wave does not decrease and there is no visible secondary solitary wave produced, which is consistent with those presented in [5]. 4.2. Two solitary waves In this example, we study the interaction of two solitary waves moving on a collision course. Here, x ∈ [−60, 100] and the initial conditions are
64
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 1. The propagation for single solitary wave.
α = 0.5, h = 0.2 and t = 0.1.
Fig. 2. (a) Error in the discrete mass for FVEM; (b) Error in the discrete mass for EFVEM. Single solitary wave.
Fig. 3. (a) Error in the discrete energy for FVEM; (b) Error in the discrete energy for EFVEM. Single solitary wave.
Table 3 Numerical results with
EFVEM M1 in [11] [5]
α = 0.5 at t = 72. h
t
L ∞ error
Energy error
0.1 0.1 0.1
0 .1 0.001 0.001
4.1403e−04 1.0064e−03 4.1998e−04
7.9492e−14
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 4. The propagation for single solitary wave.
α = 0.5, h = t = 0.1.
Fig. 5. The propagation for two solitary waves (a) and level curves (b).
2
u (x, 0) = α1 sech
ut (x, 0) = 2α1
α1 x − x01 6
α1 6
+ 2α2
6
+ α2 sech
β1
α1 x − x01
sech2
α2
6
2
sech
√
α2 x − x02 β2
α2 x − x02
β1 6
2
tanh
6
β2
6
tanh
α2 x − x02 β2
(37)
,
β1 6
α1 = α2 = 0.4, h = t = 0.1.
α1 x − x01
65
,
(38)
√
where x01 = −20, x02 = 60, β1 = 1 + 2α1 /3, β2 = − 1 + 2α2 /3. The first solitary wave is initially located at x = x01 and moving to the right with speed β1 and amplitude α1 . The other one is initially located at x = x02 and moving to the left with speed |β2 | and amplitude α2 . Firstly, we study the collision of two solitary waves with the same amplitude. Fig. 5 presents the propagation for two solitary waves with α1 = α2 = 0.4. From the figure, we can find that the solitary waves preserve their amplitudes and no displacement is observed. Additionally, there is no visible secondary solitary waves produced. Therefore, we can conclude that the collision is elastic. Fig. 6 presents the propagation for two solitary waves with α1 = α2 = 0.8. It can be observed from the figure that the secondary solitary waves are visible. Hence, the interaction is inelastic. The results are in good agreement with those reported in [5,10]. Now, we give some results for the interaction of solitary waves of different amplitudes. Fig. 7 shows interaction of two solitary waves with α1 = 0.8 and α2 = 0.4. Fig. 8 displays the collision of two solitary waves with α1 = 1.5 and α2 = 0.4. From the figures, we can find that the collisions are inelastic. According to definition used in [5,22] for the inelasticity coefficient, the coefficient K is defined as follows
K=
αm , max(α1 , α2 )
(39)
66
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 6. The propagation for two solitary waves (a) and level curves (b).
α1 = α2 = 0.8, h = t = 0.1.
Fig. 7. The propagation for two solitary waves (a) and level curves (b).
α1 = 0.8, α2 = 0.4, h = t = 0.1.
Fig. 8. The propagation for two solitary waves (a) and level curves (b).
α1 = 1.5, α2 = 0.4, h = t = 0.1.
where αm is the maximum joint amplitude at the collision time. The results for αm and K are presented in Table 4. Examination of the table shows that αm is always smaller than twice the amplitude of the max(α1 , α2 ). In addition, the coefficient K will decrease when α1 increase. The results are contrary to those presented in [23], Fig. 9 shows the error in the discrete mass for standard finite volume element scheme and energy-preserving finite volume element scheme, which well agree with the discrete mass conservation law. Fig. 10 shows the error in the discrete energy for FVEM and EFVEM. From the figure, we can make a similar conclusion. The energy-preserving finite volume
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Table 4 Numerical results for different
67
α1 and α2 at t = 80.
α1
α2
αm
K
0.4 0.8 0.8 1.5
0.4 0.8 0.4 0.4
0.7225 1.3895 1.0738 1.7188
1.8063 1.7369 1.3422 1.1458
Fig. 9. (a) Error in the discrete mass for FVEM; (b) Error in the discrete mass for EFVEM. Two solitary waves.
α1 = 0.8, α2 = 0.4, h = t = 0.1.
Fig. 10. (a) Error in the discrete energy for FVEM; (b) Error in the discrete energy for EFVEM. Two solitary waves.
α1 = 0.8, α2 = 0.4, h = t = 0.1.
element scheme can conserve the energy to machine precision. On the other hand, the standard finite volume element scheme can only conserve energy within 10−3 . 4.3. Wave break-up The initial conditions are given in the following:
u (x, 0) = α sech where
2
α x − x0 6
β
α = 1.0, x0 = 20 and β = 1 +
, 2α . 3
ut (x, 0) = 0,
(40)
The problem is solved on the region −60 x 100.
Fig. 11 shows the initial solitary wave of amplitude 1.0 breaking into two smaller diverging solitary waves. The break-up is completed at approximately t = 10 and the amplitudes of these two solitary waves are approximately equal to 0.48. It is also noted that the two solitary waves are symmetric and move in opposite directions with oscillating tail between the two waves, which are in good agreement with those presented in [9]. Fig. 12 shows the error in the discrete mass for standard finite volume element scheme and energy-preserving finite volume element scheme, which well agree with the discrete mass conservation law. In addition, we also investigate the
68
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 11. The wave break-up.
α = 1.0, h = t = 0.1.
Fig. 12. (a) Error in the discrete mass for FVEM; (b) Error in the discrete mass for EFVEM. Wave break-up.
Fig. 13. (a) Error in the discrete energy for FVEM; (b) Error in the discrete energy for EFVEM. Wave break-up.
conservation of energy for wave break-up. The associated results are presented by Fig. 13. We can find that the energypreserving finite volume element scheme can strictly conserve the energy in numerical sense, while the standard finite volume element scheme only conserves energy within 10−3 . 5. Conclusions In this paper, the energy-preserving finite volume element method has been designed for solving the initial boundary problems of the improved Boussinesq equation. We also proved the conservation of the proposed numerical schemes.
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
69
Numerical experiments are performed to validate the effectiveness of the schemes and theoretical analysis. Numerical results demonstrate that the proposed finite volume element scheme is second-order accuracy in space and time. Moreover, the new scheme can conserve mass and energy to machine accuracy. Additionally, the method can also be used to solve other similar problems. Acknowledgements This project is partially supported by the National Natural Science Foundation of China grant No. 11071123, and the Innovation Program for University Postgraduates in Jiangsu Province Nos. CXZZ12_0382, CXLX12_0388. The authors would like to thank two anonymous referees for their useful comments and suggestions which have helped to improve the paper greatly. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]
N. Zabusky, M. Kruskal, Interaction of solitons in a collisionless plasma and the recurrence of initial states, Phys. Rev. Lett. 15 (6) (1965) 240. M. Wadati, Stochastic Korteweg–de Vries equation, J. Phys. Soc. Jpn. 52 (1983) 2642–2648. G. Lin, L. Grinberg, G.E. Karniadakis, Numerical studies of the stochastic Korteweg–de Vries equation, J. Comput. Phys. 213 (2) (2006) 676–703. S. Koide, D. Furihata, Nonlinear and linear conservative finite difference schemes for regularized long wave equation, Jpn. J. Ind. Appl. Math. 26 (1) (2009) 15–40. A. Bratsos, A second order numerical scheme for the improved Boussinesq equation, Phys. Lett. A 370 (2) (2007) 145–147. A. Bratsos, A predictor–corrector scheme for the improved Boussinesq equation, Chaos Solitons Fractals 40 (5) (2009) 2083–2094. P.L. Christiansen, V. Muto, S. Rionero, Solitary wave solutions to a system of Boussinesq-like equations, Chaos Solitons Fractals 2 (1) (1992) 45–50. C. Christov, Conservative difference scheme for Boussinesq model of surface waves, in: Proceedings ICFD V, 1995, pp. 343–349. H. El-Zoheiry, Numerical study of the improved Boussinesq equation, Chaos Solitons Fractals 14 (3) (2002) 377–384. Q. Lin, Y. Wu, R. Loxton, S. Lai, Linear B-spline finite element method for the improved Boussinesq equation, J. Comput. Appl. Math. 224 (2) (2009) 658–667. D. Irk, I˙ . Da˘g, Numerical simulations of the improved Boussinesq equation, Numer. Methods Partial Differ. Equ. 26 (6) (2010) 1316–1327. Z. Zhang, F. Lu, Quadratic finite volume element method for the improved Boussinesq equation, J. Math. Phys. 53 (2012) 013505. Q. Wang, Z. Zhang, Q. Zhu, Numerical simulation of the stochastic damped improved Boussinesq equation, J. Math. Phys. 54 (2013) 013503. J. Xu, Q. Zou, Analysis of linear and quadratic simplicial finite volume methods for elliptic equations, Numer. Math. 111 (3) (2009) 469–492. Z. Zhang, Error estimates of finite volume element method for the pollution in groundwater flow, Numer. Methods Partial Differ. Equ. 25 (2) (2009) 259–274. Q. Wang, Z. Zhang, Z. Li, A Fourier finite volume element method for solving two-dimensional quasi-geostrophic equations on a sphere, Appl. Numer. Math. 71 (2013) 1–13. Q. Wang, S. Lin, Z. Zhang, Numerical methods for a fluid mixture model, Int. J. Numer. Methods Fluids 71 (1) (2013) 1–12. D. Furihata, A stable and conservative finite difference scheme for the Cahn–Hilliard equation, Numer. Math. 87 (4) (2001) 675–699. T. Matsuo, New conservative schemes with discrete variational derivatives for nonlinear wave equations, J. Comput. Appl. Math. 203 (1) (2007) 32–56. Y. Miyatake, T. Matsuo, Conservative finite difference schemes for the Degasperis–Procesi equation, J. Comput. Appl. Math. 236 (15) (2012) 3728–3740. A. Biswas, D. Milovic, A. Ranasinghe, Solitary waves of Boussinesq equation in a power law media, Commun. Nonlinear Sci. Numer. Simul. 14 (11) (2009) 3738–3742. I. Bogolubsky, Some examples of inelastic soliton interaction, Comput. Phys. Commun. 13 (3) (1977) 149–155. L. Iskandar, P.C. Jain, Numerical solutions of the improved Boussinesq equation, Proc. Indian Acad. Sci. Math. Sci. 89 (3) (1980) 171–181.
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Energy-preserving finite volume element method for the improved Boussinesq equation Quanxiang Wang a,b , Zhiyue Zhang a,∗ , Xinhua Zhang b , Quanyong Zhu a,c a b c
Jiangsu Key Laboratory for NSLSCS, School of Mathematical Sciences, Nanjing Normal University, Nanjing 210023, China College of Engineering, Nanjing Agricultural University, Nanjing 210031, China School of Science, Lishui University, Lishui 323000, China
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 2 January 2014 Received in revised form 23 March 2014 Accepted 27 March 2014 Available online 2 April 2014
In this paper, we design an energy-preserving finite volume element scheme for solving the initial boundary problems of the improved Boussinesq equation. Theoretical analysis shows that the proposed numerical schemes can conserve the energy and mass. Numerical experiments are performed to illustrate the efficiency of the scheme and theoretical analysis. While the results demonstrate that the proposed finite volume element scheme is second-order accuracy in space and time. Moreover, the new scheme can conserve mass and energy. © 2014 Elsevier Inc. All rights reserved.
Keywords: Conservative Improved Boussinesq equation Energy Finite volume element
1. Introduction The Korteweg–de Vries (KdV) equation is often used to describe nonlinear dispersive wave phenomena, such as the propagation of long waves on the surface of shallow water, the propagation of long waves in one-dimensional nonlinear lattices, and the propagation of ion-acoustic waves. In the past decades, much attention has been paid to investigate the KdV equation theoretically and numerically by many authors [1–4]. Boussinesq equation was first presented by Joseph Boussinesq in 1873, it is a fourth-order nonlinear partial differential equation and belongs to Korteweg–de Vries equation (KdV) family. The Boussinesq equation which describes motions of long waves in shallow water under gravity propagating in both directions, is given by
utt = u xx + qu xxxx + u 2
xx
(1)
,
where u = u (x, t ) is a sufficiently differentiable function, which for q = −1 gives the good Boussinesq or well-posed equation, while for q = 1 the bad or ill-posed Boussinesq equation [5,6]. In this study, we will take q = 1. If the term u xxxx in Eq. (1) is replaced with u xxtt , it gives the so-called improved Boussinesq equation
utt = u xx + u xxtt + u 2
xx
,
x ∈ [a, b].
(2)
Researchers do not observe any instabilities for the above improved Boussinesq equation [7]. In addition, the equation has in common the mass and energy conservation property [8]:
*
Corresponding author. E-mail address: [email protected] (Z. Zhang).
http://dx.doi.org/10.1016/j.jcp.2014.03.053 0021-9991/© 2014 Elsevier Inc. All rights reserved.
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
dM dt
59
b = 0,
M=
u dx,
(3)
a
dE dt
b = 0,
E=
v2 2
+
(ut )2 2
+ G (u ) dx,
(4)
a 2
3
under some suitable boundary conditions, where ut = v x , G (u ) = u2 + u3 . The improved Boussinesq equation has been the subject of extensive investigations and some numerical schemes have been designed to study it. For example, El-Zoheiry [9] designed a three-level iterative scheme based on the compact implicit method for solving the improved Boussinesq equation. Bratsos [5] proposed a second order implicit finite difference scheme by rational approximations. Recently, Bratsos [6] constructed a predictor–corrector scheme to avoid the solution of a nonlinear system. Lin et al. [10] developed a finite element scheme for solving initial boundary value problems for the improved Boussinesq equation. Irk and Da˘g [11] proposed two numerical methods based on classical five point finite difference approximation. Zhang and Lu [12] designed quadratic finite volume element schemes for the equation. Then, Wang et al. [13] numerically studied stochastic damped improved Boussinesq equation using the finite volume element schemes. However, these numerical methods are not concerned about the conservation properties of numerical schemes. The finite volume element method has been one of the most commonly used numerical methods for solving partial differential equations. One main attractive property of the method is that, the mass conservation law is maintained, which is fairly desirable for fluid and underground fluid computations. So it has been extensively used in computational fluid dynamics [14–17]. In this study, we propose an energy-preserving finite volume element method for the improved Boussinesq equation. The method conserves mass M and energy E of the improved Boussinesq equation. The numerical schemes are constructed using the discrete variational derivative method (DVDM). Furihata has used the DVDM to design a stable finite difference scheme for the Cahn–Hilliard equation [18]. Further, Matsuo et al. [19,20] have extended the method to solve nonlinear Klein–Gordon equation and Degasperis–Procesi equation. The organization of the paper is the following. In Section 2, preliminaries and notations are introduced. In Section 3, the proposed energy-preserving finite volume element schemes are presented. We provide some numerical experiments to illustrate the effectiveness and accuracy of the proposed schemes in Section 4. Conclusions are drawn in Section 5. 2. Preliminaries and notations In this section, we introduce notations and useful propositions employed in this paper. 2.1. Trial and test function spaces The region of interest [a, b] can be decomposed into a grid T h with nodes
a = x0 < x1 < x2 < · · · < x N −1 < x N = b, where xi − xi −1 = h and h = nodes
b−a . N
(5)
Denote T h = { I i : I i = [xi −1 , xi ], i = 1, 2, . . . , N }. Accordingly we place a dual grid T h∗ with
a = x0 < x1/2 < x3/2 < · · · < x N −3/2 < x N −1/2 < x N = b, where xi −1/2 = xi − h/2 (i = 1, 2, . . . , N ). Denote I 0∗
= [x0 , x1/2 ], I i∗
(6)
= [xi −1/2 , xi +1/2 ] and I ∗N
= [xN −1/2 , xN ]. Select the trial function space U h as the linear element space with respect to T h . The basis function of the node xi is
φi (x) =
1 − |x − xi |/h, 0,
x i −1 x x i +1 , elsewhere.
Then numerical solution U for Eq. (2) can be uniquely written as U = element I i , we have
U = U i −1 (1 − μ) + U i μ,
U = (U i − U i −1 )/h,
(7)
N
i =1
U i φi (x), where U i = U (xi , t ). So in the
(8) (9)
where μ = (x − xi −1 )/h. The test function space V h corresponding to T h∗ is taken as the piecewise constant function space. The test function of the nodes x j is
ψ j (x) =
1, x j −1/2 x x j +1/2 , 0, elsewhere.
(10)
60
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
The numerical solution is denoted by U kn u (kh, nt ), where t is the time mesh size. For the improved Boussinesq equation, we adopt the following periodic boundary conditions
∂ j u ∂ j u = ( j = 0, 1, 2). ∂ x j x=a ∂ x j x=b
(11)
So the discrete periodic boundary conditions are
U kn = U kn mod N
for ∀k ∈ Z.
(12)
2.2. Conserved quantities In order to see the conservation property, it is more convenient to work with the following system representation
ut = v x , v t − v xxt = ∂∂x δδGu ,
(13)
where δδGu is the Euler–Lagrange variational derivative defined by
∂G ∂ δG δG . = − δu ∂u ∂ x δux
So the numerical solution V for Eq. (13) can be uniquely written as V = Proposition 1. Let u be the solution of Eq. (13). Then the total mass
d
b a
N
i =1
V i φi (x), where V i = V (xi , t ).
u dx is independent of t. Namely,
b u dx = 0.
dt
(14)
a
The conservation of mass can be easily proved [21] and the standard finite volume element schemes can conserve the mass of IBq [13]. Proposition 2. Let u be the solution of Eq. (13). Then the total energy
d
b
dt
b a
G˜ (u , ut , v ) dx is independent of t. Namely,
G˜ (u , ut , v ) dx = 0,
(15)
a
where
v2
G˜ (u , ut , v ) =
2
+
(ut )2 2
+ G (u ).
Proof. According to the chain rule, we have
d
b
dt
G˜ dx =
a
b v v t + ut utt + a
δG ut dx δu
By the first equality of Eq. (13), we can get
b = a
δG v v t + v x v xt + v x dx δu
Using the integration by part formula, we have
b =
v v t − v xxt − a
∂ δG dx ∂ x δu
At last, we can complete the proof using the second equality of Eq. (13).
2
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
61
3. Design of schemes 3.1. Discrete variational derivative method In the following, we recall the discrete variational derivative method briefly. Firstly, we define the energy as follows:
b J (u ) =
G˜ (u , ut , v ) dx,
(16)
a
where we call G˜ the energy function throughout this paper. As the discretization of energy function, we employ the following finite volume element approximation
G˜ d (U , U t , V ) G˜ .
(17)
Then we define the total discrete energy as follows:
b J (u ) =
G˜d (U , U t , V ) dx.
(18)
a
˜ ˜ To obtain the approximation δ(Uδ G, V ) of Euler–Lagrange derivative δδGu , we need to compute the following difference:
b J (U ) − J ( V ) = a
δ G˜ (U − V ) dx. δ(U , V )
(19)
3.2. Energy-preserving finite volume element scheme Recall the definition G˜ (u , ut , v ) =
δ G δ G˜ δ u , δ ut
˜ = ut and δδGv
v2 2
+
(ut )2 2
+ G (u ). According the definition to Euler–Lagrange derivative, we have
δ G˜ δu
=
= v. So Eq. (13) can be rewritten as
δ G˜
ut = ∂∂x δ v , ˜ (1 − ∂x2 ) v t = ∂∂x δδGu ,
(20)
2 where ∂x2 = ∂∂x2 .
Using any ψ ∈ V h and ψ˜ ∈ V h to multiply both sides of Eq. (20) and integrating on the interval [a, b], we can obtain the following semi-discrete scheme
,ψ , ˜ = ∂ δ G˜ , ψ˜ , ((1 − ∂x2 ) v t , ψ) ∂ x δu (ut , ψ) =
∂ δ G˜ ∂x δv
(21)
where (·, ·) denotes the inner product of L 2 ([a, b]). Now we discretize the above semi-discrete scheme in time. To simplify the notation, we will use the difference operator 1
δn un = δ G˜ δv
un+1 −un−1 . 2t
˜ Additionally, we define the discrete version of Euler–Lagrange derivative δδGu =
u2 2
+
˜ u3 , δδuG 3 t
= ut and
= v by
δ G˜ d u n ( u n +1 + u n + u n −1 ) n u + , = 3 δ(un+1 , un , un−1 ) δ G˜ d 1 = δn un , 1 n δ(δn u ) δ G˜ d = vn. δ( v n )
(22) (23)
(24)
At last, we can get the following energy-preserving finite volume element scheme
⎧ ⎨ (δ 1 un , ψ) =
δ G˜ d ∂ ∂ x δ( v n ) , ψ , δ G˜ d ⎩ ((1 − ∂ 2 )δn1 v n , ψ) ˜ ˜ = ∂ x ∂ x δ(un+1 ,un ,un−1 ) , ψ . n
Eq. (25) has the following discrete mass and energy conservation law.
(25)
62
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Theorem 1 (Discrete mass conservation law). Under the discrete periodic boundary condition, then numerical solution by Eq. (25) conserves mass
b un dx = const.
(26)
a
Theorem 2 (Discrete energy conservation law). Under the discrete periodic boundary condition, then numerical solution by Eq. (25) conserves energy
b
G˜ nd dx = const.
(27)
a
Proof.
b
1
t
˜ 1 1 δn un+1 − δn un−1 δ G˜ d δ G d 1 n δ G˜ d 1 n + + , u v , δ , δ n n 1 2t δ( v n ) δ(un+1 , un , un−1 ) δ(δn un ) ˜ ˜ n+1 ˜ ∂ δ G˜ d δGd δGd δ G d 1 n ∂ v x − v nx −1 + , = + v , , δ n ∂ x δ( v n ) 2t δ( v n ) δ(un+1 , un , un−1 ) ∂ x δ( v n ) ˜ δ G˜ d δ G˜ d ∂ 2 1 n δ G d + δ 1 − ∂ v , , =− n x ∂ x δ(un+1 , un , un−1 ) δ( v n ) δ( v n ) = 0. 2
G˜ nd+1 − G˜ nd dx =
a
3.3. Finite volume element scheme For comparison, we employ the standard finite volume element scheme constructed based on the
u = v , t x (1 − ∂x2 ) v t = (u + u 2 )x .
(28)
Then we can obtain the following finite volume element scheme
1
(δn un , ψ) = ( v nx , ψ), 1 ˜ = ((un + (un )2 )x , ψ). ˜ ((1 − ∂x2 )δn v n , ψ)
(29)
Theorem 3 (Discrete mass conservation law). Under the discrete periodic boundary condition, then numerical solution by Eq. (29) conserves mass
b un dx = const.
(30)
a
4. Numerical experiments In this section, we present some numerical examples using the proposed energy-preserving finite volume element method. To show the efficiency of the present method, we calculate L 1 , L 2 and L ∞ error norms defined as follows:
L1 =
N uh (xi ) − u (xi ) h,
(31)
i =1
N uh (xi ) − u (xi ) 2 h, L2 = i =1
L ∞ = max u h (xi ) − u (xi ) . i
Define the rate of convergence r:
r = log2
u 2h − u
,
u h − u
(32) (33)
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
63
Table 1 The rate of convergence in space for the EFVEM. h
L ∞ error
r
L 2 error
r
L 1 error
r
0.40 0.20 0.10 0.05
9.78e−04 2.46e−04 6.23e−05 1.58e−05
– 1.99 1.98 1.98
1.95e−03 4.92e−04 1.24e−04 3.18e−05
– 1.98 1.98 1.96
6.35e−03 1.59e−03 4.04e−04 1.03e−04
– 2.00 1.98 1.97
Table 2 The rate of convergence in time for the EFVEM.
t
L ∞ error
r
L 2 error
r
L 1 error
r
0.040 0.020 0.010 0.005
1.42e−03 2.99e−04 7.47e−05 1.87e−05
– 2.24 2.00 2.00
2.53e−03 6.15e−04 1.54e−04 3.84e−05
– 2.04 2.00 2.00
8.24e−03 2.02e−03 5.06e−04 1.26e−04
– 2.03 2.00 2.01
where u h is the numerical solution with space step size h, u the analytical solution, N the number of node. The rate approaches the number 2 would indicate second-order accuracy in space and time. 4.1. Single solitary wave Eq. (2) has a special solution [22] which describes solitary wave:
u (x, t ) = α sech
2
α x − x0 − β t 6
β
,
β=
1+
2α 3
(34)
,
where α is amplitude of the solitary wave which is initially centered at x = x0 traveling with velocity β . Differentiating the above equation with respect to time and then setting t = 0, the initial conditions are obtained
u (x, 0) = α sech2
ut (x, 0) = 2α
α 6
α x − x0 6
sech2
β
α x − x0 6
(35)
,
β
tanh
α x − x0 6
β
.
(36)
To compare our results with those given in other references, we will take x0 = 0 and α = 0.5 in the following. In order to demonstrate the accuracy in space and time, we discretize the problem over the region [−60, 100] up to time t = 1. Firstly, we choose a relatively small time step t = 0.001 such that the error from the time discretization is negligible compared with those from the space. The L 1 , L 2 and L ∞ error norms for u have been computed by the energy-preserving finite volume element method (EFVEM). They are displayed in Table 1. Examination of the table shows that the error measures of the EFVEM diminish approximately quadratically as the space step is halved. In other words, the proposed energy-preserving scheme is second-order in space. To investigate the accuracy in time, the time and space step are chosen such that the ratio t /h = 0.1. We list the errors of the numerical solution computed by our new scheme in Table 2 in the L 1 , L 2 and L ∞ error norms. Second-order convergence can be observed. We perform a simulation with space step h = 0.2 and time step t = 0.1 for single solitary wave on the time interval [0, 40]. Fig. 1 presents that the solitary wave propagates towards the right without change in form. It is easy to find that the solitary wave continues to evolve with the initial velocity. We also perform a simulation with the resolution to check the conservation of schemes. Fig. 2 shows the error in the discrete mass for standard finite volume element scheme (FVEM) and energy-preserving finite volume element scheme, which well agree with the discrete mass conservation law. Fig. 3 presents the error in the discrete energy for FVEM and EFVEM. From the figure, we can find that standard finite volume scheme can conserve energy within 10−3 , which though very good is not strict in numerical sense. On the contrary, the proposed energy-preserving finite volume element scheme conserves energy within 10−14 , which is strict conservation in numerical sense and conforms the theoretical result. To make comparison with the earlier works [5,11,22], the numerical solution is given at t = 72 with α = 0.5 and x ∈ [−80, 100]. The results are presented in Table 3. We can find that our numerical method is much more efficient than those presented in the references. For example, L ∞ error = 4.1403e−04 at t = 72 is achieved with t = 0.1. Fig. 4 shows the single solitary wave at t = 72. Contrarily to [22], the amplitude of solitary wave does not decrease and there is no visible secondary solitary wave produced, which is consistent with those presented in [5]. 4.2. Two solitary waves In this example, we study the interaction of two solitary waves moving on a collision course. Here, x ∈ [−60, 100] and the initial conditions are
64
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 1. The propagation for single solitary wave.
α = 0.5, h = 0.2 and t = 0.1.
Fig. 2. (a) Error in the discrete mass for FVEM; (b) Error in the discrete mass for EFVEM. Single solitary wave.
Fig. 3. (a) Error in the discrete energy for FVEM; (b) Error in the discrete energy for EFVEM. Single solitary wave.
Table 3 Numerical results with
EFVEM M1 in [11] [5]
α = 0.5 at t = 72. h
t
L ∞ error
Energy error
0.1 0.1 0.1
0 .1 0.001 0.001
4.1403e−04 1.0064e−03 4.1998e−04
7.9492e−14
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 4. The propagation for single solitary wave.
α = 0.5, h = t = 0.1.
Fig. 5. The propagation for two solitary waves (a) and level curves (b).
2
u (x, 0) = α1 sech
ut (x, 0) = 2α1
α1 x − x01 6
α1 6
+ 2α2
6
+ α2 sech
β1
α1 x − x01
sech2
α2
6
2
sech
√
α2 x − x02 β2
α2 x − x02
β1 6
2
tanh
6
β2
6
tanh
α2 x − x02 β2
(37)
,
β1 6
α1 = α2 = 0.4, h = t = 0.1.
α1 x − x01
65
,
(38)
√
where x01 = −20, x02 = 60, β1 = 1 + 2α1 /3, β2 = − 1 + 2α2 /3. The first solitary wave is initially located at x = x01 and moving to the right with speed β1 and amplitude α1 . The other one is initially located at x = x02 and moving to the left with speed |β2 | and amplitude α2 . Firstly, we study the collision of two solitary waves with the same amplitude. Fig. 5 presents the propagation for two solitary waves with α1 = α2 = 0.4. From the figure, we can find that the solitary waves preserve their amplitudes and no displacement is observed. Additionally, there is no visible secondary solitary waves produced. Therefore, we can conclude that the collision is elastic. Fig. 6 presents the propagation for two solitary waves with α1 = α2 = 0.8. It can be observed from the figure that the secondary solitary waves are visible. Hence, the interaction is inelastic. The results are in good agreement with those reported in [5,10]. Now, we give some results for the interaction of solitary waves of different amplitudes. Fig. 7 shows interaction of two solitary waves with α1 = 0.8 and α2 = 0.4. Fig. 8 displays the collision of two solitary waves with α1 = 1.5 and α2 = 0.4. From the figures, we can find that the collisions are inelastic. According to definition used in [5,22] for the inelasticity coefficient, the coefficient K is defined as follows
K=
αm , max(α1 , α2 )
(39)
66
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 6. The propagation for two solitary waves (a) and level curves (b).
α1 = α2 = 0.8, h = t = 0.1.
Fig. 7. The propagation for two solitary waves (a) and level curves (b).
α1 = 0.8, α2 = 0.4, h = t = 0.1.
Fig. 8. The propagation for two solitary waves (a) and level curves (b).
α1 = 1.5, α2 = 0.4, h = t = 0.1.
where αm is the maximum joint amplitude at the collision time. The results for αm and K are presented in Table 4. Examination of the table shows that αm is always smaller than twice the amplitude of the max(α1 , α2 ). In addition, the coefficient K will decrease when α1 increase. The results are contrary to those presented in [23], Fig. 9 shows the error in the discrete mass for standard finite volume element scheme and energy-preserving finite volume element scheme, which well agree with the discrete mass conservation law. Fig. 10 shows the error in the discrete energy for FVEM and EFVEM. From the figure, we can make a similar conclusion. The energy-preserving finite volume
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Table 4 Numerical results for different
67
α1 and α2 at t = 80.
α1
α2
αm
K
0.4 0.8 0.8 1.5
0.4 0.8 0.4 0.4
0.7225 1.3895 1.0738 1.7188
1.8063 1.7369 1.3422 1.1458
Fig. 9. (a) Error in the discrete mass for FVEM; (b) Error in the discrete mass for EFVEM. Two solitary waves.
α1 = 0.8, α2 = 0.4, h = t = 0.1.
Fig. 10. (a) Error in the discrete energy for FVEM; (b) Error in the discrete energy for EFVEM. Two solitary waves.
α1 = 0.8, α2 = 0.4, h = t = 0.1.
element scheme can conserve the energy to machine precision. On the other hand, the standard finite volume element scheme can only conserve energy within 10−3 . 4.3. Wave break-up The initial conditions are given in the following:
u (x, 0) = α sech where
2
α x − x0 6
β
α = 1.0, x0 = 20 and β = 1 +
, 2α . 3
ut (x, 0) = 0,
(40)
The problem is solved on the region −60 x 100.
Fig. 11 shows the initial solitary wave of amplitude 1.0 breaking into two smaller diverging solitary waves. The break-up is completed at approximately t = 10 and the amplitudes of these two solitary waves are approximately equal to 0.48. It is also noted that the two solitary waves are symmetric and move in opposite directions with oscillating tail between the two waves, which are in good agreement with those presented in [9]. Fig. 12 shows the error in the discrete mass for standard finite volume element scheme and energy-preserving finite volume element scheme, which well agree with the discrete mass conservation law. In addition, we also investigate the
68
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
Fig. 11. The wave break-up.
α = 1.0, h = t = 0.1.
Fig. 12. (a) Error in the discrete mass for FVEM; (b) Error in the discrete mass for EFVEM. Wave break-up.
Fig. 13. (a) Error in the discrete energy for FVEM; (b) Error in the discrete energy for EFVEM. Wave break-up.
conservation of energy for wave break-up. The associated results are presented by Fig. 13. We can find that the energypreserving finite volume element scheme can strictly conserve the energy in numerical sense, while the standard finite volume element scheme only conserves energy within 10−3 . 5. Conclusions In this paper, the energy-preserving finite volume element method has been designed for solving the initial boundary problems of the improved Boussinesq equation. We also proved the conservation of the proposed numerical schemes.
Q. Wang et al. / Journal of Computational Physics 270 (2014) 58–69
69
Numerical experiments are performed to validate the effectiveness of the schemes and theoretical analysis. Numerical results demonstrate that the proposed finite volume element scheme is second-order accuracy in space and time. Moreover, the new scheme can conserve mass and energy to machine accuracy. Additionally, the method can also be used to solve other similar problems. Acknowledgements This project is partially supported by the National Natural Science Foundation of China grant No. 11071123, and the Innovation Program for University Postgraduates in Jiangsu Province Nos. CXZZ12_0382, CXLX12_0388. The authors would like to thank two anonymous referees for their useful comments and suggestions which have helped to improve the paper greatly. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]
N. Zabusky, M. Kruskal, Interaction of solitons in a collisionless plasma and the recurrence of initial states, Phys. Rev. Lett. 15 (6) (1965) 240. M. Wadati, Stochastic Korteweg–de Vries equation, J. Phys. Soc. Jpn. 52 (1983) 2642–2648. G. Lin, L. Grinberg, G.E. Karniadakis, Numerical studies of the stochastic Korteweg–de Vries equation, J. Comput. Phys. 213 (2) (2006) 676–703. S. Koide, D. Furihata, Nonlinear and linear conservative finite difference schemes for regularized long wave equation, Jpn. J. Ind. Appl. Math. 26 (1) (2009) 15–40. A. Bratsos, A second order numerical scheme for the improved Boussinesq equation, Phys. Lett. A 370 (2) (2007) 145–147. A. Bratsos, A predictor–corrector scheme for the improved Boussinesq equation, Chaos Solitons Fractals 40 (5) (2009) 2083–2094. P.L. Christiansen, V. Muto, S. Rionero, Solitary wave solutions to a system of Boussinesq-like equations, Chaos Solitons Fractals 2 (1) (1992) 45–50. C. Christov, Conservative difference scheme for Boussinesq model of surface waves, in: Proceedings ICFD V, 1995, pp. 343–349. H. El-Zoheiry, Numerical study of the improved Boussinesq equation, Chaos Solitons Fractals 14 (3) (2002) 377–384. Q. Lin, Y. Wu, R. Loxton, S. Lai, Linear B-spline finite element method for the improved Boussinesq equation, J. Comput. Appl. Math. 224 (2) (2009) 658–667. D. Irk, I˙ . Da˘g, Numerical simulations of the improved Boussinesq equation, Numer. Methods Partial Differ. Equ. 26 (6) (2010) 1316–1327. Z. Zhang, F. Lu, Quadratic finite volume element method for the improved Boussinesq equation, J. Math. Phys. 53 (2012) 013505. Q. Wang, Z. Zhang, Q. Zhu, Numerical simulation of the stochastic damped improved Boussinesq equation, J. Math. Phys. 54 (2013) 013503. J. Xu, Q. Zou, Analysis of linear and quadratic simplicial finite volume methods for elliptic equations, Numer. Math. 111 (3) (2009) 469–492. Z. Zhang, Error estimates of finite volume element method for the pollution in groundwater flow, Numer. Methods Partial Differ. Equ. 25 (2) (2009) 259–274. Q. Wang, Z. Zhang, Z. Li, A Fourier finite volume element method for solving two-dimensional quasi-geostrophic equations on a sphere, Appl. Numer. Math. 71 (2013) 1–13. Q. Wang, S. Lin, Z. Zhang, Numerical methods for a fluid mixture model, Int. J. Numer. Methods Fluids 71 (1) (2013) 1–12. D. Furihata, A stable and conservative finite difference scheme for the Cahn–Hilliard equation, Numer. Math. 87 (4) (2001) 675–699. T. Matsuo, New conservative schemes with discrete variational derivatives for nonlinear wave equations, J. Comput. Appl. Math. 203 (1) (2007) 32–56. Y. Miyatake, T. Matsuo, Conservative finite difference schemes for the Degasperis–Procesi equation, J. Comput. Appl. Math. 236 (15) (2012) 3728–3740. A. Biswas, D. Milovic, A. Ranasinghe, Solitary waves of Boussinesq equation in a power law media, Commun. Nonlinear Sci. Numer. Simul. 14 (11) (2009) 3738–3742. I. Bogolubsky, Some examples of inelastic soliton interaction, Comput. Phys. Commun. 13 (3) (1977) 149–155. L. Iskandar, P.C. Jain, Numerical solutions of the improved Boussinesq equation, Proc. Indian Acad. Sci. Math. Sci. 89 (3) (1980) 171–181.
