Using a cubic B-spline method in conjunction with a one-step optimized hybrid block approach to solve nonlinear partial differential equations

Abstract

In this paper, we develop an optimized hybrid block method which is combined with a modified cubic B-spline method, for solving non-linear partial differential equations. In particular, it will be applied for solving three well-known problems, namely, the Burgers equation, Buckmaster equation and FitzHugh–Nagumo equation. Most of the developed methods in the literature for non-linear partial differential equations have not focused on optimizing the time step-size and a very small value must be considered to get accurate approximations. The motivation behind the development of this work is to overcome this trade-off up to much extent using a larger time step-size without compromising accuracy. The optimized hybrid block method considered is proved to be A-stable and convergent. Furthermore, the obtained numerical approximations have been compared with exact and numerical solutions available in the literature and found to be adequate. In particular, without using quasilinearization or filtering techniques, the results for small viscosity coefficient for Burgers equation are found to be accurate. We have found that the combination of the two considered methods is computationally efficient for solving non-linear PDEs.

1 Introduction

The non-linear PDE model equations arise very frequently in various areas of physical, chemical and biological sciences. The study of these types of equations has proceeded at a very rapid pace. Finding analytical solutions of most of the non-linear problems is a challenge for researchers. Therefore, a numerical solution is the rescuer in the situation, where it is impossible to find an analytical solution of the problem. In this paper, three PDEs containing non-linear terms are considered as they are famous for their numerous practical applications:

• (i)  Burgers equation,

• (ii) Buckmaster equation,

• (iii) FitzHugh–Nagumo equation.

The Burgers equation involves both non-linear propagation effects and diffusive effects. This equation is similar to the Navier–Stokes equation without the pressure term. Therefore, it is a simpler model to analyze fluid turbulence. The Burgers equation has applications in various fields, such as statistical physics, fluid dynamics, gas dynamics, shock theory, viscous flow and turbulence, quantum fields, traffic flow, cosmology, etc. In the past years, the Burgers equation has been studied by various numerical methods based on finite element method, finite difference methods, differential transform method, the homotopy analysis method, quadrature method, Haar wavelet quasilinearization approach, splitting methods, etc. (Tamsir et al. 2016; Seydaoğlu 2018; Jiwari 2015, 2012; Kadalbajoo and Awasthi 2006; Kutluay et al. 2004; Mittal and Jain 2012; Liao 2008; Öziş and Erdoğkan 2009; Rashidi and Erfani 2009).

The Buckmaster equation is a parabolic non-linear partial differential equation with two non-linear terms. This equation is considered as a model of thin viscous sheet flow and is also used for the purpose of describing large scale and long term deformation. The Buckmaster equation has been solved numerically using finite volume method (Hussain and Alwan 2013) and B-splines collocation method (Chanthrasuwan et al. 2017).

The FitzHugh–Nagumo equation (FHN) is another important non-linear reaction–diffusion equation. The FHN equation is used in the fields of neurophysiology, flame propagation, logistic population growth, nuclear reactor theory, branching brownian motion processes and catalytic chemical reactions (Bhrawy 2013; Wazwaz and Gorguis 2004). This equation has been solved using the Haar wavelet method (Hariharan and Kannan 2010), the homotopy analysis method (Abbasbandy 2008), a polynomial quadrature method (Jiwari et al. 2014), the local meshless method (Ahmad et al. 2019) and the collocation method using cubic B-splines with SSP-RK54 (Mittal and Tripathi 2015). Analogously, there are numerous applications of their corresponding fractional models existing in the literature (Onal and Esen 2020; Li et al. 2016; Lozi et al. 2020).

In 1953, Milne (Milne and Milne 1953) proposed block methods for the first time to enhance the performance of the numerical methods. Later, to provide the essential starting values needed for the predictor schemes, Rosser (Rosser 1967) considered block methods too. Furthermore, Gragg and Stetter (Gragg and Stetter 1964) developed hybrid methods that incorporate information about the solution at off-step points. In addition, the combination of hybrid and block method has many advantages, such as the ability to evade the Dahlquist’s barrier (Dahlquist 1956) using information at off-step points, the ability to change the step size, simultaneously producing approximate solutions at several points and being self starting. Block methods have the advantage over Runge–Kutta methods that they are less expensive in terms of number of function evaluations for a given order. In recent years, work on block methods (Singh et al. 2019; Ramos and Singh 2017; Singla et al. 2021) has drawn attention of many researchers.

In this paper, we derive a one-step hybrid block method to solve first order initial value problems and is further applied to a system of first order differential equations obtained after applying a modified cubic B-splines collocation method (Mittal and Jain 2012) to PDE for space derivatives. B-splines have been chosen as basis functions, because they lead to band matrices which involve low computational cost and less storage space. The developed one-step optimized hybrid block method considers two intra-step points to improve accuracy. The method proposed in this paper has the advantages that provides accurate results by taking few time steps as compared to numerical methods present in the literature, and produces efficient results even for small values of the kinematic coefficient in the case of the Burgers equation.

2 Development of hybrid block method

Consider an ordinary differential equation of first order

$$\begin{aligned} y^\prime =f(t,y), \text {where }t\in [t_0,t_n], \end{aligned}$$

(1)

with initial condition $y(t_0)=y_{0}$. We consider the grid points $t_0< t_1<\dots <t_n$, such that the fixed step size is $k=t_{i+1}-t_{i}, \forall i=0, 1,\dots , n-1$. Let us consider the following polynomial p(t) to provide an approximate solution of (1), that is

$$\begin{aligned}&y(t)\approx p(t)=\sum _{j=0}^{4}a_jt^j , \end{aligned}$$

(2)

$$\begin{aligned}&y^\prime (t)\approx p^\prime (t)=\sum _{j=1}^{4}ja_jt^{j-1} , \end{aligned}$$

(3)

where the $a_j\in \mathbf {R}$ are unknown coefficients.

By introducing two intra-step points $t_{i+r}=t_i+rk$ and $t_{i+s}=t_i+sk$, where $0<r<s<1$, we interpolate (2) at $t_i$ and collocate (3) at $t_i, t_{i+r}, t_{i+s}$ and $t_{i+1}$, resulting into a system of equations given in the following matrix form:

