About this Journal Submit a Manuscript Table of Contents
Advances in Mathematical Physics
VolumeΒ 2011Β (2011), Article IDΒ 420608, 22 pages
http://dx.doi.org/10.1155/2011/420608
Research Article

A New High-Order Approximation for the Solution of Two-Space-Dimensional Quasilinear Hyperbolic Equations

1Department of Mathematics, Faculty of Mathematical Sciences, University of Delhi, Delhi 110 007, India
2Department of Mathematics, Aditi Mahavidyalaya, University of Delhi, Delhi 110 039, India

Received 31 March 2011; Revised 13 June 2011; Accepted 5 July 2011

Academic Editor: RicardoΒ Weder

Copyright Β© 2011 R. K. Mohanty and Suruchi Singh. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

we propose a new high-order approximation for the solution of two-space-dimensional quasilinear hyperbolic partial differential equation of the form 𝑒𝑑𝑑=𝐴(π‘₯,𝑦,𝑑,𝑒)𝑒π‘₯π‘₯+𝐡(π‘₯,𝑦,𝑑,𝑒)𝑒𝑦𝑦+𝑔(π‘₯,𝑦,𝑑,𝑒,𝑒π‘₯,𝑒𝑦,𝑒𝑑), 0<π‘₯, 𝑦<1, 𝑑>0 subject to appropriate initial and Dirichlet boundary conditions , where π‘˜>0 and β„Ž>0 are mesh sizes in time and space directions, respectively. We use only five evaluations of the function 𝑔 as compared to seven evaluations of the same function discussed by (Mohanty et al., 1996 and 2001). We describe the derivation procedure in details and also discuss how our formulation is able to handle the wave equation in polar coordinates. The proposed method when applied to a linear hyperbolic equation is also shown to be unconditionally stable. Some examples and their numerical results are provided to justify the usefulness of the proposed method.

1. Introduction

We consider the following two-space dimensional quasilinear hyperbolic partial differential equation:𝑒𝑑𝑑=𝐴(π‘₯,𝑦,𝑑,𝑒)𝑒π‘₯π‘₯+𝐡(π‘₯,𝑦,𝑑,𝑒)𝑒𝑦𝑦+𝑔π‘₯,𝑦,𝑑,𝑒,𝑒π‘₯,𝑒𝑦,𝑒𝑑,0<π‘₯,𝑦<1,𝑑>0(1.1) subject to the initial conditions 𝑒(π‘₯,𝑦,0)=πœ™(π‘₯,𝑦),𝑒𝑑(π‘₯,𝑦,0)=πœ“(π‘₯,𝑦),0≀π‘₯,𝑦≀1,(1.2) and the boundary conditions𝑒(0,𝑦,𝑑)=π‘Ž0(𝑦,𝑑),𝑒(1,𝑦,𝑑)=π‘Ž1(𝑦,𝑑),0≀𝑦≀1,𝑑β‰₯0,(1.3a)𝑒(π‘₯,0,𝑑)=𝑏0(π‘₯,𝑑),𝑒(π‘₯,1,𝑑)=𝑏1(π‘₯,𝑑),0≀π‘₯≀1,𝑑β‰₯0,(1.3b)where (1.1) is assumed to satisfy the hyperbolicity condition 𝐴(π‘₯,𝑦,𝑑,𝑒)>0 and 𝐡(π‘₯,𝑦,𝑑,𝑒)>0 in the solution region Ω≑{(π‘₯,𝑦,𝑑)∢0<π‘₯,𝑦<1,𝑑>0}. Further we assume that 𝑒(π‘₯,𝑦,𝑑)∈𝐢6, 𝐴(π‘₯,𝑦,𝑑,𝑒),𝐡(π‘₯,𝑦,𝑑,𝑒)∈𝐢4, and πœ™(π‘₯,𝑦) and πœ“(π‘₯,𝑦) are sufficiently differentiable function of as higher-order as possible.

Second-order quasilinear hyperbolic partial differential equation with appropriate initial and boundary conditions serves as models in many branches of physics and technology. Ciment and Leventhal [1, 2] have discussed operator compact implicit method to solve wave equation. Using finite volume technique, Schwartzkopff et al. [3] have studied high-order ADER schemes for first-order linear hyperbolic systems. Even much higher order ADER finite volume schemes have been presented and analyzed by Schwartzkopff et al. [4] on Cartesian meshes, and by Dumbser and KΓ€ser [5] on general two- and three-dimensional unstructured meshes. In 1996, Mohanty et al. [6] have developed a high-order numerical method for the solution of two-space dimensional second-order nonlinear hyperbolic equation. Later, Mohanty et al. [7] have extended their technique to solve second-order quasilinear hyperbolic equations. In both cases they have used seven evaluations of the function 𝑔 and 19-grid points. It has been shown that the linear schemes discussed in [6, 7] are conditionally stable. In the recent past, many researchers (see [8–16]) have developed unconditionally stable implicit finite difference methods for the solution of two-space dimensional linear hyperbolic equations with significant first derivative terms. Most recently, Mohanty and Singh [17] have derived a high accuracy numerical method based on Numerov type discretization for the solution of one space dimensional nonlinear hyperbolic equations, in which they have shown that the linear scheme is unconditionally stable.

In this paper, using nineteen grid-points, we derive a new compact three-level implicit numerical method of accuracy two in time and four in space for the solution of two-space dimensional quasilinear hyperbolic equation (1.1). In this method we require only five evaluations of the function 𝑔 as compared to seven evaluations of the same function discussed in [6, 7]. In the next section, we give formulation of the method. In Section 3, we give the complete derivation of the method. In Section 4, we discuss the application of the proposed method to two-space dimensional wave equation in polar coordinates and discuss the stability analysis. In this section, we modify our method in such a way that the solution retains its order and accuracy everywhere in the vicinity of the singularity. In Section 5, we discuss superstable method for two-space dimensional telegraphic equation. In Section 6, we examine our new method over a set of linear and nonlinear second-order hyperbolic equations whose exact solutions are known and compared the results with the results of other known methods. Concluding remarks are given in Section 7.

2. Formulation of the Numerical Method

In this section, we aim to discuss a numerical method for the solution of nonlinear wave equation: 𝑒𝑑𝑑=𝐴(π‘₯,𝑦,𝑑)𝑒π‘₯π‘₯+𝐡(π‘₯,𝑦,𝑑)𝑒𝑦𝑦+𝑔π‘₯,𝑦,𝑑,𝑒,𝑒π‘₯,𝑒𝑦,𝑒𝑑,0<π‘₯,𝑦<1,𝑑>0.(2.1)

Let β„Ž>0 and π‘˜>0 be the mesh spacing in the space and time directions, respectively. We replace the solution region Ξ©={(π‘₯,𝑦,𝑑)∣0<π‘₯,𝑦<1,𝑑>0} by a set of grid points (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), where π‘₯𝑙=π‘™β„Ž, 𝑙=0(1)𝑁+1, π‘¦π‘š=π‘šβ„Ž, π‘š=0(1)𝑁+1, and 𝑑𝑗=π‘—π‘˜, 0<𝑗<𝐽, 𝑁 and 𝐽 being positive integers and (𝑁+1)β„Ž=1. Let 𝑝=π‘˜/β„Ž>0 be the mesh ratio parameter. The exact values of 𝑒(π‘₯,𝑦,𝑑), 𝐴(π‘₯,𝑦,𝑑), and 𝐡(π‘₯,𝑦,𝑑) at the grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗) are denoted by π‘ˆπ‘—π‘™,π‘š, 𝐴𝑗𝑙,π‘š, and 𝐡𝑗𝑙,π‘š, respectively. Similarly at the grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), we denote 𝐴𝑗π‘₯𝑙,π‘š=πœ•π΄π‘—π‘™,π‘š/πœ•π‘₯,𝐴𝑗𝑦𝑙,π‘š=πœ•π΄π‘—π‘™,π‘š/πœ•π‘¦,𝐴𝑗π‘₯π‘₯𝑙,π‘š=πœ•2𝐴𝑗𝑙,π‘š/πœ•π‘₯2,… and so forth. Let 𝑒𝑗𝑙,π‘š be the approximate value of 𝑒(π‘₯,𝑦,𝑑) at the same grid point.

