A Quadratic B-Spline Galerkin Approach for Solving a Coupled KdV Equation

In this paper, a quadratic B-spline Galerkin finite element approach is applied to one-dimensional coupled KdV equation in order to obtain its numerical solutions. The performance of the method is examined on three test problems. Computed results are compared with the exact results and also other numerical results given in the literature. A Fourier stability analysis of the approach is also done.


Introduction
In this paper, we will consider the coupled KdV equation of the form is related to most types of long waves having weak dispersion. Some internal acoustic and planetary waves arising in geophysical fluid mechanics can be given as examples [6].
The coupled KdV type equations have been the most important class of nonlinear evolution equations arising in physical sciences and engineering. Therefore, the coupled KdV type equations have been analytically considered by several authors: Hirota and Satsuma [13], and Tam et al. [27] have presented a coupled Korteweg-de Vries equation and obtained soliton solutions of the equation. Tian and Gao [28] have obtained the exact solutions to the Bogoyavlenskii coupled KdV equations via an auto-Backlund transformation with the application of the Painlevė analysis and computer algebra. Fan and Zhang [10] have obtained several kinds of exact solutions for a system of coupled KdV equations using an improved homogeneous balance method. Roy [20] has examined the bi-Hamiltonian structure of two coupled KdV equations proposed by Hirota and Satsuma [13], and Ito [17]. Zhu [30] has given a difference scheme for the periodic initial-boundary problem of the coupled KdV Equation. Cao et al. [7] have obtained more exact solutions for a new coupled MKdV equations by using a direct and efficient trigonometric function transform method based on the idea of the homogeneous balance method and for KdV equations by Miura transformation. Zhou et al. [29] have obtained the periodic wave solutions to a coupled KdV equations with variable coefficients by using F-expansion method which can be thought of as an over-all generalization of Jacobi elliptic function expansion method. Karasu and Kilic [18] have studied the Painleve property of coupled, non-autonomous KdV type of systems and obtained the conditions under which the systems pass the Painleve test for integrability. Qian and Tian [24] have found the single soliton solutions for a coupled KdV equation using the non-local Lie-Backlund transformation theorem to the trivial zero solution. Inan [14] has given some exact solutions of the coupled KdV equation by using the generalized tanh function method. Ma and Zhu [21] have obtained some new exact solutions of the coupled KdV equations via the Jacobian elliptic function expansion approach and Hermite transformation. Assas [5] has applied He's variational iteration method to solve the non-linear coupled-KdV equation based on the use of Lagrange multipliers for identification of optimal value of a parameter in a functional. Abbasbandy [1] has used an analytical technique, namely the homotopy analysis method, to solve a generalized Hirota-Satsuma coupled KdV equation. Al-Khaled et al. [3] have applied both the tanh and the He's variational iteration methods for analytical study for the nonlinear coupled KdV equations. Mokhtari and Mohammadi [22] have decomposed a coupled system of nonlinear partial differential equations into a set of algebraic equations as well as an ordinary differential equation and then solved by using Exp-function method.
The coupled KdV type equations have been also considered numerically by some authors: Fan [9] has obtained numerical solutions of the coupled KdV equation by using a Riccati equation involving a parameter and its solutions to replace the tanh-function in the tanh method. Halim et al. [11] have introduced a numerical method for general coupled KdV systems valid for solving Cauchy problems for arbitrary number of equations with arbitrary constant coefficients.
Halim and Leble [12] have selected a second covariant equation to form Lax pair of a coupled KdV-MKdV system and introduced a numerical method for general KdV-MKdV system. Kaya and Inan [19] have found the explicit and numerical traveling wave solutions for a coupled KdV equation and a coupled MKdV equation by using the decomposition method with the help of symbolic computation. Alvarez-Samaniegoa and Carvajalb [4] have studied the locally well-posed conditions of coupled KdV equations for some systems. Ismail [16] has set up a numerical method for solving the coupled KdV equation based on the collocation method with quintic B-spline finite elements to simulate the solution of coupled KdV equation. Siraj-ul-Islam et al. [15] have formulated a simple classical radial basis functions collocation method for the numerical solution of the coupled KdV equations. Rady et al. [2] have considered the system of coupled KdV equations and established the transformation which turns the coupled KdV equations into the single nonlinear partial differential equation, then they obtained an auto-Backlund transformation and lax pairs using the extended homogeneous balance method. Biswas and Ismail [6] have used solitary wave ansatz to carry out the integration of the coupled KdV equation with power law nonlinearity, and then supplemented their results by numerical simulations. Recently, Rashid and Ismail [25] have obtained error estimates of spectral collocation method for the coupled KdV equations and presented their numerical solutions.
In this paper, we have applied a Galerkin quadratic B-spline finite element method to the one dimensional coupled KdV equation given by (1.1) with a set of initial and boundary conditions to obtain its numerical solutions.