$$\begin{aligned} \begin{bmatrix} 1&{}t_i&{}t_{i}^2&{}t_{i}^3&{}t_{i}^4\\ 0&{}1&{}2t_i&{}3t_{i}^2&{}4t_{i}^3\\ 0&{}1&{}2t_{i+r}&{}3t_{i+r}^2&{}4t_{i+r}^3\\ 0&{}1&{}2t_{i+s}&{}3t_{i+s}^2&{}4t_{i+s}^3\\ 0&{}1&{}2t_{i+1}&{}3t_{i+1}^2&{}4t_{i+1}^3\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ a_2\\ a_3\\ a_4\\ \end{bmatrix}= \begin{bmatrix} y_i\\ f_{i}\\ f_{i+r}\\ f_{i+s}\\ f_{i+1}\\ \end{bmatrix}, \end{aligned}$$

where $f_l=f(t_l,y_l)$, for $l=i, i+r, i+s, i+1$.

We solve this system to get the values of the coefficients $a_j$’s, then substituting the obtained values of $a_j$’s in

$$\begin{aligned} p(t_i+mk)=\sum _{j=0}^{4}a_j(t_i+mk)^j , \end{aligned}$$

after rewriting it we obtain the following

$$\begin{aligned} p(t_i +mk)=b_0y_i +k(b_1f_i+b_2f_{i+r}+b_3f_{i+s}+b_4f_{i+1}) , \end{aligned}$$

(4)

where the $b_j$’s are continuous coefficients that depend on m. Now, taking $m=1$ in (4), we get a formula that approximates the true solution at $t_{i+1}$, which is denoted as $y_{i+1}$. We expand $y(t_{i+1})$ using the Taylor series about $t_i$ to find the values of r and s for which the local truncation error $L(y(t_{i+1}), k)$ is minimized. This local truncation error is given by

$$\begin{aligned} L(y(t_{i+1}),k)=&\frac{k^5y^{(5)}(t_{i})}{5!}\left( 1+5\frac{(r+s-2rs-3)}{12}\right) \nonumber \\&+\frac{k^6y^{(6)}(t_{i})}{6!}\left( 1+\frac{-2r^2s+r^2-2rs^2-rs+r+s^2+s-3}{2}\right) \nonumber \\&+O(k^7). \end{aligned}$$

(5)

Since two more equations are needed to find unknown values of r and s, we equate the coefficients of $k^5$ and $k^6$ to zero in (5). Therefore, one gets

$$\begin{aligned} r=\frac{1}{2}-\frac{\sqrt{5}}{10},s=\frac{1}{2}+\frac{\sqrt{5}}{10}. \end{aligned}$$

(6)

Substituting the obtained values in (6) in the local truncation error formula (5), we obtain

$$\begin{aligned} L(y(t_{i+1}),k)=-\frac{k^7y^{(7)}(t_{i})}{1512000}+O(k^8). \end{aligned}$$

Now, the formulas of the block method are deduced by substituting the values of r and s from (6) in (4) for $m=r, s$ and 1, respectively, giving rise to a system of three equations as given below:

$$\begin{aligned} y_{i+r}=y_{i}+\frac{k}{12}&\Bigg [\left( \frac{11}{10}+\frac{\sqrt{5}}{10}\right) f_i+\left( \frac{5}{2}-\frac{\sqrt{5}}{10}\right) f_{i+r} +\left( \frac{5}{2}-13\frac{\sqrt{5}}{10}\right) f_{i+s}\nonumber \\&\quad +\left( \frac{-1}{10}+\frac{\sqrt{5}}{10}\right) f_{i+1}\Bigg ], \end{aligned}$$

(7a)

$$\begin{aligned} y_{i+s}=y_{i}+\frac{k}{12}&\Bigg [\left( \frac{11}{10}-\frac{\sqrt{5}}{10}\right) f_i+\left( \frac{5}{2}+13\frac{\sqrt{5}}{10}\right) f_{i+r} +\left( \frac{5}{2}-\frac{\sqrt{5}}{10}\right) f_{i+s}\nonumber \\&\quad +\left( \frac{-1}{10}-\frac{\sqrt{5}}{10}\right) f_{i+1}\Bigg ], \end{aligned}$$

(7b)

$$\begin{aligned}&y_{i+1}=y_{i}+\frac{k}{12}\left( f_i+5f_{i+r}+5f_{i+s}+f_{i+1}\right) . \end{aligned}$$

(7c)

3 Implementation of the cubic B-spline hybrid block method

Concerning the cubic B-spline collocation, let us consider $N+1$ grid points $a=x_0<x_1<x_2<\dots<x_{N-1}<x_N=b$, as uniform mesh with step size $h=x_{i+1}-x_i=\frac{b-a}{N}$ along the x-direction. Our goal is to approximate the solution u(x, t) of a given PDE problem by a B-spline function U(x, t) given in the form:

$$\begin{aligned} U(x,t)=\sum _{i=-1}^{N+1}\gamma _i(t)\mathbb {B}_i(x), \end{aligned}$$

(8)

where $\gamma _i(t)$ are time dependent coefficients to be determined by imposing boundary and collocation conditions. The cubic B-spline basis functions at the set of given grid points are defined as