At the grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), we consider the following approximations:π‘ˆπ‘—π‘‘π‘™,π‘š=ξ‚€π‘ˆπ‘—+1𝑙,π‘šβˆ’π‘ˆπ‘—βˆ’1𝑙,π‘šξ‚,(2π‘˜)(2.2a)π‘ˆπ‘—π‘‘π‘™Β±1,π‘š=ξ‚€π‘ˆπ‘—+1𝑙±1,π‘šβˆ’π‘ˆπ‘—βˆ’1𝑙±1,π‘šξ‚(2π‘˜),(2.2b)π‘ˆπ‘—π‘‘π‘™,π‘šΒ±1=ξ‚€π‘ˆπ‘—+1𝑙,π‘šΒ±1βˆ’π‘ˆπ‘—βˆ’1𝑙,π‘šΒ±1(2π‘˜),(2.2c)π‘ˆπ‘—π‘‘π‘‘π‘™,π‘š=ξ‚€π‘ˆπ‘—+1𝑙,π‘šβˆ’2π‘ˆπ‘—π‘™,π‘š+π‘ˆπ‘—βˆ’1𝑙,π‘šξ‚ξ€·π‘˜2ξ€Έ,(2.2d)π‘ˆπ‘—π‘‘π‘‘π‘™Β±1,π‘š=ξ‚€π‘ˆπ‘—+1𝑙±1,π‘šβˆ’2π‘ˆπ‘—π‘™Β±1,π‘š+π‘ˆπ‘—βˆ’1𝑙±1,π‘šξ‚ξ€·π‘˜2ξ€Έ,(2.2e)π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šΒ±1=ξ‚€π‘ˆπ‘—+1𝑙,π‘šΒ±1βˆ’2π‘ˆπ‘—π‘™,π‘šΒ±1+π‘ˆπ‘—βˆ’1𝑙,π‘šΒ±1ξ‚ξ€·π‘˜2ξ€Έ,(2.2f)β€‰π‘ˆπ‘—π‘₯𝑙,π‘š=ξ‚€π‘ˆπ‘—π‘™+1,π‘šβˆ’π‘ˆπ‘—π‘™βˆ’1,π‘šξ‚,(2β„Ž)(2.3a)π‘ˆπ‘—π‘₯𝑙±1,π‘š=ξ‚€Β±3π‘ˆπ‘—π‘™Β±1,π‘šβˆ“4π‘ˆπ‘—π‘™,π‘šΒ±π‘ˆπ‘—π‘™βˆ“1,π‘šξ‚(2β„Ž),(2.3b)π‘ˆπ‘—π‘₯𝑙,π‘šΒ±1=ξ‚€π‘ˆπ‘—π‘™+1,π‘šΒ±1βˆ’π‘ˆπ‘—π‘™βˆ’1,π‘šΒ±1(2β„Ž),(2.3c)π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š=ξ‚€π‘ˆπ‘—π‘™+1,π‘šβˆ’2π‘ˆπ‘—π‘™,π‘š+π‘ˆπ‘—π‘™βˆ’1,π‘šξ‚ξ€·β„Ž2ξ€Έ,(2.3d)π‘ˆπ‘—π‘₯π‘₯𝑙,π‘šΒ±1=ξ‚€π‘ˆπ‘—π‘™+1,π‘šΒ±1βˆ’2π‘ˆπ‘—π‘™,π‘šΒ±1+π‘ˆπ‘—π‘™βˆ’1,π‘šΒ±1ξ‚ξ€·β„Ž2ξ€Έ,(2.3e)β€‰π‘ˆπ‘—π‘¦π‘™,π‘š=ξ‚€π‘ˆπ‘—π‘™,π‘š+1βˆ’π‘ˆπ‘—π‘™,π‘šβˆ’1,(2β„Ž)(2.4a)π‘ˆπ‘—π‘¦π‘™Β±1,π‘š=ξ‚€π‘ˆπ‘—π‘™Β±1,π‘š+1βˆ’π‘ˆπ‘—π‘™Β±1,π‘šβˆ’1(2β„Ž),(2.4b)π‘ˆπ‘—π‘¦π‘™,π‘šΒ±1=ξ‚€Β±3π‘ˆπ‘—π‘™,π‘šΒ±1βˆ“4π‘ˆπ‘—π‘™,π‘šΒ±π‘ˆπ‘—π‘™,π‘šβˆ“1(2β„Ž),(2.4c)π‘ˆπ‘—π‘¦π‘¦π‘™,π‘š=ξ‚€π‘ˆπ‘—π‘™,π‘š+1βˆ’2π‘ˆπ‘—π‘™,π‘š+π‘ˆπ‘—π‘™,π‘šβˆ’1ξ‚ξ€·β„Ž2ξ€Έ,(2.4d)π‘ˆπ‘—π‘¦π‘¦π‘™Β±1,π‘š=ξ‚€π‘ˆπ‘—π‘™Β±1,π‘š+1βˆ’2π‘ˆπ‘—π‘™Β±1,π‘š+π‘ˆπ‘—π‘™Β±1,π‘šβˆ’1ξ‚ξ€·β„Ž2ξ€Έ,(2.4e) 𝐺𝑗𝑙±1,π‘šξ‚€π‘₯=𝑔𝑙±1,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™Β±1,π‘š,π‘ˆπ‘—π‘₯𝑙±1,π‘š,π‘ˆπ‘—π‘¦π‘™Β±1,π‘š,π‘ˆπ‘—π‘‘π‘™Β±1,π‘šξ‚,(2.5a)𝐺𝑗𝑙,π‘šΒ±1ξ‚€π‘₯=𝑔𝑙,π‘¦π‘šΒ±1,𝑑𝑗,π‘ˆπ‘—π‘™,π‘šΒ±1,π‘ˆπ‘—π‘₯𝑙,π‘šΒ±1,π‘ˆπ‘—π‘¦π‘™,π‘šΒ±1,π‘ˆπ‘—π‘‘π‘™,π‘šΒ±1,(2.5b)β€‰π‘ˆπ‘—π‘₯𝑙,π‘š=π‘ˆπ‘—π‘₯𝑙,π‘š+β„Ž16𝐴𝑗𝑙,π‘šξ‚€πΊπ‘—π‘™+1,π‘šβˆ’πΊπ‘—π‘™βˆ’1,π‘šξ‚βˆ’β„Ž16𝐴𝑗𝑙,π‘šξ‚€π‘ˆπ‘—π‘‘π‘‘π‘™+1,π‘šβˆ’π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1,π‘šξ‚+β„Žπ΅π‘—π‘–,π‘š16𝐴𝑗𝑙,π‘šξ‚€π‘ˆπ‘—π‘¦π‘¦π‘™+1,π‘šβˆ’π‘ˆπ‘—π‘¦π‘¦π‘™βˆ’1,π‘šξ‚+β„Ž28𝐴𝑗𝑙,π‘šξ‚€π΄π‘—π‘₯𝑙,π‘šπ‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+𝐡𝑗π‘₯𝑙,π‘šπ‘ˆπ‘—π‘¦π‘¦π‘™,π‘šξ‚,(2.6a)π‘ˆπ‘—π‘¦π‘™,π‘š=π‘ˆπ‘—π‘¦π‘™,π‘š+β„Ž16𝐡𝑗𝑙,π‘šξ‚€πΊπ‘—π‘™,π‘š+1βˆ’πΊπ‘—π‘™,π‘šβˆ’1ξ‚βˆ’β„Ž16𝐡𝑗𝑙,π‘šξ‚€π‘ˆπ‘—π‘‘π‘‘π‘™,π‘š+1βˆ’π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šβˆ’1+β„Žπ΄π‘—π‘™,π‘š16𝐡𝑗𝑙,π‘šξ‚€π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+1βˆ’π‘ˆπ‘—π‘₯π‘₯𝑙,π‘šβˆ’1+β„Ž28𝐡𝑗𝑙,π‘šξ‚€π΄π‘—π‘¦π‘™,π‘šπ‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+𝐡𝑗𝑦𝑙,π‘šπ‘ˆπ‘—π‘¦π‘¦π‘™,π‘šξ‚,(2.6b) 𝐺𝑗𝑙,π‘šξ‚΅π‘₯=𝑔𝑙,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™,π‘š,π‘ˆπ‘—π‘₯𝑙,π‘š,π‘ˆπ‘—π‘¦π‘™,π‘š,π‘ˆπ‘—π‘‘π‘™,π‘šξ‚Ά.(2.7)

Then at each grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), 𝑙,π‘š=1,2,…,𝑁, 𝑗=2,…,𝐽, an approximation with accuracy of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) for the solution of differential equation (2.1) may be written as𝐿𝑒≑𝑝2𝐴𝑗𝑙,π‘šβˆ’β„Ž26𝐴𝑗π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐴𝑗π‘₯𝑙,π‘šβˆ’β„Ž26𝐡𝑗𝑦𝑙,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐴𝑗𝑦𝑙,π‘š+β„Ž2𝐴12𝑗π‘₯π‘₯𝑙,π‘š+β„Ž2𝐴12𝑗𝑦𝑦𝑙,π‘šξƒ­π›Ώ2π‘₯π‘ˆπ‘—π‘™,π‘š+𝑝2𝐡𝑗𝑙,π‘šβˆ’β„Ž26𝐴𝑗π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐡𝑗π‘₯𝑙,π‘šβˆ’β„Ž26𝐡𝑗𝑦𝑙,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐡𝑗𝑦𝑙,π‘š+β„Ž2𝐡12𝑗π‘₯π‘₯𝑙,π‘š+β„Ž2𝐡12𝑗𝑦𝑦𝑙,π‘šξƒ­π›Ώ2π‘¦π‘ˆπ‘—π‘™,π‘š+β„Žπ‘2𝐴12𝑗𝑦𝑙,π‘šβˆ’ξƒ©π΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐴𝑗𝑙,π‘šξƒ­ξ€·π›Ώ2π‘₯2πœ‡π‘¦π›Ώπ‘¦ξ€Έπ‘ˆπ‘—π‘™,π‘š+β„Žπ‘2𝐡12𝑗π‘₯𝑙,π‘šβˆ’ξƒ©π΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐡𝑗𝑙,π‘šξƒ­ξ€·π›Ώ2𝑦2πœ‡π‘₯𝛿π‘₯ξ€Έπ‘ˆπ‘—π‘™,π‘š+𝑝2𝐴12𝑗𝑙,π‘š+𝐡𝑗𝑙,π‘šξ‚„π›Ώ2π‘₯𝛿2π‘¦π‘ˆπ‘—π‘™,π‘š=π‘˜2121βˆ’β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™+1,π‘š+1+β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1,π‘š+1βˆ’β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™,π‘š+1+1+β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™,π‘šβˆ’1+8π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šξƒ­βˆ’π‘˜2121βˆ’β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐺𝑗𝑙+1,π‘š+1+β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺπΊπ‘—π‘™βˆ’1,π‘š+1βˆ’β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐺𝑗𝑙,π‘š+1+1+β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐺𝑗𝑙,π‘šβˆ’1+8𝐺𝑗𝑙,π‘šξƒ­+𝑇𝑗𝑙,π‘š,𝑙=1,2,…,𝑁;π‘š=1,2,…,𝑁;𝑗=1,2,…,𝐽,(2.8) where 𝛿π‘₯π‘ˆπ‘—π‘™,π‘š=(π‘ˆπ‘—π‘™+12,π‘šβˆ’π‘ˆπ‘—π‘™βˆ’12,π‘š), πœ‡π‘₯π‘ˆπ‘—π‘™,π‘š=(1/2)(π‘ˆπ‘—π‘™+1/2,π‘š+π‘ˆπ‘—π‘™βˆ’(1/2),π‘š),… and so forth and 𝑇𝑗𝑙,π‘š=𝑂(π‘˜4+π‘˜4β„Ž2+π‘˜2β„Ž4).