The Finite Element Solution
The finite interval [c, d] is divided into N finite elements of equal length h by the knots x i , (i = 0, . . . , N ) such that c = So global approximations U N (x, t) and V N (x, t) to the exact solutions U (x, t) and V (x, t) are respectively sought in terms of the quadratic B-splines as where δ j and σ j are time dependent parameters to be determined from the boundary and weighted residual conditions.
In terms of the local coordinate transformation ξ = x − x m , 0 ≤ ξ ≤ h, the splines can be expressed as Since all other quadratic B-splines are zero over the element [x m , x m+1 ], the approximation (2.1) over this element can be written in terms of the basis functions (2.2) as In the above equations and throughout this paper, the prime denotes differentiation with respect to x.
which can also be written in the matrix form as A eδ e − 6aB e (δ)δ e − 2b B e (σ)σ e + aC e δ e = 0, A eσ e + 3B e (δ)σ e − C e σ e = 0. (2.5) In the above equations and throughout the paper, the dot denotes differentiation with respect to t, δ e = (δ m−1 , δ m , δ m+1 ) and σ e = (σ m−1 , σ m , σ m+1 ) are the element parameters, and A e ij , B e ikj , B e ikj and C e ij are the element matrices given by the following integrals: Assembling all contributions coming from all the elements, Eq. (2.5) yields the system of equations By substituting the Crank-Nicolson approaches δ = 0.5(δ n + δ n+1 ), σ = 0.5(σ n + σ n+1 ) and the forward finite difference approximationṡ Using the boundary conditions in (2.8) yields a 2N × 2N system of linear equations. The obtained system is easily solved by a variation of Gauss-elimination algorithm. The following inner iterations are applied two or three times at each time to improve the accuracy of the approximation

Using the relations
, since the first derivative of the approximate initial conditions shall agree with those of the exact initial conditions, initial vector δ 0 j can be obtained from the following system of linear equations: which can be solved using an appropriate algorithm. A similar system of equation with the same matrix is solved for σ 0 j . Hence, the approximate solution functions for U (x, t) and V (x, t) can be obtained from δ n and σ n using Eq. (2.8).

Stability analysis
For stability analysis, it is convenient to use the Fourier stability analysis. Since system (2.8) consists of two equations with two variables, we use the following Fourier modes [26]: where i = √ −1, q is a complex number, ϕ is a real number and P , W are harmonic amplitudes. For the stability, the condition |q| ≤ 1 must be satisfied.

Numerical Examples and Results
The importance to construct discrete schemes for which the discrete versions of conservation laws are satisfied is well known for many mathematical models. We note nonlinear Schrödinger problems, where discrete conservation laws are playing a very important role in the analysis of stability and accuracy of numerical approximations, see [8]. All numerical calculations were executed on a Pentium i7 PC in the Fortran code using double precision arithmetic. We use the L 2 and L ∞ error norms defined by and evaluate only two constants of motion to validate the conservation properties [16] To implement the performance of the method, the following three test problems have been considered: The motion of a single soliton, the interaction of two solitons and birth of solitons.
We will solve the above problem in the interval −25 ≤ x ≤ 25 by taking the following values of a, b and λ to compare our results with the earlier works. For these three cases the computed values of the invariants I 1 and I 2 with the error norms L 2 and L ∞ are presented respectively at some selected times for h = 0.2, 0.1 and ∆t = 0.01 in Tables 1-3. All tables confirm that the error norms L 2 and L ∞ are still small when the time is increased up to t = 20. The values of the invariants I 1 and I 2 remain almost constant during the computer run. For example, the relative change of the invariants I 1 and I 2 are respectively 0.093×10 −3 % and 0.415×10 −3 % for Case (a), 0.070×10 −3 % and 0.103×10 −3 % for Case (b), 0.057 × 10 −3 % and 0.140 × 10 −3 % for Case (c) with ∆t = 0.01 and h = 0.1. Table 4 displays the values of x and amplitudes of U N and V N at different values of t. It is clearly seen from the table that the computed values of the amplitudes are very close to their exact ones. It is also observed from the table that the soliton moves to the right at an almost unchanged amplitude with an increasing of time, as expected. Figure 1 shows the profiles of single solitons at t = 0, 10 and 20 with the error distributions of U N and V N . It is also clearly seen from the figure that the soliton moves to the right at an almost unchanged amplitude with an increasing of time as mentioned above and the error distributions are high around the right boundary for V N about the position x at which U N attains its highest amplitude. Table 5 displays a comparison of the values of the invariants and the error norm L ∞ obtained by the present method with those obtained by the quintic Table 2.
The numerical solutions of single soliton with a = −0.5, b = 3, λ = 0.5 and ∆t = 0.01.  Table 3. The numerical solutions of single soliton with a = −0.125, b = −3, λ = 0.5 and ∆t = 0.01.  Table 4. B-spline collocation finite element method given in Ref. [16] at different times for the values of a, b and λ given in Cases (a)-(c). It is clearly seen from the table that the invariants are in very good agreement with each other and the present method produces the error norms L ∞ at each time smaller than those given in Ref. [16].