$$\begin{aligned} \mathbb {B}_i(x)=\frac{1}{h^3}\left\{ \begin{array}{cc} (x-x_{i-2})^3 &{} x\in [x_{i-2},x_{i-1}),\\ (x-x_{i-2})^3-4(x-x_{i-1})^3 &{} x\in [x_{i-1},x_i),\\ (x_{i+2}-x)^3-4(x_{i+1}-x)^3 &{} x\in [x_{i},x_{i+1}),\\ (x_{i+2}-x)^3 &{} x\in [x_{i+1},x_{i+2}),\\ 0 &{} \text {otherwise},\\ \end{array} \right. \end{aligned}$$

where $\{\mathbb {B}_{-1}(x), \mathbb {B}_0(x),..., \mathbb {B}_{N+1}(x)\}$ form a basis over the domain $a\le x\le b$. The values of the cubic B-splines and their derivatives at the mesh points are given in Table 1.

Table 1 Values of $\mathbb {B}_i(x)$ and its derivatives

In accordance with the fact that the number of basis functions in a collocation procedure should equal to the number of collocation points, so, there is a need to use a set of modified basis functions which are defined as (Mittal and Jain 2012):

$$\begin{aligned} \tilde{\mathbb {B}_{i}}(x)=\left\{ \begin{aligned}&\mathbb {B}_0(x)+2\mathbb {B}_{-1}(x)&i=0\\&\mathbb {B}_1(x)-\mathbb {B}_{-1}(x)&i=1\\&\mathbb {B}_i(x)&i=2,3,...,N-2\\&\mathbb {B}_{N-1}(x)-\mathbb {B}_{N+1}(x)&i=N-1\\&\mathbb {B}_{N}(x)+2\mathbb {B}_{N+1}(x)&i=N, \end{aligned} \right. \end{aligned}$$

(9)

where $\{\tilde{\mathbb {B}_0}(x) ,\tilde{\mathbb {B}_1}(x),...,\tilde{\mathbb {B}_N}(x)\}$ form a basis in the domain $a\le x\le b$ and ensures that the system that discretizes the differential problem with Dirichlet boundary conditions has a diagonally dominant matrix. In fact, using this modified basis, we consider that the approximate solution may be expressed as

$$\begin{aligned} U(x,t)=\sum _{j=0}^{N}\alpha _j(t)\tilde{\mathbb {B}_j}(x). \end{aligned}$$

(10)

Thus, for a PDE problem of the form

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial u}{\partial t}=F\left( t,x, u,\frac{\partial u}{\partial x},\frac{\partial ^2 u}{\partial x^2},\dots \right) \\&u(x,0)=\psi (x),a\le x\le b\\&u(a,t)=g_1(t),u(b,t)=g_2(t),t\ge 0, \end{aligned} \right. \end{aligned}$$

(11)

one gets the approximations at the grid points in the following form:

$$\begin{aligned} \begin{aligned}&U(x_0,t)=g_1(t),\\&U(x_i,t)=\sum _{j=0}^{N}\alpha _j(t)\tilde{\mathbb {B}_j}(x_i),\quad i=1,\dots ,N-1,\\&U(x_N,t)=g_2(t), \end{aligned} \end{aligned}$$

and for the first derivatives with respect to time, we have

$$\begin{aligned}\begin{aligned}&U_t(x_0,t)=\dot{g_1}(t),\\&U_t(x_i,t)=\dot{\alpha }_{i-1}(t)+4\dot{\alpha }_{i}(t)+\dot{\alpha }_{i+1}(t),\quad i=1,\dots ,N-1,\\&U_t(x_N,t)=\dot{g_2}(t), \end{aligned}\end{aligned}$$

where the dot corresponds to the derivative with respect to the time variable t. After applying the modified cubic B-splines procedure, we get a system of first-order differential equations in the variables $\left\{ \alpha _0(t),\alpha _1(t),\dots ,\alpha _N(t) \right\}$, which will be solved using the block method (7).

Using the above strategies, we have implemented a combined B-spline hybrid block method (which will be named in short OHBCBM) on three well-known problems.

3.1 Burgers equation

Consider the one dimensional Burgers equation

$$\begin{aligned} \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu \frac{\partial ^2 u}{\partial x^2}, \end{aligned}$$

(12a)

with initial condition

$$\begin{aligned} u(x,0)=\psi (x),a\le x\le b, \end{aligned}$$

(12b)

and boundary conditions

$$\begin{aligned} u(a,t)=g_1(t),u(b,t)=g_2(t),t\ge 0, \end{aligned}$$

(12c)

where u, x, t and $\nu$ are the velocity, space variable, time and kinematic viscosity, respectively. Now, the collocation technique is applied to the Burgers equation using the modified cubic B-splines (9) and considering the $\alpha _j$’s as time-dependent coefficients (henceforth we write $\alpha _j$ instead of $\alpha _j(t)$). We obtain the following differential system

• For $j=0$, $6\dot{\alpha _0}=\dot{g_1}(t)$,

• For $j=N$, $6\dot{\alpha _N}=\dot{g_2}(t)$,

• For $j=1,2,...,N-1$,

$$\begin{aligned}&\dot{\alpha }_{j-1}+4\dot{\alpha }_{j}+\dot{\alpha }_{j+1}=\nu \frac{6}{h^2}\left( \alpha _{j-1}-2\alpha _j+\alpha _{j+1}\right) \\&\quad -\frac{3}{h}\left( \alpha _{j-1}+4\alpha _j+\alpha _{j+1}\right) \left( \alpha _{j+1}-\alpha _{j-1}\right) . \end{aligned}$$

The above system can be written in matrix notation as

$$\begin{aligned} \begin{bmatrix} 6&{}0&{}0&{} &{} &{} \\ 1&{}4&{}1&{} &{} &{} \\ &{}..&{}..&{}..&{} &{} \\ &{} &{}..&{}..&{}..&{} \\ &{} &{} &{}1&{}4&{}1 \\ &{} &{} &{} 0&{}0&{}6\\ \end{bmatrix} \begin{bmatrix} \dot{\alpha }_0\\ \dot{\alpha }_1\\ ...\\ ...\\ \dot{\alpha }_{N-1}\\ \dot{\alpha }_{N}\\ \end{bmatrix}=\frac{6\nu }{h^2} \begin{bmatrix} 0&{}0&{} 0&{} &{} &{} \\ 1&{}-2&{}1&{} &{} &{} \\ &{}..&{}..&{}..&{} &{} \\ &{} &{}..&{}..&{}..&{} \\ &{} &{} &{}1&{}-2&{}1\\ &{} &{} &{} 0&{}0&{}0\\ \end{bmatrix} \begin{bmatrix} \alpha _0\\ \alpha _1\\ ...\\ ...\\ \alpha _{N-1}\\ \alpha _{N}\\ \end{bmatrix}+C \end{aligned}$$

where $C= \begin{bmatrix} \dot{g_1}\\ -\frac{3}{h}(\alpha _0+4\alpha _1+\alpha _2)(\alpha _2-\alpha _0)\\ .....\\ .....\\ -\frac{3}{h}(\alpha _{N-2}+4\alpha _{N-1}+\alpha _{N})(\alpha _{N}-\alpha _{N-2})\\ \dot{g_2}\\ \end{bmatrix}.$

It can also be written in the following compact form:

$$\begin{aligned} A\dot{\alpha }=B\alpha +C, \end{aligned}$$

where A and B are the obvious $(N+1)\times (N+1)$ tridiagonal matrices.

3.2 Buckmaster equation

Consider the Buckmaster equation

$$\begin{aligned} \frac{\partial u}{\partial t}=\frac{\partial ^2 (u^4)}{\partial x^2}+\frac{\partial (u^3)}{\partial x}+f(x,t), \end{aligned}$$

(13a)

with initial condition

$$\begin{aligned} u(x,0)=\psi (x),a<x<b, \end{aligned}$$

(13b)

and boundary conditions

$$\begin{aligned} u(a,t)=g_1(t),u(b,t)=g_2(t),t\ge 0, \end{aligned}$$

(13c)

where x and t are the space and time variables, respectively.

Now, the B-spline collocation method is applied to the Buckmaster equation using (9) and the values in Table (1), and obtain the following system:

• For $j=0$, $6\dot{\alpha _0}=\dot{g_1}(t)$,

• For $j=N$, $6\dot{\alpha _N}=\dot{g_2}(t)$,

• For $j=1,2,...,N-1$,

$$\begin{aligned} \dot{\alpha }_{j-1}+4\dot{\alpha }_{j}+\dot{\alpha }_{j+1}=&4(\frac{6}{h^2}[\alpha _{j-1}-2\alpha _j+\alpha _{j+1}][\alpha _{j-1}+4\alpha _j+\alpha _{j+1}]^3)\\&+12(\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1})^2[\frac{3}{h}(\alpha _{j+1}-\alpha _{j-1})]^2\\&+3(\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1})^2[\frac{3}{h}(\alpha _{j+1}-\alpha _{j-1})]+f(x_j,t). \end{aligned}$$