3. Derivation Procedure of the Approximation

At the grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), let us denote 𝛼𝑗𝑙,π‘š=ξ‚΅πœ•π‘”πœ•π‘ˆπ‘₯𝑗𝑙,π‘š,𝛽𝑗𝑙,π‘š=ξ‚΅πœ•π‘”πœ•π‘ˆπ‘¦ξ‚Άπ‘—π‘™,π‘š.(3.1)

Further at the grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), the exact solution value π‘ˆπ‘—π‘™,π‘š satisfies π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šξ€·π‘₯βˆ’π΄π‘™,π‘¦π‘š,π‘‘π‘—ξ€Έπ‘ˆπ‘—π‘₯π‘₯𝑙,π‘šξ€·π‘₯βˆ’π΅π‘™,π‘¦π‘š,π‘‘π‘—ξ€Έπ‘ˆπ‘—π‘¦π‘¦π‘™,π‘šξ‚€π‘₯=𝑔𝑙,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™,π‘š,π‘ˆπ‘—π‘₯𝑙,π‘š,π‘ˆπ‘—π‘¦π‘™,π‘š,π‘ˆπ‘—π‘‘π‘™,π‘šξ‚β‰‘πΊπ‘—π‘™,π‘š(say).(3.2) Using Taylor expansion about the grid point (π‘₯𝑙,π‘¦π‘š,𝑑𝑗), first we obtain𝐿𝑒=π‘˜2121βˆ’β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™+1,π‘š+1+β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1,π‘š+1βˆ’β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™,π‘š+1+1+β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺπ‘ˆπ‘—π‘‘π‘‘π‘™,π‘šβˆ’1+8π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šξƒ­βˆ’π‘˜2121βˆ’β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐺𝑗𝑙+1,π‘š+1+β„Žπ΄π‘—π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺπΊπ‘—π‘™βˆ’1,π‘š+1βˆ’β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐺𝑗𝑙,π‘š+1+1+β„Žπ΅π‘—π‘¦π‘™,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐺𝑗𝑙,π‘šβˆ’1+8𝐺𝑗𝑙,π‘šξƒ­ξ€·π‘˜+𝑂4+π‘˜4β„Ž2+π‘˜2β„Ž4ξ€Έ.(3.3)

With the help of the approximations (2.2a)–(2.2e), and from (2.5a) and (2.5b), we obtain𝐺𝑗𝑙±1,π‘š=𝐺𝑗𝑙±1,π‘šβˆ’β„Ž23πœ•3π‘ˆπ‘—π‘™,π‘šπœ•π‘₯3𝛼𝑗𝑙,π‘š+β„Ž26πœ•3π‘ˆπ‘—π‘™,π‘šπœ•π‘¦3𝛽𝑗𝑙,π‘šξ€·π‘˜+𝑂2Β±π‘˜2β„ŽΒ±β„Ž3+β„Ž4ξ€Έ,𝐺𝑗𝑙,π‘šΒ±1=𝐺𝑗𝑙,π‘šΒ±1+β„Ž26πœ•3π‘ˆπ‘—π‘™,π‘šπœ•π‘₯3𝛼𝑗𝑙,π‘šβˆ’β„Ž23πœ•3π‘ˆπ‘—π‘™,π‘šπœ•π‘¦3𝛽𝑗𝑙,π‘šξ€·π‘˜+𝑂2Β±π‘˜2β„ŽΒ±β„Ž3+β„Ž4ξ€Έ.(3.4) Now we define the following approximations: π‘ˆπ‘—π‘₯𝑙,π‘š=π‘ˆπ‘—π‘₯𝑙,π‘š+π‘Ž11β„Žξ‚€πΊπ‘—π‘™+1,π‘šβˆ’πΊπ‘—π‘™βˆ’1,π‘šξ‚+π‘Ž12β„Žξ‚€π‘ˆπ‘—π‘‘π‘‘π‘™+1,π‘šβˆ’π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1,π‘šξ‚+π‘Ž13β„Žξ‚€π‘ˆπ‘—π‘¦π‘¦π‘™+1,π‘šβˆ’π‘ˆπ‘—π‘¦π‘¦π‘™βˆ’1,π‘šξ‚+β„Ž2ξ‚€π‘Ž14π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+π‘Ž15π‘ˆπ‘—π‘¦π‘¦π‘™,π‘šξ‚,π‘ˆπ‘—π‘¦π‘™,π‘š=π‘ˆπ‘—π‘¦π‘™,π‘š+𝑏11β„Žξ‚€πΊπ‘—π‘™,π‘š+1βˆ’πΊπ‘—π‘™,π‘šβˆ’1+𝑏12β„Žξ‚€π‘ˆπ‘—π‘‘π‘‘π‘™,π‘š+1βˆ’π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šβˆ’1+𝑏13β„Žξ‚€π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+1βˆ’π‘ˆπ‘—π‘₯π‘₯𝑙,π‘šβˆ’1+β„Ž2𝑏14π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+𝑏15π‘ˆπ‘—π‘¦π‘¦π‘™,π‘šξ‚,(3.5) where π‘Ž1𝑖 and 𝑏1𝑖, 𝑖=1,2,…,5, are free parameters to be determined.

By using the approximations (2.2a)–(2.4e) and (3.4), from (3.5), we getπ‘ˆπ‘—π‘₯𝑙,π‘š=π‘ˆπ‘—π‘₯𝑙,π‘š+β„Ž26𝑇1ξ€·π‘˜+𝑂2+π‘˜2β„Ž2+β„Ž4ξ€Έ,π‘ˆπ‘—π‘¦π‘™,π‘š=π‘ˆπ‘—π‘¦π‘™,π‘š+β„Ž26𝑇2ξ€·π‘˜+𝑂2+π‘˜2β„Ž2+β„Ž4ξ€Έ,(3.6) where 𝑇1=ξ‚€1βˆ’12π‘Ž11𝐴𝑗𝑙,π‘šξ‚π‘ˆπ‘—π‘₯π‘₯π‘₯𝑙,π‘šξ€·π‘Ž+1211+π‘Ž12ξ€Έπ‘ˆπ‘—π‘‘π‘‘π‘₯𝑙,π‘šξ‚€+12βˆ’π‘Ž11𝐡𝑗𝑙,π‘š+π‘Ž13ξ‚π‘ˆπ‘—π‘₯𝑦𝑦𝑙,π‘š+ξ‚€6π‘Ž14βˆ’12π‘Ž11𝐴𝑗π‘₯𝑙,π‘šξ‚π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+ξ‚€6π‘Ž15βˆ’12π‘Ž11𝐡𝑗π‘₯𝑙,π‘šξ‚π‘ˆπ‘—π‘¦π‘¦π‘™,π‘š,𝑇2=ξ‚€1βˆ’12𝑏11𝐡𝑗𝑙,π‘šξ‚π‘ˆπ‘—π‘¦π‘¦π‘¦π‘™,π‘šξ€·π‘+1211+𝑏12ξ€Έπ‘ˆπ‘—π‘‘π‘‘π‘¦π‘™,π‘šξ‚€+12βˆ’π‘11𝐴𝑗𝑙,π‘š+𝑏13ξ‚π‘ˆπ‘—π‘₯π‘₯𝑦𝑙,π‘š+ξ‚€6𝑏14βˆ’12𝑏11𝐴𝑗𝑦𝑙,π‘šξ‚π‘ˆπ‘—π‘₯π‘₯𝑙,π‘š+ξ‚€6𝑏15βˆ’12𝑏11𝐡𝑗𝑦𝑙,π‘šξ‚π‘ˆπ‘—π‘¦π‘¦π‘™,π‘š.(3.7) Thus from (2.7), we obtain 𝐺𝑗𝑙,π‘š=𝐺𝑗𝑙,π‘š+β„Ž26𝑇1𝛼𝑗𝑙,π‘š+β„Ž26𝑇2𝛽𝑗𝑙,π‘šξ€·π‘˜+𝑂2+π‘˜2β„Ž2+β„Ž4ξ€Έ.(3.8) Finally by the help of the approximations (3.4) and (3.8), from (2.8) and (3.3), we get 𝑇𝑗𝑙,π‘š=βˆ’π‘˜2β„Ž2ξƒ©πœ•363π‘ˆπ‘—π‘™,π‘šπœ•π‘₯3𝛼𝑗𝑙,π‘š+πœ•3π‘ˆπ‘—π‘™,π‘šπœ•π‘¦3𝛽𝑗𝑙,π‘šξ‚€π‘‡βˆ’41𝛼𝑗𝑙,π‘š+𝑇2𝛽𝑗𝑙,π‘šξ‚ξƒͺξ€·π‘˜+𝑂4+π‘˜4β„Ž2+π‘˜2β„Ž4ξ€Έ.(3.9)

