Exponential B-spline Collocation Solutions to the Gardner Equation

TL;DR

This paper introduces a collocation method using exponential B-splines to solve the Gardner equation, effectively handling its nonlinear terms and coupled PDE system with verified accuracy and conservation properties.

Contribution

The paper develops a novel exponential B-spline collocation approach for the Gardner equation, including order reduction and integration with Crank-Nicolson, demonstrating its effectiveness.

Findings

High accuracy in solving the Gardner equation

Effective handling of nonlinear quadratic and cubic terms

Preservation of conservation laws during simulations

Abstract

The exponential B-spline basis function set is used to develop a collocation method for some initial boundary value problems (IBVPs) to the Gardner equation. The Gardner equation has two nonlinear terms, namely quadratic and cubic ones. The order reduction of the equation is resulted in a coupled system of PDEs that enables the exponential B-splines to be implemented. The system is integrated in time by Crank-Nicolson implicit method. The validity of the method is investigated by calculating the discrete maximum error norm and observing the absolute relative changes of the conservation laws at the end of the simulations.

Tables

Table 1. Table 1: B m ​ ( x ) subscript 𝐵 𝑚 𝑥 B_{m}(x) and its first and second derivatives

$x$	$x_{m-2}$	$x_{m-1}$	$x_m$	$x_{m+1}$	$x_{m+2}$
$B_m$	$0$	$\frac{\sinh \left(\zeta h\right)-\zeta h}{2\left(\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)\right)}$	$1$	$\frac{\sinh \left(\zeta h\right)-\zeta h}{2\left(\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)\right)}$	$0$
$B_m^{{}^{\prime }}$	$0$	$\frac{\zeta \left(1-\cosh \left(\zeta h\right)\right)}{2\left(\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)\right)}$	$0$	$\frac{\zeta \left(\cosh \left(\zeta h\right)-1\right)}{2\left(\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)\right)}$	$0$
$B_m^{{}^{\prime \prime }}$	$0$	$\frac{{\zeta }^2\sinh \left(\zeta h\right)}{2\left(\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)\right)}$	$-\frac{{\zeta }^2\sinh \left(\zeta h\right)}{\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)}$	$\frac{{\zeta }^2\sinh \left(\zeta h\right)}{2\left(\zeta h\cosh \left(\zeta h\right)-\sinh \left(\zeta h\right)\right)}$	$0$

Table 2. Table 2: L ∞ subscript 𝐿 L_{\infty} Error norms for μ 1 = 4 , subscript 𝜇 1 4 \mu_{1}=4, μ 2 = − 3 , subscript 𝜇 2 3 \mu_{2}=-3, μ 3 = 1 , subscript 𝜇 3 1 \mu_{3}=1, Δ ​ t = 0.1 , Δ 𝑡 0.1 \Delta t=0.1, t = 2.5 𝑡 2.5 t=2.5 and 5 , 5 5, − 20 ≤ x ≤ 30 20 𝑥 30 -20\leq x\leq 30

$N$	$L_{\infty }\left(2.5\right)\left(\zeta =1\right)$	$L_{\infty }\left(2.5\right)($various $\zeta$)	$L_{\infty }\left(5\right)\left(\zeta =1\right)$	$L_{\infty }\left(5\right)($various $\zeta$)
$100$	$1.1502\times {10}^{-4}$	$\left(\zeta =0.000003\right)3.2331\times {10}^{-5}$	$2.1665\times {10}^{-4}$	$\left(\zeta =0.000003\right)5.1481\times {10}^{-5}$
$200$	$4.1696\times {10}^{-5}$	$\left(\zeta =0.000001\right)1.6622\times {10}^{-5}$	$5.7428\times {10}^{-5}$	$\left(\zeta =0.000001\right)1.8886\times {10}^{-5}$
$300$	$2.3860\times {10}^{-5}$	$\left(\zeta =0.000005\right)1.3923\times {10}^{-5}$	$2.9888\times {10}^{-5}$	$\left(\zeta =0.000004\right)1.7006\times {10}^{-5}$
$400$	$1.6985\times {10}^{-5}$	$\left(\zeta =0.000004\right)1.4470\times {10}^{-5}$	$1.8721\times {10}^{-5}$	$\left(\zeta =0.000003\right)1.5404\times {10}^{-5}$