The above system can be written as

$$\begin{aligned} \begin{bmatrix} 6&{}0&{}0&{} &{} &{} \\ 1&{}4&{}1&{} &{} &{} \\ &{}..&{}..&{}..&{} &{} \\ &{} &{}..&{}..&{}..&{} \\ &{} &{} &{}1&{}4&{}1 \\ &{} &{} &{} 0&{}0&{}6\\ \end{bmatrix} \begin{bmatrix} \dot{\alpha }_0\\ \dot{\alpha }_1\\ ...\\ ...\\ \dot{\alpha }_{N-1}\\ \dot{\alpha }_{N}\\ \end{bmatrix}= B\alpha +C, \end{aligned}$$

where

$$\begin{aligned} {\begin{matrix} C_{1,1}=&\dot{g_1}, \end{matrix}}\\ {\begin{matrix} C_{j,1}=&{}4(\frac{6}{h^2}[\alpha _{j-1}-2\alpha _j+\alpha _{j+1}][\alpha _{j-1}+4\alpha _j+\alpha _{j+1}]^3)\\ &{}+12(\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1})^2[\frac{3}{h}(\alpha _{j+1}-\alpha _{j-1})]^2\\ &{}+3(\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1})^2[\frac{3}{h}(\alpha _{j+1}-\alpha _{j-1})]+f(x_j,t), \end{matrix}}\\ {\begin{matrix} C_{N,1}=&\dot{g_2}. \end{matrix}} \end{aligned}$$

Also, the above system can be rewritten in the following compact way:

$$\begin{aligned} A\dot{\alpha }=B\alpha +C, \end{aligned}$$

where A and B are the obvious $(N+1)\times (N+1)$ tridiagonal matrices. Note that B here is a null matrix as the Buckmaster equation does not have linear terms on the right hand side.

3.3 FitzHugh–Nagumo equation

Consider the FitzHugh–Nagumo equation

$$\begin{aligned} \frac{\partial u}{\partial t}=\frac{\partial ^2 u}{\partial x^2}+u(u-\mu )(1-u), \end{aligned}$$

(14a)

with initial conditions

$$\begin{aligned} u(x,0)=\psi (x),a<x<b, \end{aligned}$$

(14b)

and boundary conditions

$$\begin{aligned} u(a,t)=g_1(t), u(b,t)=g_2(t), t\ge 0, \end{aligned}$$

(14c)

where x and t are the space and time variables, respectively.

Now, the B-spline collocation method is applied to this equation using (9) and the values from Table 1 which results in the following system:

• For $j=0$; $6\dot{\alpha _0}=\dot{g_1}(t)$,

• For $j=N$; $6\dot{\alpha _N}=\dot{g_2}(t)$,

• For $j=1,2,...,N-1$,

$$\begin{aligned} \dot{\alpha }_{j-1}+4\dot{\alpha }_{j}+\dot{\alpha }_{j+1}=&\frac{6}{h^2}[\alpha _{j-1}-2\alpha _j+\alpha _{j+1}]\\&-\mu [\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1}]+[\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1}]^2\\&+\mu [\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1}]^2-[\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1}]^3. \end{aligned}$$

The above system can be written using the matrix notation as

$$\begin{aligned}&\begin{bmatrix} 6&{}0&{}0&{} &{} &{} \\ 1&{}4&{}1&{} &{} &{} \\ &{}..&{}..&{}..&{} &{} \\ &{} &{}..&{}..&{}..&{} \\ &{} &{} &{}1&{}4&{}1 \\ &{} &{} &{} 0&{}0&{}6\\ \end{bmatrix} \begin{bmatrix} \dot{\alpha }_0\\ \dot{\alpha }_1\\ ...\\ ...\\ \dot{\alpha }_{N-1}\\ \dot{\alpha }_{N}\\ \end{bmatrix}\\&\quad = \begin{bmatrix} 0&{}0&{} 0&{} &{} &{} \\ \frac{6}{h^2}-\mu &{}\frac{-12}{h^2}-4\mu &{}\frac{6}{h^2}-\mu &{} &{} &{} \\ &{}..&{}..&{}..&{} &{} \\ &{} &{}..&{}..&{}..&{} \\ &{} &{} &{}\frac{6}{h^2}-\mu &{}\frac{-12}{h^2}-4\mu &{}\frac{6}{h^2}-\mu \\ &{} &{} &{} 0&{}0&{}0\\ \end{bmatrix} \begin{bmatrix} \alpha _0\\ \alpha _1\\ ...\\ ...\\ \alpha _{N-1}\\ \alpha _{N}\\ \end{bmatrix}+C, \end{aligned}$$

