Numerical Solutions of the Gardner Equation by Extended Form of the Cubic B-splines

TL;DR

This paper introduces an extended cubic B-spline collocation method to numerically solve the Gardner equation, aiming to improve accuracy and efficiency over classical methods, validated through analytical and stability analyses.

Contribution

The study develops an extended cubic B-spline approach for the Gardner equation, optimizing the extension parameter to enhance solution accuracy and stability.

Findings

The method accurately solves initial boundary value problems for the Gardner equation.

Conservation laws are maintained throughout the simulations.

Von Neumann analysis confirms the method's stability.

Abstract

The extended definition of the polynomial B-splines may give a chance to improve the results obtained by the classical cubic polynomial B-splines. Determination of the optimum value of the extension parameter can be achieved by scanning some intervals containing zero. This study aims to solve some initial boundary value problems con- structed for the Gardner equation by the extended cubic B-spline collocation method.The test problems are derived from some analytical studies to validate the efficiency and accuracy of the suggested method. The conservation laws are also determined to observe them remain constant as expected in theoretical aspect. The stability of the proposed method is investigated by the Von Neumann analysis.

Tables

Table 1. Table 1: Error norms for various choices of extension parameter

$N$	$L_{\infty }\left(2.5\right)\left(\lambda =0\right)$	$L_{\infty }\left(2.5\right)($various $\lambda )$	$L_{\infty }\left(5\right)\left(\lambda =0\right)$	$L_{\infty }\left(5\right)($various $\lambda )$
$100$	$3.2726\times {10}^{-5}$	$\left(\lambda =-0.00840\right)1.2330\times {10}^{-5}$	$5.22606\times {10}^{-5}$	$2.2789\times {10}^{-5}$
$200$	$2.0537\times {10}^{-5}$	$\left(\lambda =-0.00280\right)1.4819\times {10}^{-5}$	$1.91604\times {10}^{-5}$	$1.9119\times {10}^{-5}$
$300$	$1.4428\times {10}^{-5}$	$\left(\lambda =-0.00094\right)1.2509\times {10}^{-5}$	$1.70403\times {10}^{-5}$	$1.6944\times {10}^{-5}$
$400$	$1.4452\times {10}^{-5}$	$\left(\lambda =-0.00178\right)1.4440\times {10}^{-5}$	$1.61150\times {10}^{-5}$	$1.5872\times {10}^{-5}$

Table 2. Table 2: Calculated conservation laws and their absolute relative changes

$N$	$M_0$	$E_0$	$H_0$	$C\left(M_5\right)$	$C\left(E_5\right)$	$C\left(H_5\right)$
$100$	$1.0445$	$0.0601$	$0.0040$	$5.4748\times {10}^{-6}$	$3.8176\times {10}^{-8}$	$1.5233\times {10}^{-6}$
$200$	$1.0445$	$0.0601$	$0.0040$	$3.2669\times {10}^{-6}$	$5.1126\times {10}^{-8}$	$1.7003\times {10}^{-6}$
$300$	$1.0445$	$0.0601$	$0.0040$	$2.4190\times {10}^{-7}$	$2.1767\times {10}^{-8}$	$2.8351\times {10}^{-6}$
$400$	$1.0445$	$0.0601$	$0.0040$	$1.3753\times {10}^{-6}$	$2.0910\times {10}^{-10}$	$3.3939\times {10}^{-6}$

Table 3. Table 3: Error norms for optimum choices of extension parameter