In order to obtain the difference approximation of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) the coefficient of π‘˜2β„Ž2 in (3.9) must be zero, which gives the values of parameters:π‘Ž11=116𝐴𝑗𝑙,π‘š,π‘Ž12=βˆ’116𝐴𝑗𝑙,π‘š,π‘Ž13=𝐡𝑗𝑙,π‘š16𝐴𝑗𝑙,π‘š,π‘Ž14=𝐴𝑗π‘₯𝑙,π‘š8𝐴𝑗𝑙,π‘š,π‘Ž15=𝐡𝑗π‘₯𝑙,π‘š8𝐴𝑗𝑙,π‘š,𝑏11=116𝐡𝑗𝑙,π‘š,𝑏12=βˆ’116𝐡𝑗𝑙,π‘š,𝑏13=𝐴𝑗𝑙,π‘š16𝐡𝑗𝑙,π‘š,𝑏14=𝐴𝑗𝑦𝑙,π‘š8𝐡𝑗𝑙,π‘š,𝑏15=𝐡𝑗𝑦𝑙,π‘š8𝐡𝑗𝑙,π‘š.(3.10)

Thus we obtain the difference method of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) for the differential equation (2.1) and the local truncation error reduces to 𝑇𝑗𝑙,π‘š=𝑂(π‘˜4+π‘˜4β„Ž2+π‘˜2β„Ž4).

Now we consider the numerical method of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) for the solution of the quasilinear hyperbolic equation (1.1). Whenever the coefficients are 𝐴=𝐴(π‘₯,𝑦,𝑑,𝑒) and 𝐡=𝐡(π‘₯,𝑦,𝑑,𝑒), the difference scheme (2.8) needs to be modified. For this purpose, we use the following central differences: 𝐴𝑗π‘₯𝑙,π‘š=𝐴𝑗𝑙+1,π‘šβˆ’π΄π‘—π‘™βˆ’1,π‘šξ‚ξ€·β„Ž(2β„Ž)+𝑂2ξ€Έ,𝐴𝑗π‘₯π‘₯𝑙,π‘š=𝐴𝑗𝑙+1,π‘šβˆ’2𝐴𝑗𝑙,π‘š+π΄π‘—π‘™βˆ’1,π‘šξ‚ξ€·β„Ž2ξ€Έξ€·β„Ž+𝑂2ξ€Έ,𝐴𝑗𝑦𝑙,π‘š=𝐴𝑗𝑙,π‘š+1βˆ’π΄π‘—π‘™,π‘šβˆ’1ξ‚ξ€·β„Ž(2β„Ž)+𝑂2ξ€Έ,𝐴𝑗𝑦𝑦𝑙,π‘š=𝐴𝑗𝑙,π‘š+1βˆ’2𝐴𝑗𝑙,π‘š+𝐴𝑗𝑙,π‘šβˆ’1ξ‚ξ€·β„Ž2ξ€Έξ€·β„Ž+𝑂2ξ€Έ,𝐡𝑗π‘₯𝑙,π‘š=𝐡𝑗𝑙+1,π‘šβˆ’π΅π‘—π‘™βˆ’1,π‘šξ‚ξ€·β„Ž(2β„Ž)+𝑂2ξ€Έ,𝐡𝑗π‘₯π‘₯𝑙,π‘š=𝐡𝑗𝑙+1,π‘šβˆ’2𝐡𝑗𝑙,π‘š+π΅π‘—π‘™βˆ’1,π‘šξ‚ξ€·β„Ž2ξ€Έξ€·β„Ž+𝑂2ξ€Έ,𝐡𝑗𝑦𝑙,π‘š=𝐡𝑗𝑙,π‘š+1βˆ’π΅π‘—π‘™,π‘šβˆ’1ξ‚ξ€·β„Ž(2β„Ž)+𝑂2ξ€Έ,𝐡𝑗𝑦𝑦𝑙,π‘š=𝐡𝑗𝑙,π‘š+1βˆ’2𝐡𝑗𝑙,π‘š+𝐡𝑗𝑙,π‘šβˆ’1ξ‚ξ€·β„Ž2ξ€Έξ€·β„Ž+𝑂2ξ€Έ,(3.11) where𝐴𝑗𝑙,π‘šξ‚€π‘₯=𝐴𝑙,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™,π‘šξ‚,𝐴𝑗𝑙±1,π‘šξ‚€π‘₯=𝐴𝑙±1,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™Β±1,π‘šξ‚,𝐴𝑗𝑙,π‘šΒ±1ξ‚€π‘₯=𝐴𝑙,π‘¦π‘šΒ±1,𝑑𝑗,π‘ˆπ‘—π‘™,π‘šΒ±1,𝐡𝑗𝑙,π‘šξ‚€π‘₯=𝐡𝑙,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™,π‘šξ‚,𝐡𝑗𝑙±1,π‘šξ‚€π‘₯=𝐡𝑙±1,π‘¦π‘š,𝑑𝑗,π‘ˆπ‘—π‘™Β±1,π‘šξ‚,𝐡𝑗𝑙,π‘šΒ±1ξ‚€π‘₯=𝐡𝑙,π‘¦π‘šΒ±1,𝑑𝑗,π‘ˆπ‘—π‘™,π‘šΒ±1.(3.12)

By the help of the approximations (3.11), it is easy to verify that 𝐴𝑗𝑙,π‘šβˆ’β„Ž26𝐴𝑗π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐴𝑗π‘₯𝑙,π‘šβˆ’β„Ž26𝐡𝑗𝑦𝑙,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐴𝑗𝑦𝑙,π‘š+β„Ž2𝐴12𝑗π‘₯π‘₯𝑙,π‘š+β„Ž2𝐴12𝑗𝑦𝑦𝑙,π‘š=𝐴𝑗𝑙,π‘šβˆ’124𝐴𝑗𝑙,π‘šξ‚€π΄π‘—π‘™+1,π‘šβˆ’π΄π‘—π‘™βˆ’1,π‘šξ‚2βˆ’124𝐡𝑗𝑙,π‘šξ‚€π΅π‘—π‘™,π‘š+1βˆ’π΅π‘—π‘™,π‘šβˆ’1𝐴𝑗𝑙,π‘š+1βˆ’π΄π‘—π‘™,π‘šβˆ’1+1𝐴12𝑗𝑙+1,π‘šβˆ’2𝐴𝑗𝑙,π‘š+π΄π‘—π‘™βˆ’1,π‘šξ‚+1𝐴12𝑗𝑙,π‘š+1βˆ’2𝐴𝑗𝑙,π‘š+𝐴𝑗𝑙,π‘šβˆ’1ξ‚ξ€·β„Ž+𝑂4ξ€Έ,𝐡𝑗𝑙,π‘šβˆ’β„Ž26𝐴𝑗π‘₯𝑙,π‘šπ΄π‘—π‘™,π‘šξƒͺ𝐡𝑗π‘₯𝑙,π‘šβˆ’β„Ž26𝐡𝑗𝑦𝑙,π‘šπ΅π‘—π‘™,π‘šξƒͺ𝐡𝑗𝑦𝑙,π‘š+β„Ž2𝐡12𝑗π‘₯π‘₯𝑙,π‘š+β„Ž2𝐡12𝑗𝑦𝑦𝑙,π‘š=𝐡𝑗𝑙,π‘šβˆ’124𝐡𝑗𝑙,π‘šξ‚€π΅π‘—π‘™,π‘š+1βˆ’π΅π‘—π‘™,π‘šβˆ’12βˆ’124𝐴𝑗𝑙,π‘šξ‚€π΄π‘—π‘™+1,π‘šβˆ’π΄π‘—π‘™βˆ’1,π‘šπ΅ξ‚ξ‚€π‘—π‘™+1,π‘šβˆ’π΅π‘—π‘™βˆ’1,π‘šξ‚+1𝐡12𝑗𝑙+1,π‘šβˆ’2𝐡𝑗𝑙,π‘š+π΅π‘—π‘™βˆ’1,π‘šξ‚+1𝐡12𝑗𝑙,π‘š+1βˆ’2𝐡𝑗𝑙,π‘š+𝐡𝑗𝑙,π‘šβˆ’1ξ‚ξ€·β„Ž+𝑂4ξ€Έ,…,andsoforth.(3.13)

Thus substituting the central difference approximations (3.11) into (2.8), we obtain the required numerical method of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) for the solution of quasilinear hyperbolic partial differential equation (1.1).