Table 3. Table 3: Absolute relative changes of conservation laws for for ζ = 1 𝜁 1 \zeta=1

$N$	$M_0$	$E_0$	$H_0$	$C\left(M_5\right)$	$C\left(E_5\right)$	$C\left(H_5\right)$
$100$	$1.04458$	$0.06013453$	$0.00407022$	$5.5668\times {10}^{-6}$	$2.6168\times {10}^{-8}$	$1.2174\times {10}^{-5}$
$200$	$1.04458$	$0.06013453$	$0.00407022$	$2.9640\times {10}^{-6}$	$5.0740\times {10}^{-8}$	$1.0597\times {10}^{-6}$
$300$	$1.04458$	$0.06013453$	$0.00407022$	$2.3326\times {10}^{-7}$	$2.2152\times {10}^{-8}$	$2.7126\times {10}^{-6}$
$400$	$1.04458$	$0.06013453$	$0.00407022$	$1.1862\times {10}^{-6}$	$8.8551\times {10}^{-10}$	$3.3555\times {10}^{-6}$

Table 4. Table 4: L ∞ subscript 𝐿 L_{\infty} Error norms for μ 1 = 1 , subscript 𝜇 1 1 \mu_{1}=1, μ 2 = − 5 , subscript 𝜇 2 5 \mu_{2}=-5, μ 3 = 1 , subscript 𝜇 3 1 \mu_{3}=1, Δ ​ t = 0.1 , Δ 𝑡 0.1 \Delta t=0.1, t = 12 , 𝑡 12 t=12, − 80 ≤ x ≤ 80 80 𝑥 80 -80\leq x\leq 80

$N$	$L_{\infty }\left(12\right)\left(\zeta =1\right)$	$L_{\infty }\left(12\right)($various $\zeta$)
$100$	$3.8436\times {10}^{-4}$	$\left(\zeta =0.000001\right)2.3022\times {10}^{-5}$
$200$	$1.0016\times {10}^{-4}$	$\left(\zeta =0.000002\right)5.8623\times {10}^{-6}$
$400$	$2.5327\times {10}^{-5}$	$\left(\zeta =0.000004\right)1.3684\times {10}^{-6}$
$600$	$1.1280\times {10}^{-5}$	$\left(\zeta =0.000006\right)5.3420\times {10}^{-7}$
$800$	$6.3476\times {10}^{-6}$	$\left(\zeta =0.000008\right)2.3800\times {10}^{-7}$

Table 5. Table 5: Absolute relative changes of lowest three conservation laws for ζ = 1 𝜁 1 \zeta=1

$N$	$M_0$	$E_0$	$H_0$	$C\left(M_{12}\right)$	$C\left(E_{12}\right)$	$C\left(H_{12}\right)$
$100$	$16.1599$	$3.0129$	$0.0979$	$4.9493\times {10}^{-4}$	$5.3092\times {10}^{-4}$	$5.4405\times {10}^{-4}$
$200$	$16.0799$	$2.9969$	$0.0974$	$4.9750\times {10}^{-4}$	$5.3387\times {10}^{-4}$	$5.4720\times {10}^{-4}$
$400$	$16.0399$	$2.9889$	$0.0971$	$4.9875\times {10}^{-4}$	$5.3531\times {10}^{-4}$	$5.4871\times {10}^{-4}$
$600$	$16.0266$	$2.9862$	$0.0971$	$4.9917\times {10}^{-4}$	$5.3578\times {10}^{-4}$	$5.4922\times {10}^{-4}$
$800$	$16.0199$	$2.9849$	$0.0970$	$4.9937\times {10}^{-4}$	$5.3602\times {10}^{-4}$	$5.4947\times {10}^{-4}$