Interaction of two solitons
Secondly, the coupled KdV equation (1.1) has been considered with the boundary conditions given by (1.2) and the initial conditions , For this problem, all calculations are done in the range −10 ≤ x ≤ 120 for λ 1 = 1, λ 2 = 0.6, γ 1 = 10, γ 2 = 30, a = 0.5 and b = −3. Tables 6 Table 5. Comparison of the numerical solutions of the single soliton with results from [16] with ∆t = 0.01 and h = 0.   Table 8 displays a comparison of the invariants obtained by the present method with those obtained by the quintic B-spline collocation finite element method given in Ref. [16]. From Table 8, it is clearly seen that both results are in good agreement with each other.  The experiment was run from t = 0 to t = 90 to allow the interaction of two solitons to take place. Figure 2 shows the interaction of two solitons for U N . As it is seen from Figure 2, the amplitudes of the big and small waves at t = 0 are 1.996731 at the point x = 10.2 and 0.719566 at the point x = 53.1, respectively. As time increases, both of the waves move forward the right and then the big wave catches up the smaller one. At t = 90, it is observed that the big wave leaves from the small one and they both move forward. At t = 90, the amplitude of the big wave is 1.987528 at the point x = 103.1 whereas the amplitude of the smaller one is 0.719579 at the point x = 80.9. Figure 3 illustrates the behaviour of the interaction of two solitons for V N at times from t = 0 to t = 90. At t = 0, the amplitudes of the big and small waves are 1.413057 at the point x = 10.2 and 0.508964 at the point x = 53.1, respectively. As time increases, both of the waves move forward the right and the big wave catches up the smaller one. It is observed that at t = 90 the big wave leaves from the small one and then they both move forward. At t = 90, the amplitude of small wave is 0.508966 with having the position x = 80.9 while the amplitude of the big one is 1.409834 with having the position x = 103.1.   As a result, Figures 2 and 3 show the interaction of two positive amplitude waves having the same position at the point x and at the same time t having different amplitudes.  Table 9. The table also presents a comparison of the invariants obtained by the present method with those given in Ref. [16]. It is obviously seen from the table that the results are in good agreement with each other.

Birth of solitons
It is clearly seen from the table that the variations of the invariants I 1 and I 2 are satisfactorily well. For example, the relative change of the invariants I 1 and I 2 are respectively 0.039% and 0.140% with ∆t = 0.02 and h = 0.2; 0.014% and 0.130% with ∆t = 0.02 and h = 0.1; 0.031% and 0.015% with ∆t = 0.01 and h = 0.2; 3.946 × 10 −3 % and 6.687 × 10 −3 % with ∆t = 0.01 and h = 0.1.
The solutions U N and V N obtained for a = 0.5, b = −3 and t = 0, 30 and 50 are illustrated in Figure 4 and 5, respectively. It is seen from the figures that at t = 0, there is only a single wave for each of U N and V N with the amplitude 1.0 at x = 0, and then each single wave for U N and V N moves to the right when time increases and then each one creates a large number of waves behind    itself with different amplitudes at t = 50. The positions and amplitudes of the successive waves at time t = 50 are given in Table 10.

Conclusions
In this paper, the numerical solutions of the coupled KdV equation with various initial and boundary conditions were obtained using the Galerkin quadratic Bspline finite element method. The efficiency of the method was tested on three numerical experiments of wave propagation: The motion of a single soliton, the interaction of two solitons and the birth of solitons. The accuracy of the method was examined by the error norms L 2 and L ∞ . The obtained results show that the error norms are reasonably small and the conservation properties are all very good. As a conclusion, the method can also be used efficiently for solving a large number of coupled nonlinear partial differential equations.