Note that the initial and Dirichlet boundary conditions are given by (1.2), (1.3a) and (1.3b), respectively. Incorporating the initial and boundary conditions, we can write the method (2.8) in a block tridiagonal matrix form. If the differential equation (1.1) is linear, we can solve the linear system using operator splitting method or, alternating direction implicit method (see [18, 19]); in the nonlinear case, we can use Newton-Raphson iterative method to solve the non-linear system (see [20–24]).

4. Stability Analysis

In this section, we aim to discuss a stable difference scheme for the two-space dimensional linear hyperbolic equation with singular coefficients and ensure that the numerical method developed here retains its order and accuracy.

Let us consider the equation of the following form:𝑒𝑑𝑑=𝑒π‘₯π‘₯+𝑒𝑦𝑦+𝐷(π‘₯)𝑒π‘₯+𝑓(π‘₯,𝑦,𝑑),0<π‘₯,𝑦<1,𝑑>0.(4.1) Applying the approximation (2.8) to the differential equation (4.1), we obtain𝑝2𝛿2π‘₯+𝛿2𝑦+16𝛿2π‘₯𝛿2π‘¦ξ‚„π‘ˆπ‘—π‘™,π‘š=π‘˜2ξ‚Έ12π‘ˆπ‘—π‘‘π‘‘π‘™+1,π‘š+π‘ˆπ‘—π‘‘π‘‘π‘™βˆ’1,π‘š+π‘ˆπ‘—π‘‘π‘‘π‘™,π‘š+1+π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šβˆ’1+8π‘ˆπ‘—π‘‘π‘‘π‘™,π‘šξ‚Ήβˆ’π‘˜2ξ‚Έ12𝐺𝑗𝑙+1,π‘š+πΊπ‘—π‘™βˆ’1,π‘š+𝐺𝑗𝑙,π‘š+1+𝐺𝑗𝑙,π‘šβˆ’1+8𝐺𝑗𝑙,π‘šξ‚Ήξ€·π‘˜+𝑂4+π‘˜4β„Ž2+π‘˜2β„Ž4𝑙=1,2,…,𝑁;π‘š=1,2,…,𝑁;𝑗=2,3,…,𝐽,(4.2) where we denote𝐷𝑙π‘₯=𝐷𝑙,𝑓𝑗𝑙,π‘šξ€·π‘₯=𝑓𝑙,π‘¦π‘š,𝑑𝑗,𝐺𝑗𝑙±1,π‘š=𝐷𝑙±1π‘ˆπ‘—π‘₯𝑙±1,π‘š+𝑓𝑗𝑙±1,π‘š,𝐺𝑗𝑙,π‘šΒ±1=π·π‘™π‘ˆπ‘—π‘₯𝑙,π‘šΒ±1+𝑓𝑗𝑙,π‘šΒ±1,𝐺𝑗𝑙,π‘š=π·π‘™π‘ˆπ‘—π‘₯𝑙,π‘š+𝑓𝑗𝑙,π‘š.(4.3)

Note that the linear scheme (4.2) is of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) for the solution of differential equation (4.1). However, the scheme (4.2) fails to compute at 𝑙=1, if the coefficient 𝐷(π‘₯) and/or the function 𝑓(π‘₯,𝑦,𝑑) contains the singular terms like 1/π‘₯,1/π‘₯2,…, and so forth. We overcome this difficulty by modifying the method (4.2) in such a way that the solution retains its order and accuracy everywhere in the vicinity of singularity π‘₯=0.

We need the following approximations: 𝐷𝑙±1=π·π‘™Β±β„Žπ·π‘₯𝑙+β„Ž22𝐷𝑙π‘₯π‘₯ξ€·β„ŽΒ±π‘‚3ξ€Έ,𝑓𝑗𝑙±1,π‘š=𝑓𝑗𝑙,π‘šΒ±β„Žπ‘“π‘—π‘₯𝑙,π‘š+β„Ž22𝑓𝑗π‘₯π‘₯𝑙,π‘šξ€·β„ŽΒ±π‘‚3ξ€Έ,𝑓𝑗𝑙,π‘šΒ±1=𝑓𝑗𝑙,π‘šΒ±β„Žπ‘“π‘—π‘¦π‘™,π‘š+β„Ž22𝑓𝑗𝑦𝑦𝑙,π‘šξ€·β„ŽΒ±π‘‚3ξ€Έ.(4.4)

Using the approximations (4.4) in the scheme (4.2) and neglecting local truncation error, we get11+𝛿122π‘₯+𝑅0ξ€·2πœ‡π‘₯𝛿π‘₯+𝛿2𝑦𝛿122𝑑𝑒𝑗𝑙,π‘š=𝑝2𝑅1𝛿2π‘₯+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯ξ€Έ+𝛿2𝑦+𝑅06𝛿2𝑦2πœ‡π‘₯𝛿π‘₯ξ€Έ+16𝛿2π‘₯𝛿2𝑦𝑒𝑗𝑙,π‘š+𝑓,(4.5) where𝑅0=β„Žπ·π‘™2,𝑅1β„Ž=1+2ξ‚€122𝐷π‘₯𝑙+𝐷𝑙2,𝑅2=𝑅0+β„Ž3𝐷24π‘₯π‘₯𝑙+𝐷𝑙𝐷π‘₯𝑙,ξ“π‘˜π‘“=2ξ‚€1212𝑓𝑗𝑙,π‘š+β„Ž2𝑓𝑗π‘₯π‘₯𝑙,π‘š+𝑓𝑗𝑦𝑦𝑙,π‘š+𝐷𝑙𝑓𝑗π‘₯𝑙,π‘š.(4.6) The previous scheme in product form may be written as11+𝛿122π‘₯+𝑅0ξ€·2πœ‡π‘₯𝛿π‘₯𝛿1+2𝑦𝛿122𝑑𝑒𝑗𝑙,π‘š=𝑝2𝑅1𝛿2π‘₯+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯ξ€Έ+𝛿2𝑦+𝑅06𝛿2𝑦2πœ‡π‘₯𝛿π‘₯ξ€Έ+16𝛿2π‘₯𝛿2𝑦𝑒𝑗𝑙,π‘š+𝑓.(4.7)

It is difficult to find the stability region for the scheme (4.7). In order to obtain a valid stability region, we may modify the scheme (4.7) (see [23]) as11+𝑅121𝛿2π‘₯+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯𝛿1+2𝑦𝛿122𝑑𝑒𝑗𝑙,π‘š=𝑝2𝑅1𝛿2π‘₯+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯ξ€Έ+𝛿2𝑦+16𝑅1𝛿2π‘₯+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯𝛿2𝑦𝑒𝑗𝑙,π‘š+𝑓.(4.8)

The additional terms added in (4.7) and (4.8) are of high orders and do not affect the accuracy of the scheme.

For stability, we put 𝑒𝑗𝑙,π‘š=π΄π‘™π΅π‘šπœ‰π‘—π‘’π‘–π›½π‘™π‘’π‘–π›Ύπ‘š (where πœ‰=π‘’π‘–πœ™ such that |πœ‰|=1) in the homogeneous part of the scheme (4.8), and we getπœ‰βˆ’2+πœ‰βˆ’1=βˆ’4sin2πœ™2=𝑝2𝑀1+𝑀2+(1/6)𝑀1𝑀2ξ€»ξ€Ί1+(1/12)𝑀1ξ€»ξ€Ί1+(1/12)𝑀2ξ€»,(4.9) where𝑀1=𝑅1𝐴+π΄βˆ’1ξ€Έξ€·cosπ›½βˆ’2+π‘–π΄βˆ’π΄βˆ’1ξ€Έξ€Ύsin𝛽+𝑅2ξ€½ξ€·π΄βˆ’π΄βˆ’1ξ€Έξ€·cos𝛽+𝑖𝐴+π΄βˆ’1ξ€Έξ€Ύ,𝑀sin𝛽2=𝐡+π΅βˆ’1ξ€Έξ€·cosπ›Ύβˆ’2+π‘–π΅βˆ’π΅βˆ’1ξ€Έsin𝛾,(4.10) where 𝐴 and 𝐡 are nonzero real parameters to be determined. Left-hand side of (4.9) is a real quantity. Thus the imaginary part of right-hand side of (4.9) must be zero. Thus we obtain 𝑅1+𝑅2ξ€Έξ€·π‘…π΄βˆ’1βˆ’π‘…2ξ€Έπ΄βˆ’1=0,π΅βˆ’π΅βˆ’1=0.(4.11) On solving, we get √𝐴=(𝑅1βˆ’π‘…2)/(𝑅1+𝑅2), π΄βˆ’1=√(𝑅1+𝑅2)/(𝑅1βˆ’π‘…2), and 𝐡=π΅βˆ’1=1, provided that 𝑅1±𝑅2>0,

Substituting the values of 𝐴, π΄βˆ’1, 𝐡 and π΅βˆ’1 in (4.10), we obtain𝑀1𝑅=βˆ’21+𝑅21βˆ’π‘…22ξ‚΅2sin2𝛽2βˆ’1ξ‚Άξ‚Ή,𝑀2=βˆ’4sin2𝛾2.(4.12)

