Differential Quadrature Solution of Hyperbolic Telegraph Equation

Abstract

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𝛼𝑢_𝑡(𝑥,𝑡)+𝛽²𝑢(𝑥,𝑡)=𝑢_𝑥𝑥(𝑥,𝑡)+𝑓(𝑥,𝑡),inΩ₁,(1.1) and the second-order two-dimensional hyperbolic telegraph equation: 𝑢_𝑡𝑡(𝑥,𝑦,𝑡)+2𝛼𝑢_𝑡(𝑥,𝑦,𝑡)+𝛽²𝑢(𝑥,𝑦,𝑡)−𝑔(𝑥,𝑦,𝑡)=𝑢_𝑥𝑥(𝑥,𝑦,𝑡)+𝑢_𝑦𝑦(𝑥,𝑦,𝑡)inΩ₂,(1.2) where Ω₁=[𝑎<𝑥<𝑏]×[𝑡>0] and Ω₂=[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)=𝑣₁(𝑥),𝑢_𝑡(𝑥,0)=𝑣₂(𝑥),𝑎<𝑥<𝑏(1.3) with Dirichlet type boundary conditions 𝑢(𝑎,𝑡)=ℎ₁(𝑡),𝑢(𝑏,𝑡)=ℎ₂(𝑡),𝑡≥0,(1.4) or Neumann type boundary conditions 𝑢_𝑥(𝑎,𝑡)=𝑘₁(𝑡),𝑢_𝑥(𝑏,𝑡)=𝑘₂(𝑡),𝑡≥0,(1.5) for (1.1), and 𝑢(𝑥,𝑦,0)=𝜙₁(𝑥,𝑦),𝑢_𝑡(𝑥,𝑦,0)=𝜙₂(𝑥,𝑦),0<𝑥,𝑦<1(1.6) with the boundary conditions 𝑢(𝑥,𝑦,𝑡)=𝑝₁(𝑥,𝑦,𝑡)at𝑥,𝑦=±1,𝑡≥0,(1.7) or 𝑢𝑢_𝑥(𝑥,𝑦,𝑡)=𝑝₂(𝑥,𝑦,𝑡)at𝑥=±1,𝑡≥0,_𝑦(𝑥,𝑦,𝑡)=𝑝₃(𝑥,𝑦,𝑡)at𝑦=±1,𝑡≥0(1.8) for (1.2).

The functions 𝑓(𝑥,𝑡),ℎ₁(𝑡),ℎ₂(𝑡),𝑣₁(𝑥),𝑣₂(𝑥),𝑘₁(𝑡),𝑘₂(𝑡) are continuous functions defined on Ω₁, and similarly 𝑔(𝑥,𝑦,𝑡),𝜙₁(𝑥,𝑦),𝜙₂(𝑥,𝑦),𝑝₁(𝑥,𝑦,𝑡),𝑝₂(𝑥,𝑦,𝑡),𝑝₃(𝑥,𝑦,𝑡) are defined and continuous on Ω₂.

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)𝑥_𝑖=𝑀⁽¹⁾𝑥_𝑖𝑤𝑥_𝑖−𝑥_𝑗𝑀⁽¹⁾𝑥_𝑗,𝑖≠𝑗,(2.2a)_𝑖(1)𝑥_𝑖=−^𝑁_{𝑗=1,𝑗≠𝑖}𝑤^𝑗(1)𝑥_𝑖,(2.2b)where 𝑀⁽¹⁾𝑥_𝑗=^𝑁_𝑘=_{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 𝑓(𝑥)=𝑐₀+^𝑁/2_𝑘=1𝑐^𝑘cos𝑘𝜋𝑥𝑏−𝑎+𝑑^𝑘sin𝑘𝜋𝑥𝑏−𝑎,(2.5) the Lagrange interpolated trigonometric polynomials are taken as 𝑤_𝑘(𝑥)=𝑆(𝑥)𝑃𝑥_𝑘sin𝜋𝑥−𝑥_𝑘/2(𝑏−𝑎),𝑘=0,1,…,𝑁,(2.6) where 𝑆(𝑥)=^𝑁_𝑘=0sin𝜋𝑥−𝑥^𝑘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)𝑤_𝑖(2)𝑥_𝑖=−^𝑁_{𝑗=1,𝑖≠𝑗}𝑤^𝑗(1)𝑥_𝑖,𝑥_𝑖=−^𝑁_{𝑗=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𝑎_𝑙𝑘𝑢_𝑖𝑘+𝛽²𝑢_𝑖𝑙+^𝑁_𝑘=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)=𝑣₁(𝑥) 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 𝑣₂𝑥_𝑖=𝑢_𝑡𝑥_𝑖,0=^𝐿_𝑘=1𝑎_1𝑘𝑢_𝑖𝑘,𝑖=1,2,…,𝑁.(2.16) Equation (2.16) will also be a system which can be described as 𝐵₀{𝑢}=𝑣₂,(2.17) where the matrix [𝐵₀] and the vector {𝑣₂} contain coefficients 𝑎_1𝑘 and 𝑣₂(𝑥_𝑖) 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 [𝐵₀] 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 [𝐵^neu₀] in (2.17) will be (𝐿−1)𝑁×(𝐿−1)𝑁 and 𝑁×(𝐿−1)𝑁, respectively. Neumann boundary conditions are also discretized using PDQ as 𝑘₁𝑡_𝑙=𝑢_𝑥𝑎,𝑡_𝑙=^𝑁_𝑘=1𝑘₂𝑎_1𝑘𝑢_𝑘𝑙,𝑡_𝑙=𝑢_𝑥𝑏,𝑡_𝑙=^𝑁_𝑘=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 𝑘₁(𝑡_𝑖),𝑘₂(𝑡_𝑖) as entries. Then, the whole system (2.15) with [𝐵^neu], (2.17) with [𝐵^neu₀] 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𝑎_𝑙𝑘𝑢_𝑖𝑗𝑘+𝛽²𝑢_𝑖𝑗𝑙−^𝑁_𝑘=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 𝑥_𝑖=1𝑦21−cos𝑁𝑖−−11𝜋,𝑖=1,2,…,𝑁,_𝑗=121−cos𝑀𝑗−−11𝜋,𝑗=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𝐿𝑙−−11𝜋,𝑙=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)²nod^1/2Relativeerror=∑_nod𝑖=1,(𝑢_exact−𝑢_computed)²∑^nod_𝑖=1𝑢²_exact^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𝛼+𝛽²)𝑒^−𝑡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⁻⁵ accuracy is achieved even with a coarse mesh (Δ𝑡=0.5) in time direction. It is improved with a finer mesh (Δ𝑡=0.25).