where

$$\begin{aligned} {\begin{matrix} C_{1,1}=&\dot{g_1}, \end{matrix}}\\ {\begin{matrix} C_{j,1}=&[\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1}]^2(1+\mu )-[\alpha _{j-1}+4\alpha _{j}+\alpha _{j+1}]^3, \end{matrix}}\\ {\begin{matrix} C_{N,1}=&\dot{g_2}. \end{matrix}} \end{aligned}$$

Also, the above system can be rewritten in the following compact way:

$$\begin{aligned} A\dot{\alpha }=B\alpha +C, \end{aligned}$$

where A and B are the obvious $(N+1)\times (N+1)$ tridiagonal matrices.

Note that in all the above examples, the form of the differential system is $A\dot{\alpha }=B\alpha +C,$ with the matrix A being the same, which is invertible, and C a $N+1$-vector containing boundary values and nonlinear terms. Thus, the system of first-order differential equations may be written as

$$\begin{aligned} \dot{\alpha }=A^{-1}B\alpha +A^{-1}C. \end{aligned}$$

(15)

Note that the necessary initial values to get a unique solution of the system can be readily obtained through the initial condition in (11), from which we have

$$\begin{aligned}\begin{aligned}&U(x_0,0)=6{\alpha }_{0}(0)=\psi (x_0),\\&U(x_i,0)={\alpha }_{i-1}(0)+4{\alpha }_{i}(0)+{\alpha }_{i+1}(0)=\psi (x_i),\quad i=1,\dots ,N-1,\\&U(x_N,0)=6{\alpha }_{N}(0)=\psi (x_N), \end{aligned}\end{aligned}$$

which in fact results in an algebraic tridiagonal system of $(N+1)\times (N+1)$ equations that may be written as

$$\begin{aligned} A \alpha (0)=b, \end{aligned}$$

where

$$\begin{aligned}\begin{aligned}&\alpha (0)=\left( {\alpha }_{0}(0),{\alpha }_{1}(0),\dots ,{\alpha }_{N}(0) \right) ^T,\\&b=\left( \psi (x_0),\psi (x_1),\dots ,\psi (x_N) \right) ^T. \end{aligned}\end{aligned}$$

This system can be solved using the Thomas’ algorithm and provides the initial values for the system in (15).

By applying the hybrid block method in (7) to the system (15), we get

$$\begin{aligned} \alpha ^{i+r}=\alpha ^{i}+\frac{k}{12}&\Bigg [\left( \frac{11}{10}+\frac{\sqrt{5}}{10}\right) \bar{f}_i+\left( \frac{5}{2}-\frac{\sqrt{5}}{10}\right) \bar{f}_{i+r}\nonumber \\&\quad +\left( \frac{5}{2}-13\frac{\sqrt{5}}{10}\right) \bar{f}_{i+s}+\left( \frac{-1}{10}+\frac{\sqrt{5}}{10}\right) \bar{f}_{i+1}\Bigg ], \end{aligned}$$

(16a)

$$\begin{aligned} \alpha ^{i+s}=\alpha ^{i}+\frac{k}{12}&\Bigg [\left( \frac{11}{10}-\frac{\sqrt{5}}{10}\right) \bar{f}_i+\left( \frac{5}{2}+13\frac{\sqrt{5}}{10}\right) \bar{f}_{i+r}\nonumber \\&\quad +\left( \frac{5}{2}-\frac{\sqrt{5}}{10}\right) \bar{f}_{i+s}+\left( \frac{-1}{10}-\frac{\sqrt{5}}{10}\right) \bar{f}_{i+1}\Bigg ], \end{aligned}$$

(16b)

$$\begin{aligned}&\alpha ^{i+1}=\alpha ^{i}+\frac{k}{12}\left( \bar{f_i}+5\bar{f}_{i+r}+5\bar{f}_{i+s}+\bar{f}_{i+1}\right) , \end{aligned}$$

(16c)

where

$$\begin{aligned} \bar{f}_i=A^{-1}B\alpha ^i+A^{-1}C_i, \end{aligned}$$

$$\begin{aligned} \bar{f}_{i+r}=A^{-1}B\alpha ^{i+r}+A^{-1}C_{i+r}, \end{aligned}$$

$$\begin{aligned} \bar{f}_{i+s}=A^{-1}B\alpha ^{i+s}+A^{-1}C_{i+s}, \end{aligned}$$

$$\begin{aligned} \bar{f}_{i+1}=A^{-1}B\alpha ^{i+1}+A^{-1}C_{i+1}. \end{aligned}$$

Here, $\alpha ^m$ and $C_m$ are values of $\alpha$ and C at time $t_m$, respectively. For each time step, the system (16) has been solved.

4 Characteristics of the hybrid block method and stability analysis

4.1 Zero-stability

The hybrid block method (7) is said to be zero-stable if the roots of its first characteristic equation $\rho (R)=0$ have modulus $|R_m|\le 1$ and if $|R_m|$=1 then the multiplicity does not exceed 1. By taking limit $k\rightarrow 0$, the proposed method is written in matrix form as

$$\begin{aligned} IY_{\eta }-BY_{\eta -1}=0 \end{aligned}$$

where

$$\begin{aligned} Y_{\eta }=(y_{i+1},y_{i+r},y_{i+s})^T ,Y_{\eta -1}=(y_{i},y_{i+r-1},y_{i+s-1})^T, \end{aligned}$$

$$\begin{aligned} B= \begin{bmatrix} 1&{}0&{}0\\ 1&{}0&{}0\\ 1&{}0&{}0\\ \end{bmatrix}, \end{aligned}$$

and I is the identity matrix of order three. Therefore, the characteristic equation is $\rho (R)=det[RI-B]=R^2(R-1)=0$. The roots are $\left\{ 0, 0 , 1\right\}$, and thus, the proposed method is zero-stable.

4.2 Order of the hybrid block method

Let z(t) be a sufficiently differentiable function and $\tau$ be the linear difference operator associated with the hybrid block method formula in (7):

$$\begin{aligned} \tau [z(t);k]=E_1Z-kE_2F, \end{aligned}$$

