A novel B-spline collocation method for Hyperbolic Telegraph equation

1.   Introduction

Hyperbolic partial differential equations (HPDEs) have drawn much attention in recent years. This type of equations play a key role in understanding physical phenomena such as vibrations of structures and several atomic physics fields. The hyperbolic TE, used in the modelling of signal analysis for transmission and propagation of electric signals, is one of the basic type of HPDEs. In this paper, we deal with the approximate solution of the following TE

which satisfies the following boundary conditions (BCs)

and initial conditions (ICs)

where $\alpha$ and $\lambda$ are positive constants and $g$ is adequately differentiable forcing function of $x$ and $t$. Suppose further that $\Gamma _{0}, \Gamma _{1}, \Gamma _{2}, \Gamma _{3}$ and their derivatives are continuous functions of $t$ and similarly $\Psi _{0}, \Psi _{1}$ and their derivatives are continuous function of $x$. A variety of the numerical methods have been proposed in literature to establish the solutions of the TE, such as the Rothe-wavelet approximation ^[1], legendre multiwavelet Galerkin technique ^[2], the Rothe-wavelet Galerkin method ^[3], the reproducing Kernel Hilbert space method ^[4], the Chebyshev wavelets approach ^[5], the meshfree procedure with the radial basis functions ^[6], the computational method based on the polynomial scaling functions ^[7], the numerical approach associated Hermite orthogonal functions ^[8], the Bessel collocation method ^[9], the Galerkin approach ^[10], the implicit three-level difference scheme ^[11], the differential quadrature method ^[12,13,14]. In addition, in the past few years, various B-spline collocation techniques have been implemented to obtain the approximate solution of the TE. Alshomrani et al. ^[15] designed a new algorithm based on modified cubic trigonometric B-spline functions. Dehghand and Shokri ^[16] proposed a numerical scheme structured on thin plate splines. The collocation approach based on quartic, septic and cubic B-splines were presented to compute the approximate solution of the TE in the studies ^[17,18] and ^[19,20,21] respectively. Nazir et al. ^[22] applied the cubic trigonometric B-spline approach for the approximate solution of the TE. Sharifi and Rashidinia ^[23] proposed a collocation approach with extended cubic B-spline functions to solve numerically the TE. Singh et al. ^[24] presented a method obtained by using exponential B-spline collocation procedure in space and second-order Runge-Kutta scheme in time. Later, Singh et al. ^[25] derived a fourth-order cubic spline technique to solve numerically the TE. A high-order new numerical scheme is introduced in the study ^[26] in which the spatial integration of Eq 1.1 is managed through collocation technique and fourth-order implicit scheme is used for temporal variables.

In the present study, our aim is to produce numerical results with high order accurate by enhancing the accuracy in both time and space. For this purpose, a novel optimal B-spline collocation method based on QBS basis functions is developed to discretize the spatial domain and fourth-order implicit scheme is derived for time integration. The advantage of the suggested technique over the existing ones in the literature is that the suggested technique has an accuracy of $O\left(h^{6}\right)$ in the spatial direction and an accuracy of $O\left(\Delta t^{4}\right)$ in temporal direction and produces excellent results even with less points of the temporal and spatial domains. The structure of this paper is as follows: In section 2, the time discretization is carried out and a novel optimal B-spline collocation method is constructed. In section 3, the application of the proposed method to (1.1) is explained. In section 4, the numerical experiments are provided and the results are reported in the table form. The conclusion of the study is given in the section 5 with the remarks about the main observations.

2.   Derivation of the proposed method

In this section, we derive a new high-order approach to solve the problem (1.1). Firstly, by introducing an auxiliary variable $v = u_{t}$, the original equation (1.1) is transformed to the following system of first-order (in time) equations

The relevant BCs and ICs are rewritten as

and

Consider the uniform partition of the spatial and temporal domains represented by $\Omega = \left[ a, b\right] \; \times \; \left[ 0, T\right]$ with the grid points $(x_{j}, t_{n})$, where $x_{j} = a+jh, \; j = 0, 1, 2, ..., N, \; t_{n} = n\Delta t, \; n = 0, 1, 2, ...$.

2.1.   Temporal discretization

Utilizing finite difference method which is fourth-order implicit scheme, the time discretization of Eqs (2.1) and (2.2) is derived as follows

and

where

obtained by taking partial derivative of both sides Eq (2.2) with respect to $t$ and using Eq (2.1) and $\theta _{1}$, $\theta _{2}$, $\theta _{3}$, $\theta _{4}$ are unknown parameters to be defined later.