$N$	$L_{\infty }\left(4\right)\left(\lambda =0\right)$	$L_{\infty }\left(4\right)($various $\lambda )$	$L_{\infty }\left(12\right)\left(\lambda =0\right)$	$L_{\infty }\left(12\right)($various $\lambda )$
$100$	$8.4150\times {10}^{-6}$	$\left(\lambda =-0.01850\right)3.8974\times {10}^{-6}$	$2.3158\times {10}^{-5}$	$\left(\lambda =-0.01850\right)1.2330\times {10}^{-5}$
$200$	$2.1207\times {10}^{-6}$	$\left(\lambda =-0.00574\right)1.0194\times {10}^{-6}$	$5.9956\times {10}^{-6}$	$\left(\lambda =-0.00574\right)2.9662\times {10}^{-6}$
$400$	$5.3296\times {10}^{-7}$	$\left(\lambda =-0.00115\right)2.5440\times {10}^{-7}$	$1.5016\times {10}^{-6}$	$\left(\lambda =-0.00115\right)7.7413\times {10}^{-7}$
$600$	$2.2377\times {10}^{-7}$	$\left(\lambda =-0.00057\right)1.1335\times {10}^{-7}$	$6.6655\times {10}^{-7}$	$\left(\lambda =-0.00057\right)3.3921\times {10}^{-7}$
$800$	$1.4601\times {10}^{-7}$	$\left(\lambda =-0.00024\right)6.3749\times {10}^{-8}$	$5.2835\times {10}^{-6}$	$\left(\lambda =-0.00024\right)5.2779\times {10}^{-6}$

Table 4. Table 4: Calculated conservation laws and their absolute relative changes

$100$	$16.1599$	$3.0129$	$0.0979$	$4.9504\times {10}^{-3}$	$5.3104\times {10}^{-3}$	$5.4423\times {10}^{-3}$
$200$	$16.0799$	$2.9969$	$0.0974$	$4.9751\times {10}^{-3}$	$5.3388\times {10}^{-3}$	$5.4721\times {10}^{-3}$
$400$	$16.0399$	$2.9889$	$0.0972$	$4.9875\times {10}^{-3}$	$5.3531\times {10}^{-3}$	$5.4871\times {10}^{-3}$
$600$	$16.0266$	$2.9862$	$0.0971$	$4.9916\times {10}^{-3}$	$5.3578\times {10}^{-3}$	$5.4922\times {10}^{-3}$
$800$	$16.0199$	$2.9849$	$0.0970$	$4.9938\times {10}^{-3}$	$5.3603\times {10}^{-3}$	$4.9481\times {10}^{-3}$

Table 5. Table 5: Calculated conservation laws and their absolute relative changes

$t$	$M_0$	$E_0$	$H_0$	$C\left(M_{1}5\right)$	$C\left(E_{1}5\right)$	$C\left(H_{1}5\right)$
$5$	$5.2255$	$1.5033$	$1.5994$	$8.0719\times {10}^{-7}$	$3.0588\times {10}^{-5}$	$1.2886\times {10}^{-3}$
$10$	$5.2255$	$1.5033$	$1.5994$	$2.7652\times {10}^{-6}$	$4.1342\times {10}^{-5}$	$1.8485\times {10}^{-3}$
$15$	$5.2255$	$1.5033$	$1.5994$	$7.0380\times {10}^{-6}$	$6.1132\times {10}^{-4}$	$2.1571\times {10}^{-3}$

Equations

$$u_t+{\mu }_{1}u{u}_x+{\mu }_2{u}^2{u}_x+{\mu }_3{u}_{xxx}=0$$

$$u_t+\left({\mu }_{1}u+{\mu }_2{u}^2\right)u_x+{\mu }_3{v}_{xx}$$

$$v-u_x$$

$$u(x,0)$$

$$v(x,0)$$

$$u_x(a,t)$$

$$u_{xx}(a,t)$$

$$v_x(a,t)$$

$$v_{xx}(a,t)$$

$$\pi :a=x_0<x_1<\dots <x_N=b$$