Let 𝑀1=βˆ’2𝑀3 and 𝑀2=βˆ’4𝑀4, where 𝑀3=𝑅1+𝑅21βˆ’π‘…22ξ‚΅2sin2𝛽2ξ‚Άβˆ’1,𝑀4=sin2𝛾2.(4.13)

Substituting the values of 𝑀1 and 𝑀2 in (4.9), we obtain sin2πœ™2=𝑝2𝑀3+2𝑀4βˆ’(2/3)𝑀3𝑀4ξ€»2ξ€Ί1βˆ’(1/6)𝑀3ξ€»ξ€Ί1βˆ’(1/3)𝑀4ξ€».(4.14)

Since sin2(πœ™/2)≀1, hence the method (4.8) is stable as long as𝑝max2𝑀3+2𝑀4βˆ’23𝑀3𝑀4ξ‚€21≀min1βˆ’6𝑀311βˆ’3𝑀4.(4.15) this implies 0<𝑝2≀𝑅4βˆ’2+𝑅22βˆ’π‘…32ξ‚Ά3ξ‚΅2+𝑅2+𝑅22βˆ’π‘…32ξ‚Ά,(4.16) which is the required stability interval for the scheme (4.8).

In order to facilitate the computation, we may rewrite (4.8) in operator split form (see [18, 19]) as 𝛿1+𝑦2𝑒12βˆ—π‘™,π‘š=𝑝2𝑅1𝛿π‘₯2+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯ξ€Έ+𝛿𝑦2+16𝑅1𝛿π‘₯2+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯𝛿𝑦2𝑒𝑗𝑙,π‘š+1𝑓,1+𝑅121𝛿π‘₯2+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯𝛿𝑑2𝑒𝑗𝑙,π‘š=π‘’βˆ—π‘™,π‘š(4.17) or,11+𝑅121𝛿π‘₯2+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯𝑒𝑗+1𝑙,π‘š=π‘’βˆ—π‘™,π‘š+11+𝑅121𝛿2π‘₯+𝑅2ξ€·2πœ‡π‘₯𝛿π‘₯ξ€Έξ€Έξ‚„ξ‚€2𝑒𝑗𝑙,π‘šβˆ’π‘’π‘—βˆ’1𝑙,π‘šξ‚,(4.18) where π‘’βˆ—π‘™,π‘š is an intermediate value. The intermediate boundary values required for solving (4.17) can be obtained from (4.18).

The left-hand sides of (4.17) and (4.18) are factorizations into 𝑦- and π‘₯-differences, respectively, which allows us to solve by sweeping first (4.17) in the 𝑦- and then (4.18) in the π‘₯-direction. It will be seen that these sweeps require only the solution of tri-diagonal systems.

5. Super Stable Method for Two-Dimensional Linear Hyperbolic Equation

Let us consider the equation of the following form:𝑒𝑑𝑑+2𝛼𝑒𝑑+𝛽2𝑒=𝑒π‘₯π‘₯+𝑒𝑦𝑦+𝑓(π‘₯,𝑦,𝑑),0<π‘₯,𝑦<1,𝑑>0,𝛼>0,𝛽β‰₯0(5.1) subject to appropriate initial and Dirichlet boundary conditions that are prescribed. The equation above represents two-space dimensional linear telegraphic hyperbolic equation. Applying the difference method (2.8) to the differential equation (5.1), we obtain𝛿2𝑑𝑒𝑗𝑙,π‘š+βˆšπ‘Žξ€·2πœ‡π‘‘π›Ώπ‘‘ξ€Έπ‘’π‘—π‘™,π‘š+βˆšπ‘Žξ€·π›Ώ122π‘₯+𝛿2𝑦2πœ‡π‘‘π›Ώπ‘‘ξ€Έπ‘’π‘—π‘™,π‘š+𝑏𝑒𝑗𝑙,π‘š+𝑏12βˆ’π‘2𝛿2π‘₯+𝛿2𝑦𝑒𝑗𝑙,π‘šβˆ’π‘26𝛿2π‘₯𝛿2π‘¦π‘ˆπ‘—π‘™,π‘š+1𝛿122π‘₯+𝛿2𝑦𝛿2π‘‘π‘ˆπ‘—π‘™,π‘š=𝐹𝐹,(5.2) whereπ‘Ž=𝛼2π‘˜2,𝑏=𝛽2π‘˜2π‘˜,𝐹𝐹=2𝑓12𝑗𝑙+1,π‘š+π‘“π‘—π‘™βˆ’1,π‘š+𝑓𝑗𝑙,π‘š+1+𝑓𝑗𝑙,π‘šβˆ’1+8𝑓𝑗𝑙,π‘šξ‚.(5.3)

The previous scheme is conditionally stable (see [7]). In order to obtain a superstable method, we simply follow the ideas given by Chawla [24] and the scheme (5.2) may be re-written as1+πœ‚π‘2+ξ‚€112βˆ’π›Ύπ‘2𝛿2π‘₯+𝛿2𝑦𝛿2𝑑𝑒𝑗𝑙,π‘š+βˆšπ‘Žξƒ©π›Ώ1+2π‘₯+𝛿122𝑦ξƒͺξ€·122πœ‡π‘‘π›Ώπ‘‘ξ€Έπ‘’π‘—π‘™,π‘š=𝑝2βˆ’π‘ξ‚ξ€·π›Ώ122π‘₯+𝛿2𝑦+𝑝26𝛿2π‘₯𝛿2π‘¦ξ‚Ήπ‘’βˆ’π‘π‘—π‘™,π‘š+𝐹𝐹,(5.4) where πœ‚ and 𝛾 are free parameters to be determined. The additional terms are of high orders and do not affect the accuracy of the scheme.

To study stability for the scheme (5.4), we substitute 𝑒𝑗𝑙,π‘š=πœ‰π‘—π‘’π‘–πœƒπ‘™π‘’π‘–πœ™π‘š in the homogeneous part of the (5.4), and we obtain the characteristic equation:π΄πœ‰2+π΅πœ‰+𝐢=0,(5.5) where𝐴=1+πœ‚π‘2+βˆšβˆšπ‘Žβˆ’π‘Ž3ξ‚΅sin2πœƒ2+sin2πœ™2ξ‚Ά+ξ‚€4𝛾𝑝2βˆ’13sin2πœƒ2+sin2πœ™2ξ‚Ά,𝐡=βˆ’2βˆ’2πœ‚π‘2𝑝+42βˆ’π‘ξ‚ξ‚΅12sin2πœƒ2+sin2πœ™2ξ‚Ά+ξ‚€2+𝑏3βˆ’8𝛾𝑝2sin2πœƒ2+sin2πœ™2ξ‚Άβˆ’8𝑝23sin2πœƒ2sin2πœ™2,𝐢=1+πœ‚π‘2βˆ’βˆšβˆšπ‘Ž+π‘Ž3ξ‚΅sin2πœƒ2+sin2πœ™2ξ‚Ά+ξ‚€4𝛾𝑝2βˆ’13sin2πœƒ2+sin2πœ™2ξ‚Ά.(5.6)

Using the transformation πœ‰=(1+𝑧)/(1βˆ’π‘§) in (5.5), we obtain the transformed characteristic equation:(π΄βˆ’π΅+𝐢)𝑧2+2(π΄βˆ’πΆ)𝑧+(𝐴+𝐡+𝐢)=0.(5.7) For stability of the scheme (5.4), we must have the following conditions: 𝐴+𝐡+𝐢>0,π΄βˆ’πΆ>0,π΄βˆ’π΅+𝐢>0.(5.8) It is easy to verify that the coefficient is as follows:𝑏𝐴+𝐡+𝐢=6ξ‚΅sin2πœƒ2+sin2πœ™2ξ‚Ά+𝑏2ξ‚΅cos2πœƒ2+cos2πœ™2ξ‚Ά+4𝑝2ξ‚΅sin2πœƒ2+sin2πœ™2βˆ’2sin2πœƒ2sin2πœ™2ξ‚Ά+16𝑝23sin2πœƒ2sin2πœ™2>0for𝛼>0,𝛽β‰₯0,andallπœƒ,πœ™,exceptπœƒ=πœ™=0or2πœ‹.(5.9) We can treat this case separately.

The coefficient is as follows:βˆšπ΄βˆ’πΆ=π‘Žξ‚΅13ξ‚΅sin2πœƒ2+sin2πœ™2ξ‚Ά+cos2πœƒ2+cos2πœ™2ξ‚Ά>0for𝛼>0,𝛽β‰₯0,andallπœƒ,πœ™.(5.10) For stability, it is required that the third coefficient is as follows:π΄βˆ’π΅+𝐢=4+4πœ‚π‘2π‘βˆ’π‘+3ξ‚΅sin2πœƒ2+sin2πœ™2ξ‚Άξ‚€+4(4π›Ύβˆ’1)𝑝2βˆ’13sin2πœƒ2+sin2πœ™2ξ‚Ά+8𝑝23sin2πœƒ2sin2πœ™2>0.(5.11) Multiplying throughout of (5.11) by 16πœ‚, we get (64πœ‚βˆ’1)+(8πœ‚π‘βˆ’1)2+16π‘πœ‚3ξ‚΅sin2πœƒ2+sin2πœ™2+64πœ‚(4π›Ύβˆ’1)𝑝2βˆ’13ξ‚„ξ‚΅sin2πœƒ2+sin2πœ™2ξ‚Ά+128𝑝23sin2πœƒ2sin2πœ™2>0.(5.12) Thus the scheme is stable if1πœ‚β‰₯64,𝛾β‰₯1+3𝑝212𝑝2.(5.13) For πœƒ=πœ™=0 or 2πœ‹ and 𝛽=0, we have the characteristic equation:ξ‚€βˆš1+π‘Žξ‚πœ‰2ξ‚€βˆšβˆ’2πœ‰+1βˆ’π‘Žξ‚=0(5.14) whose roots are πœ‰1,2=1, √(1βˆ’βˆšπ‘Ž)/(1+π‘Ž)=(1βˆ’π›Όπ‘˜)/(1+π›Όπ‘˜).