Using Eqs (2.1) and (2.8) in Eq (2.6) yields

Equation (2.10) can be rearranged to the form given as Eq (2.11)

where

In a similar manner, using Eqs (2.2) and (2.9) in Eq (2.7) gives

After simplifying the above equation, we get,

where

Lemma 1. Suppose that $u, v, f\in C^{6}\left(\Omega \right)$ Then, when

the suggested scheme is consistent and has fourth-order accurarcy in temporal direction for the norm $\left \Vert \cdot \right \Vert _{\infty }$.

Proof. For the proof, see ^[10].

It can be also easily observed that by selecting the following parameters in Eqs (2.6) and (2.7)

we achieve Crank-Nicolson scheme having an accurarcy of $O(\Delta t^{2})$ in time.

2.2.   Space discretization

In this subsection, we describe a fully discrete scheme by means of a novel optimal B-spline collocation (NOBSC) method based on QBS basis functions in the spatial discretization. To this end, we first provide some properties of quintic spline interpolant (QSI) which will be used later in the formulation of NOBSCM.

2.2.1.   Properties of QBS interpolant

Let $\pi = \left \{ a = x_{0} < x_{1} < ... < x_{N} = b\right \}$ be a uniform partition over the interval $\left[ a, b\right]$, where $x_{j}-x_{j-1} = h$ is the step size. In order to construct the QBS functions, ten additional mesh points are required outside the interval $\left[ a, b\right]$, which are positioned as $x_{-5} < x_{-4} < x_{-3} < x_{-2} < x_{-1}$ and $x_{N+1} < x_{N+2} < x_{N+3} < x_{N+4} < x_{N+5}$. By making use of the results in ^[27,28], we can define QBS basis functions $\phi _{i}(x)$ for $i = -2, -1, ..., N+1, N+2$ as follows

where $f_{j} = x-x_{j}$. The set of QBS $\left \{ \phi _{-2}, \phi _{-1}, \ldots, \phi _{N+1}, \phi _{N+2}\right \}$ generates a basis over the solution interval$\ \left[ a, b\right]$. Let $S(x, t)$ and $R(x, t)$ be QBS approximations of the exact solutions $u\left(x, t\right)$ and $v\left(x, t\right)$, respectively Then, $S(x, t)$ and $R(x, t)$ are expressed by the sum of QBS basis functions as:

where $c_{j}$ and $d_{j}$, $j = -2, -1, 0, ..., N+2$, are unknown time dependent quantities to be determined via collocation technique with the convenient BCs. Throughout this study, we denote

Let $S(x, t)$ and $R(x, t)$ be the QSI of $u(x, t)$ and $v(x, t)$, satisfying the following interpolation conditions

Theorem 1. Let $S(x, t)$ and $R(x, t)$ be the QSI for $u(x, t)$ and $v(x, t)\in C^{8}\left(\Omega \right)$, respectively, and satisfy the interpolation conditions (2.15–2.18). Then, the below given results hold at the knot points $x_{j}$, $j = 0, 1, ..., N$ :

and

where $K^{\left(l\right) }, \; u^{\left(l\right) }, \; P^{\left(l\right) }$ and $v^{\left(l\right) }$ denote the $l$ th derivative w.r.t. '$x$'.

Proof. For the proof see ^[29].

2.2.2.   Construction of NOBSC method

Here, we construct a NOBSC method by developing a new approximation for $S^{\prime \prime }(x, t)$ and $R^{\prime \prime }(x, t)$. To do this, we first define the discrete operator $\xi$ as follows

for any function $M$ defined at points of spatial discretization.

Lemma 2. Let $S(x, t)$ and $R(x, t)$ be the QSI for $u(x, t)$ and $v(x, t) \; \in C^{8}\left(\Omega \right)$ and satisfy the interpolation conditions (2.15–2.18). Then we get,

Proof. We prove relation (2.24). From Eq (2.19), we have the following equality

The application of the operator $\xi$ defined by (2.23) on (2.26) gives

By using Taylor expansion and finite differences for the $u$ terms at $x_{j}$ on the right hand side of (2.27), we obtain,

This proves relation (2.24). In a similar manner, one can prove the relation (2.25).

Corollary 1. If $u(x, t)$ and $v(x, t)\in C^{8}\left(\Omega \right)$, then for $\ j = 2, 3, ..., N-2$, the below given approximations hold.

Proof. By Lemma 2 and Theorem 1, one can prove this corollary.