Table 6. Table 6: Absolute relative changes of lowest three conservation laws for for ζ = 1 𝜁 1 \zeta=1

$t$	$M_0$	$E_0$	$H_0$	$C\left(M_t\right)\left(p=1\right)$	$C\left(E_t\right)\left(p=1\right)$	$C\left(H_t\right)\left(p=1\right)$
$5$	$5.2255$	$1.5033$	$1.5994$	$1.2040\times {10}^{-6}$	$3.7180\times {10}^{-5}$	$21608\times {10}^{-3}$
$10$	$5.2255$	$1.5033$	$1.5994$	$3.6819\times {10}^{-6}$	$5.2527\times {10}^{-5}$	$3.1907\times {10}^{-3}$
$15$	$5.2255$	$1.5033$	$1.5994$	$8.9144\times {10}^{-6}$	$5.8526\times {10}^{-4}$	$3.5478\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)$$

$$\pi :x_m=a+mh,m=0,1,\dots N$$

$$B_{m}(x)=\left\{\begin{array}[]{lcc}b_{2}\left((x_{m-2}-x)-\frac{1}{\zeta}\sinh(\zeta(x_{m-2}-x))\right)&,&[x_{m-2},x_{m-1}]\\ a_{1}+b_{1}(x_{m}-x)+c_{1}\exp(\zeta(x_{m}-x))+d_{1}\exp(-\zeta(x_{m}-x))&,&[x_{m-1},x_{m}]\\ a_{1}+b_{1}(x-x_{m})+c_{1}\exp(\zeta(x-x_{m}))+d_{1}\exp(-\zeta(x-x_{m}))&,&[x_{m},x_{m+1}]\\ b_{2}\left((x-x_{m+2})-\frac{1}{\zeta}\sinh(\zeta(x-x_{m+2}))\right)&,&[x_{m+1},x_{m+2}]\\ 0&,&otherwise\end{array}\right.$$

$$\begin{array}[]{l}a_{1}=\dfrac{\zeta h\cosh(\zeta h)}{\zeta h\cosh(\zeta h)-\sinh(\zeta h)},\\ b_{1}=\dfrac{\zeta}{2}\dfrac{\cosh(\zeta h)(\cosh(\zeta h)-1)+\sinh^{2}(\zeta h)}{(\zeta h\cosh(\zeta h)-\sinh(\zeta h))(1-\cosh(\zeta h))},\\ b_{2}=\dfrac{\zeta}{2(\zeta h\cosh(\zeta h)-\sinh(\zeta h))},\\ c_{1}=\dfrac{1}{4}\dfrac{\exp(-\zeta h)(1-\cosh(\zeta h))+\sinh(\zeta h)(\exp(-\zeta h)-1))}{(\zeta h\cosh(\zeta h)-\sinh(\zeta h))(1-\cosh(\zeta h))},\\ d_{1}=\dfrac{1}{4}\dfrac{\exp(\zeta h)(\cosh(\zeta h)-1)+\sinh(\zeta h)(\exp(\zeta h)-1))}{(\zeta h\cosh(\zeta h)-\sinh(\zeta h))(1-\cosh(\zeta h))},\end{array}$$

$$U(x,t)=\sum _{m=-1}^{N+1}{\delta }_m{B}_m(x),V(x,t)=\sum _{m=-1}^{N+1}{\varphi }_m{B}_m(x).$$

$$\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_{m-1}^{n}+\delta_{m}^{n}+\alpha_{1}\delta_{m+1}^{n}\\ L=\alpha_{1}\phi_{m-1}^{n}+\phi_{m}^{n}+\alpha_{1}\phi_{m+1}^{n}\end{array}$$

$${\alpha }_1$$

$${\gamma }_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},).$$

$$U_x(a,t)$$

$$U_x(b,t)$$

$${\delta }_{-1}$$

$${\varphi }_{-1}$$

$${\delta }_{N-1}$$

$${\varphi }_{N-1}$$