In this case also |πœ‰|≀1 and the scheme (5.4) is stable.

Hence for 𝛼>0, 𝛽β‰₯0, πœ‚β‰₯1/64, 𝛾β‰₯(1+3𝑝2)/12𝑝2 (which are independent of β„Ž and π‘˜), the scheme (5.4) is superstable.

Further, the scheme (5.4) in product form may be written as11+12βˆ’π›Ύπ‘2𝛿2𝑦1+πœ‚π‘2+ξ‚€112βˆ’π›Ύπ‘2𝛿2π‘₯𝛿2𝑑𝑒𝑗𝑙,π‘š+βˆšπ‘Žξƒ©π›Ώ1+2π‘₯ξƒͺξ€·122πœ‡π‘‘π›Ώπ‘‘ξ€Έπ‘’π‘—π‘™,π‘šξƒ°=𝑝2βˆ’π‘ξ‚ξ€·π›Ώ122π‘₯+𝛿2𝑦+𝑝26𝛿2π‘₯𝛿2π‘¦ξ‚Ήπ‘’βˆ’π‘π‘—π‘™,π‘š+𝐹𝐹.(5.15)

In this case also additional terms are of high orders, which do not affect the accuracy of the scheme. In order to facilitate the computation, we may rewrite the scheme (5.15) in two-step operator split form:11+12βˆ’π›Ύπ‘2𝛿2π‘¦ξ‚„π‘’βˆ—π‘™,π‘š=𝑝2βˆ’π‘ξ‚ξ€·π›Ώ122π‘₯+𝛿2𝑦+𝑝26𝛿2π‘₯𝛿2π‘¦ξ‚Ήπ‘’βˆ’π‘π‘—π‘™,π‘šξ‚ƒ+𝐹𝐹,(5.16a)1+πœ‚π‘2+ξ‚€112βˆ’π›Ύπ‘2𝛿2π‘₯𝛿2𝑑𝑒𝑗𝑙,π‘š+βˆšπ‘Žξƒ©π›Ώ1+2π‘₯ξƒͺξ€·122πœ‡π‘‘π›Ώπ‘‘ξ€Έπ‘’π‘—π‘™,π‘š=π‘’βˆ—π‘™,π‘š,(5.16b)where π‘’βˆ—π‘™,π‘š is any intermediate value, and the intermediate boundary conditions required for the solution of π‘’βˆ—π‘™,π‘š may be obtained from (5.16b). The left-hand side matrices represented by (5.16a) and (5.16b) are tridiagonal, thus very easily solved in the region 0<π‘₯, 𝑦<1, 𝑑>0 using a tridiagonal solver.

6. Numerical Illustrations

In this section, we have solved some benchmark problems using the method described by (2.8) and compared our results with the results of the fourth-order numerical method discussed in [6, 7]. The exact solutions are provided in each case. The initial and boundary conditions may be obtained using the exact solution as a test procedure. The linear difference equations have been solved using a tridiagonal solver, whereas non-linear difference equations have been solved using the Newton-Raphson method. While using the Newton-Raphson method, the iterations were stopped when absolute error tolerance ≀10βˆ’12 was achieved. In order to demonstrate the fourth-order convergence of the proposed method, for the computation of Example 6.1, we have chosen the fixed value of the parameter 𝜎=π‘˜/β„Ž2=3.2 for Example 6.1, and for other examples we have chosen the value of 𝜎=1.6. For this choice, our method behaves like a fourth-order method in space. All computations were carried out using double precision arithmetic.

Note that the proposed method (2.8) for second-order hyperbolic equations is a three-level scheme. The value of 𝑒 at 𝑑=0 is known from the initial condition. To start any computation, it is necessary to know the numerical value of 𝑒 of required accuracy at 𝑑=π‘˜. In this section, we discuss an explicit scheme of 𝑂(π‘˜2) for 𝑒 at first time level, that is, at 𝑑=π‘˜ in order to solve the differential equation (1.1) using the method (2.8), which is applicable to problems in both Cartesian and polar coordinates.

Since the values of 𝑒 and 𝑒𝑑 are known explicitly at 𝑑=0, this implies that all their successive tangential derivatives are known at 𝑑=0; that is, the values of 𝑒,𝑒π‘₯,𝑒π‘₯π‘₯,𝑒𝑦,𝑒𝑦𝑦,…,𝑒𝑑,𝑒𝑑π‘₯,𝑒𝑑𝑦,…, and so forth are known at 𝑑=0.

An approximation for 𝑒 of 𝑂(π‘˜2) at 𝑑=π‘˜ may be written as𝑒1𝑙,π‘š=𝑒0𝑙,π‘š+π‘˜π‘’0𝑑𝑙,π‘š+π‘˜22𝑒𝑑𝑑0𝑙,π‘šξ€·π‘˜+𝑂3ξ€Έ.(6.1)

From (1.1), we have𝑒𝑑𝑑0𝑙,π‘š=𝐴(π‘₯,𝑦,𝑑,𝑒)𝑒π‘₯π‘₯+𝐡(π‘₯,𝑦,𝑑,𝑒)𝑒𝑦𝑦+𝐺π‘₯,𝑦,𝑑,𝑒,𝑒π‘₯,𝑒𝑦,𝑒𝑑0𝑙,π‘š.(6.2)

Thus using the initial values and their successive tangential derivative values, from (6.2) we can obtain the value of (𝑒𝑑𝑑)0𝑙,π‘š, and then ultimately, from (6.1) we can compute the value of 𝑒 at first time level, that is, at 𝑑=π‘˜.

Example 6.1 (the telegraphic equation with forcing function). 𝑒𝑑𝑑+2𝛼𝑒𝑑+𝛽2𝑒=𝑒π‘₯π‘₯+𝑒𝑦𝑦+ξ€·2βˆ’4𝛼+𝛽2ξ€Έπ‘’βˆ’2𝑑sinhπ‘₯sinh𝑦,0<π‘₯,𝑦<1,𝑑>0,(6.3) where 𝛼>0 and 𝛽β‰₯0 are real parameters. The exact solution is 𝑒=π‘’βˆ’2𝑑sinhπ‘₯sinh𝑦. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 1 at 𝑑=5 for 𝛼>0 and 𝛽β‰₯0.

tab1
Table 1: Example 6.1: the maximum absolute errors (using proposed 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4)-method) (with CPU time).

Example 6.2 (wave equation in cylindrical polar coordinates). 𝑒𝑑𝑑=π‘’π‘Ÿπ‘Ÿ+𝑒𝑧𝑧+1π‘Ÿπ‘’π‘Ÿβˆ’ξ‚€3coshπ‘Ÿ+sinhπ‘Ÿπ‘Ÿξ‚cosh𝑧sin𝑑,0<π‘Ÿ,𝑧<1,𝑑>0.(6.4) The aforementioned equation represents the two-space dimensional wave equation in cylindrical polar coordinates. The exact solution is 𝑒=coshπ‘Ÿcosh𝑧sin𝑑. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 2 at 𝑑=1 and 𝑑=2.

tab2
Table 2: Example 6.2: the maximum absolute errors (with CPU time).

Example 6.3 (Van der Pol type nonlinear wave equation). 𝑒𝑑𝑑=𝑒π‘₯π‘₯+𝑒𝑦𝑦𝑒+𝛾2ξ€Έπ‘’βˆ’1𝑑+ξ€·2πœ‹2+𝛾2π‘’βˆ’2𝛾𝑑sin2(πœ‹π‘₯)sin2𝑒(πœ‹π‘¦)βˆ’π›Ύπ‘‘sin(πœ‹π‘₯)sin(πœ‹π‘¦),0<π‘₯,𝑦<1,𝑑>0,(6.5) with exact solution 𝑒=π‘’βˆ’π›Ύπ‘‘sin(πœ‹π‘₯)sin(πœ‹π‘¦). The maximum absolute errors and CPU time (in seconds) are tabulated in Table 3 at 𝑑=2 for 𝛾=1, 2, and 3.

tab3
Table 3: Example 6.3: the maximum absolute errors (with CPU time).

Example 6.4 (dissipative nonlinear wave equation). 𝑒𝑑𝑑=𝑒π‘₯π‘₯+π‘’π‘¦π‘¦βˆ’2𝑒𝑒𝑑+ξ€·2sin(πœ‹π‘₯)sin(πœ‹π‘¦)cos𝑑+2πœ‹2ξ€Έβˆ’1sin(πœ‹π‘₯)sin(πœ‹π‘¦)sin𝑑,0<π‘₯,𝑦<1,𝑑>0,(6.6) with exact solution 𝑒=sin(πœ‹π‘₯)sin(πœ‹π‘¦)sin𝑑. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 4 at 𝑑=1.