Lemma 3. Suppose that $S(x, t)$ and $R(x, t)$ are the QSI for $u(x, t)$ and $v(x, t) \; \in C^{8}\left(\Omega \right)$ and satisfy the interpolation conditions (2.15–2.18). Then, we have

For $j = 0, 1:$

For $\left(j, k\right) = \left(N-1, 1\right), \left(N, 0\right) :$

For $j = 0, 1:$

For $\left(j, k\right) = \left(N-1, 1\right), \left(N, 0\right) :$

Proof. We first prove relation (2.28) for $j = 1$. Consider the following approximation for $u^{(6)}(x_{1}, t)$,

Making use relation (2.24) for $j = 2, 3$, from Eq (2.32), we obtain

Hence, this proves relation (2.28) for $j = 1$. Now, we prove relation (2.28) for $j = 0$. For this purpose, we take into consideration the following approximation for $u^{(6)}(x_{0}, t)$,

Considering relation (2.28) for $j = 1$ and relation (2.24) for $j = 2$, it follows from Eq (2.33) that

Hence, the relation (2.28) for $j = 0$ is obtained. In a similar way, the remaining relations in Lemma 3 can be proved. Thus, the proof of Lemma 3 is completed.

Corollary 2. If $u(x, t)$ and $v(x, t) \; \in C^{8}\left(\Omega \right)$, then for $j = 0, 1, N-1, N$, the following relations are obtained

Proof. With the help of Lemma 3 and Theorem 1, one can prove corollary 2.

Let $S(x, t)$ and $R(x, t)$ be optimal QBS approximate solutions of $u(x, t)$ and $v(x, t)$, respectively. Then, the values of $S(x, t)$ and $R(x, t)$ and their first two higher-order derivatives are as follows:

and

3.   Implementation of the proposed method

Substituting Eq (2.14) into Eqs (2.11) and (2.12), using Eqs (2.34) and (2.35), we get the following systems:

For $j = 0$,

For $j = 1$,

For $2\leq j\leq N-2$

For $j = N-1$,

For $j = N,$

And for $j = 0$,

For $j = 1$,

For $2\leq j\leq N-2$

For $j = N-1$,

For $j = N,$

where

and for $0\leq j\leq N$

The final system obtained by combining Eqs (3.1)–(3.5) with Eqs (3.6)–(3.10) consists of $2N+2$ linear equations involving $2N+10$ unknown parameters, $\mathbf{c} = \left(c_{-2}^{n+1}, c_{-1}^{n+1}, ..., c_{N+2}^{n+1}\right)$ and $\mathbf{d} = \left(d_{-2}^{n+1}, d_{-1}^{n+1}, ..., d_{N+2}^{n+1}\right)$. To set the number of unknown parameters equal to the number of equations, BCs are implemented to eliminate unknown parameters $c_{-2}^{n+1}$, $d_{-2}^{n+1}$, $c_{-1}^{n+1}$, $d_{-1}^{n+1}$, $c_{N+1}^{n+1}$, $d_{N+1}^{n+1}$, $c_{N+2}^{n+1}$, $d_{N+2}^{n+1}$ from the system. Therefore, the resultant system can be written in the following form

where $A$ is $\left(2N+2\right) \times \left(2N+2\right)$ matrix, $B$ and $X$ are $\left(2N+2\right) \times 1$ column matrices. In order to initiate the time evolution process, the initial vector

is first determined by utilizing ICs and BCs of the problem. After getting the values of the initial parameters $\mathbf{c}^{0}$ and $\mathbf{d}^{0}$, the unknown vector

for $n = 0, 1, 2, ...$ is iteratively calculated at the required any time level. Pseudo code for the numerical algorithm is given below.

 

4.   Numerical results

In this section, three numerical examples are presented to validate the applicability and efficiency of the suggested technique. We compute the $L_{\infty }$ error norm defined by the following formula:

where $u$ and $U_{N}$ denote analytical and approximate solutions, respectively. The convergence orders in time and space are calculated by the following formulas:

where $\left(L_{\infty }\right) _{h_{i}}$and $\left(L_{\infty }\right) _{\Delta t_{i}}$ are the error norm $L_{\infty }$ for space step $h_{i}$ and time step $\Delta t_{i}$, respectively.

4.1.   Test problem 1

Consider the following TE

with BCs and ICs

where

The analytical solution of this problem is also taken as