where $E_1$ and $E_2$ are matrices of coefficients of dimensions $3\times 4$, and $Z=[z(t),z(t+rk),z(t+sk),z(t+k)]^T$, $F=[z'(t),z'(t+rk),z'(t+sk),z'(t+k)]^T$. Using Taylor series, we expand $\tau [z(t);k]$ resulting that

$$\begin{aligned} \tau [z(t);k]=C_0z(t)+C_1kz^\prime (t)+C_2k^2z^{\prime \prime }(t)+...+C_pk^pz^{(p)}(t)+... , \end{aligned}$$

with $C_0=C_1=C_2=...=C_p=0$ and $C_{p+1}\ne 0$. For the proposed method it is found that $C_0=C_1=...=C_4=0$ and

$$\begin{aligned} C_5=[-\frac{\sqrt{5}}{30000},\frac{\sqrt{5}}{30000},0]^T. \end{aligned}$$

Thus, the derived hybrid block method (7) has at least order $p=4$ which implies that the method is consistent.

4.3 Convergence

A linear multi-step method is convergent if it is both consistent and zero stable (PK 1962; Lambert 1973). Therefore, the proposed hybrid block method (7) is convergent.

4.4 Stability analysis

First, to analyze the linear stability of the hybrid block method we consider the Dahlquist’s test equation

$$\begin{aligned} y^\prime (t)=\lambda y(t),\text {where }Re(\lambda )<0. \end{aligned}$$

(17)

The general solution of the above equation is $y(t)=e^{\lambda t}$, which will be damped out as $t\rightarrow \infty$. The numerical method is said to be linearly stable if the numerical solution it provides has a similar qualitative behavior as the theoretical solution. To determine the region in which the numerical method exhibits a similar behavior, the proposed block method in (7) is applied to (17), thereby resulting in

$$\begin{aligned} A \begin{bmatrix} y_{i+r}\\ y_{i+s}\\ y_{i+1}\\ \end{bmatrix} =B \begin{bmatrix} y_{i+r-1}\\ y_{i+s-1}\\ y_{i}\\ \end{bmatrix}, \end{aligned}$$

where

$$\begin{aligned} A= \begin{bmatrix} \frac{\bar{k}\sqrt{5}-25\bar{k}+120}{120}&{}\frac{13\sqrt{5}\bar{k}-25\bar{k}}{120}&{}\frac{-\sqrt{5}\bar{k}+\bar{k}}{120}\\ \frac{-13\sqrt{5}\bar{k}-25\bar{k}}{120}&{}\frac{-\bar{k}\sqrt{5}-25\bar{k}+120}{120}&{}\frac{\sqrt{5}\bar{k}+\bar{k}}{120}\\ \frac{-5\bar{k}}{12}&{}\frac{-5\bar{k}}{12}&{}\frac{-\bar{k}+12}{12}\\ \end{bmatrix} , \end{aligned}$$

$$\begin{aligned} B=\begin{bmatrix} 0&{}0&{}\frac{11\sqrt{5} \bar{k}+120\bar{k}}{120}\\ 0&{}0&{}\frac{-\sqrt{5} \bar{k}+120+11\bar{k}}{120}\\ 0&{}0&{}\frac{\bar{k}+12}{12}\\ \end{bmatrix} . \end{aligned}$$

The above system can be written as

$$\begin{aligned} \begin{bmatrix} y_{i+r}\\ y_{i+s}\\ y_{i+1}\\ \end{bmatrix} =M(\bar{k}) \begin{bmatrix} y_{i+r-1}\\ y_{i+s-1}\\ y_{i}\\ \end{bmatrix} , \end{aligned}$$

where $M(\bar{k})=A^{-1}B$ with $\bar{k}=\lambda k$. The eigenvalues of the stability matrix $M(\bar{k})$ are found to be $\{0,0,P(\bar{k})\}$, where

$$\begin{aligned} P(\bar{k})=\frac{-\bar{k}^3-12\bar{k}^2-60\bar{k}-120}{\bar{k}^3-12\bar{k}^2+60\bar{k}-120}. \end{aligned}$$

Now, consider a set S which is the stability region of the method (7) and is defined as

$$\begin{aligned} S=\{\bar{k}\in \mathbf {C} : \vert P(\bar{k})\vert <1\}. \end{aligned}$$

If the portion of the complex plane with real part $Re(z)<0$ is contained in S then the numerical method is said to be A-stable. The graph showing the stability region of the proposed method (7) is shown in Fig. 1. This figure shows that the left half plane is contained in S. Hence, the proposed method in (7) is A-stable.

Fig. 1

Stability region is shaded left Half Plane

Now, we address the stability of the differential system in (15). To investigate the stability of the proposed B-spline approximation, we linearize the non-linear terms $uu_x$, $4u^3u_{xx}+12u^2{u_x}^2+3u^2u_x$ and $u(u-\mu )(1-u)$ in the differential equations (12), (13) and (14), respectively, by assuming that u is locally constant, and thus we take $u(x,t)=U_{ij}$ for (x, t) in a ball of the node $(x_i,t_j)$. Thus, the application of the modified cubic B-splines approach in (9) produces a system of ordinary differential equations of the form:

$$\begin{aligned} \dot{a}=La+M, \end{aligned}$$

where $a=[\alpha _1,\alpha _2,...,\alpha _{N-1}]^T$ and M is a $(N-1)$-vector containing the boundary conditions and non homogeneous parts. The matrix L is a square matrix of order $N-1$, which depending on the problem considered is given, respectively, by

• $L=(\nu L_2-U_{ij}L_1)$ for the Burgers equation,

• $L=(4U_{ij}^3L_2+3U_{ij}^2L_1)$ for the Buckmaster equation,

• $L=(L_2-\mu I)$ for the FitzHugh–Nagumo equation,

where K, $L_1$ and $L_2$ are matrices of order $N-1$ given by

$$\begin{aligned} K= \begin{bmatrix} 4&{}1&{}0&{} &{} &{} \\ 1&{}4&{}1&{} &{} &{}\\ &{} &{}..&{}..&{}..&{}\\ &{} &{}..&{}..&{}..&{}..\\ &{} &{} &{} 0&{}1&{}4\\ \end{bmatrix}, L_1=\frac{3}{h}K^{-1} \begin{bmatrix} 0&{}1&{}0&{} &{} &{}\\ -1&{}0&{}1&{} &{} &{}\\ &{}..&{}..&{}..&{} &{}\\ &{} &{}..&{}..&{}..&{}\\ &{} &{} &{}0&{}-1&{}0\\ \end{bmatrix}, \end{aligned}$$