tab4
Table 4: Example 6.4: the maximum absolute errors (with CPU time).

Example 6.5 (nonlinear wave equation with variable coefficients). 𝑒𝑑𝑑=ξ€·1+π‘₯2𝑒π‘₯π‘₯+ξ€·1+𝑦2𝑒𝑦𝑦𝑒+𝛾𝑒π‘₯+𝑒𝑦+𝑒𝑑+ξ€Ίπ›Ύπ‘’βˆ’π‘‘(coshπ‘₯coshπ‘¦βˆ’sinh(π‘₯+𝑦))βˆ’1βˆ’π‘₯2βˆ’π‘¦2ξ€»π‘’βˆ’π‘‘coshπ‘₯cosh𝑦,0<π‘₯,𝑦<1,𝑑>0,(6.7) with exact solution 𝑒=π‘’βˆ’π‘‘coshπ‘₯cosh𝑦. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 5 at 𝑑=1 for 𝛾=1, 2 and 5.

tab5
Table 5: Example 6.5: The maximum absolute errors (with CPU time).

Example 6.6 (quasilinear hyperbolic equation). 𝑒𝑑𝑑=ξ€·1+𝑒2𝑒π‘₯π‘₯+𝑒𝑦𝑦𝑒+𝛾𝑒π‘₯+𝑒𝑦+𝑒𝑑+ξ€Ίπ›Ύπ‘’βˆ’π‘‘ξ€»π‘’(sinhπ‘₯coshπ‘¦βˆ’cosh(π‘₯+𝑦))βˆ’1βˆ’π‘‘ξ€·π‘’sinhπ‘₯coshπ‘¦βˆ’2βˆ’π‘‘ξ€Έsinhπ‘₯cosh𝑦3,0<π‘₯,𝑦<1,𝑑>0,(6.8) with exact solution 𝑒=π‘’βˆ’π‘‘sinhπ‘₯cosh𝑦. The maximum absolute errors and CPU time (in seconds) are tabulated in Table 6 at 𝑑=1 for 𝛾=1, 2, and 5.
The order of convergence may be obtained by using the following formula: 𝑒logβ„Ž1ξ€Έξ€·π‘’βˆ’logβ„Ž2ξ€Έξ€·β„Žlog1ξ€Έξ€·β„Žβˆ’log2ξ€Έ,(6.9) where π‘’β„Ž1and π‘’β„Ž2 are maximum absolute errors for two uniform mesh widths β„Ž1 and β„Ž2, respectively. For computation of order of convergence of the proposed method, we have considered errors for last two values of β„Ž, that is, β„Ž1=1/32, β„Ž2=1/64 for two linear problems, and β„Ž1=1/16, β„Ž2=1/32 for all non-linear and quasilinear problems, and results are reported in Table 7.

tab6
Table 6: Example 6.6: The maximum absolute errors (with CPU time).
tab7
Table 7: Order of convergence.

7. Concluding Remarks

Available numerical methods for the numerical solution of second-order quasilinear wave equations are of order four, which require 19-grid points. In this article, using the same number of grid points and five evaluations of the function 𝑔 (as compared to seven evaluations of the function 𝑔 discussed in [6, 7]), we have derived a new stable numerical method of 𝑂(π‘˜2+π‘˜2β„Ž2+β„Ž4) accuracy for the solution of quasilinear wave equation (1.1). Ultimately, we use less algebra for computation, and for a fixed parameter value 𝜎=π‘˜/β„Ž2, the proposed method behaves like a fourth-order method, which is exhibited from the computed results. The proposed numerical method is applicable to wave equation in polar coordinates, and for the damped wave equation and telegraphic equation the method is shown to be unconditionally stable.

Acknowledgments

The authors thank the anonymous reviewer for his valuable suggestions, which substantially improved the standard of the paper. This research was supported by The University of Delhi under research Grant no. Dean (R)/R & D/2010/1311.

References

  1. M. Ciment and S. H. Leventhal, β€œHigher order compact implicit schemes for the wave equation,” Mathematics of Computation, vol. 29, no. 132, pp. 985–994, 1975. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  2. M. Ciment and S. H. Leventhal, β€œA note on the operator compact implicit method for the wave equation,” Mathematics of Computation, vol. 32, no. 141, pp. 143–147, 1978. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  3. T. Schwartzkopff, C. D. Munz, and E. F. Toro, β€œAder: a High-Order Approach for Linear Hyperbolic Systems in 2D,” Journal of Scientific Computing, vol. 17, pp. 231–240, 2002.
  4. T. Schwartzkopff, M. Dumbser, and C. -D. Munz, β€œFast high order ADER schemes for linear hyperbolic equations,” Journal of Computational Physics, vol. 197, no. 2, pp. 532–539, 2004. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  5. M. Dumbser and M. Käser, β€œArbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems,” Journal of Computational Physics, vol. 221, no. 2, pp. 693–723, 2007. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  6. R. K. Mohanty, M. K. Jain, and K. George, β€œHigh order difference schemes for the system of two space second order nonlinear hyperbolic equations with variable coefficients,” Journal of Computational and Applied Mathematics, vol. 70, no. 2, pp. 231–243, 1996. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  7. R. K. Mohanty, U. Arora, and M. K. Jain, β€œLinear stability analysis and fourth-order approximations at first time level for the two space dimensional mildly quasi-linear hyperbolic equations,” Numerical Methods for Partial Differential Equations, vol. 17, no. 6, pp. 607–618, 2001. View at Publisher Β· View at Google Scholar Β· View at MathSciNet
  8. R. K. Mohanty and M. K. Jain, β€œAn unconditionally stable alternating direction implicit scheme for the two space dimensional linear hyperbolic equation,” Numerical Methods for Partial Differential Equations, vol. 22, no. 6, pp. 983–993, 2001. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  9. R. K. Mohanty, β€œAn operator splitting method for an unconditionally stable difference scheme for a linear hyperbolic equation with variable coefficients in two space dimensions,” Applied Mathematics and Computation, vol. 152, no. 3, pp. 799–806, 2004. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  10. M. Dehghan and A. Shokri, β€œA numerical method for solving the hyperbolic telegraph equation,” Numerical Methods for Partial Differential Equations, vol. 24, no. 4, pp. 1080–1093, 2008. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  11. M. Dehghan and A. Mohebbi, β€œHigh order implicit collocation method for the solution of two-dimensional linear hyperbolic equation,” Numerical Methods for Partial Differential Equations, vol. 25, no. 1, pp. 232–243, 2009. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  12. M. Dehghan and A. Shokri, β€œA meshless method for numerical solution of a linear hyperbolic equation with variable coefficients in two space dimensions,” Numerical Methods for Partial Differential Equations, vol. 25, no. 2, pp. 494–506, 2009. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  13. R. K. Mohanty, β€œNew unconditionally stable difference schemes for the solution of multi-dimensional telegraphic equations,” International Journal of Computer Mathematics, vol. 86, no. 12, pp. 2061–2071, 2009. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  14. J. Liu and K. Tang, β€œA new unconditionally stable ADI compact scheme for the two-space-dimensional linear hyperbolic equation,” International Journal of Computer Mathematics, vol. 87, no. 10, pp. 2259–2267, 2010. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  15. Y. Y. Hu and H. W. Liu, β€œAn unconditionally stable spline difference scheme for solving the second 2D linear hyperbolic equation,” in Proceedings of the 2nd International Conference on Computer Modeling and Simulation (ICCMS '10), vol. 4, pp. 375–378, 2010. View at Publisher Β· View at Google Scholar
  16. S. Karaa, β€œUnconditionally stable ADI scheme of higher-order for linear hyperbolic equations,” International Journal of Computer Mathematics, vol. 87, no. 13, pp. 3030–3038, 2010. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  17. R. K. Mohanty and S. Singh, β€œHigh accuracy Numerov type discretization for the solution of one-space dimensional non-linear wave equations with variable coefficients,” Journal of Advanced Research in Computer Science, vol. 3, pp. 53–66, 2011.
  18. M. Lees, β€œAlternating direction methods for hyperbolic differential equations,” vol. 10, pp. 610–616, 1962. View at Zentralblatt MATH
  19. A. R. Gourlay and A. R. Mitchell, β€œA classification of split difference methods for hyperbolic equations in several space dimensions,” SIAM Journal on Numerical Analysis, vol. 6, pp. 62–71, 1969. View at Publisher Β· View at Google Scholar
  20. C. T. Kelly, Iterative Methods for Linear and Non-Linear Equations, SIAM Publications, Philadelphia, Pa, USA, 1995.
  21. Y. Saad, Iterative Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2nd edition, 2003.
  22. L. A. Hageman and D. M. Young, Applied Iterative Methods, Dover Publications Inc., Mineola, NY, USA, 2004.
  23. R. K. Mohanty, β€œStability interval for explicit difference schemes for multi-dimensional second-order hyperbolic equations with significant first-order space derivative terms,” Applied Mathematics and Computation, vol. 190, no. 2, pp. 1683–1690, 2007. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet
  24. M. M. Chawla, β€œSuperstable two-step methods for the numerical integration of general second order initial value problems,” Journal of Computational and Applied Mathematics, vol. 12, pp. 217–220, 1985. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet