A Collocation Method for Numerical Solution of Hyperbolic Telegraph Equation with Neumann Boundary Conditions

We present a technique based on collocation of cubic B-spline basis functions to solve second order one-dimensional hyperbolic telegraph equation with Neumann boundary conditions. The use of cubic B-spline basis functions for spatial variable and its derivatives reduces the problem into system of first order ordinary differential equations. The resulting system subsequently has been solved by SSP-RK54 scheme.The accuracy of the proposed approach has been confirmed with numerical experiments, which shows that the results obtained are acceptable and in good agreement with the exact solution.


Introduction
The hyperbolic partial differential equations model the vibrations of structures (e.g., building, beams, and machines) and are the basis for fundamental equations of atomic physics.In this paper we consider nonlinear second order onedimensional hyperbolic telegraph equation   (, ) + 2  (, ) +  2  (, ) =   (, ) +  (, ) +  () ,  ≤  ≤ ,  ≥ 0 (1) with initial conditions  (, 0) =  1 () ,   (, 0) =  2 () and Neumann boundary conditions   (, ) =  1 () ,   (, ) =  2 () ,  ≥ 0, where  and  are known real constants.If () = 0, (1) represents a linear hyperbolic telegraph equation in which both the electric voltage and current, in a double conductor, satisfy the equation, where  is distance and  is time.For  > 0,  = 0 it represents a damped wave equation and for  >  > 0 it is called a telegraph equation.Linear second order hyperbolic telegraph equation arise in the study of propagation of electric signal in a cable of transmission line and wave phenomena.Interaction between convection and diffusion or reciprocal action of reaction and diffusion describes a number of nonlinear phenomena in physical and biological process.In fact telegraph equation is more suitable than ordinary diffusion equation in modeling reaction diffusion for such branches of sciences.Numerical solutions of one-dimensional linear hyperbolic equation with Dirichlet boundary conditions have been studied by many authors.El-Azab and El-Ghamel [1] have used Routhe-wavelet method for the numerical solution of telegraph equation.Dehghan and Shokri [2] presented a meshless method based on collocation with radial basis functions.Spline solutions of hyperbolic telegraph equation have been studied in [3][4][5][6].L. B. Liu and H. W. Liu [3] used the quartic spline method.In [4], H.-W. Liu and L.-B.Liu applied an unconditionally stable spline difference method; Dosti and Nazemi [5] presented a quartic -spline collocation method; modified cubic -spline collocation method has been proposed by Mittal and Bhatia [6] to find the numerical solution of one-dimensional linear hyperbolic telegraph equation.Numerical solution of linear and nonlinear onedimensional hyperbolic telegraph equation with variable coefficient has been presented in [7,8].Unconditionally stable finite difference schemes have been proposed in [9,10].Dehghan and Lakestani [11] used the Chebyshev cardinal functions, whereas Saadatmandi and Dehghan [12] used the Chebyshev tau method for expanding the approximate solution of one-dimensional telegraph equation.Mohebbi and Dehghan [13] reported a higher order compact finite difference approximation of fourth order in space and used collocation method for time direction.Other techniques used for numerical solutions of one-dimensional hyperbolic telegraph equation with Dirichlet boundary conditions are interpolating scaling function technique [14] and radial basis function technique [15].Thus much work has been done to solve (1) with Dirichlet boundary conditions.Little has been done for numerical solution of one-dimensional hyperbolic telegraph equation with Neumann boundary conditions.In [16], Dehghan and Ghesmati reported a dual reciprocity boundary integral equation (DRBIE) method, in which three different types of radial basis functions have been used to approximate the solution of one-dimensional linear hyperbolic telegraph equation.Recently, L. B. Liu and H. W. Liu [17] developed an unconditionally stable compact difference scheme for one-dimensional telegraphic equations with Neumann boundary conditions.
In the literature the -spline collocation method has been successfully applied to find the numerical solutions of various linear and nonlinear partial differential equations [6,18,19].-splines have some special properties like local support, smoothness, and capability of handling local phenomena, which make them suitable to solve linear and nonlinear partial differential equations easily and elegantly.The use of cubic -spline basis functions leads to a global solution for which not only the functions but also their first and second derivatives are continuous.The present collocation method has two great advantages: the set-up procedure does not involve integrations and the resulting matrix system is banded with small band width.Hence, -spline when combined with the collocation provides a simple solution procedure of linear and nonlinear partial differential equations.In this paper, we have developed a collocation method based on cubic -spline basis function to solve hyperbolic telegraph equation with Neumann boundary conditions.Equation ( 1) is converted into a system of equations using the transformation   = V.Then for approximating the solution we use the collocation of cubic -spline basis functions.Finally we get a system of first order systems of ordinary differential equations, which have been solved by five-stage, fourth-order strong stability preserving time stepping Runge Kutta (SSP-RK54) [20] scheme.The method needs less storage space and causes less accumulation of numerical errors.
The outline of the paper is as follows.In Section 2, cubic -spline basis functions are introduced.In Section 3, numerical scheme is presented for telegraph equation using cubic -spline functions based collocation method.In Section 4, initial vector is computed for given initial conditions.In Section 5, stability of the scheme is discussed and found to be unconditionally stable.In Section 6 computational results for telegraph equations (1)-(3) for some test problems are illustrated and finally conclusions are included in Section 6.

