Differential quadrature method (DQM) is proposed for the numerical solution of one- and two-space dimensional hyperbolic telegraph equation subject to appropriate initial and boundary conditions. Both polynomial-based differential quadrature (PDQ) and Fourier-based differential quadrature (FDQ) are used in space directions while PDQ is made use of in time direction. Numerical solution is obtained by using Gauss-Chebyshev-Lobatto grid points in space intervals and equally spaced and/or GCL grid points for the time interval. DQM in time direction gives the solution directly at a required time level or steady state without the need of iteration. DQM also has the advantage of giving quite good accuracy with considerably small number of discretization points both in space and time direction.

1. Introduction

Consider the second-order linear hyperbolic telegraph equation in one-dimensional space: 𝑢𝑡𝑡(𝑥,𝑡)+2𝛼𝑢𝑡(𝑥,𝑡)+𝛽2𝑢(𝑥,𝑡)=𝑢𝑥𝑥(𝑥,𝑡)+𝑓(𝑥,𝑡),inΩ1,(1.1) and the second-order two-dimensional hyperbolic telegraph equation: 𝑢𝑡𝑡(𝑥,𝑦,𝑡)+2𝛼𝑢𝑡(𝑥,𝑦,𝑡)+𝛽2𝑢(𝑥,𝑦,𝑡)−𝑔(𝑥,𝑦,𝑡)=𝑢𝑥𝑥(𝑥,𝑦,𝑡)+𝑢𝑦𝑦(𝑥,𝑦,𝑡)inΩ2,(1.2) where Ω1=[ğ‘Ž<𝑥<𝑏]×[𝑡>0] and Ω2=[0<𝑥,𝑦<1]×[𝑡>0].

The hyperbolic differential equations are the basis for fundamental equations of atomic physics and commonly used in signal analysis for transmission and propagation of electrical signals. Equations (1.1)-(1.2) are damped wave equations for 𝛼>0,𝛽=0 of which the solution is of great importance in wave phenomena and telegraph equation for 𝛼≥𝛽>0. The initial conditions are assumed to be 𝑢(𝑥,0)=𝑣1(𝑥),𝑢𝑡(𝑥,0)=𝑣2(𝑥),ğ‘Ž<𝑥<𝑏(1.3) with Dirichlet type boundary conditions 𝑢(ğ‘Ž,𝑡)=ℎ1(𝑡),𝑢(𝑏,𝑡)=ℎ2(𝑡),𝑡≥0,(1.4) or Neumann type boundary conditions 𝑢𝑥(ğ‘Ž,𝑡)=𝑘1(𝑡),𝑢𝑥(𝑏,𝑡)=𝑘2(𝑡),𝑡≥0,(1.5) for (1.1), and 𝑢(𝑥,𝑦,0)=𝜙1(𝑥,𝑦),𝑢𝑡(𝑥,𝑦,0)=𝜙2(𝑥,𝑦),0<𝑥,𝑦<1(1.6) with the boundary conditions 𝑢(𝑥,𝑦,𝑡)=𝑝1(𝑥,𝑦,𝑡)at𝑥,𝑦=±1,𝑡≥0,(1.7) or 𝑢𝑥(𝑥,𝑦,𝑡)=𝑝2𝑢(𝑥,𝑦,𝑡)at𝑥=±1,𝑡≥0,𝑦(𝑥,𝑦,𝑡)=𝑝3(𝑥,𝑦,𝑡)at𝑦=±1,𝑡≥0(1.8) for (1.2).

The functions 𝑓(𝑥,𝑡),ℎ1(𝑡),ℎ2(𝑡),𝑣1(𝑥),𝑣2(𝑥),𝑘1(𝑡),𝑘2(𝑡) are continuous functions defined on Ω1, and similarly 𝑔(𝑥,𝑦,𝑡),𝜙1(𝑥,𝑦),𝜙2(𝑥,𝑦),𝑝1(𝑥,𝑦,𝑡),𝑝2(𝑥,𝑦,𝑡),𝑝3(𝑥,𝑦,𝑡) are defined and continuous on Ω2.

For one-dimensional telegraph equation there are several studies which give numerical solutions. Mohebbi and Dehghan [1] give a compact finite difference approximation of fourth order in space and use collocation method for the time direction. For large values of time level they are forced to increase the degree of polynomial in collocation approach to obtain accurate results. In [2], a scheme similar to finite difference method is proposed using collocation points and approximating the solution with thin plate splines radial basis functions. Dehghan and Lakestani [3] and Saadatmandi and Dehghan [4] make use of Chebyshev cardinal functions and shifted Chebyshev polynomials, respectively, for expanding the approximate solution of one-dimensional hyperbolic telegraph equation. In [4], the advantage is obtaining the closed form of the approximate solution. Dehghan and Ghesmati [5] have applied also dual reciprocity boundary element method (DRBEM) for solving second-order one-dimensional hyperbolic telegraph equation. For time discretization, Crank-Nicolson finite difference has been used, and the solution has been obtained iteratively.

Dehghan and Mohebbi [6] have used the same idea given in [1] for solving the two-dimensional linear hyperbolic telegraph equation which is the combination of finite difference in space and collocation in time directions. In [7], Mohanty and Jain introduced a new unconditionally stable alternating direction implicit (ADI) scheme of second order accurate for two-dimensional telegraph equation. The solution is progressed in time direction by splitting the systems in x- and y-directions and solving two systems for each time level. Mohanty [8] and Mohanty et al. [9] extended their studies for linear hyperbolic equations with variable coefficients in two-space dimensions. Mohanty [10] and Mohanty et al. [11] extended also ADI method solution procedure to two- and three-space dimensional telegraph equations. Ding and Zhang [12] proposed a three-level compact difference scheme of fourth order for the solution of two-dimensional, second-order, nonhomogenous linear hyperbolic equation for positive coefficients. Meshless method has also been used by Dehghan and Ghesmati [13] and Dehghan and Shokri [14] for solving two-dimensional telegraph equation. The conventional moving least squares approximation is exploited in order to interpolate the solution by using monomials from the Pascal triangle in [13], and thin plate splines radial basis functions are used for the approximation of the solution in [14]. In both of these studies, another time integration scheme has been used (finite difference) for time derivatives, and the solutions are obtained iteratively.

The usage of DQM in both space and time directions is encountered in [15]. The authors proposed to solve time-dependent problems (first order time derivatives) by a block-marching methodology in time direction. In each time block, DQM is applied both in space and time directions. The novelty of this approach is in the higher order of accuracy and less computational effort compared to 4-stage Runge-Kutta method.

In this paper, we employ the differential quadrature method in both time and space directions to obtain numerical solutions of one- and two-dimensional linear hyperbolic telegraph equations which contain second-order time derivatives. The use of DQM in time direction, which also discretizes the initial condition 𝑢𝑡, automatically results in an overdetermined system. The numerical scheme then provides the solution at any time level without an iteration. This makes the main difference from the conventional time integration methods. The solution is obtained directly at all required time levels by solving one overdetermined system which contains the solution at the grid points in space and time directions. The Gauss-Chebyshev-Lobatto (GCL) points are used in space direction whereas either equally spaced or GCL grid points can be taken in time direction. The numerical procedure requires very small number of grid points in space directions and appropriate number of time grid points for reaching a certain time level.

The organization of the paper is as follows. In Section 2, we introduce the polynomial and Fourier-based DQM and discuss the formulation of the method both in time and space directions for (1.1)-(1.2). In Section 3, we discuss the accuracy and efficiency of the proposed method applying to some test problems in one- and two-dimensional cases. Section 4 gives the conclusion about the findings.

2. Method of Solution

2.1. Differential Quadrature Method

The differential quadrature method approximates the derivative of a function at a grid point by a linear summation of all the functional values in the whole computational domain.

When the function 𝑓(𝑥) is approximated by a high-order polynomial, Shu [16] presented some explicit formulations to compute the weighting coefficients. The weighting coefficients in that weighted sum are determined using the Lagrange interpolation polynomial which has no limitation on the choice of the grid points. This leads to the polynomial-based differential quadrature method (PDQ). When Fourier series expansion is used, Fourier expansion-based differential quadrature (FDQ) is presented.

Suppose the degree of the polynomial is 𝑁−1, then a single variable function 𝑓(𝑥) and its first- and second-order derivatives can be approximated at a grid point 𝑥𝑖 by PDQ approach as 𝑓𝑥𝑖=𝑁𝑗=1𝑤𝑗𝑥𝑖𝑓𝑥𝑗,𝑓𝑥𝑥𝑖=𝑁𝑗=1𝑤𝑗(1)𝑥𝑖𝑓𝑥𝑗,𝑓𝑥𝑥𝑥𝑖=𝑁𝑗=1𝑤𝑗(2)𝑥𝑖𝑓𝑥𝑗,(2.1) where 𝑖,𝑗=1,2,…,𝑁 and 𝑁 is the number of grid points in the whole (one-dimensional) domain, and 𝑤𝑗(𝑥) are the Lagrange interpolated polynomials. The weighting coefficients for the first-order derivative are [16]𝑤𝑗(1)𝑥𝑖=𝑀(1)𝑥𝑖𝑥𝑖−𝑥𝑗𝑀(1)𝑥𝑗𝑤,𝑖≠𝑗,(2.2a)𝑖(1)𝑥𝑖=−𝑁𝑗=1,𝑗≠𝑖𝑤𝑗(1)𝑥𝑖,(2.2b)where 𝑀(1)𝑥𝑗=𝑁𝑘=1,𝑘≠𝑗𝑥𝑗−𝑥𝑘.(2.3)

The weighting coefficients for the second-order derivative are given in terms of weighting coefficients for the first-order derivatives as𝑤𝑗(2)𝑥𝑖=2𝑤𝑗(1)𝑥𝑖𝑤𝑖(1)𝑥𝑖−1𝑥𝑖−𝑥𝑗𝑤,𝑖≠𝑗,(2.4a)𝑖(2)𝑥𝑖=−𝑁𝑗=1,𝑗≠𝑖𝑤𝑗(2)𝑥𝑖.(2.4b)As it can be seen from these equations, there is a recurrence relationship for the second and so higher-order derivatives.

When the function 𝑓(𝑥) on any interval [ğ‘Ž,𝑏] is approximated by a Fourier series expansion in the form 𝑓(𝑥)=𝑐0+𝑁/2𝑘=1𝑐𝑘cosğ‘˜ğœ‹ğ‘¥ğ‘âˆ’ğ‘Ž+𝑑𝑘sin𝑘𝜋𝑥,ğ‘âˆ’ğ‘Ž(2.5) the Lagrange interpolated trigonometric polynomials are taken as 𝑤𝑘(𝑥)=𝑆(𝑥)𝑃𝑥𝑘𝜋sin𝑥−𝑥𝑘/2(ğ‘âˆ’ğ‘Ž),𝑘=0,1,…,𝑁,(2.6) where 𝑆(𝑥)=𝑁𝑘=0𝜋sin𝑥−𝑥𝑘𝑥2(ğ‘âˆ’ğ‘Ž),𝑃𝑖=𝑁𝑗=0,𝑖≠𝑗𝜋𝑥sin𝑖−𝑥𝑗.2(ğ‘âˆ’ğ‘Ž)(2.7)

Thus, we can write 𝑓𝑥𝑖=𝑁𝑘=1𝑤𝑘𝑥𝑖𝑓𝑥𝑘,𝑖=1,2,…,𝑁.(2.8)

In FDQ approach for nonperiodic problems, weighting coefficients for the first- and second-order derivatives are, respectively: 𝑤𝑗(1)𝑥𝑖=𝜋𝑃𝑥2(ğ‘âˆ’ğ‘Ž)𝑖𝑃𝑥𝑗𝜋𝑥sin𝑖−𝑥𝑗𝑤/2(ğ‘âˆ’ğ‘Ž),𝑖≠𝑗,𝑗(2)𝑥𝑖=𝑤𝑗(1)𝑥𝑖2𝑤𝑖(1)𝑥𝑖−𝜋(ğœ‹î€·ğ‘¥ğ‘âˆ’ğ‘Ž)cot𝑖−𝑥𝑗𝑤2(ğ‘âˆ’ğ‘Ž),𝑖≠𝑗,𝑖(1)𝑥𝑖=−𝑁𝑗=1,𝑖≠𝑗𝑤𝑗(1)𝑥𝑖,𝑤𝑖(2)𝑥𝑖=−𝑁𝑗=1,𝑗≠𝑖𝑤𝑗(2)𝑥𝑖,(2.9) where 𝑖,𝑗=1,2,…,𝑁.

Different grid distributions may be chosen. In this study, we choose the Gauss-Chebyshev-Lobatto (GCL) grid points in space direction which are clustered through the end points and equally spaced (ES) or GCL grid points in time direction. GCL points are the extremal values of Chebyshev polynomials −1≤𝑇𝑘(𝑥)≤1 and are defined as 𝑥𝑗=cos𝑗𝜋𝑘,𝑗=0,1,…,𝑘,(2.10) in [−1,1] where 𝑘 is the degree of the Chebyshev polynomial.

2.2. Application of DQM To Time and Space Directions for Hyperbolic Telegraph Equation

Application to one- and two-dimensional hyperbolic telegraph equations is going to be shown using PDQ method. FDQ application will be similar, the only difference being in the weighting coefficients.

2.2.1. One-Dimensional Hyperbolic Telegraph Equation

PDQ approximations for the derivatives in (1.1) can be written as 𝑢𝑥=𝑁𝑘=1ğ‘Žğ‘–ğ‘˜ğ‘¢ğ‘˜ğ‘™,𝑢𝑥𝑥=𝑁𝑘=1𝑏𝑖𝑘𝑢𝑘𝑙,𝑢𝑡=𝐿𝑘=1ğ‘Žğ‘™ğ‘˜ğ‘¢ğ‘–ğ‘˜,𝑢𝑡𝑡=𝑁𝑘=1𝑏𝑙𝑘𝑢𝑖𝑘,(2.11) where 𝑖=1,2,…,𝑁;𝑙=1,2,…,𝐿, 𝑁 and 𝐿 are the number of discretization points, and ğ‘Žğ‘–ğ‘˜,𝑏𝑖𝑘, and ğ‘Žğ‘™ğ‘˜, 𝑏𝑙𝑘 are the weighting coefficients for the first- and second-order derivatives, in space and time directions, respectively. These coefficients are computed by (2.2a)–(2.4b) when PDQ is used. Therefore, (1.1) will be discretized as 𝐿𝑘=1𝑏𝑙𝑘𝑢𝑖𝑘+2𝛼𝐿𝑘=1ğ‘Žğ‘™ğ‘˜ğ‘¢ğ‘–ğ‘˜+𝛽2𝑢𝑖𝑙+𝑁𝑘=1𝑏𝑖𝑘𝑢𝑘𝑙𝑥=𝑓𝑖,𝑡𝑙,(2.12) where 𝑖=1,2,…,𝑁 and 𝑙=1,2,…,𝐿.

GCL grids in 𝑥-space on an interval [ğ‘Ž,𝑏] are computed as 𝑥𝑖=(𝑏+ğ‘Ž)2−(ğ‘âˆ’ğ‘Ž)2cos((𝑖−1)𝜋)(𝑁−1),𝑖=1,2,…,𝑁.(2.13) In time direction, equally spaced (ES) grid points on an interval [0,𝑇] as 𝑡𝑙=(𝑙−1)𝑇(𝐿−1),𝑙=1,2,…,𝐿,(2.14) are used.

The main system equation (2.12) may be denoted as an algebraic system [𝐵]{𝑢}={𝑓},(2.15) where 𝑢 is the unknown vector to be determined with the entries at the grid points (𝑥𝑖,𝑡𝑙), and the matrix [𝐵] consists of the known weighting coefficients. The known vector 𝑓 contains the function values 𝑓(𝑥𝑖,𝑡𝑙) as entries. The initial condition 𝑢(𝑥,0)=𝑣1(𝑥) is inserted to the system (2.15) directly modifying the matrix [𝐵] and the vector 𝑓. The other initial condition 𝑢𝑡(𝑥,0) is expanded by DQM formulation as 𝑣2𝑥𝑖=𝑢𝑡𝑥𝑖=,0𝐿𝑘=1ğ‘Ž1𝑘𝑢𝑖𝑘,𝑖=1,2,…,𝑁.(2.16) Equation (2.16) will also be a system which can be described as 𝐵0𝑣{𝑢}=2,(2.17) where the matrix [𝐵0] and the vector {𝑣2} contain coefficients ğ‘Ž1𝑘 and 𝑣2(𝑥𝑖) as entries.

If the Dirichlet type boundary conditions are given, these conditions are directly inserted to (2.15) and (2.17) with the given initial conditions 𝑢(𝑥,0). So, the coefficient matrices [𝐵] and [𝐵0] will be of size (𝑁−2)(𝐿−1)×(𝑁−2)(𝐿−1) and (𝑁−2)×(𝑁−2)(𝐿−1), respectively.

The systems (2.15) and (2.17) form an overdetermined system. Therefore, least square method or QR factorization will be made use of for obtaining the solution vector {𝑢}.

If the Neumann type boundary conditions are given, the sizes of the new coefficient matrix [𝐵neu] in (2.15) and [𝐵neu0] in (2.17) will be (𝐿−1)𝑁×(𝐿−1)𝑁 and 𝑁×(𝐿−1)𝑁, respectively. Neumann boundary conditions are also discretized using PDQ as 𝑘1𝑡𝑙=ğ‘¢ğ‘¥î€·ğ‘Ž,𝑡𝑙=𝑁𝑘=1ğ‘Ž1𝑘𝑢𝑘𝑙,𝑘2𝑡𝑙=𝑢𝑥𝑏,𝑡𝑙=𝑁𝑘=1ğ‘Žğ‘ğ‘˜ğ‘¢ğ‘˜ğ‘™,𝑙=1,…,𝐿(2.18) which can be formed as system 𝐵𝜕Ω1{𝑢}={𝑘},(2.19) where [𝐵𝜕Ω1] is a matrix of size 2(𝐿−1)×(𝐿−1)𝑁 containing the coefficients ğ‘Ž1𝑘,ğ‘Žğ‘ğ‘˜ and the vector {𝑘} contains 𝑘1(𝑡𝑖),𝑘2(𝑡𝑖) as entries. Then, the whole system (2.15) with [𝐵neu], (2.17) with [𝐵neu0] and (2.19) again is an overdetermined system which will be solved either by QR or least square method.

Also, the ordering of the unknown vector is important since the structure of [𝐵] or [𝐵neu] depends on this ordering. To get a well-conditioned matrix [𝐵], the order of the unknown vector {𝑢} is arranged as 𝑢{𝑢}=𝑖1,𝑢𝑖2,𝑢𝑖3,…,𝑢𝑖𝐿𝑇,(2.20) where 𝑖=2,3,…,𝑁−1 if boundary conditions are Dirichlet type and 𝑖=1,2,…,𝑁 if Neumann type boundary conditions are given.

2.2.2. Two-Dimensional Hyperbolic Telegraph Equation

The space and time derivatives in (1.2) can be discretized by using PDQ as 𝑢𝑥=𝑁𝑘=1ğ‘Žğ‘–ğ‘˜ğ‘¢ğ‘˜ğ‘—ğ‘™,𝑢𝑥𝑥=𝑁𝑘=1𝑏𝑖𝑘𝑢𝑘𝑗𝑙,𝑢𝑦=𝑀−𝑘=1ğ‘Žğ‘—ğ‘˜ğ‘¢ğ‘–ğ‘˜ğ‘™,𝑢𝑦𝑦=𝑀−𝑘=1𝑏𝑗𝑘𝑢𝑖𝑘𝑙,𝑢𝑡=𝐿𝑘=1ğ‘Žğ‘™ğ‘˜ğ‘¢ğ‘–ğ‘—ğ‘˜,𝑢𝑡𝑡=𝐿𝑘=1𝑏𝑙𝑘𝑢𝑖𝑗𝑘,(2.21) and the discretized form of (1.2) takes the form 𝐿𝑘=1𝑏𝑙𝑘𝑢𝑖𝑗𝑘+2𝛼𝐿𝑘=1ğ‘Žğ‘™ğ‘˜ğ‘¢ğ‘–ğ‘—ğ‘˜+𝛽2𝑢𝑖𝑗𝑙−𝑁𝑘=1𝑏𝑖𝑘𝑢𝑘𝑗𝑙−𝑀−𝑘=1𝑏𝑗𝑘𝑢𝑖𝑘𝑙𝑥=𝑔𝑖,𝑦𝑗,𝑡𝑙,(2.22) where 𝑁,𝑀,𝐿 are the number of discretization points in 𝑥,𝑦,𝑡 directions and 𝑖=1,2,…,𝑁;𝑗=1,2,…,𝑀;𝑙=1,2,…,𝐿, respectively.

GCL grids for 𝑥- and 𝑦-spaces are taken as 𝑥𝑖=121−cos𝑖−1𝜋𝑦𝑁−1,𝑖=1,2,…,𝑁,𝑗=121−cos𝑗−1𝜋𝑀−1,𝑗=1,2,…,𝑀,(2.23) when 0≤𝑥,𝑦≤1. Furthermore, time direction is divided equally as in (2.14) or by GCL grid points on an interval [0,𝑇] as 𝑡𝑙=𝑇21−cos𝑙−1𝜋𝐿−1,𝑙=1,2,…,𝐿.(2.24)

The initial condition 𝑢𝑡(𝑥,𝑦,0) and the Neumann type boundary conditions are added to the system equation (2.22) discretizing them using PDQ approximations as follows: 𝑢𝑡(𝑥,𝑦,0)=𝐿𝑘=1ğ‘Ž1𝑘𝑢𝑖𝑗𝑘𝑢,𝑙=1,𝑖=1,…,𝑁,𝑗=1,…,𝑀,𝑥(0,𝑦,𝑡)=𝑁𝑘=1ğ‘Ž1𝑘𝑢𝑘𝑗𝑙𝑢,𝑖=1,𝑥(1,𝑦,𝑡)=𝑁𝑘=1ğ‘Žğ‘ğ‘˜ğ‘¢ğ‘˜ğ‘—ğ‘™ğ‘¢,𝑖=𝑁,𝑗=1,…,𝑀,𝑙=1,…,𝐿,𝑦(𝑥,0,𝑡)=𝑀−𝑘=1ğ‘Ž1𝑘𝑢𝑖𝑘𝑙𝑢,𝑗=1,𝑦(𝑥,1,𝑡)=𝑀−𝑘=1ğ‘Žğ‘€ğ‘˜ğ‘¢ğ‘–ğ‘˜ğ‘™,𝑗=𝑀;𝑖=1,…,𝑁,𝑙=1,…,𝐿.(2.25)

As in one-dimensional telegraph equation, the DQM discretized system (2.22) will be solved together with the initial and boundary conditions. If the boundary conditions are of Dirichlet type, they will be inserted to the overdetermined system (2.22) combined with the system resulting from initial conditions. For Neumann boundary conditions, the system (2.22) and (2.25) will be solved together.

One of the difficulties in two-dimensional hyperbolic telegraph equation is that the system to be solved becomes larger as 𝑁,𝑀,𝐿 are increased. This causes more memory and CPU usage. To overcome this problem, the system is reduced by removing the entries on the coefficient matrix of the system which correspond to known information (e.g., initial condition 𝑢(𝑥,𝑦,0) and the Dirichlet type boundary conditions). Meanwhile, the right hand side of the reduced system is also modified taking into account removed known entries in the coefficient matrix.

The order of the unknown vector 𝑢 in this case is organized to get a well-conditioned system. Consider the unknown vector as a matrix 𝑈 whose each row entry corresponds to a time level. Notationally, âŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢1𝑗1𝑢2𝑗1…𝑢𝑁𝑗1𝑢1𝑗2𝑢2𝑗2…𝑢𝑁𝑗2𝑢⋮⋮⋱⋮1𝑗𝐿𝑢2ğ‘—ğ¿â€¦ğ‘¢ğ‘ğ‘—ğ¿âŽ¤âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¦,(2.26) where 𝑗=1,2,…,𝑀, and the matrix of size 𝐿×𝑁𝑀 will be rewritten as a vector writing columns consecutively.

The solution will be obtained by solving only one system with the initial and boundary conditions being all inserted, and the aforementioned reduction of known entries is performed to reduce the size of the system. The solution vector contains all required time level values in it.

The solvability of the overdetermined system of equations depends on the column rank of the coefficient matrix which is 𝑁𝑀𝐿 in this case. When the initial and/or Neumann type boundary conditions are discretized using DQM, and added to the system, the row size is certainly greater than 𝑁𝑀𝐿 which makes the system overdetermined. The choice of the grid points in both space and time domains affects the stability of the system. As mentioned in the Shu’s book [16], the solution with GCL grid points becomes more stable than equally spaced grid points in both space and time directions. Moreover, appropriate choice of 𝑁,𝑀, and 𝐿 makes the final coefficient matrix full column rank.

3. Numerical Results

The discretization is performed by taking 𝑁,𝑀 GCL grid points in x- and y-directions. When equally spaced points are taken in time direction on [0,𝑇], 𝑇𝐿=Δ𝑡+1,(3.1) in which 𝑇 is the required uptime level.

To measure the accuracy of the method, we use the following errors which can be defined as [3, 13] ∑RMSerror=nod𝑖=1(𝑢exact−𝑢computed)2nod1/2,∑Relativeerror=nod𝑖=1(𝑢exact−𝑢computed)2∑nod𝑖=1𝑢2exact1/2,(3.2) where nod is the total number of all grid points in space direction, and the errors can be computed at any time level.

We present some numerical results of one- and two-dimensional hyperbolic telegraph equation for different 𝛼 and 𝛽 values to observe the accuracy and efficiency of the DQ method. Exact solutions of the test problems are available. Thus, the variations of RMS and relative errors can be obtained for discussing accuracy of DQM solutions.

3.1. Problem 1

We consider the one-dimensional hyperbolic telegraph equation (1.1) in the interval 0≤𝑥≤2𝜋, 0<𝑡≤3 with 𝛼=4,𝛽=2. The exact solution is taken as [2] 𝑢(𝑥,𝑡)=𝑒−𝑡sin(𝑥).(3.3) In this case, 𝑓(𝑥,𝑡)=(2−2𝛼+𝛽2)𝑒−𝑡sin(𝑥). The initial conditions are given as 𝑢𝑢(𝑥,0)=sin(𝑥),𝑡(𝑥,0)=−sin(𝑥).(3.4)

We solve the equation first by taking Dirichlet type boundary conditions and then with the Neumann type boundary conditions. The Dirichlet type boundary conditions from the exact solution are [].𝑢(0,𝑡)=0,𝑢(2𝜋,𝑡)=0,𝑡∈0,3(3.5) The Neumann type boundary conditions are 𝑢𝑥(0,𝑡)=𝑒−𝑡,𝑢𝑥(2𝜋,𝑡)=𝑒−𝑡[].,𝑡∈0,3(3.6)

Table 1 shows RMS errors obtained by using Dirichlet type boundary conditions with different Δ𝑡 (equally spaced time grid is used) values. It is noticed that 10−5 accuracy is achieved even with a coarse mesh (Δ𝑡=0.5) in time direction. It is improved with a finer mesh (Δ𝑡=0.25).

Table 2 reports RMS errors obtained by using Neumann boundary conditions with different number of discretization points in [0,2𝜋]. As can be seen from Table 2, increasing 𝑁 from 11 to 21 improves the results in terms of accuracy, which was dropped with Neumann boundary conditions being used, with Δ𝑡=0.25. Table 3 gives RMS errors when FDQ weighting coefficients are used in space direction and PDQ weighting coefficients in time direction for the same Neumann problem. For this particular problem, the accuracy of the results by using FDQ even with 𝑁=11 doubled from the accuracy of the results obtained by PDQ.

Figure 1 depicts the very well agreement of DQM and exact solutions at different time levels; even equally spaced time grid points are used. The number of GCL grid points in space is small (𝑁=21), and Δ𝑡=0.25 is considerably large compared to other time integration schemes.

3.2. Problem 2

We consider now two-dimensional telegraph equation (1.2) in the region 0≤𝑥,𝑦≤1,𝑡>0, with 𝛼=𝛽=1. The analytical solution given by [13] is 𝑢(𝑥,𝑦,𝑡)=cos(𝑡)sin(𝑥)sin(𝑦),(3.7) with the initial conditions 𝑢𝑢(𝑥,𝑦,0)=sin(𝑥)sin(𝑦),𝑡(𝑥,𝑦,0)=0,(3.8) and the boundary conditions are Dirichlet type 𝑢𝑢(𝑥,0,𝑡)=0,0≤𝑥≤1,𝑦=0,(0,𝑦,𝑡)=0,0≤𝑦≤1,𝑥=0,𝑢(1,𝑦,𝑡)=cos(𝑡)sin(1)sin(𝑦),0≤𝑦≤1,𝑡≥0,𝑢(𝑥,1,𝑡)=cos(𝑡)sin(𝑥)sin(1),0≤𝑥≤1,𝑡≥0.(3.9)

We extract 𝑔(𝑥,𝑦,𝑡) from the analytical solution as [].𝑔(𝑥,𝑦,𝑡)=2sin(𝑥)sin(𝑦)cos(𝑡)−sin(𝑡)(3.10)

Tables 4 and 5 present the comparison of the DQM solution and the exact solution in terms of relative errors for different Δ𝑡 values when equally spaced time grid points are used. In this two-dimensional problem, FDQ approximation for space derivatives, keeping PDQ approximation for time derivatives, is also studied. As can be seen from the Tables 4 and 5, both PDQ and FDQ approximations in space directions give almost the same accuracy even with coarse grids (Δ𝑡=0.5 and Δ𝑡=0.25) taking a little more space grid points in FDQ. Furthermore, as 𝐿 is increased (Δ𝑡 is decreased), the accuracy is increased.

3.3. Problem 3

In this case, (1.2) is taken with 𝛼=𝛽=1 and 𝛼=5,𝛽=1 in the region 0≤𝑥,𝑦≤1,𝑡>0. The analytical solution is given by [6] 𝑢(𝑥,𝑦,𝑡)=ln(1+𝑥+𝑦+𝑡),(3.11) with the initial conditions 𝑢𝑢(𝑥,𝑦,0)=ln(1+𝑥+𝑦),𝑡1(𝑥,𝑦,0)=,1+𝑥+𝑦(3.12) and the boundary conditions are of Dirichlet and Neumann type on the walls parallel to each other: 𝑢𝑦1(𝑥,0,𝑡)=𝑢1+𝑥+𝑡,0≤𝑥≤1,𝑥1(1,𝑦,𝑡)=2+𝑦+𝑡,0≤𝑦≤1,𝑡≥0,𝑢(𝑥,1,𝑡)=ln(2+𝑥+𝑡),0≤𝑥≤1,𝑢(0,𝑦,𝑡)=ln(1+𝑦+𝑡),0≤𝑦≤1.(3.13)

The inhomogeneity is extracted from the exact solution as 𝑔(𝑥,𝑦,𝑡)=2𝛼1+𝑥+𝑦+𝑡+𝛽21ln(1+𝑥+𝑦+𝑡)+(1+𝑥+𝑦+𝑡)2.(3.14)

Relative errors between numerical solution (PDQ) and exact solution for different Δ𝑡 values are presented in Table 6. Even with a large Δ𝑡=1, about 10−4 accuracy is reached. As Δ𝑡 is decreased, accuracy is increased as expected. In Table 7, the error decreases with a higher 𝛼 value. Again, it is noticed that with considerably small number of grid points both in space direction (𝑁=𝑀=11) and time direction (Δ𝑡=0.5  or  Δ𝑡=0.25), at least 10−6 accuracy is obtained.

3.4. Problem 4

Consider the homogenous (𝑔(𝑥,𝑦,𝑡)=0) (1.2) in the region 0≤𝑥,𝑦≤1,𝑡>0, with 𝛼=1+𝜋2,𝛽=1. The exact solution is given by [12] 𝑢(𝑥,𝑦,𝑡)=𝑒−𝑡sin(𝜋𝑥)sin(𝜋𝑦),(3.15) with the initial conditions 𝑢𝑢(𝑥,𝑦,0)=sin(𝜋𝑥)sin(𝜋𝑦),𝑡(𝑥,𝑦,0)=−sin(𝜋𝑥)sin(𝜋𝑦),(3.16) and Dirichlet and Neumann boundary conditions on parallel walls are 𝑢𝑢𝑢(𝑥,0,𝑡)=0,0≤𝑥≤1,(1,𝑦,𝑡)=0,0≤𝑦≤1,𝑡≥0,𝑦(𝑥,1,𝑡)=−𝜋𝑒−𝑡𝑢sin(𝜋𝑥),0≤𝑥≤1,𝑥(0,𝑦,𝑡)=𝜋𝑒−𝑡sin(𝜋𝑦),0≤𝑦≤1.(3.17)

Table 8 compares the numerical solution (both PDQ and FDQ) with the exact solution for different time levels in terms of RMS errors. FDQ weighting coefficients are used in space directions remaining PDQ weighting coefficients in time direction (FDQ-PDQ) and PDQ approximation is used in both directions (PDQ-PDQ). In this problem, we use GCL grid points in both space and time directions. FDQ-PDQ approximation gives better accuracy than PDQ-PDQ approximation for this problem.

To emphasize the importance of the GCL grid points in time direction, we present RMS errors versus the number of grid points in time direction in Figures 2 and 3. The number of GCL grid points in space direction (𝑁=𝑀=11) at 𝑇=3 is fixed, and only PDQ-PDQ approximation is considered. When the number of equally spaced points in time exceeds 28, the system becomes rank deficient. This means that large number of ES grid points in time direction causes unstable solution. On the other hand, GCL grid points in time direction still give very good accuracy for 11≤𝐿≤41 as can be seen in Figure 3. Moreover, RMS error using 11 to 23 ES grid points increases faster than RMS error using the same number of GCL grid points in time. The accuracy with GCL grid points in time direction remains in a scale between 0.2×10−9 and 1.8×10−9.

This observation is important for physical problems which require the solution at a high time level. In that case, the number of grid points in time direction should be large and endure oscillations of the solution. Thus, GCL grid points in time are more preferable than ES grid points.

4. Conclusion

The polynomial-based or Fourier expansion-based DQM in space domain and polynomial-based DQM in time domain have been proposed to solve one- and two-dimensional hyperbolic telegraph equations. Differential quadrature method has the capability of producing highly accurate results using considerably small number of grid points, and thus resulting in minimal computational effort. We have used Gauss-Chebyshev-Lobatto grid in space directions which is a nonuniform grid distribution clustering through the end points. For the time direction, equally spaced points may be used for achieving the solution at a specific time level. Considering stability, GCL grids in time direction are more reliable than ES grids when the solution is required at a high time level. The accuracy and efficiency of the results have been demonstrated by root mean square (RMS) and relative errors. The use of FDQ increases the accuracy for one-dimensional telegraph equation especially when the exact solution contains trigonometric functions. For the two-dimensional hyperbolic telegraph equations, the combination of FDQ-PDQ in space and time directions is preferable.

The advantage of the use of DQM both in space and time directions lies in the fact that the solution can be obtained at a required time level by solving one system. There is no need to employ an iteration between the time levels, and large time steps can be used.