$$\begin{array}[]{l}U_{x}(a,0)=0=\delta_{-1}^{0}-\delta_{1}^{0},\\ U(x_{m},0)=\alpha_{1}\delta_{m-1}^{0}+\alpha_{2}\delta_{m}^{0}+\alpha_{1}\delta_{m+1}^{0}=u(x_{m},0),m=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_{m},0)=\alpha_{1}\phi_{m-1}^{0}+\alpha_{2}\phi_{m}^{0}+\alpha_{1}\phi_{m+1}^{0}=v(x_{m},0),m=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})$$

Full text

1 Introduction

The Gardner equation

[TABLE]

where $\mu_{i},i=1,2,3$ are constant parameters has two nonlinear terms besides the third order dissipative term. When these two nonlinear terms are considered, the equation is a combined form of the Korteweg-de Vries and the modified Korteweg-de Vries equations. The Gardner equation is a significant model for the motion of negative ion-acoustic plasma waves [1]. The equation can also describe internal waves with large amplitudes and weakly nonlinear dispersive waves [2]. Moreover, undular bore type solutions for the Gardner equation and the relation of these solutions with the sign of the cubic nonlinear term are studied in details by Kamchatnov et al. in the same paper. Various patterns covering undular bores in the bright and dark cnoidals, trigonometric bores, kinks and some other types constructed from them are developed in [2]. Occurrence of large ocean waves in unexpected meaning can be modeled by in the Gardner equation and in some cases the modulational instability can be observed [3]. An analytical study deals with transcritical flow of a fluid passing a local topographical obstacle [4].

Interaction of various types of solitons and other wave profiles are examined in details in different studies. Collision of a large amplitude soliton with a limiting soliton is derived by the aid of the Darboux transformation [5]. A large class of interaction solutions of solitons with cnoidal and periodic waves is constructed by the consistent $\tanh$-expansion method [6]. Some collusion models of cnoidal wave to soliton are derived by the consistent Riccati expansion approach [7]. Some exact traveling wave type solutions constructed in terms of $\tanh$ and $\coth$ are constructed by implementing the extended form of the $\tanh$ method [8]. Projective Riccati equations also enable to generate some hyperbolic type solitary wave and periodic solutions expressed in the finite series form [9]. These type solutions can also be constructed by generalized form of the $(G^{\prime}/G)-$expansion [10, 11, 12], $\exp-$ function [13] methods. A bunch of traveling wave type exact solutions to the Gardner equation are determined by using various hyperbolic ansatzes[14]. Zayed and Abdelaziz [15] focus on new exact traveling wave solutions constructed by $(G^{\prime}/G,1/G)$ approach. Mapping method can also be used to integrate the Gardner equation [16]. Lie group and tan-cot methods are also effective methods to construct some exact solutions to the Gardner equation [17, 18].

The three conservation laws representing various physical quantities such as momentum, energy and etc. are given for the generalized Gardner equation power law nonlinearities by using some algebraic and derivative manipulations [19]. Some finite difference and restrictive Taylor’s approximation are used to determine the numerical solutions to the Gardner equation [20, 21]. In the present study, we construct exponential B-spline collocation algorithms for the numerical solutions to some IBVPs for the Gardner equation. Even though the cubic degree of the basis does not allow to approach the third order derivative, the reduction of the derivative order permits to implement of the method to the resultant nonlinear PDE system.

Let $v=u_{x}$. Then, the Gardner equation (1) is reduced to a coupled system nonlinear PDEs

[TABLE]

The initial data

[TABLE]

and the zero Neumann conditions at both end of the artificial finite problem interval $[a,b]$ are combined with the system (2) for the mathematical representation of the IBVPs.

2 Exponential B-spline Approach

Let $\pi$ be a uniformly distributed grids of the finite interval $[a,b]$, such as,

[TABLE]

where $h=\frac{b-a}{N}$. Then, the exponential cubic B-splines set

[TABLE]

where $m=-1,0,\ldots N+1$,

[TABLE]