Cubic B-Spline Collocation Method
The solution domain  ≤  ≤  is partitioned into a mesh of uniform length ℎ =  +1 −   by the knot   , where  = 0, 1, 2, . . .,  − 1,  such that In the cubic -spline collocation method, the approximate solution can be written as the linear combination of cubic -spline basis functions for the approximation space under consideration.
Our numerical treatment for solving (1)-( 3) using the collocation method with cubic -spline is to find an approximate solution (, ) to the exact solution (, ) in the form where   () are the time dependent quantities to be determined from boundary conditions and collocation from the differential equation.
The cubic -spline   () at the knots is given by where the set of functions { −1 ,  0 ,  1 , . . .,  −1 ,   ,  +1 } forms a basis for the function defined over the region  ≤  ≤  with the obvious adjustment of the boundary base functions to avoid undefined knots.Each cubic -spline covers four elements so that an element is covered by four cubic splines.
The values of   () and its derivatives are tabulated in Table 1.Since   () = 0 outside the interval ( −2 ,  +2 ), so there is no need to tabulate   () for other values of .
Then, using approximate function ( 4) and Table 1, the approximate values of (, ) and its two derivatives at the knots are determined in terms of the time parameters   as follows:

Numerical Scheme
The telegraph equation ( 1) is decomposed into a system of partial differential equations using the following transformation.
Then the transformed form of ( 1) is For solving the coupled equations ( 7) we assume our approximate solution as the linear combination of cubic -spline basis function Using approximate solution (8), approximate values of   (, ) can be written as follows: where ċ () is the derivative of   () with respect to time .
Once the vector  = [ 0 ,  1 , . . .,   ]  is computed at a specific time level, the approximate solution (, ) can be computed at the required knots.
System (16) represents a system of (2 + 2) first order ordinary differential equations, which is solved using SSP-RK54 [20] scheme, with a variant of Thomas algorithm and consequently the approximate solution   (, ) is completely known.

Initial Vectors
To solve the resulting systems of first order ordinary differential equations, we need initial vectors  0 and V 0 , which are computed from the initial and boundary conditions as follows.
4.1.Initial Vector  0 .Initial vector  0 is computed from the initial condition and boundary values of the derivatives as the following expressions: Using ( 10)- (12), in (18) we get a ( + 1) × ( + 1) system of equations of the form where The matrix  is tridiagonal matrix.System ( 19) is solved using Thomas algorithm and hence initial vector  0 is computed.

Initial
Vector V 0 .Initial vector V 0 can be computed using the initial condition (2) as We have

Stability of Scheme
The stability of system ( 16) is very important since it is related to the stability of numerical scheme for solving it.If the system of ordinary differential equations ( 16) is unstable, then stable numerical scheme for temporal discretization may not generate converged solution.The stability of this system depends on the eigenvalues of coefficient matrix, since its exact solution can be directly found using its eigenvalues.For linear case, that is, if () = 0, we consider where is the unknown vector and () is a column vector of order 2(−1).Consider where  1 and  2 are symmetric tridiagonal matrix of order  − 1, given by International Journal of Computational Mathematics 5 where  and  are identity and null matrix of order  − 1, respectively.
Stability of system (23) depends on the eigenvalues of the coefficient matrix  −1 .If all the eigenvalues are having negative real part, then the system (23) will be stable.
Let  be any eigenvalue of  −1  and  1 and  2 are two components each of order ( − 1), of eigenvector corresponding to eigenvalue .We have From (26), we have From the above system of equations, we get ⇒ ( + 2) is an eigenvalue of  2  −1 1 .Since  1 and  2 are tridiagonal matrices, so their eigenvalues are 1 and  2 are symmetric matrices, so eigenvalues of matrix  2  −1 1 are which is negative and real.Let  =  + , where  and  are real numbers.
We have that ( + )( +  + 2) is real and negative, which implies From the above set of equations, we get the solutions as follows: (i)  +  = 0 and  is arbitrary real number; ⇒  is negative real number, since  is real and positive; (ii)  = 0 ⇒ ( + 2) < 0; ⇒  is negative and real, since  is real and positive.
Hence, the real part of eigenvalues of  −1  should always be negative for stability.When () ̸ = 0, that is, when (1) is nonlinear, we have the following system: The stability can be discussed by finding the eigenvalues of the matrix  −1 , where  = / is the Jacobian matrix.For stability eigenvalues of matrix  −1  should be negative.