$$\begin{aligned} L_2=\frac{6}{h^2}K^{-1} \begin{bmatrix} -2&{}1&{}0&{} &{} &{}\\ 1&{}-2&{}1&{} &{} &{}\\ &{} &{}..&{}..&{}..&{}\\ &{} &{} &{}..&{}..&{}..\\ &{} &{} &{} 0&{}1&{}-2\\ \end{bmatrix}. \end{aligned}$$

The stability of the above systems of first order differential equations of the linearized problems will imply the stability of the non-linear systems discretizing the partial differential equations (12), (13) and (14), which in turn is related to the eigenvalues of matrix L. We have that the real part of the eigenvalues of the coefficient matrix L, which further depend upon eigenvalues of $L_1$ and $L_2$, are either zero or negative as verified for different values of h analogous to (Tamsir et al. 2016) as it is shown in Figs. 2 and 3, respectively, for different values of h. For the Burgers equation and the Buckmaster equation, the system is stable. The system concerned with the FitzHugh–Nagumo equation has a bifurcation parameter $\mu$ deciding the character of the equilibrium point. Therefore, for $\mu >0$, the system is found to be stable considering 0 as an equilibrium point. Therefore, the proposed OHBCBM is stable and there is no restriction on the step-size due to stability issues.

Fig. 2

Eigen values of $L_1$

Fig. 3

Eigen values of $L_2$

5 Numerical experiments

This section comprises some numerical results obtained by applying OHBCBM to solve two well-known examples of the Burgers’ equation, two examples of the Buckmaster equation and one example of the FitzHugh–Nagumo equation. The results are compared to the exact solutions and results of some other methods appeared in the literature. To measure the accuracy of the discussed algorithm, the following formulas are used to calculate root mean square error, $L_{rms}$ and maximum absolute error, $L_{\infty }$:

$$\begin{aligned} L_{rms}=\left( \sum _{i=1}^{N}\frac{e_i^2}{N}\right) ^{1/2}; L_{\infty }=\max _{1\le i\le N}\mid {e_i}\mid , \end{aligned}$$

where $e_i=(u(x_i,t)-U(x_i,t))$, being t a specified value in the grid time interval, and $u(x,t), U(x_i,t)$ the exact and the numerical solutions, respectively. The MatlabR2017a software has been used for computations.

5.1 Burgers equation

Example 1

In this example, we consider the initial and boundary conditions for the Burgers equation in (12) given by

$$\begin{aligned} u(x,0)=\sin (\pi x ), 0\le x\le 1, \end{aligned}$$

$$\begin{aligned} u(0,t)=u(1,t)=0, t\ge 0. \end{aligned}$$

Using the Hopf-Cole transformation, the analytical solution of the problem is given by

$$\begin{aligned} u(x,t)=2\pi \nu \frac{\sum _{n=1}^{\infty }c_n \exp (-n^2\pi ^2\nu t)n \sin (n\pi x)}{\sum _{n=1}^{\infty }c_n\exp (-n^2\pi ^2\nu t)\cos (n\pi x)}, \end{aligned}$$

where

$$\begin{aligned} c_0=\int _0^1 \exp \{-(2\pi \nu )^{-1}[1-\cos (\pi x)]\}dx,\\ c_n=2\int _0^1 \exp \{-(2\pi \nu )^{-1}[1-\cos (\pi x)]\}\cos (n\pi x)dx,\\ (n=1,2,3...). \end{aligned}$$

The approximate solutions provided by our method are compared to the ones in Seydaoğlu (2018); Jiwari (2015); Kadalbajoo and Awasthi (2006); Kutluay et al. (2004) for $\nu =0.01, 0.003, 0.004,$ and $N=80,200$, as shown in Tables (2) and (3), respectively, at different grid points. We have taken $k=0.1$ with the OHBCBM method, while the methods used for comparisons used $k=0.01, 0.001$ as given in the tables. The last rows show the maximum absolute errors of each column, been the proposed method the most efficient. Figure 4 shows the physical behavior of the numerical solution for $\nu =0.001$ at different values of t. The results of OHBCBM for $\nu =0.0001$ are compared with the exact solution and the results in Seydaoğlu (2018), as shown in Table 4, and found to be more accurate without using a filtering technique as was done in Seydaoğlu (2018). These tables show that the numerical results are in good agreement with exact solution even though it is performed with a larger value of k, whereas the methods used for comparisons have chosen much smaller time step-sizes.

Table 2 Comparison of numerical solutions and exact solution for Example 1 with values of $\nu =0.01$, $N=80$

Fig. 4

Numerical solution of Example 1 for $\nu =0.001$

Fig. 5

Numerical solution of Example 2 for $\nu =0.001$.

Table 3 Comparison of numerical and exact solutions for Example 1 with values of $\nu =0.003$ and 0.004, $N=200$

Table 4 Comparison of approximate and exact solutions for Example 1 for $\nu =0.0001$, $N=19$, at time $t=1$

Example 2

Consider the initial and boundary conditions for the Burgers equation in (12) as follows

$$\begin{aligned} u(x,0)=4x(1 - x) , 0\le x\le 1, \end{aligned}$$

$$\begin{aligned} u(0,t)=u(1,t)=0, t>0. \end{aligned}$$

The exact solution of this example is

$$\begin{aligned} u(x,t)=2\pi \nu \frac{\sum _{n=1}^{\infty }c_n\exp (-n^2\pi ^2\nu t)n\sin (n\pi x)}{\sum _{n=1}^{\infty }c_n\exp (-n^2\pi ^2\nu t)\cos (n\pi x)} , \end{aligned}$$

where

$$\begin{aligned} c_0=\int _0^1 \exp \{\frac{-1}{3\nu }(3x^2-2x^3)\}dx,\\ c_n=2\int _0^1 \exp \{\frac{-1}{3\nu }(3x^2-2x^3)\}\cos (n\pi x)dx,\\ (n=1,2,3...). \end{aligned}$$