$$E_{j}(x)=\frac{1}{24h^{4}}\left\{\begin{array}[]{ll}4h(1-\lambda)(x-x_{j-2})^{3}+3\lambda(x-x_{j-2})^{4},&\left[x_{j-2},x_{j-1}\right],\\ \begin{array}[]{l}(4-\lambda)h^{4}+12h^{3}(x-x_{j-1})+6h^{2}(2+\lambda)(x-x_{j-1})^{2}\\ -12h(x-x_{j-1})^{3}-3\lambda(x-x_{j-1})^{4}\end{array}&\left[x_{j-1},x_{i}\right],\\ \begin{array}[]{l}(4-\lambda)h^{4}-12h^{3}(x-x_{j+1})+6h^{2}(2+\lambda)(x-x_{j+1})^{2}\\ +12h(x-x_{j+1})^{3}-3\lambda(x-x_{j+1})^{4}\end{array}&\left[x_{i},x_{j+1}\right],\\ 4h(\lambda-1)(x-x_{j+2})^{3}+3\lambda(x-x_{j+2})^{4},&\left[x_{j+1},x_{j+2}\right],\\ 0&\text{otherwise.}\end{array}\right.$$

$$U(x,t)=\sum _{j=-1}^{N+1}{\delta }_j{E}_j(x),V(x,t)=\sum _{j=-1}^{N+1}{\varphi }_j{E}_j(x)$$

$$U_j=U(x_j,t)=\frac{4-\lambda }{24}{\delta }_{j-1}+\frac{8+\lambda }{12}{\delta }_j+\frac{4-\lambda }{24}{\delta }_{j+1},$$

$$\begin{array}[]{r}\dfrac{U^{n+1}-U^{n}}{\Delta t}+\mu_{1}\dfrac{(UU_{x})^{n+1}+(UU_{x})^{n}}{2}+\mu_{2}\dfrac{(U^{2}U_{x})^{n+1}+(U^{2}U_{x})^{n}}{2}+\mu_{3}\dfrac{V_{xx}^{n+1}+V_{xx}^{n}}{2}=0\\ \\ \dfrac{V^{n+1}+V^{n}}{2}-\dfrac{U_{x}^{n+1}+U_{x}^{n}}{2}=0\end{array}$$

$$(U{U}_x{)}^{n+1}=U^{n+1}{U}_x^n+U^n{U}_x^{n+1}-U^n{U}_x^n$$

$$(U^2{U}_x{)}^{n+1}=2U^{n+1}{U}^n{U}_x^n+(U^n{)}^2{U}_x^{n+1}-2(U^n{)}^2{U}_x^n$$

$$\left[\left(\frac{2}{\Delta t}+{\mu }_{1}L+2{\mu }_{2}KL\right){\alpha }_1+\left({\mu }_{1}K+{\mu }_2{K}^2\right){\beta }_1\right]{\delta }_{j-1}^{n+1}+\left[{\mu }_3{\gamma }_1\right]{\varphi }_{j-1}^{n+1}+\left[\left(\frac{2}{\Delta t}+{\mu }_{1}L+2{\mu }_{2}KL\right){\alpha }_2\right]{\delta }_j^{n+1}$$

$$+\left[{\mu }_3{\gamma }_2\right]{\varphi }_j^{n+1}\left[\left(\frac{2}{\Delta t}+{\mu }_{1}L+2{\mu }_{2}KL\right){\alpha }_1-\left({\mu }_{1}K+{\mu }_2{K}^2\right){\beta }_1\right]{\delta }_{j+1}^{n+1}+\left[{\mu }_3{\gamma }_1\right]{\varphi }_{j+1}^{n+1}$$

$$=\left[\left(\frac{2}{\Delta t}+{\mu }_{2}KL\right){\alpha }_1\right]{\delta }_{j-1}^n-{\mu }_3{\gamma }_1{\varphi }_{j-1}^n+\left[\left(\frac{2}{\Delta t}+{\mu }_{2}KL\right){\alpha }_2\right]{\delta }_j^n-{\mu }_3{\gamma }_2{\varphi }_j^n+\left[\left(\frac{2}{\Delta t}+{\mu }_{2}KL\right){\alpha }_1\right]{\delta }_{j+1}^n$$

$$-{\beta }_1{\delta }_{j-1}^{n+1}+{\alpha }_1{\varphi }_{j-1}^{n+1}+{\alpha }_2{\varphi }_j^{n+1}+{\beta }_1{\delta }_{j+1}^{n+1}+{\alpha }_1{\varphi }_{j+1}^{n+1}={\beta }_1{\delta }_{j-1}^n-{\alpha }_1{\varphi }_{j-1}^n-{\alpha }_2{\varphi }_j^n-{\beta }_1{\delta }_{j+1}^n-{\alpha }_1{\varphi }_{j+1}^n$$

$$m=0,...,N,n=0,\mathrm{1...},$$

$$\begin{array}[]{l}K=\alpha_{1}\delta_{j-1}^{n}+\alpha_{2}\delta_{j}^{n}+\alpha_{1}\delta_{j+1}^{n}\\ L=\alpha_{1}\phi_{j-1}^{n}+\alpha_{2}\phi_{j}^{n}+\alpha_{1}\phi_{j+1}^{n}\end{array}$$

$${\alpha }_1$$

$${\beta }_1$$

$${\mathbf{Ax}}^{n+1}={\mathbf{Bx}}^n$$

$$\mathbf{A}=\left[\begin{array}{cccccccccc}{\nu }_{m1}& {\nu }_{m2}& {\nu }_{m3}& {\nu }_{m4}& {\nu }_{m5}& {\nu }_{m2}& & & & \\ -{\beta }_1& {\alpha }_1& 0& {\alpha }_2& {\beta }_1& {\alpha }_1& & & & \\ & & {\nu }_{m1}& {\nu }_{m2}& {\nu }_{m3}& {\nu }_{m4}& {\nu }_{m5}& {\nu }_{m2}& & \\ & & -{\beta }_1& {\alpha }_1& 0& {\alpha }_2& {\beta }_1& {\alpha }_1& & \\ & & & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & & & {\nu }_{m1}& {\nu }_{m2}& {\nu }_{m3}& {\nu }_{m4}& {\nu }_{m5}& {\nu }_{m2}\\ & & & & -{\beta }_1& {\alpha }_1& 0& {\alpha }_2& {\beta }_1& {\alpha }_1\end{array}\right]$$

$$\mathbf{B}=\left[\begin{array}{cccccccccc}{\nu }_{m6}& -{\nu }_{m2}& {\nu }_{m7}& -{\nu }_{m4}& {\nu }_{m6}& -{\nu }_{m2}& & & & \\ {\beta }_1& -{\alpha }_1& 0& -{\alpha }_2& -{\beta }_1& -{\alpha }_1& & & & \\ & & {\nu }_{m6}& -{\nu }_{m2}& {\nu }_{m7}& -{\nu }_{m4}& {\nu }_{m6}& -{\nu }_{m2}& & \\ & & {\beta }_1& -{\alpha }_1& 0& -{\alpha }_2& -{\beta }_1& -{\alpha }_1& & \\ & & & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & & & {\nu }_{m6}& -{\nu }_{m2}& {\nu }_{m7}& -{\nu }_{m4}& {\nu }_{m6}& -{\nu }_{m2}\\ & & & & {\beta }_1& -{\alpha }_1& 0& -{\alpha }_2& -{\beta }_1& -{\alpha }_1\end{array}\right]$$

$$\begin{array}[]{ll}\nu_{m1}=\left(\frac{2}{\Delta t}+\mu_{1}L+2\mu_{2}KL\right)\alpha_{1}+\left(\mu_{1}K+\mu_{2}K^{2}\right)\beta_{1}&\nu_{m6}=\left(\frac{2}{\Delta t}+\mu_{2}KL\right)\alpha_{1}\\ \nu_{m2}=\mu_{3}\gamma_{1}&\nu_{m5}=\left(\frac{2}{\Delta t}+\mu_{1}L+2\mu_{2}KL\right)\alpha_{2}-\left(\mu_{1}K+\mu_{2}K^{2}\right)\beta_{1}\\ \nu_{m3}=\left(\frac{2}{\Delta t}+\mu_{1}L+2\mu_{2}KL\right)\alpha_{2}-\left(\mu_{1}K+\mu_{2}K^{2}\right)\beta_{1}&\nu_{m7}=\left(\frac{2}{\Delta t}+\mu_{2}KL\right)\alpha_{2}\\ \nu_{m4}=\mu_{3}\gamma_{2}&\end{array}$$

$${\mathbf{x}}^{n+1}=({\delta }_{-1}^{n+1},{\varphi }_{-1}^{n+1},{\delta }_0^{n+1},{\varphi }_0^{n+1},\dots ,{\delta }_{n+1}^{n+1},{\varphi }_{n+1}^{n+1},).$$

$$\begin{array}[]{l}\delta_{-1}=\delta_{1}\\ \phi_{-1}=\phi_{1}\\ \delta_{N-1}=\delta_{N+1}\\ \phi_{N-1}=\phi_{N+1}\end{array}$$

$$\begin{array}[]{l}U_{x}(a,0)=0=\delta_{-1}^{0}-\delta_{1}^{0},\\ U(x_{j},0)=\alpha_{1}\delta_{j-1}^{0}+\alpha_{2}\delta_{j}^{0}+\alpha_{1}\delta_{j+1}^{0}=u(x_{j},0),j=1,...,N-1\\ U_{x}(b,0)=0=\delta_{N-1}^{0}-\delta_{N+1}^{0},\\ V_{x}(a,0)=0=\phi_{-1}^{0}-\phi_{1}^{0}\\ V(x_{j},0)=\alpha_{1}\phi_{j-1}^{0}+\alpha_{2}\phi_{j}^{0}+\alpha_{1}\phi_{j+1}^{0}=v(x_{j},0),j=1,...,N-1\\ V_{x}(a,0)=\phi_{N-1}^{0}-\phi_{N+1}^{0}\end{array}$$

$${\delta }_j^n$$

$${\varphi }_j^n$$

$$\rho =\frac{{\xi }^{n+1}}{{\xi }^n}$$

$$a_1{\delta }_{j-1}^{n+1}+a_2{\delta }_j^{n+1}+a_1{\delta }_{j+1}^{n+1}+\frac{\lambda k\epsilon }{2}(a_3{\delta }_{j-1}^{n+1}-a_3{\delta }_{j+1}^{n+1})+\frac{k{\mu }_3}{2}(a_4{\varphi }_{j-1}^{n+1}+a_5{\varphi }_j^{n+1}+a_4{\varphi }_{j+1}^{n+1})$$

$$a_3{\delta }_{j-1}^{n+1}-a_3{\delta }_{j+1}^{n+1}-a_1{\varphi }_{j-1}^{n+1}-a_2{\varphi }_j^{n+1}-a_1{\varphi }_{j+1}^{n+1}$$

Full text

Keywords: Gardner equation; soliton; perturbation; wave generation; extended B-spline.

1 Introduction

Consider the Gardner equation (or combined KdV-mKdV) of the form

[TABLE]

where $u=u(x,t)$ and $\mu_{1}$, $\mu_{2}$ and $\mu_{3}$ are constants. The Gardner equation has two nonlinear terms in the quadratic and cubic forms and the dissipative term is of third order. The Gardner equation is an integrable system and Miura transformation connects it to the KdV equation[1].

The Gardner equation is a useful model to understand the propagation of negative ion acoustic plasma waves[2]. The equation can be derived from the system of plasma motion equations in one dimension with arbitrarily charged cold ions and inertia neglected isothermal electrons.

The Gardner equation can also be a good description of internal waves with large amplitudes[3]. The modulation system in the Riemann invariant form and classified solutions for large $t$ values subject to some particular initial conditions are studied deeply to explain the undular bore formation. The structure of the nonlinear terms and the positions of the initial step are significant for the solution classification. The sign of $\mu_{2}$ has also critical role for the structures of the solutions. The patterns of the solutions can be of the forms bright or dark cnoidal or trigonometric bores, kinks, rarefaction waves or combinations of them.

Some particular forms of the internal ocean rogue waves in the coasts occur in the Gardner equation when the $\mu_{2}$ and $\mu_{3}$ have the same sign which enables modulational instability[4]. The Gardner equation can also be a useful model to study stratified fluid transcritical flow passing a topographic obstacle when a forcing term is added to the equation[5].

Interacting two-soliton type solution is derived by Darboux transformations in Slyunyaev and Pelinovskii’s study[6]. The consistent tanh method is also capable of generating interacting solutions for the Gardner equation[7]. Some interacting of two wave solutions such as soliton-cnoidal wave or soliton-periodic wave are presented in this study. The consistent Riccati expansion is another method that is capable to obtain soliton-cnoidal wave interaction type solutions[8].

Dynamics of solitons, in the adiabatic parameter case, of the perturbed Gardner equation is studied in details by Biswas and Zerrad[9]. The motion integrals are also determined for the perturbed form of the equation. The effects of perturbation to a single soliton type solution are discussed by considering the $\tanh$ type initial data[10].

Exact solutions to the Gardner equation are set up by using various methods. some solutions containing $\tanh$ and $\coth$ functions are proposed by the extended form of the $\tanh$ method[11]. Solitary wave and periodic solutions are constructed by aid of the projective Riccati equations[12]. These solutions have various terms including trigonometric or hyperbolic functions in rational forms. $G^{\prime}/G$ is another expansion method to determine the exact solutions of the Gardner equation. Some solitary wave, periodic, exponential, rational and complex-type traveling wave solutions are found by this method[13]. Some other exact solutions in terms of trigonometric functions [14, 15, 16], hyperbolic functions[14, 15, 17], kink solutions[18] are determined by using different expansion or ansatz methods.

Numerical solutions to the Gardner equations are also concerned in various studies. The conservative finite difference schemes are developed to determine propagation of one soliton and collusion of two soliton solutions numerically[19]. Restrictive Taylor’s technique has also been implemented to simulate the propagation of some solutions numerically[20].

In the present study, we develop an extended cubic B-spline collocation method to the solutions of the Gardner equation. Having only two continuous derivatives of the cubic B-splines force us to reduce the order of the third order derivative term. Setting $v=u_{x}$ reduces the order of the Gardner equation (1) in the resultant coupled system of nonlinear PDEs of the form

[TABLE]

In order to complete the mathematical notation of the initial-boundary value problems (IVPs), we use the initial data

[TABLE]

and the homogeneous Neumann conditions

[TABLE]

at both ends of the problem interval $[a,b]$.

2 Numerical Approximate

Consider the equal grid distribution

[TABLE]

of the finite interval $[a,b]$ where $h=(b-a)/N$ is the equal mesh size. An extended B-spline function is defined as

[TABLE]

with the extension parameter $\lambda$. The set of extended B-spline functions $\{E_{-1}(x),E_{0}(x),\dots,E_{N+1}(x)\}$ constitutes a basis function set for the functions defined in this interval. The approximate solutions $U(x,t)$ and $V(x,t)$ to $u(x,t)$ and $v(x,t)$, respectively, can be written in terms of the extended B-splines as

[TABLE]

where $\delta_{j}$ and $\phi_{j}$ are time dependent parameters. These parameters are determined after implementation of the collocation method and complementary data. Since each extended B-spline has lowest two derivatives, the nodal derivative values of both $U$ and $V$ can be summarized in terms of extended B-splines as

[TABLE]

The time integration of the space-splitted system (2) is performed by the Crank-Nicolson method as

[TABLE]

where the superscript p represents the solution at the $p$th time level with equal time step size $\Delta t$ satisfying $t^{n+1}=t^{n}+\Delta t$.

The nonlinear terms $(UU_{x})^{n+1}$ and $\left(U^{2}U_{x}\right)^{n+1}$ in Eq. (8) are converted to linear forms by using

[TABLE]

and

[TABLE]

defined in [21]. The resultant linear system is discretized in time by using Crank-Nicolson method as

[TABLE]

where

[TABLE]

[TABLE]

This system can be written in the matrix notation as

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The system (10) consists of $2N+2$ linear equations and $2N+6$ unknowns

[TABLE]

The unique solution of this system requires additional four constraints. The boundary data $U_{x}(a,t)=0,$ $V_{x}(a,t)=0$ and $U_{x}(b,t)=0,$ $V_{x}(b,t)=0$ can be written in terms of parameters as the following equations:

[TABLE]

The parameters $\delta_{-1},\phi_{-1},\delta_{N+1},\phi_{N+1}$ in Eq.(10) are eliminated from the system by using the boundary data to determine a solvable system. In order to initialize the iteration algorithm, initial parameters $\delta_{j}^{0},\phi_{j}^{0},$ $j=-1,\ldots,N+1$ are determined by using the data obtained from the complementary data as

[TABLE]

3 Stability Analysis

The stability of the method is investigated by performing the Von-Neumann analysis where

[TABLE]

[TABLE]

Here, $A_{1}$ and $A_{2}$ represent the harmonics amplitude. $k$ is the mode number, $\rho$ is the amplification factor and $\varphi=kh$. The term $U+U^{2}$ is assumed as locally constant and replaced $\varepsilon$. Substituting 11 into the system

[TABLE]

[TABLE]

gives

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It can be concluded from both (14) and (15) that $\left|\rho\right|$ is less than or equal to $1$. Thus, the proposed method method is unconditionally stable.

4 Numerical Illustrations

The numerical solutions for some IBVPs set up with the Gardner equation are summarized in this section. The accuracy of the results determined by the extended B-spline collocation method is discussed by examining graphical representations, measuring the error between the numerical and the analytical solutions and the preservation of conservation laws. The error of the numerical solution is measured by using the discrete maximum norm defined as

[TABLE]

where $U_{j}^{n}$ and $u(x_{j},t)$ are numerical and analytical solutions at the discrete time $t$.

The conservation laws can also be indicators of the validity of the proposed algorithms even when the analytical solution does not exist. The conservation laws of the Gardner equation

[TABLE]

are expected to keep their initial quantities during numerical simulations [22]. The relative changes of these quantities at a discrete time $t>0$ are measured by using $C(M_{t})$, $C(E_{t})$ and $C(H_{t})$ defined as

[TABLE]

where $M_{t}$, $E_{t}$, $H_{t}$, ($t\geq 0$) are the measured quantities at the time $t$.

4.1 Propagation of Initial Single Positive Pulse

In the first numerical illustration, we study propagation of an initial single pulse with positive amplitude. The equation parameters are chosen as $\mu_{1}=4$, $\mu_{2}=-3$ and $\mu_{3}=1$. The initial data are determined from the exact solution [18]

[TABLE]

by assuming $t=0$. Since the exact solution approaches zero as $x$ approaches infinity or minus infinity, the choice of homogeneous boundary conditions is compatible with the solution. The artificial interval $[-20,30]$ is chosen for the numerical simulation and the designed algorithm is run up the simulation ending time $t=5$ with the fixed time discretization parameter $\Delta t=0.1$ and various spatial discretization numbers $N$. The propagation of a positive single solitary wave is depicted in Fig 1.

The error distribution for the extension parameter $\lambda=0$ and the discretization parameters $\Delta t=0.1$ and $N=100$ at the simulation ending time $t=5$ is depicted in Fig 2(a). Usage of the same discretization parameters for optimum extension parameter $\lambda=-0.00840$ results in the error distribution at the simulation ending time as given in Fig 2(b). A simple comparison shows that the results are improved approximately two times when the optimum extension parameter is used in the algorithm.

The discretized maximum error norms for various values of grid size are summarized in Table 1. In the meanwhile, the algorithm seeks for the optimal value of the extension parameter $\lambda$ by comparing the discrete maximum error norm at each $\lambda$. Even though we do not observe an improvement in decimal digits in the results, using the optimum extension parameter $\lambda$ generates improved results for all choices of $N$.

The conservation laws are required to preserve their initial values as time proceeds during the simulation. The initial values of these laws are calculated by using Maple, Table 2. The absolute relative changes of conservation laws are obtained at least six decimal digits at the simulation ending time $t=5$. These preservation rates can be accepted as indicators of a valid algorithm.

4.2 Propagation of Kink-like Wave

Kink-like wave solution of the Gardner equation is

[TABLE]

for the equation parameters $\mu_{1}=1$, $\mu_{2}=-5$ and $\mu_{3}=1$[23]. This wave propagates along the $x-$axis with the speed $1/30$. The initial data are generated from the analytical solution by assuming $t=0$. Since the analytical solution disappears as $x\rightarrow\infty$ and $u(x,t)\rightarrow 0.2$ as $x\rightarrow-\infty$, the homogeneous Neumann data are compatible. The artificial problem interval is chosen as $[-80,80]$ and the designed algorithm is run up to the ending time $t=12$ with $\Delta t=0.1$ and various grid numbers used in the space domain. The plot Fig 3 is a summary of the propagation in this finite interval.

The designed algorithm is scanned the extension parameter in $[-1,1]$ with step size $\Delta\lambda=0.000001$ for an optimum value to improve the results. Optimum extension parameters are determined for all choices of the number of grids, Table 3. Determination of the optimum extension parameter reduces the maximum absolute error to the half for all choices of $N$.

The initial values of the conservation laws are calculated by integration the quantities by substituting the initial data of the IBVP. The absolute relative changes of the conservation laws indicate a reliable solution in Table 4.

4.3 Wave Generation from an Initial Pulse

The perturbed Gardner equation of the form

[TABLE]

for some nonzero real $\epsilon$ can be useful to study the wave generation from an initial positive pulse. The initial data are chosen as

[TABLE]

with the parameter choice $\mu_{1}=10$, $\mu_{2}=-3$ and $\mu_{3}=1$ in the Gardner Equation (19). The designed program is run in the finite interval $[-40,60]$ with the parameters $N=400$ and $\Delta t=0.01$ up to the $t=15$. The initial positive pulse of height $0.4305$ positioned at $x=5$ propagates along the $x$-axis, Fig 5(a). When the propagation time reaches $t=5$, a frontier wave of height $0.6568$ is positioned at $x=18.25$, Fig 5(b). The first follower solitary wave is clearly observable at this time. This follower is of height $0.3318$ and it is positioned at $x=12$. Even though it is not clearly observable, a bulge at the left of the first follower can be evaluated as an indicator of a second follower wave. The first three waves are clearly observable at $t=10$, Fig 5(c). The height of the frontier reaches $0.6871$ and its peak point is positioned at $x=28.75$. The peak of the first follower wave of height $0.3913$ is positioned at $x=18.25$. At the left of the first follower, almost completely formed second follower of height $0.1736$ is positioned at $x=9.75$. The height of the frontier is measured as $0.6941$ at the time $t=15$, Fig 5(d). Its peak is positioned at $x=39$ at this time. The height of the first follower wave reaches $0.3998$ and its peak position is measured as $x=24.75$. The second follower of height $0.1910$ is positioned at $x=13$. The bulge appearing at the left of the second follower is the indicator of formation of one more solitary.

Even though we can not compare the results with the analytical solutions, the absolute relative changes of the conservation laws is a good indicator to evaluate the efficiency of the proposed method. The initial values of the conservation laws and their absolute relative changes are summarized in Table 5. The absolute relative change of the first conservation law is only in six decimal digits, as the second one in four decimal digits and the third one in three decimal digits at the simulation ending time $t=15$.

5 Conclusion

For the sake of the improve the results obtained by the cubic B-splines, the extended cubic B-spline based collocation method is constructed for some initial boundary value problems for the Gardner Equation. The scan of the optimum extension parameter in the interval $[-1,1]$ gives opportunity to obtain more reliable results when compare classical cubic B-splines. Since the extended cubic b-splines have only first and second order derivatives, the reduction of the order of the third order derivative is required. Thus, the coupled system of nonlinear PDEs is obtained. The extended B-spline function are used to approximate the solutions of this system. Following the spatial discretization, the linearization procedure is followed. At the end, the time integration is done Crank-Nicolson method. Von-Neumann stability analysis shows that the suggested algorithm is unconditionally stable.

The first two examples give opportunity to measure the error between the analytical and the numerical solutions by calculating maximum error norms for various choices of the discretization parameters. Both graphical representations and the absolute relative changes of the conservation laws are indicators of a reliable and valid method. The results are improved by determining the optimum value of the extension parameter.

In the third example, a non analytical problem simulating wave generation from an initial single solitary wave is studied. The proposed algorithm simulates the expected results successfully. The absolute relative changes of the conservation laws confirm the valid results.

6 Acknowledgements

This study is a part of the project with number 2016/19052 supported by Eskisehir Osmangazi University Scientific Research Projects Committee and was partially presented at 3rd International Conference on Pure and Applied Sciences, Dubai, 2017.