with real parameter $\zeta$ constitutes a basis functions set for the functions defined over the problem interval $[a,b]$ [22]. Each exponential cubic B-spline $B_{m}(x)$ has two continuous first and second order derivatives defined in the interval $[x_{m-2},x_{m+2}]$. Exponential B-splines are used as basis functions in some methods to solve problems appearing in various fields [23, 24, 25, 26, 27, 28]. The functional and derivative values of the exponential B-splines are summarized in Table 1 .

Let $U$ and $V$ be approximate solutions to $u$ and $v$ respectively. Then,

[TABLE]

where $\delta_{m}$ and $\phi_{m}$ are time dependent parameters to be determined from the collocation points $x_{m},m=0,...,N$ and the complementary data. The nodal values $U$ and its first and second derivatives at the knots can be found from the (5) as

[TABLE]

When Gardner equation is space-splitted as (2), the resultant system has maximum second order derivatives that enables to construct a smooth cubic B-spline approximation with the exponential B-splines. Implementation of the Crank-Nicolson method to the space-slitted system (2) yields

[TABLE]

where $t^{n+1}=t^{n}+\Delta t$, and the superscripts $n$ and $n+1$ denote $(n)$th and $(n+1)$th time levels, respectively.

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

[TABLE]

and

[TABLE]

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

[TABLE]

where

[TABLE]

[TABLE]

The system (15) can be rewritten in the matrix form as

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The system (16) consists of $2N+2$ linear equations and $2N+6$ unknown parameters

[TABLE]

Additional four constraints are required to have the unique solution to this system. These constraints can be determined by the imposition of the boundary data

[TABLE]

to give

[TABLE]

When the unknowns $\delta_{-1},\phi_{-1},\delta_{N+1},\phi_{N+1},$ are eliminated from the system (16), the resultant system is a solvable system of $2N+2$ linear equations with $2N+2$ unknowns. The algorithm derived from the classical Thomas algorithm is used to solve (16) in each time level.

In order to start the iteration algorithm, the parameters $\delta_{m}^{0},\phi_{m}^{0},$ $m=-1,\ldots,N+1$ are required to be determined from