For this example, the performance of OHBCBM has been compared for $\nu =0.01$, $N=80$, in Table 5 and for $\nu = 0.003, 0.004$, $N=200$, in Table 6 with (Seydaoğlu 2018; Jiwari 2015; Kadalbajoo and Awasthi 2006; Kutluay et al. 2004) and found to be the most efficient. Figure 5 shows the physical behavior of the numerical solution for $\nu =0.001$ which indicates the efficient performance of OHBCBM for a small value of the viscosity coefficient.

Table 5 Comparison of approximate and exact solutions for Example 2 for $\nu =0.01$, $N=80$.

Table 6 Comparison of numerical solutions and exact solution for Example 2 for $\nu =0.003,0.004$ and $N=200$.

5.2 Buckmaster equation

Example 3

In this example we consider the Buckmaster equation in (13) for $x\in [0,1]$ with

$$\begin{aligned} f(x,t)=-12x^2e^{4t}-3x^2e^{3t}+xe^{t}, \end{aligned}$$

boundary conditions

$$\begin{aligned} u(0,t)=0,u(1,t)=e^{t},t>0 , \end{aligned}$$

and initial condition $u(x,0)=x,x\in (0,1).$

The exact solution of this problem is $u(x,t)=xe^{t}$. Numerical results with OHBCBM calculated using h=0.2 and $k=0.05$ at time $t=0.05$ and comparisons with the results in Chanthrasuwan et al. (2017) for $k=0.01$ using the Crank–Nicolson scheme and a fully implicit scheme are given in Table 7. In Hussain and Alwan (2013), numerical results have been calculated using $h=0.01$ and $k=0.01$ at $t=0.03$, resulting

$$\begin{aligned} \max _{1\le i\le N}\mid u(x_i,t)-U(x_i,t)\mid =0.022533, \end{aligned}$$

while for the OHBCBM method, using $h=0.01$ and $k=0.03$, at $t=0.03$ we found that

$$\max _{1\le i\le N}\mid u(x_i,t)-U(x_i,t)\mid =4.8210\times {10}^{-11},$$

showing a great performance.

Table 7 Comparison of exact and approximate solutions for Example 3 with $h=0.2$ at $t=0.05$

Example 4

We consider the Buckmaster equation in (13) for $x\in [0,1]$ with

$$\begin{aligned} f(x,t)=-12x^2\cos ^4(t)-3x^2\cos ^3(t)+x\sin (t), \end{aligned}$$

with boundary conditions

$$\begin{aligned} u(0,t)=0,u(1,t)=\cos (t),t>0, \end{aligned}$$

and initial condition

$$\begin{aligned} u(x,0)=x,x\in (0,1). \end{aligned}$$

The exact solution of this problem is $u(x,t)=x\cos (t)$. The approximate solution is calculated using $h=0.2$ and $k=0.05$ at time $t=0.05$ and comparison with the results in Chanthrasuwan et al. (2017) using the Crank-Nicolson scheme and a fully implicit scheme are given in Table (8), which indicate that the results with the OHBCBM method fully agree with the exact solutions.

Table 8 Comparison of approximate and exact solution for Example 4 with $h=0.2$ at $t=0.05$

5.3 FitzHugh–Nagumo equation

Example 5

Consider the FHN equation in (14) for $x\in [-10,10]$ with boundary conditions

$$\begin{aligned} u(-10,t)=\frac{1}{2}+\frac{1}{2}\tanh \Big (\frac{1}{2\sqrt{2}}\Big (-10-\frac{2\mu -1}{\sqrt{2}}t\Big )\Big ), \end{aligned}$$

$$\begin{aligned} u(10,t)=\frac{1}{2}+\frac{1}{2}\tanh \Big (\frac{1}{2\sqrt{2}}\Big (10-\frac{2\mu -1}{\sqrt{2}}t\Big )\Big ), \end{aligned}$$

and initial condition

$$\begin{aligned} u(x,t)=\frac{1}{2}+\frac{1}{2}\tanh \Big (\frac{x}{2\sqrt{2}}\Big ),x\in [-10,10]. \end{aligned}$$

The exact solution is

$$\begin{aligned} u(x,t)=\frac{1}{2}+\frac{1}{2}\tanh \Big (\frac{1}{2\sqrt{2}}\Big (x-\frac{2\mu -1}{\sqrt{2}}t\Big )\Big ). \end{aligned}$$

To measure the accuracy, $L_{rms}$ and $L_{\infty }$ have been calculated for $\mu =0.75$ with $h=0.2$ and $k=0.1$ and the results compared with those in Jiwari et al. (2014); Ahmad et al. (2019) and presented in Table (9). These results indicate that OHBCBM performs better than the method in Jiwari et al. (2014) and competes favorably with the one in Ahmad et al. (2019) with major differences for the values of k considered. Also, a 3-D plot of the exact and the numerical solutions for $\mu =0.75$, $t\in [0,5]$ are shown in Fig. 6.

Table 9 Comparison of numerical results by OHBCBM with results from Jiwari et al. (2014); Ahmad et al. (2019) for Example 5 with $\mu =0.75$, $N=100$

Fig. 6

Exact solution on left and numerical solution of Example 5 for $\mu =0.75$ with $h=0.2$ and $k=0.1$ upto time $t=5$

6 Conclusions

In this article, a combined scheme OHBCBM using an optimized hybrid block method and a modified cubic B-spline based collocation approach has been proposed. It has been implemented on three non-linear problems, Burgers equation, Buckmaster equation and FitzHugh–Nagumo equation, to prove the efficiency of the method. The suggested scheme has the following advantages:

• 1. Accuracy is achieved with larger time step-sizes than those used by of other methods, as can be seen in the numerical experiments. That means that k (or $\Delta t$) needs not to be very small for the suggested method.

• 2. The proposed numerical method does not require quasilinearization or any filtering technique to deal with non-linear terms.

• 3. It shows a correct physical behavior for small values of the viscosity coefficient $\nu =0.001$ as shown in the two examples discussed of the Burgers equation.

• 4. The proposed method is stable, simple and accurate.

Numerical results are compared with the results of existing methods and are in good agreement with exact solutions, which shows that the method is efficient and can be considered accurate for a variety of non-linear equations.