Numerical Experiments
In this section, three experiments including linear and nonlinear form of problem (1)-( 3) are considered to demonstrate the accuracy and applicability of the proposed method.The numerical efficiency is shown by calculating the  ∞ ,  2 , and root mean square error norms.
Example 1.We consider the following linear telegraph equation: with initial conditions and boundary conditions The exact solution is given [16] by 2 ,  ∞ , and RMS errors are reported in Table 2 with ℎ = .02,ℎ = .05,and Δ = .001.RMS errors with different space step sizes are reported in Table 3 and compared with the results given by Dehghan and Ghesmati [16].It is noticed that our results are in good agreement with the results of [16].
From Figure 1, it is clear that approximate solution coincides  with the exact solution at  = 1, 2, 3, with ℎ = .05and Δ = .001. Figure 2 depicts the space-time graph of approximate solution up to time  = 3 with ℎ = .05and Δ = .001.Similar figures have been depicted in [16,17].It is clear that the obtained numerical results are accurate and comparable favorably with the exact solutions and earlier studies.
The exact solution is given [16] by 2 ,  ∞ , and RMS errors are reported in Table 4 with ℎ = .05and Δ = .001.In Table 5, RMS errors are reported with different space step sizes and compared with the results given by Dehghan and Ghesmati [16].Figure 3 shows that the approximate solution and exact solution coincide for  = 1, 2, 3 with ℎ = .05and Δ = .001.Space-time graph of approximate solution is depicted in Figure 4 with ℎ = .05and Δ = .001up to  = 3.Similar figures have been depicted in [16,17].with initial conditions and boundary conditions The exact solution is given [17] by 2 ,  ∞ , and RMS errors are computed for different values of  and reported in Table 6 with ℎ = .05and Δ = .001.In Table 7, error norms are reported for ℎ = .1 and Δ = .01and compared with the results given by L. B. Liu and H. W. Liu [17].From Figure 5, we observe that approximate solution and exact solution coincide for  = 1, 2, 3 with ℎ = .05and Δ = .001.Space-time graph of approximate solution is depicted in Figure 6 with ℎ = .05and Δ = .001,which shows the approximate solution profile up to  = 3.Similar figure has been depicted in [17].
It should be noticed that in this case we obtain a system of nonlinear first order ordinary differential equations.Unlike the numerical methods [17], which use finite difference methods to solve nonlinear hyperbolic wave equation to get a system of nonlinear algebraic equations and finally solved using the Newton-Raphson method, we do not need to do any extra effort to handle the nonlinear term and solutions are easily computed by solving the obtained system of ODEs using SSP-RK54 scheme.Hence, the present scheme reduces the computational cost.

Conclusions
The main idea of the present paper is to first convert the given problem into a coupled system of partial differential equations and then using cubic -spline basis functions for spatial variable and derivatives, the coupled system of equations is converted into system of first order ordinary differential equations.The resulting system of first order ordinary differential equations has been solved by SSP-RK54 scheme.The stability of the scheme is discussed using matrix stability analysis and found to be unconditionally stable.The method is also capable of solving nonlinear telegraph equation with Neumann boundary conditions.It has been noticed that obtained numerical solutions are in good agreement with the exact solution and earlier studies.These facts illustrate that the proposed method is a reliable, valid, and powerful tool for solving telegraph equation with Neumann boundary conditions.
Numerical solution at t = 1 Exact solution at t = 1 Numerical solution at t = 2 Exact solution at t = 2 Numerical solution at t = 3 Exact solution at t = 3

Table 1 :
Values of cubic -spline and its derivatives at different knots.

Table 2 :
The errors of numerical solutions with Δ = .001.

Table 3 :
Comparison of RMS errors at  = 3.

Table 4 :
The errors of numerical solutions with Δ = .001.

Table 5 :
Comparison of RMS errors at  = 2 for Example 2.

Table 6 :
The errors of numerical solutions with Δ = .001.

Table 7 :
Comparison of RMS errors at different time levels for Example 3.