[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 18 into the system

[TABLE]

[TABLE]

gives

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

4 Illustrations

This section contains implementation of the proposed algorithm to some IVPs to validate its accuracy and efficiency. The accuracy of the numerical results is measured by using the discrete maximum norm defined as

[TABLE]

at the time $t$ when the analytical solution exists.

The conservation laws also validates the accuracy of the proposed algorithms if they keep their initial values even in the nonexistence of the analytical solutions case. The lowest three conservation laws

[TABLE]

are expected to keep their initial values as time proceeds[19]. In order to measure the absolute relative changes of these quantities at any time $t>0$, $C(M_{t})$, $C(E_{t})$ and $C(H_{t})$ are defined as

[TABLE]

where $M_{0}$, $E_{0}$ and $H_{0}$ are initial, $M_{t}$, $E_{t}$ and $H_{t}$ are the quantities at the time $t>0$.

4.1 Example 1

The solution representing the propagation of an initial positive pulse is demonstrated by using the initial condition

[TABLE]

and the homogeneous Neumann conditions at both ends of the problem interval $[-20,30]$. The analytical solution can be written as

[TABLE]

when the compatible parameters are chosen $\mu_{1}=4,$ $\mu_{2}=-3$ and $\mu_{3}=1$ in the Gardner equation (1). This solution represents propagation of a positive initial pulse along the $x-$axis as time proceeds. The simulation of the solution is depicted in Fig 1 in the time domain $[0,5]$

Various discretization parameters $h$ and $\Delta t$ are used to illustrate the numerical solutions. The optimum value of exponential spline parameter $\zeta$ is scanned by comparing the discrete maximum error norms at the simulation terminating time. The maximum values of the errors are measured near the peak points as expected for both $\zeta=1$ and $\zeta=0.000003$. The maximum error distribution for $\zeta=1$ and $\zeta=0.000003$ are depicted in Fig 2(a) and Fig 2(b), respectively, for $h=0.5$ and $\Delta t=0.1$.

The discrete maximum error norms at some distinct times are reported in Table 2 for various values of the discretization parameter $h$ and fixed $\Delta t$. Even though the solutions are improved when the spatial discretization parameter size is reduced when $\zeta=1$ at both $t=2.5$ and $t=5$. For various values of the spatial discretization parameter, the optimum exponential B-spline parameter is scanned to reduce the maximum error. When $h$ is larger, the optimum $\zeta$ improves the results by decreasing the error approximately one third or one fourth. Even though the better results are obtained for smaller $h$ values, the optimum choice of $\zeta$ increases the accuracy two times.

The initial values of the conservation laws are computed with computer algebra tools and reported in Table 3. A stable numerical method is expected to preserve these values as time goes. The absolute relative changes of all related conservation laws reported in Table 3. They all are preserved at least six decimal digits for various values of $h$ and $\zeta=1$. The preservation of conservation laws is an indicator of a reliable numerical approach.

4.2 Example 2

The solution of the Gardner equation (1) representing the motion of kink like wave is of the form

[TABLE]

when the equation parameters are chosen as $\mu_{1}=1,$ $\mu_{2}=-5$ and $\mu_{3}=1$. The initial condition required for the numerical simulation is derived from the analytical solution by substituting $t=0$ in the analytical solution. The homogeneous Neumann boundary conditions are used to complement the problem statement. The designed algorithm is run up to the terminating time $t=12$ in the finite interval $[-80,80]$. The results are graphed in Fig 3.

The maximum error distributions for $\zeta=1$ and $\zeta=0.000001$ are depicted in Fig 4(a) and Fig 4(b), respectively, for $N=100$. The errors are larger around $x=0$ as expected due to the shape of the kink like wave.

The maximum error norms are calculated for various values of $N$ for $\zeta=1$. In order to improve the results, the optimum value of the exponential B-spline parameter are investigated by scanning. Determination of optimum values of $rho$ different from $1$ gives at least one decimal digit better results, Table 4. Even though the results improve when $N$ increases, optimum values of $\zeta$ improves results independent on the discretization parameters.

The initial values and absolute relative changes of the conservation laws are tabulated in Table 5. Absolute relative changes of all the conservation laws are calculated in four decimal digits independently on number of spatial discretization point number, Table 5.

4.3 Example 3

The wave generation from an initial positive pulse is studied as the last example. The initial condition used in the first example is perturbed carefully to generate new waves. Thus, an initial condition is produced for $\mu_{1}=10,$ $\mu_{2}=-3$ and $\mu_{3}=1$ as

[TABLE]

This initial pulse is expected to generate new pulses behind as propagating towards to the right along the horizontal axis. The numerical simulation is accomplished with the discretization parameters $h=0.5$ and $\Delta t=0.1$ in the interval $[-40,60]$. The simulation is depicted in Fig 5(a) - 5(d).

The initial values and the absolute relative changes of the conservation laws are recorded during the simulation process duration in Table 6. It is clearly observable that all the conservation laws are preserved successfully at least three decimal digits in the simulation period.

5 Conclusion

The exponential B-spline collocation method is implemented for IBVPs with analytical and non analytical solutions. Having no continuous third order derivative of exponential B-splines forces us to reduce the order of third order derivative term. Thus, the Gardner equation reduces to a coupled nonlinear PDE system with maximum second order derivative term. The exponential B-spline approach to the solution of this system is substituted into the system following the Crank-Nicolson implicit time integration method. Linearization of both non linear terms with quadratic and cubic nonlinearities is followed by implementation of boundary data to obtain a solvable system of equation having equal numbers of equations and unknowns. The algorithm gets ready to run after arranging the initial state by the aid of initial and boundary data.

The validity of the method is checked by computing the error between the analytical and numerical solutions in the case of having an analytical solution. The absolute relative changes of the conservation laws can also be a good indicator to observe the validity and accuracy of the proposed method even when there exits no analytical solution. The scan to determine the optimum exponential B-spline parameter value shows that the accuracy of the solutions can be improved when compared the choice $\zeta=1$.

6 Acknowledgements

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