The computational process is carried out at different levels of time for the step sizes $h = 0.001$, $h = 0.01$, $h = 0.02$, $\Delta t = 0.001$ and $\Delta t = 0.01$. The error $L_{\infty }$ is reported in Table 1 along with the results obtained by the methods in ^{[8,14,20,23,25,26]}. From Table 1, we can observe that our method provides far better results than the methods in ^{[8,14,20,23,25,26]}. Moreover, Table 1 demonstrates the improvement in results even with fewer collocation points and number of divisions of the temporal domain, which makes the proposed method computationally efficient. The spatial and temporal orders of the convergence along with the error norm $L_{\infty }$ are listed in Tables 2 and 3, respectively. As expected from the theoretical analysis, the used technique approaches fourth-order accuracy in time and nearly sixth-order accuracy in space. The propagation of the approximate solution obtained for $\Delta t = 0.01$ and $h = 0.01$ at different times up to $t = 1$ is plotted in Figure 1. Furthermore, the absolute error at $t = 1$ is given in Figure 2.

4.2.   Test problem 2

Consider the TE given by

with $g\left(x, t\right) = 23\exp \left(-2t\right) \sinh \left(x\right)$ over interval $\left[ 0, 1\right]$ for space variable. The ICs and BCs are considered as:

The analytical solution of the above test problem is

The computations are performed at various times by taking $h = 0.05$, $h = 0.01$, $\Delta t = 0.01$ and $\Delta t = 0.001$. The error $L_{\infty }$ is presented in Table 4. It is seen from Table 4 that the obtained error $L_{\infty }$ is satisfactorily better and smaller than those of the methods given in ^[10,12] even with less number of divisions of the temporal domain. We have listed the spatial and temporal order of convergence along with errors in Tables 5 and 6, respectively. From the Tables 5 and 6, it can be said that the proposed scheme has fourth-order temporal accuracy and sixth-order rate of convergence in space. The simulation of the approximate solution for $h = 0.05, \Delta t = 0.01$ at various values of $t$ is shown graphically in Figure 3. The absolute error of the presented method for $h = 0.05$ and $\Delta t = 0.01$ is drawn in Figure 4.

4.3.   Test problem 3

Consider the following TE

with ICs and BCs

where

The analytical solution of this problem is considered as

This problem has been solved by choosing $\Delta t = 0.01$, $\Delta t = 0.001$, $h = 0.05$ and $h = 0.005$. The error $L_{\infty }$ at various values of $t$ is reported in Table 7 which confirms that the present method produces far better results than the methods in ^{[14,17,20,24,25]}. In addition, Table 7 shows the improvement in results even with small partitions in space and time domains. Table 8 displays the error $L_{\infty }$ and order of spatial convergence with different spatial steps $h$ at $t = 1$. In Table 9, the error $L_{\infty }$ and the order of temporal convergence are calculated with different time steps $\Delta t$ at $t = 1$. When examined the results in Tables 8 and 9, it can be seen that our method has almost accuracy of $O\left(h^{6}\right)$ in the spatial direction and an accuracy of $O\left(\Delta t^{4}\right)$ in temporal direction. This shows that the orders of the temporal and spatial convergence obtained by our method are compatible with the theoretical analysis. The simulation of the approximate solution obtained for $h = 0.05$ and $\Delta t = 0.01$ at various values of $t$ is given in Figure 5. Figure 6 displays the absolute error at $t = 1.$

5.   Conclusions

In this study, we have constructed a new high-order numerical scheme with fourth-order accuracy in time and sixth-order accuracy in space. A new approximation for the second-order spatial derivative is developed. The local truncation error analysis is discussed for time discretization. The error estimation of optimal QBS technique is carried out. Three test problems are provided to justify the effectiveness and improvement in the obtained numerical results and to check the theoretical rate of convergence. It is obviously observed that theoretical prediction is consistent with the numerical calculation. Moreover, the obtained approximate results are compared with other some existing methods applied to find the numerical solution of the TE. Comparison confirms that the present scheme provides far better results than the methods given in ^{[8,10,12,14,17,20,23,24,25,26]}. A major improvement of the proposed technique is that satisfactory results are provided even with fewer points in space and time domains. To conclude, the proposed method approximates very well solution of the TE and is computationally efficient for solving the TE. Also, as a future direction, the technique suggested for the spatial discretization can also be applied to two-dimensional problems in space ^[30] and time fractional-type problems ^[31].

Acknowledgments

The author would like to thank the anonymous referees and the editor for very helpful suggestions and comments.

Conflict of interest

The author declares that he has no competing interest.