Table 1 

RMS errors for Problem 1 with 𝑁=11 using Dirichlet boundary conditions.

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.

Table 2 

RMS errors for Problem 1 with Δ𝑡=0.25 using PDQ with Neumann boundary conditions.

Table 3 

RMS errors for Problem 1 with 𝑁=11 using FDQ with Neumann boundary conditions.

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.

Figure 1 

Numerical and exact solutions at different times for Problem 1, with Δ𝑡=0.25, 𝑁=17 using Neumann boundary conditions.

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 𝑢𝑢((𝑥1,,1𝑦,,𝑡𝑡))==𝑢𝑢((cc𝑥0oo,,ss𝑦0,,((𝑡𝑡𝑡𝑡))))ss==iinn00((,,𝑥1))00ssiinn≤≤((𝑦𝑥𝑦1))≤≤,,11,,00𝑥𝑦≤≤==𝑦𝑥00≤≤,,11,,𝑡𝑡≥≥00,.(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.

Table 4 

Relative errors of Problem 2 with 𝑁=𝑀=17 using PDQ.

Table 5 

Relative errors of Problem 2 using 𝑁=𝑀=21 using FDQ.

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+𝑥+𝑦),_𝑡(𝑥,𝑦,0)=11+𝑥+𝑦,(3.12) and the boundary conditions are of Dirichlet and Neumann type on the walls parallel to each other: 𝑢_𝑦(𝑥,0,𝑡)=1𝑢1+𝑥+𝑡,0≤𝑥≤1,_𝑥(1,𝑦,𝑡)=1𝑢𝑢((𝑥0,,1𝑦,,𝑡𝑡))==2llnn+((𝑦21+++𝑡𝑥𝑦,++𝑡𝑡0)),,≤𝑦00≤≤≤1𝑦𝑥,≤≤𝑡≥11.,0,(3.13)

The inhomogeneity is extracted from the exact solution as 𝑔(𝑥,𝑦,𝑡)=2𝛼1+𝑥+𝑦+𝑡+𝛽²ln(1+𝑥+𝑦+𝑡)+1(1+𝑥+𝑦+𝑡)².(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⁻⁴ 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⁻⁶ accuracy is obtained.

Table 6 

Relative errors of Problem 3 with 𝑁=𝑀=11,𝛼=𝛽=1.

Table 7 

Relative errors of Problem 3 with 𝑁=𝑀=11,𝛼=5,𝛽=1.

3.4. Problem 4

Consider the homogenous (𝑔(𝑥,𝑦,𝑡)=0) (1.2) in the region 0≤𝑥,𝑦≤1,𝑡>0, with 𝛼=1+𝜋²,𝛽=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 𝑢𝑢(1𝑢,𝑦(𝑥,𝑡,)0=,𝑡)0=,00,≤0𝑦≤≤𝑥1≤,𝑡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.

Table 8 

RMS errors for 𝑢(𝑥,𝑡) of Problem 5 using 𝑀=𝑁=11,𝐿=13.

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⁻⁹ and 1.8×10⁻⁹.

(a)

(b)

Figure 2 

RMS errors at 𝑇 versus ES grid points in time direction for Problem 4.

Figure 3 

RMS errors at 𝑇 versus GCL grid points in time direction for Problem 4.

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.